5款流行比对工具大比拼

之所以有这个需求，是因为在做一个单细胞转录组测序数据的分析，发现一些非常诡异的比对情况。比如：尽管是全局比对率达到80%以上，但是过半数居然左右失衡，也就是说只有双端测序reads的其中一条成功比对了，另一条莫名其妙的比对失败。这种情况可以发生，只是超过50%就太夸张了，非常值得探究，究竟是reads本身的特性呢，还是比对工具的选择不够合适。

5个比对软件

我这里探究了hisat2,subjunc,star,bwa,bowtie2这5个比对工具，它们的安装方式文件不赘述了，使用代码如下:

1. hisat='/home/jianmingzeng/biosoft/HISAT/current/hisat2';

2. star="/home/jianmingzeng/biosoft/STAR/STAR-2.5.3a/bin/Linux_x86_64/STAR";

3. subjunc="/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/subjunc";

4. tophat2="/home/jianmingzeng/biosoft/TopHat/current/tophat2";

5. bowtie="/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2";

6. hisat2_mm10_index='/home/jianmingzeng/reference/index/hisat/mm10/genome';

7. star_mm10_index='/home/jianmingzeng/reference/index/star/mm10/';

8. subjunc_mm10_index='/home/jianmingzeng/reference/index/subread/mm10';

9. bowtie2_mm10_index='/home/jianmingzeng/reference/index/bowtie/mm10';

10. bwa_mm10_index='/home/jianmingzeng/reference/index/bwa/mm10';

11. fq1=SC2-01_S1_L002_R1_001.fastq.gz

12. fq2=SC2-01_S1_L002_R2_001.fastq.gz

13. sample=SC2-01_S1_L002

14. $hisat -p 5 -x $hisat2_mm10_index -1 $fq1 -2 $fq2 -S $sample.sam 2>$sample.hisat.log

15. samtools sort -O bam -@ 5  -o ${sample}_hisat.bam $sample.sam

16. $subjunc -T 5  -i $subjunc_mm10_index -r $fq1  -R $fq2 -o ${sample}_subjunc.bam

17. ## 很有趣，虽然subjunc直接输出bam，但是是按照reads名字排序，不是按照基因组坐标排序。

18. bwa mem -t 5 -M  $bwa_mm10_index $fq1 $fq2 1>$sample.sam 2>/dev/null

19. samtools sort -O bam -@ 5  -o ${sample}_bwa.bam $sample.sam

20. $bowtie -p 5 -x $bowtie2_mm10_index -1 $fq1  -2 $fq2 | samtools sort  -O bam  -@ 5 -o - >${sample}_bowtie.bam

21. ## star软件载入参考基因组非常耗时，约10分钟，也比较耗费内存，但是比对非常快，5M的序列就两分钟即可

22. $star --runThreadN  5 --genomeDir $star_mm10_index --readFilesCommand zcat --readFilesIn  $fq1 $fq2 --outFileNamePrefix  ${sample}_star

23. ## --outSAMtype BAM  可以用这个参数设置直接输出排序好的bam文件

24. samtools sort -O bam -@ 5  -o ${sample}_star.bam ${sample}_starAligned.out.sam

把比对后的bam文件统一进行 samtools flagstat统计，整理一下比对情况如下表：

bowtie2	bwa	subjunc	hisat	star
4,457,852	5,028,723	4,457,852	4,636,906	3,732,166	in total (QC-passed reads + QC-failed reads)
0	570,871	0	179,054	327,880	secondary
0	0	0	0	0	supplementary
0	0	0	0	0	duplicates
2,846,896	4,898,291	3,694,805	3,418,483	3,732,166	mapped(比对率)
4,457,852	4,457,852	4,457,852	4,457,852	3,404,286	paired in sequencing
2,228,926	2,228,926	2,228,926	2,228,926	1,702,143	read1
2,228,926	2,228,926	2,228,926	2,228,926	1,702,143	read2
1,977,832	3,741,658	2,284,400(太低了)	3,014,912	3,404,286	properly paired（双端reads比对率）
2,568,660	4,320,816	3,531,396	3,150,876	3,404,286	with itself and mate mapped
278,236	6,604	163,409	88,553	0	singletons
65,730	102,160	41,864	8,860	0	with mate mapped to a different chr
25,582	77,062	41,864	6,474	0	with mate mapped to a different chr

让我更不能理解的是，这个是RNA-seq数据，可是居然用BWA可以成功比对率是最高的，虽然比对文件里面有着大量的H/S情况，但是不影响它仍然可以比对成功。高达97.41% 的比对率，还有83.93% 左右两端reads匹配的比对率，遥遥领先与其它比对工具。完全颠覆了我以前的想象，一直以为对转录组数据不能用bwa来比对。(大家可以思考一下其中的原理) 可以看到bwa默认是容许多比对情况的，一条reads可以输出多个比对记录。

1. ## samtools flagstat SC2-01_S1_L002_bwa.bam

2. 5028723 + 0 in total (QC-passed reads + QC-failed reads)

3. 570871 + 0 secondary

4. 0 + 0 supplementary

5. 0 + 0 duplicates

6. 4898291 + 0 mapped (97.41% : N/A)

7. 4457852 + 0 paired in sequencing

8. 2228926 + 0 read1

9. 2228926 + 0 read2

10. 3741658 + 0 properly paired (83.93% : N/A)

11. 4320816 + 0 with itself and mate mapped

12. 6604 + 0 singletons (0.15% : N/A)

13. 102160 + 0 with mate mapped to a different chr

14. 77062 + 0 with mate mapped to a different chr (mapQ>=5)

而subjunc的比对结果里面的双端reads比对率与全局reads的比对率始终是太诡异了，如下，明明还有82.88%的比对率，但是左右两端匹配的比对率下降到51.24% ，理论上我应该好好探究一下各种原因。不过，本文先放过这个细节。(欢迎大家留言提出自己的见解)

1. ## samtools flagstat SC2-01_S1_L002_subjunc.bam

2. 4457852 + 0 in total (QC-passed reads + QC-failed reads)

3. 0 + 0 secondary

4. 0 + 0 supplementary

5. 0 + 0 duplicates

6. 3694805 + 0 mapped (82.88% : N/A)

7. 4457852 + 0 paired in sequencing

8. 2228926 + 0 read1

9. 2228926 + 0 read2

10. 2284400 + 0 properly paired (51.24% : N/A)

11. 3531396 + 0 with itself and mate mapped

12. 163409 + 0 singletons (3.67% : N/A)

13. 41864 + 0 with mate mapped to a different chr

14. 41864 + 0 with mate mapped to a different chr (mapQ>=5)

而且这里star软件取了一个巧，只在bam文件里面输出那些成功比对的reads，这样比对率就都是100%啦，不过也可能是我对star软件的说明书没有吃透，参数不恰当。

而且我还应该再探究一下那些被subjunc比对不上的序列却能被bwa软件成功比对到参考基因组的是写什么情况。这个也以后再分享。

表达量探索

因为转录组测序最重要是得到表达量，我们并不关心这些reads具体比对到参考基因组是被切掉了还是跨越了内含子的成功比对 (如果要考虑转录本定量，或者可变剪切，就需要真正的跨越外显子的比对)

这里我们还是选用最开始便捷的featureCounts来进行基于基因的定量，代码如下；

1. featureCounts='/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts';

2. mm10_gtf='/home/jianmingzeng/reference/gtf/gencode/gencode.vM12.annotation.gtf';

3. # Summarize paired-end reads and count fragments (instead of reads):

4. $featureCounts -T 5 -p -t exon -g gene_id \

5. -a $mm10_gtf -o counts.txt   *.bam

这里只有两个reads都成功比对到参考基因组才计数，所以数的是fragment数量，比对总结如下：

Status	bowtie.bam	bwa.bam	hisat.bam	star.bam	subjunc.bam
Assigned	1,180,215	2,183,516	1,295,795	1,311,956	1,518,380
Unassigned_Unmapped	666,360	61,914	564,935	0	299,819
Unassigned_MappingQuality	0	0	0	0	0
Unassigned_Chimera	0	0	0	0	0
Unassigned_FragmentLength	0	0	0	0	0
Unassigned_Duplicate	0	0	0	0	0
Unassigned_MultiMapping	0	0	165,425	254,209	0
Unassigned_Secondary	0	0	0	0	0
Unassigned_Nonjunction	0	0	0	0	0
Unassigned_NoFeatures	282,738	408,715	228,180	232,006	307,639
UnassignedOverlappingLength	0	0	0	0	0
Unassigned_Ambiguity	99,613	145,544	73,526	67,912	103,088

可以看到bwa软件把几乎所有的RNA-seq测序的reads都成功的比对到了基因组，并且成功的注释到了基因，计数成为了基因的表达量。从这一点来说，BWA最大限度的保留了这个RNA-seq测序得到的基因表达信息，现在只需要探究它们的表达量的相关性，看看BWA那种截取reads的部分片段进行比对，是否会出现大幅度的基因定量的偏差。

简单的R代码探究如下：

1. a=read.table("tmp.txt",header=T) ## 49585 genes in gencode

2. head(a)

3. cor(a)

4. write.table(cor(a),file="tmp.sum",quote=F,sep="\t")

5. cor(a[order(apply(a,1,mad), decreasing = T)[1:50],])

6. cor(a[order(apply(a,1,mad), decreasing = T)[1:500],])

7. dim(a[rowSums(a)>10,]) ## 5823 expressed genes in single cells

8. cor(a[rowSums(a)>10,])

相关性结果如下：

	bowtie.bam	bwa.bam	hisat.bam	star.bam	subjunc.bam
bowtie.bam	1	0.885440896	0.879128734	0.862654678	0.910767788
bwa.bam	0.885440896	1	0.921680804	0.923104636	0.954310013
hisat.bam	0.879128734	0.921680804	1	0.996418718	0.971852796
star.bam	0.862654678	0.923104636	0.996418718	1	0.970620816
subjunc.bam	0.910767788	0.954310013	0.971852796	0.970620816	1

发现只有bowtie2这个软件的比对跟其它工具得到的表达量相关性比较差。hisat,star,subjunc这3款软件都是设计来专门做转录组数据分析，所以它们之间的相关性非常高。

虽然相关性非常高，也还是有些基因区域被这5款软件计数后的表达量差异却非常大，如下：

上表中最后五列是不同比对得到的表达量，分别是bowie2,bwa,hisat2,star,subjunc。 可以看到对于这种单外显子基因来说，bwa和bowtie2比较类似，但是跟hisat2,star,subjunc这3款转录组比对工具完全不一样。真奇怪。

这里我提出了3个问题，但是代码都给你了，你只需找个fastq文件就可以开始玩啦，希望你可以跟我讨论这3个问题背后的逻辑。

还有，markdown表格可能需要在浏览器打开，效果好一点。