使用Tophat+cufflinks分析差异表达

使用Tophat+cufflinks分析差异表达

 2017-06-15 19:09:43     522     0     0

                           

使用TopHat+Cufflinks的流程图

 

 

    序列的比对是RNA分析流程中核心的一步。序列的比对，或者说是字符串的比对本身就是计算机科学中的一个经典问题，在生物信息学中更加频繁的出现。序列比对中的错配，插入、缺失可以识别出样本和基因组之间的多态性，甚至可以找出肿瘤样本中的gene fusion。而map到没有注释的基因可能是新的编码基因，或者是非编码RNA。同时RNA-seq的序列比对可以揭示新的选择性剪接和同工型（isoform）。

    此外，序列的比对也可以用作精确定量基因或者转录本的表达量，因为显然表达丰度与产生的reads是直接成比例的（需要标准化）

    Tophat使用Bowtie作为比对引擎，Bowtie使用了一种极其紧凑的数据结构FM index 来储存参考基因组的序列，并且能够迅速的查找（ tens of millions per CPU hour）。但是Bowtie的比对不允许gap的出现（Bowtie2中已经可以允许了），所以Bowtie不能比对跨越了内含子的reads。

    Tophat将Bowtie不能align的reads打断成更小的片段称作segments。一般情况下，当单独处理时，这些segmens可以map到基因组，当一个read的几个segmens可以map到基因组时,Tophat推断这个read跨越了剪接位点同时推测剪接位点的位置。通过处理每个‘initially unmappable’ read，Tophat在没有剪接位点注释的信息下能够在转录组上建立起一个剪接位点的索引。但是这个剪接位点图谱的构建还不够完整，即使是在一些深入研究的模式植物中，每个转录组测序中都能发现新的剪接位点。

Transcript assembly with Cufflinks

    计算基因表达量的另一个问题是，因为选择性剪接的原因，几个不同的转录本（isoform）可能拥有相同的外显子，此时难以确定reads到底来自其中哪个转录本（isoform）。所以能否确定所有的splice variants（isoform）决定着表达量计算的准确性。而这个又很难确定，所以cufflinks通过map到基因组的reads组装起一个简陋的转录组，用reads拼接成含有重叠部分但是长度不同的转录本（称作”transfrags“，作为splice variants的推测。拼接以后，Cufflinks计算使用严格的统计模型来计算每个transfrag的表达量。

    当有多个样本的时候，一种方法是将所有样本的reads合起来，拼接成一个转录组。这种方法的缺点是：

1. 大量reads带来的计算不便，需要更高配置的服务器和大量的时间（可以参考 de no vo Trinity，动辄上百G的内存需求）
2. 多个样本重叠使得确定所有的splice variants（isoform）更加困难

所以Cufflinks采用的策略是先单独拼接每一个样本的reads，然后使用cuffmerge来综合所有样本的拼接结果

    

 Differential analysis with Cuffdiff

Cufflinks还包括一个cuffdiff程序，用于计算基因表达丰度和统计推断。

除了基本的差异表达分析外，cuffdiff还

---

---

运行环境需求：

• Bowtie2
• SAM tools         

要求Bowtie2和SAM tools 都已经安装且添加到系统环境变量中

• Tophat
• Cufflinks
• Cummerbund
• 64-bit linux or Mac os x，最低要求是4G内存，推荐16G内存以上

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125031#references

数据来自于plos one上的一篇文献，测序策略是一小部分测序数据用来拼接转录组（双端），因为当时茶树的基因组还没有公布，只能deno vo拼接。剩下的做定量的单端测序，两个重复。

今年发表了一篇茶树基因组的文章，上传了一个拼接的基因组fasta文件，我们可以用做参考。

TopHat

安装

Tophat在运行中会调用Bowtie1 或者Bowtie2，所以首先要确定你的系统中安装了Bowtie，并且已经添加到了环境变量中

Tophat提供了已经编译好的包可以直接下载，也可以直接下载源码编译，为了方便我们当然选择前一种。网页下载好压缩包或者直接

1. wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.1.1.Linux_x86_64.tar.gz
2. tar -xvf tophat-2.1.1.Linux_x86_64.tar.gz
3. export PATH="${PATH}:/usr/local/biosoft/tophat-2.1.1.Linux_x86_64/"
4. #添加快捷方式
5. ln -s ~/tophat-2.0.0.Linux_x86_64/tophat2 .

解压后将目录添加到环境变量中，或者添加一个快捷链接

为了确保安装正确，最好进行一下测试，尤其是使用编译源码的方式安装的时候。下载test_data

1.     
2.     tar -xvf test_data.tar.gz
3.     cd test_data
4.     tophat -r 20 test_ref reads_1.fq reads_2.fq

结果也没让我们失望，果然失败了。。

研究了一下原因，发现是因为bowtie2的版本太新的原因，换成比较旧的版本就ok了

tophat 官网上最新更新的版本还是2016/02/23的，而且官网上也强调了

