之前写过博文《从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
或者
[…]
Recent Comments