Tophat2比对原理及命令

2018-09-18

相对速度快并且占用内存小的TopHat在RNA-Seq中常常用于跨内含子比对，这里我们来关注一下TopHat2，TopHat2的安装是依赖Bowtie2的(当然Bowtie1也是可行的)，TopHat2适合于75bp以上Read的比对。TopHat2是一个多步骤比对的程序，如果基因组的注释文件存在的话，那么它会先将Read比对到转录组上。这就大大提高了比对的准确性，还可以避免序列匹配到假基因上，加快了整个比对的进程。TopHat2容错率很低，并不会因为没有比对上就而截短Read从而去比对上，所以低质量的碱基比对的效果可能不是那么好。此外，TopHat2可以察觉基因组的易位，将Read比对到潜在的融合转录子上。

TopHat2比对有三个主要的过程，转录组的比对(步骤一)，基因组的比对(步骤二)，剪接比对(步骤三到六，如图4.1显示的一样)，双端测序的Read首先每端数据要分别单独进行比对，然后考虑比对的片段长度以及方向将单独比对结果合并在一起形成双端比对结果。

图4.1Tophat2剪接比对的程序

(a)那些没有比对到基因组或者转录组的Read将会片段化，形成更小的片段，并且再一次比对到基因组上，如果TopHat2发现Read片段化后的左右两个片段的位于用户定义的最大内含子长度范围内，TopHat2就会将Read映射到整个基因组的区域以便于去寻找那些含有剪接信号的剪接位点(b)潜在的剪接位点侧翼基因组序列被串联起来，并构建索引，未映射的Read片段(在这里由星号标记)用Bowtie2比对串联的侧翼区域。(c)被分割的片段重新连接形成完整的Read。

1. 如果注释信息是已知的，那么TopHat2就会将Read比对到转录组上，它是通过一个GTF或是GFF格式的文件从bowtie2建立的索引文件里提取转录本序列。

2. 没有完全比对到转录组的Read再次通过Bowtie2比对到基因组，在这一阶段，那些能够连续比对到单一外显子的Read将会映射到基因组上，而比对到多个外显子的Read将不会映射到基因组上。

3. 那些没有被映射上的Read将会被片段化，通常默认大小为25bp，之后再次比对到基因组上，如果TopHat2发现左右的片段化的Read映射到自定义的最大内含子区长度内，TopHat2就会将Read映射到整个基因组的区域以便于去寻找那些含有剪接信号(GT-AG, GC-AG, or AT-AC))的剪接位点，TopHat2在这一步也能用于寻找缺失以及融合转录子。

4.潜在的剪接位点的侧翼基因组序列被串联起来，并构建索引，未映射的Read片段用Bowtie2比再次比对串联的侧翼序列。

5. 将步骤3和步骤4的小片段拼接在一起完成整个Read的比对。

6. 在步骤2中Read有一部分碱基比对到内含子的将会用新的剪接位点的信息进行重新比对。

7. 在多重比对中为了确定哪一个位置比对是最合适的，TopHat2将会依据匹配到剪接位点的或者indel的Read数量信息重新判分以选择最优的比对。

构建索引

如果你想使用TopHat2，那么你就要像前文描述Bowtie2一样建立基因组的索引，TopHat2一样需要相应的基因组的FASTA格式的文件，所以索引建立好后用不着删除，如果FASTA格式的基因组不与索引在同一个文件夹下，那么Tophat2在运行过程中就会自动生成，所以TopHat2的使用是一个高耗时的。

如果GTF或是GFF格式的基因组注释文件存在的话，那么读长将会优先比对到注释文件的转录组，GTF的文件可在http://www.ensembl.org/info/data/ftp/index.html网站找到，选择对应的生物和GTF选项即可。

如果可行的话，应该事先准备好转录组索引以节约后续比对的时间。

tophat2 -G GRCh37.74.gtf --transcriptome-index= GRCh37.74.tr GRCh37.74

在这里我们使用了注释文件GRCh37.74.gtf和Bowtie2生成的基因组索引GRCh37.74来建立名为GRCh37.74.tr转录组索引，这里需要注意的是，基因组索引里染色体的名称一定要和GTF文件里的相同。在这里Bowtie2生成的基因组索引文件必须用到，所以TopHat2的使用是依赖于Bowtie2的。运行完毕后，生成的文件如：

GRCh37.74.tr.1.bt2

GRCh37.74.tr.2.bt2

GRCh37.74.tr.3.bt2

GRCh37.74.tr.4.bt2

