【直播】我的基因组(15):提取未比对的测序数据

之前我们说了比对上的数据,那么会有人想到有没有没有比对上的数据呢?

既然是从我的血液里面提取到的DNA进行测序的,那么理论上测序仪出来的所有测序reads都应该是我的,也应该都可以比对到人类的参考基因组,但是实际过程中的确存在着未必对上的数据。

现有人类基因组毕竟只是个参考,也许我有某些独特的DNA序列呢?而且也不一定测序数据就都是人类的,也许我血液里面会有那么些微不纯粹呢?也有可能是某些片段变异的太多了,超过了常见的比对软件的承受能力,可以试用SHRIMP这个软件来提取一下。当然,这不是本讲的重点。

前面我们已经详细讲解了sam文件的格式,就是为了给这个做铺垫,如果还不清楚的,可以回过头再仔细阅读(http://genome.sph.umich.edu/wiki/SAM)。sam格式文件的第3和第7列,可以用来判断某条reads是否比对成功到了基因组的染色体,左右两条reads是否比对到同一条染色体。有两个方法可以提取未比对成功的测序数据,sam文件的第3列是*的(如果是PE数据,需要考虑第6,7列),或者sam文件的flag标签包含0x4的,代码如下:

samtools view -f4 sample.bam > sample.unmapped.samsamtools view sample.bam |perl -alne '{print if $F[2] eq "*" or $F[5] eq "*" }' > sample.unmapped.sam

虽然上面两个方法得到的结果是一模一样的,但是这个perl脚本运行速度远远比不上上面的samtools自带的参数。

sam文件的说明书里面有这样一句话;https://samtools.github.io/hts-specs/SAMv1.pdf

An unmapped segment without coordinate has a ‘*’ at this field.  However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting. If RNAME is ‘*’, no assumptions can be made about POS and CIGAR.

(其实也不一定要自己写脚本,我们前面讲到的用来把巨大的bam文件按照染色体分割的小软件bamtools也可以完成这个需求,用 bamtools -split -in my.bam -mapped 即可!)

小写的f是提取,大写的F是过滤。因为我们测序数据的双端的,那么sam文件的第3列是reads1的比对情况,第6列是reads2的比对情况。所以未比对成功的测序数据可以分成3类,仅reads1,仅reads2,和两端reads都没有比对成功。

也可以用下面的代码分步提取这3类未比对成功的reads:

samtools view -u  -f 4 -F264 alignments.bam  > tmps1.bam samtools view -u -f 8 -F 260 alignments.bam  > tmps2.bam samtools view -u -f 12 -F 256 alignments.bam > tmps3.bamsamtools merge -u - tmps[123].bam | samtools sort -n - unmapped

bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq

可以简单的统计一下未比对成功的reads有多少:

cut -f 3,6 P_jmzeng.unmapped.sam |sort |uniq -c >unmapped.counts

结果如下

如果对bamtools软件的结果来统计:

samtools view P_jmzeng.final.REF_unmapped.bam |cut -f 3,7 |sort |uniq -c >unmapped.counts

得到的结果只有7492014 * * 说明它只考虑了PE reads均为比对成功的情况。

很奇怪,看起来我的未比对成功的测序数据里面竟然没有右端成功,而左端失败的情况,这个我没办法解释。 我也还需要学习才能搞明白这件事。

(其实之前我也搞错了,如果PE reads的左右两端均没有比对成功,那么第3,6,7列都是*,4,5,8,9都是0,第2列flag只有77,141这两种情况。(77代表PE,而且PE的两条reads都是unmanned的,141跟77一样,只是它们分别指代unmanned的的PE的reads的两端,结合https://broadinstitute.github.io/picard/explain-flags.html来理解)

如果是左右两端reads只有一个比对成功,另一个reads没有比对上,如果是read1比对了,read2失败了,那么第3列应该是read1的染色体,第7列应该是*号表明read2没有比对成功。同理,如果read2比对成功,read1失败,按照道理,我们应该看得第7列有染色体,第3列是*号,但是我们在提取的unmapped文件里面,没有发现这种情况。

但其实不管是左端还是右端,第3列都是有染色体的,第7列是=号,但这并不能说明左端跟右端有着同样的比对结果。而第6列CIGAR是*,这个才是判断左右端是否匹配失败的标准。

sam文件的说明书里面有这样一句话;https://samtools.github.io/hts-specs/SAMv1.pdf

对于第6列CIGAR来说,An unmapped segment without coordinate has a ‘*’ at this field. However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting. If RNAME is ‘*’, no assumptions can be made about POS and CIGAR.

这也就是我为什么没有发现第7列有染色体,第3列是*号的reads。即使PE reads的右端匹配,左端未匹配,它只会把这个read比对的染色体写在第3列,而不是第7列!所以说要想探究它是左端还是右端未比对成功,得看flag。

这样就提取出来了未比对的测序数据,但是还需要做进一步分析,看看这些reads究竟是何方大神!具体要在第25讲之后了,敬请期待!

请扫描以下二维码关注我们,获取直播系列的所有帖子!

菜鸟团公众号二维码

菜鸟团公众号二维码

Comments are closed.