2018-12-04 Hisat2 map 结果与 samtools flagstat 结果不一致

如果用Hista2将原始的reads 比对到参考基因组上，会产生一个log，这个log 会显示你比对的结果，我的结果如下：

[will@200server Arb-hisat]$ hisat2 -t -x E_coli -1 SRR1770413_1.fastq.gz -2 SRR1770413_2.fastq.gz -S RNA.sam
Time loading forward index: 00:00:00
Time loading reference: 00:00:00
Multiseed full-index search: 00:01:14
643253 reads; of these:
  643253 (100.00%) were paired; of these:
    240418 (37.38%) aligned concordantly 0 times
    400101 (62.20%) aligned concordantly exactly 1 time
    2734 (0.43%) aligned concordantly >1 times
    ----
    240418 pairs aligned concordantly 0 times; of these:
      1872 (0.78%) aligned discordantly 1 time
    ----
    238546 pairs aligned 0 times concordantly or discordantly; of these:
      477092 mates make up the pairs; of these:
        432923 (90.74%) aligned 0 times
        43549 (9.13%) aligned exactly 1 time
        620 (0.13%) aligned >1 times
66.35% overall alignment rate      #比对率
Time searching: 00:01:14
Overall time: 00:01:14

我们发现比对率是 66.35%

接着使用 samtools 的 flagstat 进行统计，也能得到比对率

[will@200server Arb-hisat]$ samtools flagstat RNA.bam
1299533 + 0 in total (QC-passed reads + QC-failed reads)
13027 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
866610 + 0 mapped (66.69% : N/A)
1286506 + 0 paired in sequencing
643253 + 0 read1
643253 + 0 read2
805670 + 0 properly paired (62.62% : N/A)
812960 + 0 with itself and mate mapped
40623 + 0 singletons (3.16% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

我们发现这里的比对率是 66.69%

差别不是很大，但是为什么同样的文件，统计的比对率是不一样的呢？

经过研究发现，是两者的计算方式不同
我们先来看第一个hisat2的结果

643253 reads; of these:
  643253 (100.00%) were paired; of these:
    240418 (37.38%) aligned concordantly 0 times
    400101 (62.20%) aligned concordantly exactly 1 time
    2734 (0.43%) aligned concordantly >1 times
    ----
    240418 pairs aligned concordantly 0 times; of these:
      1872 (0.78%) aligned discordantly 1 time
    ----
    238546 pairs aligned 0 times concordantly or discordantly; of these:
      477092 mates make up the pairs; of these:
        432923 (90.74%) aligned 0 times
        43549 (9.13%) aligned exactly 1 time
        620 (0.13%) aligned >1 times
66.35% overall alignment rate

它是以pair作为基数的总共有 643253 pair
所以
((400101+2734+1872)X2 + 43549 + 620) / (643253*2) = 66.35%
hista2 计算正确！

然后我们来看 samtools 的统计结果

1299533 + 0 in total (QC-passed reads + QC-failed reads)
13027 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
866610 + 0 mapped (66.69% : N/A)
1286506 + 0 paired in sequencing
643253 + 0 read1
643253 + 0 read2
805670 + 0 properly paired (62.62% : N/A)
812960 + 0 with itself and mate mapped
40623 + 0 singletons (3.16% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

它的map 率计算就是以map 上的reads 除以total 的reads

866610 / 1299533 = 66.69%

计算也完全正确

那么大家可能觉得奇怪，怎么计算都一样的呢

其实差别的根源在于， alignment 有两次，一次是 primarily 一个是 secondary 。

samtools output
1299533 + 0 in total (QC-passed reads + QC-failed reads)
13027 + 0 secondary

如果我们把总的reads 减去 secondary 的

1299533 - 13027 = 1286506
我们发现，这个数字就是 hista2结果pair的两倍

hisat2 result
643253 reads; of these:
  643253 (100.00%) were paired; of these:
    240418 (37.38%) aligned concordantly 0 times
    400101 (62.20%) aligned concordantly exactly 1 time
    2734 (0.43%) aligned concordantly >1 times

643253 X 2 = 1286506

1299533 - 13027 = 1286506 # samtools

643253 X 2 = 1286506 # hisat2

这样两个结果就可以联系起来了

两个结果都没有错，只是统计方法和策略不一样

samtools 用了 secondary 的结果，所以比对率会稍微高一点

PS. 有疑问真的应该多去国外论坛搜索，我中文找了一圈也没找到答案，多亏洲更发我samtools的源码。

最后编辑于：2018.12.04 11:20:56