GRCh37.74.tr.fa

GRCh37.74.tr.fa.tlst

GRCh37.74.tr.gff

GRCh37.74.tr.rev.1.bt2

GRCh37.74.tr.rev.2.bt2

GRCh37.74.tr.ver

比对Read

FASTA和FASTQ类型的文件都可以作为TopHat2的输入文件，.gz结尾压缩的Read是可以直接使用的，但是.tgz 或者 .tar.gz结尾的压缩文件是需要被解压后才能使用。下面的例子分别展示对了单端测序和双端测序的序列比对的命令，如果有必要的话，Tophat2也可以将单端测序的Read用于双端测序Read的比对中。

下面两个命令都可以作为单端比对方式，在例子中，序列被比对到人类参考基因组(GRCh37.74)，但是第一个命令用了先前就有的转录组的索引，而第二个命令的索引则是在执行过程中生成的。

tophat2 -o outputFolder --transcriptomeindex=GRCh37.74.tr -p 8 --phred64-quals GRCh37.74 reads1.fastq.gz

tophat2 -o outputFolder -G GRCh37.74.gtf -p 8 --phred64-quals GRCh37.74 reads1.fastq.gz

如果是近些年测序产生的数据，那么TopHat2则默认碱基质量得分为质量字符的ASCII值-33，但是如果测序较早，则应该使用参数(--phred64-quals)，则碱基质量得分为质量字符的ASCII值-64，参数(-p 8)可以设置8个线程加快比对速度。此外TopHat2处理链特异性测序数据，一般Illumina数据的默认library-type为fr-unstranded。TopHat2还有许多比对时可以使用的参数，参数(-T)使得Read仅比对到转录组上，或者通过参数(-g)改变每个Read最大比对到参考基因组的次数，一般情况默认是20次。

The align_summary.txt 显示了 79.3% 的Read映射上了:

Reads:

Input : 34232081

Mapped : 27140089(79.3% of input)

of these: 1612317 (5.9%)have multiple

alignments (2771 have >20)

79.3% overall read mapping rate.

双端测序的序列比对的方式如下，需要注意的是两个文件的Read的顺序应当一致以保证TopHat2的正常运行。如果有多个文件，那么仅需要用逗号将按顺序将文件隔开就行。

tophat2 –o outputFolder ––transcriptomeindex=GRCh37.74.tr –p 8 --phred64-quals GRCh37.74 reads1.fastq.gz reads2.fastq.gz

TopHat2有些参数特意针对应用双端测序得到的序列的比对，，参数(-r)可以依据你的数据设置需要的两比对末端的距离，通常默认是50(样品的插入片段大小为200bp，读长一般为75bp，200-2*75=50，译者注目前大多数为双端150bp)，参数(--no-discordant)可以对于成对的Read，并且只报告方向以及插入片段一致的映射；如果Tophat2不能将成对Read一块映射到基因组，那么默认条件下它将单独映射相应的Read，参数(--no-mixed)可以改变这一默认条件。

Left reads:

Input : 34232081

Mapped : 27143093(79.3% of input)

of these: 1014796 (3.7%)have multiple

alignments (3621 have >20)

Right reads:

Input : 34232081

Mapped : 22600062 (66.0% of input)

of these: 759539 (3.4%)have multiple

alignments (3193 have >20)

72.7% overall read mapping rate.

Aligned pairs: 21229613

of these: 702920 (3.3%)have multiple alignments

336032 (1.6%)are discordant alignments

61.0% concordant pair alignment rate.

TopHat2可以产生多个输出文件

accepted_hits.bam 以BAM文件的形式储存着比对信息，这些比对信息是依据染色体的顺序排列储存的。

junctions.bed 是以Bed格式储存着发现的剪接位点，每个位点都含有左右两个部分，每个部分长度是Read能够跨越剪接位点匹配到基因组的最远距离，最终的得分是跨越剪接位点的Read数量。

insertions.bed 包含着查询到的插入位点信息，chromLeft指的是插入到基因组的位置。

deletions.bed 包含着查询到的缺失位点信息，chromLeft指的是缺失片段的第一个碱基。

align_summary.txt能够显示比对效率以及多少Read能够多处比对。

免责声明：本文仅代表文章作者的个人观点，与本站无关。其原创性、真实性以及文中陈述文字和内容未经本站证实，对本文以及其中全部或者部分内容文字的真实性、完整性和原创性本站不作任何保证或承诺，请读者仅作参考，并自行核实相关内容。

http://www.pinlue.com/style/images/nopic.gif