RNA-seq差异表达分析工作流程

RNA-seq差异表达分析工作流程  

原文地址： http://blog.sina.com.cn/s/blog_74cbb8e80101gqkc.html

 

需要使用到的工具：
TopHat
Cufflinks
Samtools

参考：http://vallandingham.me/RNA_seq_differential_expression.html

首先安装的tophat需要事先安装好bowtie。至于安装方面的问题，这里不至赘述。

整个pipeline非常明确：Sequences → TopHat → Manual Check → Cufflinks → Analysis

第一个问题，是否需要做duplicate removal，如果要做，什么时候做？在回答这两个问题之前，我们还是先来看看什么是duplicate。我们将deep sequence中完全相同的序列统称为duplicate。通常这种重复会有几个来源，一，测序模板中存在一模一样的片断；二，测序过程中PCR产生的重复；三，信号读取过程中读到了同一pcr产物。按照这里的 讨论，对于 copy number detection, SV detection, ChIP-seq, and RNA-seq都应该做duplicate removal。去除的优点是可以大量的减少计算，降低假阳性。但是去除的话也有造成数据大量损失的风险，也就是说会降低真阳性结果。有文章对 相同的library做了两次测序，一次是single end, 一次paired end。比较发现，SE的duplicate高达28%，而PE的duplicate只有8%。当把PE的结果当成SE结果来处理时，duplicate 又升至28%。还有些私下的讨论认为，实际的duplicate应该只有1%左右。这里强调了去除duplicate对于数据完整性的影响。那么为什么人们在做CN/SV/ChIP-seq/RNA-seq的时候倾向于做duplicate removal呢？这主要的理论依据是在准备library的步骤中，所有模板小片段都是由超声波震断的，而相同的mRNA分子在同一地方被打断的可能性几乎为零。另一方面，当测序深度过深时，不可避免的，同一模板会被多次测序。这时候更应该去除duplicate，可以消除饱和。对于一些由酶切产生的片段，比如clip-seq, REDseq (Restriction Enzyme digestion sequence)等，就不需要做去除duplicate。在做去除duplicate之前，首先要在genome browser中观察一下mapped好的序列，看看其duplicate的存在的程度。肉眼观察这种事情，因为没有一定的标尺，所以非常不好总结。做这件事情的唯一好处就是，看得多了，就明白什么是好的测序结果。

那么duplicate removal什么时候做呢？现在的观点一般都认为是在map之后做。这样的好处是不依据序列一致就去除它，因为同一段序列可能map至不同的位点。在 map之后，使用samtools rmdup或者Picard MarkDuplicates。这里需要注意的是，无论是samtools还是Picard，在duplicate removal时，所有的mapping结果对于每个read，应该只保留一个位置，而不是多个位置。2，对于PE的结果，LR的名字应该一致，否则程序可能无法识别。这些工具的出发点都是PE的，如果是SE的测序，可能需要指定参数。

java -Xmx2g -jar /path/to/MarkDuplicates.jar INPUT=accepted_hits_sorted.bam OUTPUT=duplicated.removed.bam METRICS_FILE=picard_info.txt REMOVE_DUPLICATES=true ASSUME_SORTED=true VALIDATION_STRINGENCY=LENIENT

或者

 samtools rmdup

第二个问题，bowtie的index文件哪里下载？在bowtie2的主站上提供了一个很有用的链接：iGenome, 这里集中了目前大部分的index文件以及相关的注释文件，可以很方便的下载。本教程就是使用从这里获取的打包文件。

tophat

tophat是针对mRNA-seq对bowtie的map结果进行了优化，它表现在两个方面，第一，从基因组中提取出mRNA junction的可能组合，对没有map结果的reads进行二次比对。第二，对于pe测序的结果进行依照mRNA测序的特点进行调整。它在参数设置方面，非常简单，只需要搞清楚几个重要参数即可。

Average Mate-Pair Inner Distance： -r/–mate-inner-dist

这个参数就是设置mate的两个测序reads之间的平均距离应该是多少。通常PE测序制库时通常片段大小的平均值为300bp左右，这个 300bp包括了两端的adapter, barcodes, 以及序列本身，而-r参数需要设置的是两个测序结果之间的距离，所以它应该是总长-2*(adpter+barcode+reads)。如下图所示：

