使用TOPHAT對轉錄組數據進行比對

編輯

Tophat (http://ccb.jhu.edu/soflware/tophat/index.shtml)主要运用于将 RNA-seq 產生 的reads比對到參考基因組,從而來尋找基因的剪切位點(splicejunction)。該軟件通過調 用Bowtie,或Bowtie2來將reads比對到參考基因組上,分析比對結果,從而尋找出外顯子之間的結合位點。

Tophat相比於Bowtie的優勢在於能對跨越intron的reads進行剪切性比對。

  • 建立在 Bowtie 上並使用相同的基因組索引
  • 用於將 RNA-Seq read與基因組進行比對
  • 優化雙末端,Illumina 序列讀數 >70bp

Tophat的下載和安裝

編輯

推薦下載適合於Linux x86_64的二進制文件壓縮包,解壓縮即可使用。使用Tophat前需要安裝 Bowtie、Bowtie2、Samtool 和 Boost C++ libraries (http://www.boost.org/) 等。其中Boost C++ libraries主要用來幫助編譯Tophat源碼包。如果是直接下載Tophat的二進制文件壓縮包,則可以不需要Boost C++ libraries。

# 下载并安装Tophat:

# installing Boost C++ libraries
tar jxf ~/software/boost_1_57_0.tar.gz
cd boost_1_57_0/
./bootstrap.sh
sudo ./b2 install
cd ..
sudo rm boost_1_57_0 -rf

# installing tophat
wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.1.1.Linux_x86_64.tar.gz -P ~/software/
tar zxf ~/software/tophat-2.1.1.Linux_x86_64.tar.gz -C /opt/biosoft/
echo 'PATH=$PATH:/opt/biosoft/tophat-2.1.1.Linux_x86_64/' >> ~/.bashrc
source ~/.bashrc

# 下载并安装samtools:
wget https://github.com/samtools/samtools/releases/download/1.10/samtools-1.10.tar.bz2 -P 〜/software/
tar jxf 〜/software/samtools-1.10.tar.bz2
cd samtools-1.10/
./configure --prefix=/opt/biosoft/samtools-1.10/
make -j 4
make install
echo 'PATH=$PATH:/opt/biosoft/samtools-1.10/bin/' >> ~/.bashrc
source ~/.bashrc
cd ../ && rm samtools-1.10 -rf

Tophat的使用

編輯

tophat使用的注意事項

編輯

tophat使用的一些注意事項:

1.使用Tohat时,bowtie2(或bowtie,下同), bowtie2-align, bowtie2-inspect, bowtie2-build  samtools,这些命令必须要在系统环境变量PATH路径中。
2.Tophat能比对的最大reads为1024bp。
3.不能将多种不同类型的reads混合起来进行比对,这样会给出不好的结果。比如:单端测序reads 不能和双端测序reads混合;不同插入片段长度的双端测序reads不能混合。
4.如果有多种不同类型的reads进行比对,则需要:首先,对第一种类型的reads使用合适的参数运行tophat;然后,使用Tophat软件包中的程序bed_to_juncs将前一次的运行结果junctions.bed 转换成下一次运行topha所的-j参数所需的junction文件;最后,再次运行tophat对下一种类型的reads使用-j参数进行运行。

tophat的用法和示例

編輯

運行Tophat,需要先使用Bowtie創建index歡據庫。多組插入片段長度一致的paired-end數據一起進行比對的時候,數據文件路徑之間以逗號分開。

tophat的使用
$ tophat [options]* <index_base> <reads1_1[,...,readsN_1]> <reads1_2[,...,readsN_2]>
可以看出,tophat必须要的条件是比对的index数据库,以及要比对的reads。可以为多个paired-end reads数据以逗号分开。

tophat常用例子:
$ tophat -N 5 —read-edit-dist 5 -r 50 --mate-std-dev 20 -p 4 -a 10 -i 20 -1 4000 --min-segment-intron 20 --max-segment-intron 4000 --min-coverage-intron 20 --max-coverage-intron 4000 --coverage-search --microexon-search -G genome_annot.gff3 --library-type fr-firststrand genome reads.1.fq reads.2.fq

tophat的參數與算法

編輯

常用參數

編輯
-h | --help
-v | --version

-N | --read-mismatches  default: 2
    丢弃不匹配碱基数超过该数目的比对结果。最终比对结果中,若不匹配碱基数超过该数目,则舍弃。假若测序错误率为1%,双倍体杂合率为 1%,则read长度为100bp时,不匹配碱基数目为2的比对结果会出现很多。因此,该默认值比较严谨,个人建议将此值设置为5,来增加匹配率。
--read-gap-length  default: 2
    丢弃gap总长度超过该数目的比对结果。最终比对结果中,若gap总长度超过该数目,则舍弃。
--read-edit-dist  default: 2
    丢弃read的edit distance大于该值的比对结果。最终比对结果中,若edit distance大于该值,则舍弃。该值必须大于或等于-N参数的值。
--read-realign-edit-dist <int> default: "read-edit-dist" + 1
    一些跨越多个exons的reads可能会被错误地比对到geneome上。Tophat有多个比对步骤,每个比对步骤过后,比对结果中包含了edit distance的值。该参数能让Tophat对那些edit distance的值 >= 该参数的reads重新进行比对。若设置该参数值为0,则每个read在多个比对步骤中每次都要进行比对。这样会加大地增加比对精确性和运行时间。
    默认下该参数比上一个参数的值大,则表示对reads进行重新比对。
--bowtie1  default: bowtie2
    使用Bowtie1来代替Bowtie2进行比对。特别是使用colorspace reads时,因为只有Bowtie1支持,而Bowtie2不支持。
-o | --output <string>  default: ./tophat_out
    输出的文件夹路径
-r | --mate-inner-dist <int>  default: 50
    成对的reads之间的平均inner距离。例如:ragments长度300bp,reads长度50bp, 则其inner距离为200bp,该值该设为200。
--mate-std-dev <int>  default:20
    inner距离的标准偏差。
-a | --min-anchor-length <int>  default: 8
    read的锚定长度:该参数能设定的最小值为3;锚定在junction两边的reads长度只有都大于此值,才能用于junction的验证。现在Illumina测序的读长越来越长,建议将该值设大一些,比如设为read长度的10%。
-m | --splice-mismatches <int>  default: 0
    对于一个剪切比对,其在锚定区能出现的最大的不匹配碱基数。当-a参数设置较大时候,则需要加大该值了。
-i | --min-intron-length <int>  default:70
    最小的intron长度。Tophat会忽略比该长度要小的donor/acceptor pairs,认为该区属于exon。真核生物中intron的长度最小也要大于10bp,—般情况下,intron长度大 于30bp。个人建议将此参数设置为20。
--I | --max-intron-length <int>  default:500000
    最大的intron长度。Tophat会忽略长度大于该值的donor/acceptor pairs,除非有long read支持。真核生物绝大部分intron的长度小于2000bp。个人建议根据物种的特性将 此参数设置为3000-20000。
--max-insertion-length <int>  defautl: 3
    最大的插入长度
--max-deletion-length <int>  default: 3
    最大的缺失长度
--solexa-quals
    fastq文件使用Solexa的碱基质量格式。不设置此参数和下一个参数,则默认为Phred33。
--solexa1.3-quals | --phred64-quals
    使用Illumina GA pipeline version 1.3的碱基质量格式,即Phred64。
-Q | --quals
    说明是使用单独的碱基质量文件
--inter-quals
    有空格隔开的整数值来代表碱基质量。当使用 -C 参数时,该参数为默认参数。
-C | --color
    Colorspace reads。使用这一种reads的时候命令如下:
$ tophat --color --quals --bowtie1 [other options]* <colorspac
e_index_base> <reads1_1[,...,readsN_1]> <reads2_1[,...,readN_2]>
 <quals1_1[,...,qualsN_1]> <quals1_2[,...,qualsN_2]>
-p | --num-threads <int>  default: 1
    比对reads的线程数
-g | --max-multihits <int>  default: 20
    对于一个reads,可能会有多个比对结果,但tophat根据比对得分,最多保留的比对结果数目。如果没有 --report-secondary-alignments 参数,则只会报告出最佳的比对结果。若最佳比对结果数目超过该参数值,则只随机报告出该数目的最佳比对结果;若有 --report-secondary-alignments 参数,则按得分顺序报告出比对结果,直至达到默认的数目为止。将RNA-seq数据比 对到基因组后,当使用cufflinks计算表达量时,若出现read比对到N个位置,则每个位置的 count数计算为1/N。因此,个人不推荐设置"-g 1"。
--report-secondary-alignments
    是否报告additional or secondary alignments(基于比对分值AS来确定的)。
--no-discordant
    对于paired reads,仅仅报告concordant mappings。
--no-mixed
    对于paired reads,只报告concordant mappings  discordant mappings。默认上,是所有的比对结果都报告。
--no-coverage-search
    取消以覆盖度为基础来搜寻junctions,和下一个参数对立,该参数为默认参数。
--coverage-search
    确定以覆盖度为基础来搜寻junctions。该参数能增大敏感性。
--microexon-search
    使用该参数,pipeline会尝试寻找micro-exons。仅仅在reads长度>=50bp时有效。
--library-type default: fr-unstranded
    Tophat处理的reads具有链特异性。比对结果中将会有XS标签。一般Illumina数据的library-type为 fr-unstranded。
    设置是否为链特异性测序或其种类。非链特异性的RNA-Seq,此值要设为默认的fr-unstranded; 链特异性的RNA-Seq中需要设置为fr-firststrand或fr-secondstrand。链特异性测序的比对结 果中将会有个XS标签进行说明。
    fr-firststrand: right reads (reversed reads)是先被测序的,即测序结果文件中 reads1.fq 中的序列都比对到基因的负义链上。一般Illumina的链特异性测序是这种。
    fr-secondstrand: left reads (forward reads)是先被测序的,即测序结果文件中 reads1.fq  的序列都比对到基因的正义链上。
当有-G参数指定了基因结构时,设置正确此参数则更好。

高級參數

編輯
--segment-length default:25
	将read切成片段,片段长度最小为此长度。这些片段将独立地比对到参考序列。
--segment-mismatches default: 2
	将切碎得到的小片段比对到参考序列上,允许的最小mismach数。
--min-segment-intron default:50
	Split-segment search方法进行intron鉴定,所允许的最小intron长度。建议将此值设置成与-i参数一致。
--max-segment-intron default:500000
	Split-segment search方法进行intron鉴定,所允许的最大intron长度。建议将此值设置成与-I参数一致。
--min-coverage-intron default:50
	Coverage search方法进行intron鉴定,所允许的最小intron长度。建议将此值设置成与-i参数一致。
--max-coverage-intron default:500000
	Split-segment search方法进行intron鉴定,所允许的最大intron长度。建议将此值设置成与-I参数一致。
--keep-tmp
    保留中间文件和临时文件,对于debug有用
--keep-fasta-order
    对比对结果按基因组fasta文件进行排序。该参数会使输出的SAM/BAM文件和tophat的1.41或以前版本不兼容
--no-sort-bam
    输出的BAM文件不是坐标排序coordinate-sorted.
--no-convert-bam
    不要转换成bam格式。输出结果为sam格式。
-R | --resume <string>
    设置该参数输入上个Tophat命令未完成的结果文件夹,程序会从最末尾的成功完成点处,接着运行Tophat。使用方法为:$ tophat -R tophat_out
-z | --zpacker  default:gzip
    用来对临时文件进行压缩的的压缩程序

使用tophat2的時候,可能需要用一些參數來設置bowtie2的運行。這些參數都以』b2』 開頭。其實,這些參數使用默認的即可。

Bowtie2的特別參數

編輯
end-to-end模式(Tophat2不能使用local alignment):
--b2-very-fast
--b2-fast
--b2-sensitive
--b2-very-sensitive

比对参数:
--b2-N  default: 0
--b2-L  default: 20
--b2-i  default: S,1,1.25
--b2-n-ceil  default: L,0,0.15
--b2-gbar  defaut: 4

得分参数:
--b2-mp  default: 6,2
--b2-np  default: 1
--b2-rdg  default: 5,3
--b2-rfg  default: 5,3
--b2-score-min  default: L,-0.6,-0.6

Effort参数
--b2-D  default: 15
--b2-R  default: 2

Tophat進行融合轉錄本搜尋的參數

編輯

如果設定開啟 –fusion-search 參數,則有些reads能比對到潛在的融合轉錄本(fusion transcripts)上。這些融合比對結果使用SAM文件中的XF和XP標籤來標記。額外融合信息保存在 fusions.out 中。比對結束後,可以使用Tophat軟件包中附帶的程序tophat-fusion-post來提取融合轉錄本。

--fusion-search
    开启融合转录子的比对
--fusion-anchor-length  default: 20
    read比对到融合子的两边,每以边至少匹配的碱基数。

提供的轉錄結構註釋數據:值得注意的提供的GTF文中中的染色體名稱和Bowtie index中的一致。這些名稱是區分大小寫的。

-j | --raw-juncs <.juncs file>
    提供junctions文件。该文件可以使用tophat同一目录下的程序bed_to_juncs程序来处理tophat的结果文件junctions.bed生成。其命令为:$ bed_to_juncs <junctions.bed>
    junctions文件是tab分隔的文件,内容为:
<chrom> <leftgt> <rigth> <+/->
其中left和right数值是0-based的junction两端的值。
--no-novel-juncs
		只搜寻和GFF或junctions文件中提供的junctions想匹配的reads。如果没有 -G  -j 参数,则该参数无效。
-G | --GTF <GTF/GFF3 file>
		提供基因模型的注释文件,GTF 2.2 或者 GFF 3 格式的文件。如果设置了该参数,Tophat则先提取出转录子序列,然后使用Bowtie2将reads比对到提取的转录组中;只有不能比对上的reads再比对到genome;将比对上的reads再打断转变成genomic mappings;再融合新的mappings和junctions作为最后的输出。将比对上到虚拟转录组的比对结构转换成genomic mapping 的比对结果;将前两者的比对结果进行融合作为最后的输出。
值得注意的是GTF/GFF文件代表chromosome和contig的第一列要和bowtie index中的参考序列名一致。 使用命令 $ bowtie-inspect --names your_index命令可以获得bowtie的index数据库中的序列名。
--transcriptome-index <dir/prefix>
是使用了 -G 参数后,Tophat提取转录子序列,然后使用bowtie2-build来建立index,这个过程会消耗不少时间。于是,使用该参数,会将index文件生成到指定文件夹。则后续的运用同样的index则不再需要额外耗时了。

使用RNA-seq數據對提供的插入。插入/缺失(insertions/deletions)信息進行驗證:

--insertions | --deletions <.juncs file>
    juncs文件中每行代表一个INDEL。数字代表的位置都是0-based类型。
	缺失的写法:<chrom> <left> <right>从left到right之间的碱基发生缺失。
	缺失的例子:chr1 20564 20567表示第20565号碱基到第20566号碱基(共2个碱基)缺失。
	插入的写法:<chrom> <left> <dummy>〈inserted sequence〉dummy  left 数值一致,表示在此位置前插入了的碱基序列为〈inserted sequence〉。
	插入的例子:chrl 17491 17491 CA表示在参考序列的第17490号碱基和17291号碱基之间插入了2个碱基‘CA’。
    juncs文件例子:
chr1 20564 20567
0-based数值。表示有20565和20566这2个碱基缺失
chr1 17491 17491 CA
表示在17491处插入了2个碱基CA
--no-novel-indels
    仅仅只搜寻在已给的位点的reads。仅仅只搜寻在已给的INDEL位点的reads。如果没有上述参数提供INDELs文件,则不进行INDEL检测。

Tophat的輸出結果

編輯
1. accepted_hits.bam	比对结果文件,默认格式下为sorted的结果。
2. junctions.bed	剪接位点信息,UCSC BED格式。
3. insertions.bed 和 deletions.bed INDELs 信息结果。
4. align_summary.txt 比对结果的统计信息
5. unmapped.bam 未比对上的序列