RNA-seq差异表达分析工作流程 (转帖)

之前写过博文《从RNA-seq结果到差异表达》给出了一些关于RNA-seq分析的描述,这篇博文的目的是给出一个示例性质的工作流程。

需要使用到的工具: 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

或者

[…]