如果用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的源码。
网友评论