软件7 —— bowtie和bowtie2

一、 基本介绍

(1) DNA-seq和RNA-seq

基因组测序（DNA-seq）和转录组测序（mRNA-seq）的区别在于，真核生物的RNA经过剪接后，序列与DNA序列并不完全相同。如果是mRNA测序，那map工作就会在DNA测序map的基础上再多一步，map到转录组上去。最为流行的做法是，使用BWA来比对ChIP-seq测序，使用bowtie来map DNA测序，使用tophat来map RNA测序。

(2) bowtie和bowtie2

bowtie1出现的早，所以对于测序长度在50bp以下的序列效果不错，而bowtie2主要针对的是长度在50bp以上的测序的。Bowtie 2对最长序列没有要求，但是Bowtie 1最长不能超过1000bp（bowtie最长能读取1024个碱基的片段，模板最小尺寸不能小于1024碱基，而短序列最长而不能超过1024碱基）。Bowtie 2支持有空位的比对，支持局部比对，也可以全局比对，但是Bowtie2不能比对colorspace reads。

Bowtie2是将测序reads比对到参考序列上的工具，适用于长度大约为50到1000bp，是ChIP-seq、BS-seq等分析的第一步。Bowtie2使用FM索引（基于Burrows-Wheeler Transform，BWT）对基因组进行索引，以此来保持其占用较小内存。Bowtie2支持间隔、局部和双端比对模式。可以同时使用多个处理器来极大的提升比对速度。

(3) bowtie和tophat

Tophat/Tophat2本身不能进行比对，它是通过调用bowtie/bowtie2进行比对的。tophat1和tophat2的差别最主要的就是调用了bowtie1还是bowtie2。tophat调用bowtie，tophat2调用bowtie2，当然如果你只安装了bowtie1的话，tophat2也可以调用它。bowtie2不是bowtie的升级版，但是Tophat2是Tophat的升级版。因此Tophat只可以调用bowtie，而Tophat2不仅可以调用bowtie2（默认）还可以更改设置调用bowtie。

bowtie和bowtie2各自有自己的官网，无论是bowtie还是bowtie2在2019年仍然是有更新的。Tophat的网站只有一个，仅仅是在2012年4月9日Tophat发布了2.0.0版本，宣布支持bowtie2的比对，而我们通常也将支持bowtie2版本的Tophat称之为Tophat2。Tophat到了2016年2月便停止了更新，其原作者们不再更新Tophat2，转而开发了一个新的比对工具HISAT2，推荐人们使用HISAT2，声称其速度更快，内存占用率更小，准确率更高。此外，HISAT2不仅支持RNA-seq的比对还支持DNA-seq比对，唯一需要做的就是加上一个参数--no-spliced-alignment。但是就目前来看，大部分人都是使用HISAT2做RNA-seq，没人使用它做DNA-seq。

Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis文中给出了一个它认为比较好的RNA-seq组合方式：比对用HISAT2；reads count方法定量，采用FeatureCounts；差异表达选用DESeq2。

二、 背景知识

(1) 比对时需要解决RNA剪接的问题

真核生物的基因为断裂基因，转录得到的pre-mRNA需要经过剪接将内含子序列移除，形成仅保留外显子序列的成熟mRNA。并且，在高等生物中常常存在选择性剪接，同一基因转录产生的前体mRNA通过不同的剪接方式可以形成不同的成熟mRNA。因此RNA-Seq的比对需要解决外显子、内含子和可变剪接的问题。

(2) 算法比较

• TopHat采用外显子优先（exon-first）的方法，第一步使用BWT方法将完整的匹配外显子的reads优先回贴到参考基因组，第二步将没有比对到基因组的reads切割成小片段，分别比对到基因组后，在比对区域附近寻找有没有可能的剪接位点（GT-AG、GC-AG、AT-AC）而判断该read是否跨越了内含子。这种策略会受到假基因的干扰，会有很多基因比对到假基因上。
• Hisat2使用了叠加的比对算法，全局比对整个基因组建立index，辅助成千上万个小的局部index。这两种算法也是基于BWT/FM体系建立起来的。相比于Tophat2，速度提升了50倍。
• STAR采用了种子和扩展（seed and extend）的方法，首先将所有reads分割成k-mer片段并建立查找表，然后将k-mers片段映射到基因组，在种子位置附近扩展，直至完整测序片段。与外显子优先的方式相比，种子和扩展的方法在计算机资源和运行速度方面不占优势，但却不依赖外显子定位，能够更准确地处理跨越剪接位点的reads，更好地识别可能存在的新的剪接位点。
• 另外，已经开发出了计算高效的免比对（alignment-free）工具，例如Sailfish、Kallisto与Salmon，这些工具可以直接将测序读长与转录本进行关联，从而无需单独的定量步骤。这些工具在那些表征更高丰度（以及更长的）转录本方面表现得非常良好；然后它们在那些定量低丰度或短转录本方面表现不佳。

(3) 特点比较

• Tophat：早期转录组数据分析的经典策略是tophat + cufflinks + cuffdiff。但是，tophat的速度很慢，基于FPKM值得到的差异结果和实验手段如qPCR验证的一致性较差。
• HISAT2：是tophat的升级版本，HISAT2 + StringTie + ballgown。利用改进的BWT算法进行序列比对。HISAT2比STAR快2.5倍，比TopHat快大约100倍。HISAT2找到junction正确率最高，但是在总数上却比TopHat和STAR少。从这里可以看出HISAT2的二类错误（纳伪）比较少，但是一类错误（弃真）就高起来。
• STAR：就唯一比对而言，STAR是三者最佳的，主要是因为它不会像TopHat和HISAT2一样在PE比对不上的情况还强行把SE也比对到基因组上。而且在处理较长的read和较短read的不同情况，STAR的稳定性也是最佳的。STAR有处理较长reads的能力。非常适合于大量数据的并行计算，速度非常快。对于同时有参考基因组和参考转录组的物种，比对的准确率很高，不过index很大，至少需要30G以上内存才能运行。

三、 使用方法

(1) 建立索引

bowtie-build  GENOME.fa  GENOME 
bowtie2-build  species.fa  species 

产生六个索引文件。注意，bowtie和bowtie2的索引不通用。
索引文件也可以直接从官网下载：ftp://ftp.ccb.jhu.edu/pub/data/bowtie_indexes/或ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/。
bowtie的索引是：GENOME.1~4.ebwt和GENOME.rev.1~2.ebwt
bowtie2的索引是：species.1~4.bt2和species.rev.1~2.bt2

bowtie-build  --threads 4  -f  genome.fa  index/bowtie_dm/bowtie_dm  2>>log/bowtie_build_log.txt  &
bowtie2-build  --threads 4  -f  genome.fa index/bowtie2_dm/bowtie2_dm  2>>log/bowtie2_build_log.txt  & 

--threads：线程
-f：参考基因组

(2) 序列比对

bowtie  [options]*  <ebwt>  {-1 <m1> -2 <m2> | --12 <r> | <s>}  [<hit>]
bowtie  -p 4  -x GENOME  -1 input_1.fq  -2 input_2.fq  -S output 
bowtie  -p 4  -x GENOME  -U input.fq  -S output 

bowtie  -p 4  -x  index/bowtie_dm/bowtie_dm  SRR1449909.fastq  -S SRR1449909.sam  2>>log/bowtie_log.txt  &

<ebwt>：由bowtie2-build所生成的索引文件的前缀。
-1 <m1> -2 <m2>：双末端测序的两个fq文件。
<s>：单末端测序的一个fq文件。
<hit>：File to write hits to (default: stdout)
-q：query input files are FASTQ .fq/.fastq (default)
-f：query input files are (multi-)FASTA .fa/.mfa
-p：线程数
--phred33-quals：input quals are Phred+33 (default)

bowtie2  [options]*  -x <bt2-idx>  {-1 <m1> -2 <m2> | -U <r>}  [-S <sam>]
bowtie2  -p 4  -x species  -f  -1 reads_1.fa  -2 reads_2.fa  -S bowtie2_output.sam
bowtie2  -p 4  -x species  (-q)  -1 input_1.fq  -2 input_2.fq  |  samtools sort  -@ 4  (-O bam)  -o -  >  output.bam

-x <bt2-idx>：由bowtie2-build所生成的索引文件的前缀（去掉了.X.bt2）。
-1 <m1>：双末端测序对应的文件1。可以为多个文件，并用逗号分开，中间没有空白；多个文件必须和 -2 <m2> 中制定的文件一一对应。比如 -1 flyA_1.fq,flyB_1.fq -2 flyA_2.fq,flyB_2.fq。 测序文件中的reads的长度可以不一样。可以是gz、bz2压缩文件。
-2 <m2>：双末端测序对应的文件2。
-U <r>：单末端测序对应的文件。可以为多个文件，并用逗号分开。测序文件中的reads的长度可以不一样。
-S <sam>：File for SAM output (default: stdout)
-q：query input files are FASTQ .fq/.fastq (default)
-f：query input files are (multi-)FASTA .fa/.mfa。当设置-f时，--ignore-quals也被设置了。
-p：线程数
--phred33：输入的碱基质量等于ASCII码值加上33。
--local：local alignment; ends might be soft clipped (off)。在这种模式下，Bowtie 2不要求整个read从一端到另一端对齐。

(3) colorspace选项

./shellscript/bowtie-0.12.3/bowtie-build  -C  -f  genome.fa  index/bowtie_c/bowtie_c  &
./shellscript/bowtie-0.12.3/bowtie  -p 4  -C  index/bowtie_c/bowtie_c  SRR1449909.fastq  -S SRR1449909.sam  2>>log/bowtie_c_log.txt  &

-C/--color：旧版本bowtie-0.12.3支持colorspace数据的比对，参数稍有不同，索引也是需要单独构建的。

(4) tophat使用方法

tophat  [options]*  <bowtie_index>  <reads1_1[,…,readsN_1]>  [reads1_2,…readsN_2]
tophat  -o ./tophat_out  -p 4  GENOME  A.reads1_1.fastq,B.reads1_1.fastq  A.reads1_2.fastq,B.reads1_2fastq 

主要的结果文件是：accepted_hits.bam、junctions.bed、insertions.bed和deletions.bed。
<bowtie_index>：参考基因组的index文件的具体目录，参考基因组应该和index文件放在同一目录中。index数据文件需要给出目录以及目录文件的共同前缀，例如，index文件存放在当前目录下的index文件夹，文件的名字是hg19.*.*, index数据的文件应该是：./index/hg19。
reads：PE reads必须放在不同的两个文件中，文件名必须按照*_1,*_2的规范成对出现。如：A.reads1_1.fastq,B.reads1_1.fastq A.reads1_2.fastq,B.reads1_2fastq。
-o/--output-dir：输出文件夹，默认./tophat_out。
-p：线程数
-G：提供基因模型的注释文件GTF2或者GFF3。如果设置了该参数，Tophat则先提取出转录本序列，然后使用Bowtie将reads比对到提取的转录组中；只有不能比对上的reads再比对到genome；比对上的reads再打断转变成genomic mappings；再融合新的mappings和junctions作为最后的输出。注释文件的第一列必须与所建索引的参考序列名称一样。可以使用bowtie-inspect --names your_index进行查看。
--bowtie1：对于tophat2来说，设置使用bowtie1比对，默认使用bowtie2。