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

之前写过博文《从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

或者

 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)。如下图所示:
mateDist
我们需要注意的是,很多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来分析过滤结果了。这里就不做介绍。

333 thoughts on “RNA-seq差异表达分析工作流程

  1. Reply 汤冬 7月 14,2017 1:58 上午

    博主你好 我用cummeRbund 读取cuffdiff结果是 录入信息全为0 是什么原因,ps diff_out有结果?
    > library(“cummeRbund”)
    > cuff cuff
    CuffSet instance with:
    0 samples
    0 genes
    0 isoforms
    0 TSS
    0 CDS
    0 promoters
    0 splicing
    0 relCDS

  2. Reply 汤冬 7月 21,2017 3:28 上午

    > cuff cuff
    重新跑了一遍cuffdiff,再次导入cummeRbund 出现上面错误
    附上genes.fpkm_tracking文件前两行
    tracking_id class_code nearest_ref_id gene_id gene_short_name tss_id locus length coverage G-39_FPKM G-39_conf_lo G-39_conf_hi G-39_status J-41_FPKM J-41_conf_lo J-41_conf_hi J-41_status P-40_FPKM P-40_conf_lo P-40_conf_hi P-40_status W-42_FPKM W-42_conf_lo W-42_conf_hi W-42_status X-38_FPKM X-38_conf_lo X-38_conf_hi X-38_status Xue-43_FPKM Xue-43_conf_lo Xue-43_conf_hi Xue-43_status
    XLOC_000001 – – XLOC_000001 ENSGALG00000045540 TSS1 1:29-2665 – – 0 0 0 OK 1.23274 0 4.24099 OK 0.272951 0 1.07974 OK 0.894853 0 3.14556 OK 0.430304 0 1.66738 OK 0.508284 0 1.91057 OK
    六个处理无生物学重复。求指导,万分感谢!!!

  3. Reply 汤冬 7月 21,2017 3:29 上午

    > cuff <- readCufflinks()
    Creating database C:/Users/Administrator/Desktop/diff_out/cuffData.db
    Reading Run Info File C:/Users/Administrator/Desktop/diff_out/run.info
    Writing runInfo Table
    Reading Read Group Info C:/Users/Administrator/Desktop/diff_out/read_groups.info
    Writing replicates Table
    Reading Var Model Info C:/Users/Administrator/Desktop/diff_out/var_model.info
    Writing varModel Table
    Reading C:/Users/Administrator/Desktop/diff_out/genes.fpkm_tracking
    Checking samples table…
    Populating samples table…
    Error: Column name mismatch.
    In addition: There were 50 or more warnings (use warnings() to see the first 50) 错误提示

    • Reply admin 7月 21,2017 8:43 上午

      升级你的RSQLite至最新版,或者从source重新安装,然后再重装一下cummeRbund试一下。如果不行,我就无能为力了。

Leave a Reply

  

  

  

%d 博主赞过: