之前我们说了比对上的数据，那么会有人想到有没有没有比对上的数据呢？
既然是从我的血液里面提取到的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的，代码如下：
虽然上面两个方法得到的结果是一模一样的，但是这个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:
bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq
可以简单的统计一下未比对成功的reads有多少：
结果如下
如果对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讲之后了，敬请期待！
请扫描以下二维码关注我们，获取直播系列的所有帖子！
菜鸟团公众号二维码