我们需要注意的是，很多adapter可能比我们想象的要长，所以需要搞清楚具体的实验时的adapter长度。

Gene Model Annotations： -G/–GTF

通常这个参数我们并不设置。设置了之后可能会提高效率，但是也可能会产生倾向性。对于参照基因模型，通常的做法是在后面的步骤中再传入。

Threads：-p/–num-threads

线程数。现在多核处理器非常普遍，所以如果有四核的话，我们不妨设置-p 3，使用其中三个核，留下一核用于其它任务。但有一点要非常注意，如果你使用MPI的话，这个程序并不是MPI书写的程序，所以在cluster上运行时，需要申请独占模式，也就是使用-pe single参数而不是openmpi参数，并且使用-l mem_free=16G申请16G以上的内存以独占计算机，而不是多个内核。否则，多线程不安全。

Output：-o/–output-dir

输出目录。

–library-type

库类型，有三种，fr-unstranded Standard Illumina， fr-firststrand dUTP, NSR, NNSR， 以及fr-secondstrand Ligation, Standard SOLiD。默认为fr-unstranded。

下面是一个运行tophat2的例子：

export PATH=/share/bin/samtools:/share/bin/bowtie2/:/share/bin/tophat2/:/share/bin/python:$PATH python /share/bin/tophat2/tophat --library-type fr-unstranded --mate-inner-dist 70 -p 8 \ --output-dir tophat_output \ /path/to/Genomes/UCSC/mm10/Sequence/Bowtie2Index/genome \ R1_reads.fastq R2_reads.fastq

运行完tophat之后，我们需要检查一下map的质量，这时可以使用fastqc或者samtools来查看。但其实在run log中就已经记录了所有这些信息，可以从日志文件中查找我们需要的质量信息。

此进可以使用前文所述的脚本进行duplicates removal动作。

Cufflinks

对于Cufflinks，无论我们是否需要发现新的基因，最好还是走完三步。否则的话，结果可能会有一定倾向性。这三步是cufflinks, cuffmerge以及cuffdiff

对于每一个bam文件，我们都需要运行一遍cufflinks。

/share/bin/cufflinks/cufflinks -p 8 -o cufflinks_outputs accepted_hits.bam

然后我们将所有生成的transcripts.gtf文件的路径写入一个assemblies.txt文件中。我们知道，cufflinks的目的是通过重组mRNA来确定所有的transcripts，所以它会生成gtf文件来标记哪些是exon,这些exon组成了哪些transcripts。而下一步运行cuffmerge就是为了将所有发现的transcripts与已知的基因模型进行合并。
这里是一个典型的assemblies.txt的范例：

~/scratch/RNAseq/cufflinks_sample1/transcripts.gtf ~/scratch/RNAseq/cufflinks_sample2/transcripts.gtf

运行完cufflinks之后运行cuffmerge， 这一步需要把注释文件传入。

/share/bin/cufflinks/cuffmerge -g /path/to/Genomes/UCSC/mm10/Annotation/Genes/genes.gtf \ -s /path/to/Genomes/UCSC/mm10/Sequence/WholeGenomeFasta/genome.fa \ -o merged_asm -p 8 \ /path/to/assemblies.txt

之后就是运行cuffdiff了。这一步将需要比较的两组的bam文件都传来。需要介绍几个参数，其中，–min-reps-for-js- test是指每组中至少多少个样品，默认为3，可以依据具体的实验调整。–labels之后跟着的是两个组标，用于标记传入cuffdiff计算的两个组的名字，注意不要把顺序标错。而之后的gtf文件就是第二步cuffmerge之后的merged文件。

/share/bin/cufflinks/cuffdiff -o cuffdiff_outputs -p 8 \ --labels group1label,group2label --min-reps-for-js-test 3 \ merged_asm/merged.gtf \ gp1rep1/accepted_hits.bam,gp1rep2/accepted_hits.bam,gp1rep3/accepted_hits.bam \ gp2rep1/accepted_hits.bam,gp2rep2/accepted_hits.bam,gp2rep3/accepted_hits.bam

结过了以上的步骤，我们就可以得到诸如gene_exp.diff, isoform之类的结果。具体的，可以在输出目录中去一一查看。

之后就可以使用R来分析过滤结果了。这里就不做介绍。