使用BWA比对环境转录组至参考基因组时SAM文件损坏问题

背景与问题 #

我有经过质控的双端环境转录组数据，需要比对到由8个蓝藻基因组fasta文件合并而成的参考数据库。使用BWA v0.7.3a执行比对后，生成的SAM文件出现损坏，导致用htseq-count计数映射到各基因的reads时出错。

所用命令 #

BWA比对命令 #

bwa mem -k 10 -aM -t 128 Merged_reference Input_file1 Input_file2 > Out.sam

HTSeq-count计数命令 #

htseq-count -f sam -r pos -a 0 -m intersection-strict -s no Out.sam Merged_ref2.gtf2 > HTSeq_out

报错信息 #

[W::sam_parse1] mapped query cannot have zero coordinate; treated as unmapped
[E::sam_parse1] CIGAR and query sequence are of different length
[W::sam_read1_sam] Parse error at line 8839134
Error occured when processing input (record #8809785 in file Out.sam):
truncated file
[Exception type: OSError, raised in libcalignmentfile.pyx:1876]

SAM文件中的异常行 #

出错的第8839134和8839135行内容如下：

A00521:380:HFKMLDSX5:4:1163:32081:25332 65      lcl|NC_013771.1_cds_WP_012953431.1_1    0       0       268435455S89S   lcl|NC_013771.1_cds_WP_012954293.1_845  87      0       Sequence   Quality     NM:i:8388607    AS:i:73 XS:i:72

A00521:380:HFKMLDSX5:4:1163:32081:25332 129     lcl|NC_013771.1_cds_WP_012954293.1_845  87      0       2S15M72S        lcl|NC_013771.1_cds_WP_012953431.1_1    0       0      Sequence  Quality  NM:i:0  AS:i:15 XS:i:15

临时修复方案 #

将上述异常行修改为未比对标记后，对数据子集运行htseq-count可消除错误：

A00521:380:HFKMLDSX5:4:1163:32081:25332 65      *    0       0       *   *      0     0 Sequence   Quality     NM:i:8388607    AS:i:73 XS:i:72

A00521:380:HFKMLDSX5:4:1163:32081:25332 129     *    0       0       *   *      0     0      Sequence  Quality  NM:i:0  AS:i:15 XS:i:15

求助 #

虽然标记为未比对reads能解决问题，但我的SAM文件约有1.15亿行，且有60多个这类文件，错误每400万行出现一次，手动修改完全不现实，希望找到根本原因来高效解决问题。每次错误均由参考文件中的第一个序列lcl|NC_013771.1_cds_WP_012954293.1_845触发，有没有人遇到过这个问题？实在搞不懂BWA为什么会出现这种错误！

内容的提问来源于stack exchange，提问作者Clayton Tracey