Please note that TopHat has entered a low maintenance, low support stage as it is now largely superseded by HISAT2 which provides the same core functionality (i.e. spliced alignment of RNA-Seq reads), in a more accurate and much more efficient way.

 建立基因组索引

对于一些研究比较深入的物种，可以下载已经建立好的索引 using a pre-built index

或者根据你自己的fasta参考文件使用Bowtie2建立，

1. bowtie2-build /home/shady/RNA-seq/reference_genome_tea/Teatree_Assembly.fas Tea_genome

运行成功后产生后缀为.1.bt2, .2.bt2, .3.bt2, .4.bt2, .rev.1.bt2, and .rev.2.bt2  这些文件组成了索引

用我自己的笔记本，大概3G的基因组，build了将近一个小时。。

准备reads

Tophat 可以接受FASTQ or FASTA 格式，推荐fastq

可以用fastqc看一下数据的质量，主要是看看有没有特别差的数据以及接头有没有去除

fastqc有图形界面，也可以使用命令行模式

具体的参数和结果输出请看这篇文章，讲的比较详细

因为是已经发表了的数据，这里只是简单的看一下

1. #直接安装
2. sudo apt-get fastqc
3. fastqc *.gz --outdir=/home/shady/RNA-seq/quality_control --extrct
4. #结果对每个文件都生成一个html文件和一个压缩包 
5. #参数--extrcat 会自动解压压缩包，进去看一下summary.txt
6. cd /home/shady/RNA-seq/quality_control/SRR1747100-25du-4h-1.sra_fastqc
7. cat summary.txt

运行 Tophat 

Tophat首先会调用Bowtie2使用全局比对的方法map你的reads到参考基因组，因为你的reads来自于经过剪接的转录本，所以reads  map到基因组位置堆积起来就像一个参考基因组上的"islands"，许多个这样的"island" 就组成了所谓的外显子

Tophat然后运行一个程序，使用没有map到的reads来寻找剪接位点

一个典型的使用方法为

使用如下脚本处理我们的6个数据

1. #参数 -p 表示运行的线程数
2. #参数 -G 可以选择基因组的注释文件（可选）
3. ls `pwd` |grep "gz"| while read filename
4. do 
5. tophat -p 4 -G Teatree.gff3 -o ${filename%%.*} Tea_genome ${filename}
6. done

参数详解  一般用默认参数

输出

Tophat运行结束后，将在当前目录产生以下几个文件，可以直接在基因组浏览器中查看如IGV等

1. accepted_hits.bam   一个SAM格式的read alignment 列表，SAM是一种紧凑的短序列比对格式文件，里面包含可以map到基因组上的比对信息
2. junctions.bed       一个UCSC BED 文件。是ucsc 的genome browser的一个格式，报告剪接位点 
3. insertions.bed and deletions.bed 

---

参数解释

 常用

o 输出目录，默认值为 “./tophat_out”。
–solexa-quals/solexa1.3-quals 质量编码，关于质量编码格式请参考《FastQ格式介绍》
-p 线程数，默认值为单线程1.，可以使用多线程
-G/–GTFSupply TopHat with a set of gene model annotations and/or known transcripts, as a GTF 2.2 or GFF3 formatted file.指定已有转录本信息
–no-novel-juncs 不查找新的可变剪切
-r 比对时两成对引物间的距离中值。比如说，如果你的插入片段有300bp，而每个引物有50bp，那么r值就应该是200=(300+50*2)/2。没有默认值，如果是末端配对比对时这个值是必须的。
–mate-std-dev 末端配对时中间插入片段的长度的标准差，默认值为20bp

 参考：http://ccb.jhu.edu/software/tophat/manual.shtml#output

 

使用Cufflinks计算差异表达

官方手册 

Cufflinks包含了以下程序：

• Cufflinks 使用经过Tophat比对好的SAM格式来拼接转录组（每一个样本都进行拼接），如果有基因组注释文件的话，效果会更好。这个步骤的目的在前面已经讲过，包括后面的cuffmerge，即为了更精确地确定转录组的所有转录本以及其所有的剪接变体，从而更精确的进行定量分析。

        如果已经有了参考基因组和基因组注释文件，这个步骤和cuffmerge都可以跳过，也就是说      

        cuffquant/ cuffdiff 可以直接接受Tophat输出的文件进行定量分析，差别就是选项中的基因

        组注释文件用你之前下载的基因组注释文件，而不是由cufflinks、cuffmerge生成的注释文

        件，一般不推荐这么做。因为即使是研究比较深入的模式植物，每次转录组测序都会发现新

        的转录本。

• cuffmerge 将Cufflinks拼接的所有转录组整合起来
• cuffquant v2.2.0 后出现的新功能，在之前的cuffmerge 和 cuffdiff 之间插入了一步，作用是计算每个样本的表达量（应该是raw count），产生一个中间文件，这个文件可以用cuffdiff  和cuffnorm分析，是一个可选步骤。
• cuffdiff 可接受Tophat的输出/cuffquant的输出/或者跳过cuffquant

         对每个处理都两两进行差异分析，分析哪些基因上调或者下调，显著性如何等，后面介绍。

• cuffnorm 也是v2.2.0 后新出来的功能，不进行差异分析，只计算FPKM等值，进行标准化

下面脚本进行了从Cufflinks到cuffdiff的一个流程，绝大部分参数使用默认参数即可，其他参数请看官方手册，因为中间涉及几个不同程序，需要要搞清楚每个程序需要的输入文件和输出文件。

1. #运行Cufflinks，分别拼接每一个样本的转录组
2. ls `pwd` |grep 'gz'|while read filename
3. do 
4. cufflinks -p 4 -o "./${filename%%.*}_clout" "${filename%%.*}/accepted_hits.bam"
5. done
6. #创建一个 assemblies.txt
7. #运行cuffmerge，将所有样本拼接好的转录组合成一个
8. #参数 -g 添加基因组注释文件，可选
9. #参数 -s 为参考基因组 fasta格式
10. cuffmerge -g Teatree.gff3 -s Teatree_Assembly.fas -o ./cuffmerge_out/-p 4 assemblies.txt
11. #运行 cuffquant
12.  
13. ls `pwd` |grep 'gz'|while read filename
14. do
15. cuffquant -p 4 -o "./cuffquant_result/${filename%%.*}_cqout/" ./cuffmerge_out/merged.gtf  "${filename%%.*}/accepted_hits.bam"
16. done
17.  
18. #运行cuffdiff
19. cuffdiff -p 4 -o ./cuffdiff_out -L T25_4h,T4_4h,T4_8h ./cuffmerge_out/merged.gtf
20.  ./cuffquant_result/SRR1747100-25du-4h-1_cqout/abundances.cxb
21. ,./cuffquant_result/SRR1747101-25du-4h-2_cqout/abundances.cxb
22.  ./cuffquant_result/SRR1747102-4du-4h-1_cqout/abundances.cxb
23. ,./cuffquant_result/SRR1747103-4du-4h-2_cqout/abundances.cxb
24.  ./cuffquant_result/SRR1747105-4du-8h-1_cqout/abundances.cxb
25. ,./cuffquant_result/SRR1747107-4du-8h-2_cqout/abundances.cxb
26. ~                                                              

最终cuffdiff运行结束后产生的结果如下

产生了一大堆文件，大致分为以下几类

本次分析有3个样本每个样本两个重复，共6个文件

FPKM trcking files 计算的每个样本的FPKM值

其中又分为如下四种，下面的类似

isoforms.fpkm_tracking	Transcript FPKMs
genes.fpkm_tracking	Gene FPKMs. 因为一个基因可以有多个转录本，所以gene FPKM值是具有相同基因id的转录本的Fpkm值的总和
cds.fpkm_tracking	Coding sequence FPKMs. 具有相同蛋白id转录本fpkm值的总和
tss_groups.fpkm_tracking	Primary transcript FPKMs. 具有相同转录起始位点的转录本FPKM值的总和

Count tracking files  计算的每个样本的raw ount数，依旧分为

       

isoforms.count_tracking	Transcript counts
genes.count_tracking	Gene counts. Tracks the summed counts of transcripts sharing each gene_id
cds.count_tracking	Coding sequence counts. Tracks the summed counts of transcripts sharing each p_id, independent of tss_id
tss_groups.count_tracking	Primary transcript counts. Tracks the summed counts of transcripts sharing each tss_id

Read group tracking files  计算的样本里每个重复的表达数 和count数

isoforms.read_group_tracking	Transcript read group tracking
genes.read_group_tracking	Gene read group tracking. Tracks the summed expression and counts of transcripts sharing each gene_id in each replicate
cds.read_group_tracking	Coding sequence FPKMs. Tracks the summed expression and counts of transcripts sharing each p_id, independent of tss_id in each replicate
tss_groups.read_group_tracking	Primary transcript FPKMs. Tracks the summed expression and counts of transcripts sharing each tss_id in each replicate

Differential expression tests

样本之间的两两比较的差异表达和统计推断

isoform_exp.diff	Transcript-level differential expression.
gene_exp.diff	Gene-level differential expression. Tests differences in the summed FPKM of transcripts sharing each gene_id
tss_group_exp.diff	Primary transcript differential expression. Tests differences in the summed FPKM of transcripts sharing each tss_id
cds_exp.diff	Coding sequence differential expression. Tests differences in the summed FPKM of transcripts sharing each p_id independent of tss_id

header的具体信息

Differential splicing tests - splicing.diff

对每个转录前体产生的可变剪接变体的数目

Differential coding output - cds.diff

 在这个文件中只记录了能产生多个CDS的基因的差异表达数

Differential promoter use - promoters.diff

在这个文件中只记录了能产生不同转录前体（i.e. muti-promoter genes）的基因的差异表达

Read group info - read_groups.info

记录了样本重复之间计算定量的一些信息

使用cummeRbund 做一些常见的统计图