从RNA-seq结果到差异表达

翻译自:From RNA-seq reads to differential expression, Oshlack et al. Genome Biology 2010, 11:220

高通量测序技术,也就是下一代测序技术已经成为现代生物学研究的一个较为常规的实验手段了。这一技术的发展极大地推动了基因组学,表观基因组学以及翻译组学的研究。RNA-seq通过测定稳定状态下的RNA样品的序列来对RNA样品进行研究,从而避免了许多之前研究手段的不足,比如象基因芯片或者PCR就需要背景知识。而且RNA-seq还可以触及以前无法研究的领域,比如复杂结构的转录体。RNA-seq可以应用于以下几个方面的研究,1. SNPs;2. novel transcripts;3. alternative splicing;4. RNA editing。无论如何,使用RNA-seq最多的还是比较两组样品基因水平表达差异,比如野生型与突变型,用药组与对照组,不同组织之间,癌细胞与正常细胞,等等。我们把这种基因水平差异表达,简称为DE (differential expression,注,不是ED啊〜〜〜)。

常用的RNA-seq操作平台有Illumina GA/ HiSeq, SOLiD 还有Roche 454。它们都是提取RNA后,纯化,打碎,逆转录成cDNA,然后测序。测序的结果被称为short reads,短序。通常一个短序的长度为25-300bp之间。如果测序只测一端可能会带来比对时的困难,于是这些操作平台提供了两端都测的办法,这样的结果成对出现,中间有一定的间隔,但是因为测序长度一下子提高了一倍,所以比对会精准很多。人们把这种测序结果称为’paired-end’ reads,成对短序。一般来讲,测序结果会直接转换成一行一行的由字母组成的短序列,可能是fasta,fastq等等不同格式。

然而,这一技术产生的海量数据分析却给生物学家带来了难题。一个测序的结果文件少则几Gb,多则几十Gb,单独对比拼接,就会用去几个小时,而后再得出差异表达的结果,其耗时耗力,并非实验生物学家可以应付得了的。于是生物信息学的研究人员努力做出一些软件,以降低结果分析的难度。但是,即使这样,还是必须对分析过程有个较为细致地了解,才能正确地使用这些软件,从而得到比较接近事实的结果。

一般的来讲,RNA-seq后DE的工作流程是这样的(图1),首先,将短序映射到基因组相应的位置上去,其次,对映射的结果进行基因水平,外显子水平,以及转录水平的拼接,而后对结果进行数据统计,标准化之后生成表达水平报告文件,最后由生物学者依据系统生物学相关知识,来对数据结果进行分析。

RNA-seq分析工作流程
RNA-seq分析工作流程

不同步骤涉汲的软件和方法:

分析步骤 方法
mapping General aligner GMAP/GSNAP
BFAST
BOWTIE
CloudBurst
GNUmap
MAQ/BWA
PerM
RzaerS
Mrfast/mrsfast
SOAP/SOAP2
SHRiMP
De novo annotator QPALMA/GenomeMapper/PALMapper
SpliceMap
SOAPals
G-Mo.R-Se
TopHat
SplitSeek
De novo transcript assembler Qases
MIRA
Summarization Isoform-based Cufflinks
ALEXA-seq
Gene-based Count exons only
Exon junction libraries
Normalization library size
RPKM: reads per kilobase of exon model per million mapped reads ERANGE
TMM: trimmed mean of M-values edgeR
Upper quartile Myrna
Differential expression Poisson GLM (generalized linear model) DEGseq
Myrna
Negative binomial edgeR
DESeq
baySeq
Systems biology Gene Ontology analysis GOseq

映射至基因组(Mapping)

第一步的工作是比对(alignment)。对于RNA-seq的比对,从来都不是一件容易的事情。其难点如下:

  1. 没有很好的比对模板。现在的比对模板都是基因组模板,而不是真正的转录组模板,也就是说,这对本来就不是很长的短序来说,它很有可能是界于两个exon之间。我们在比对junction的时候,一般还是假设它如果没能在基因组模板中找到合适的位置的时候,才考虑它是否是界于junction上。这种人为的假设可能并不准确。
  2. SNPs,碱基插入,删除,错配,或者质量不高的测序结果,从模板至比对序列本身,都存在着比基因比对更为复杂的问题。
  3. 短序可能会有多个100%的匹配位点。
  4. 有些基因组可能需要庞大的内存空间。

为了解决最后一个问题,人们使用了很多办法,但基本上都会基于事先建立的引索库。即所谓“启发式”比对(heuristic match)。首先使用一定长度的(通常是11个碱基)的序列做为索引用的关键字,在匹配这一索引字之后,就很大程度地缩小了其需要匹配的模板范围。但是这一办法的问题在于不容易解决问题2中的空格,错配问题。所以在很多软件使用时,会要求人工确认高保真区,以及最高允许2〜3个错配。

现在比较快的“启发式”比对主要有两种算法,一种是哈希表(hash table),一种是BW压缩转换(Burrows Wheeler transform, BWT)。前者速度快,但是对内存要求比后者要高。

对于问题3,一般而言,大部分软件使用的办法是只保留一个匹配位点,其中,有些是只保留第一个匹配位点,有些是按照概率分布选取保留的位点。当然,前面已经提到过,可以使用paired-end read来尽量避免问题3的出现。

对于问题1,可以使用外显子库来确定junction reads。有两种办法,一种是依靠已知的外显子库来构建,另一种办法就是依据已经匹配好的短序来构建外显子库(de novo assembly of transcriptome)。后者的不足是运算量大,对测序覆盖范围要求高,最好是使用paired-end reads。

还有人发现,对于ploy(A)的处理会减少不能映身的短序数。比如,Pickrell et al.就发现,对于46bp的Illumina reads,87%的短序可以映射至模板,7%可以映射至junction library。如果对那些不能映射的短序,将在头或者尾含有的超过连续4个的A或者T去除,就可以得到约0.005%的映射。

综合评价(Summarizing mapped reads)

这一步,主要是基本于不同水平(外显子水平,转录水平,或者基因水平)进行统计。最简单的办法就是统计落在每个外显上的短序数。但是有研究表明,很多(可能超过15%)的短序会落在外显子两侧,这会影响统计的结果。另一种办法就是统会落在内显子区域的短序数。

无论如何,即使是基因水平的综合评价,也还是有其它的一些问题。比如overlapping的基因的统计。比如junction的统计。

标准化(Normalization)

标准化对于样品内及样品间的比较而言是非常重要的。标准化被分为两类,样品内及样品间(between- and within-library)。

样品内标准化使得在同一样品内不得基因之间的表达差异变得有意义。最常用到的一个办法就是使用落在同一基因内的短序数除以单位基因长度。比较常用的单位是RPKM (reads per kilobase of exon model per million mapped reads)。但是这一方法也受到样品制备和测序方法的干扰。

而对于样品间标准化,最简单而直接的办法使用短序总数来平衡表达量。然而短序总数受测序深度的干扰,而且单个基因的短序数与实际的表达量并不一定会呈线性比较关系。人们又使用四分位(quantile normlization)标准化的办法。但是有研究说这一办法并没有实际的价值。还有提出使用对数分布法则(power law distributions)来进行样品间标准化。但没有研究对这一处理方式进行验证。

差异表达(Differential expression)

差异表达分析的最终目的是将那些差异表达的基因(外显子等等)从海量数据中提取出来。最终的结果显示一般来说是表格化的,这一表格按照一定的规则排序,让人们能够尽可能简单地拿到想要的结果。

由于RNA-seq结果的离散性,人们一般都会使用统计模型来拟合实验得到的结果。一般而言,RNA-seq的结果是比较附合伯松分布(poisson distribution)的。这一结果得到了单通道Illumina GA测序结果的实验验证。但是,伯松分布分析结果常常在多组重复的样品间带来较高的假阳性,因为它低估了生物取样的样品间误差。所以RNA-seq如何设置重复是一个很重要的问题。为了平衡重复样品所带来的误差,人们使用了serial analysis of gene expression (SAGE) data。

现有的软件一般都是针对较为简单的实验设计的。而对于复杂的实验设计,比如说成对样品,时间依赖样品等等,还没有专门的,较好的解决方案。大多数都使用edgeR的线性模型来进行分析。

后期系统生物学分析

简单地讲,前景是广阔的,但目前为止手段还是比较有限的,基本上就是GO分析。

48 thoughts on “从RNA-seq结果到差异表达

  1. 你好,请问你对goseq分析是否熟悉。我现在再做两个细菌基因组的RNAseq分析,现在做到goseq分析了,但是由于基因组不是goseq支持的,再准备gene lengthdata和category mapping 文件时出现困难,是否可以请教相关问题。

  2. 你好!请问一下,我现在在利用cufflinks做表达谱的差异表达分析,需要一个transcripts.gtf作为参考,但是我没有这个文件。我的基因组和基因都具备,请问我该如何利用这些数据生成一个transcripts.gtf?

  3. 非常谢谢您的耐心解答!您那篇文章我也早已拜读,收获很大。问题就在这是一个我们自己拼接的基因组。我刚刚接手工作,而他预测基因的原始数据我并没有,导致我手上只有相关的基因序列(含intron和不含intron),并没有成熟的gene.gtf等相关注释。那我该如何制作呢。实在抱歉,我查阅了很多资料,但是由于无人指导,所以真的不清楚!真的谢谢您。

  4. 请教一个问题,我从ECODE里面找到几个TF的CHIP-SEQ数据,我如何查看TF结合的哪些基因呢?需要用什么软件处理吗?谢谢

    1. 谢谢您的留言。
      通常ChIP-seq的数据分析有几步,1,mapping,这里使用到的是bowtie/bowtie2或者bwa。2,peak detection,峰检测一般使用macs。3,annotation,通常使用R/bioconductor中的ChIPpeakAnno。

  5. 您好,我想問一下,我將兩個樣品各別的Reads對回Combined在一起進行組裝完成之序列,並進行分析各樣品在不同Unigenes之表現量,我不懂的是這些 reads 對回unigene 的原理是甚麼??
    比如:————————-這是一個unigene
    read: — —- —- —
    — — —- —- —- 那如何判別他的表現量依據?這些reads全都列入計算
    嗎??還是說會類似接併的方法將重疊的部分不列入計算之
    類的?? 抱歉,真的不懂這部分!!

    1. 感谢您的留言。简体回复请见谅。
      reads计数的过程,通常不以碱基数为单位,而是以个为单位。如果有重叠,无论重叠多少,都不影响计数。计数的说明,基本上现在以HTSEQ的htseq-count说明文档为准。而关于表达量(您的说法是表现量)的计数方法,可以参见我的如何计算cuffdiff中的FPKM值以及cufflinks

      在之后的peaks统计时,可能会小涉及重叠的部分,比如:http://tiling.mbb.yale.edu/frame1.html以及http://info.gersteinlab.org/ACT_IntegratedExample。

      希望这可以解决您的疑问。

  6. 您好,我现在按照一般的流程tophat-cufflinks-cuffmerge得到merged.gtf文件了,我想筛选lncRNA,我现在才筛选了length和exon符合要求的,我想问下,后续您知道怎么再按照FPKM筛选吗?FPKM的阈值要怎么定?有说是画ROC曲线来定?

    1. 在筛选lncRNA时,你需要注意的是你要自己选择合适的包含lncRNA注释的gtf文件。之后,你可以依据FPKM值来时行比较。至于你后面的问题,我不太明白,能再解释一下吗?

  7. 你好,我想问下,用Cuffdiff出来的差异表达结果可以用GOseq或者GOstats做富集吗?大概流程可以讲下不?谢谢

  8. 你好,我有个问题如下:
    比如我有野生型W和突变w’,处理T和不处理t两种状态,也就是会有2×2组样品,Wt,WT,w’t,w’T,四组样品,如果我要对这四组样品进行每一个基因的在这四组里面的差异表达,该如何实现呢?

    DESeq和edegR里面提到的多因子貌似都是两两比较,我看到一个可能比较类似的例子,“mRNA and Small RNA Transcriptomes Reveal Insights into
    Dynamic Homoeolog Regulation of Allopolyploid Heterosis
    in Nascent Hexaploid Wheat ”这篇文章里的Figure 7A,对几组样品进行了差异表达的比较,不知道他们的倍数是如何计算出来的。

    1. DESeq和edegR里的两两比较只是一个简单的例子。你可以参与limma的教程,对design和contrast.matrix稍做修改就可以得到多组比较的结果了。

      1. 感谢回复。
        嗯,原来是我理解有误,我以为glm模型里最后还是基于两两比较,没办法多组比较的。试着用deseq做了下,发现它只要估计了sizeFactors之后就可以把各组标准化到同一水平,那这样多组比较应该是可能的。
        再对比了下,deseq的formula似乎不太好理解?如我之前留言的例子,genotypextreatment 2×2的设计,如果看作四组,在deseq里面用
        fit0=fitNbinomGLMs(countdata,count~1)
        fit1=fitNbinomGLMs(countdata,count~group)
        再做test,无假设是没有任何一组具有差异,
        如果是设计成
        fit0=fitNbinomGLMs(countdata,count~genotype+treatment+genotype)
        fit1=fitNbinomGLMs(countdata,count~genotype+treatment+genotype:treatment)
        我不太理解这个无假设是什么,具体比较的是哪些组,能请问下你是怎么理解的吗?

        1. 那我再贴一次公式(formula)的理解。希望它能帮助到你。你的后两个公式的差别就在于是否与它们的相互作用有关。我个人觉着你使用第一个fit0就好了。

          公式

          模型

          意义

          Y ∼ A Y = β0 + β1 A 一次线性关系, 并且允许在 Y 方向上有截矩
          Y ∼ −1 + A Y = β1 A 过原点的一次线性关系
          Y ∼ A + I (A2 ) Y = β0 + β1 A + β2 A2 多项式模型; 其中 I( ) 函数可以引入普通的数学公式.
          Y ∼ A + B Y = β0 + β1 A + β2 B 二元一次方程式.
          Y ∼ A : B Y = β0 + β1 AB 与 A 及 B 的相互关系成一次线性关系
          Y ∼ A * B Y = β0 + β1 A + β2 B +β3 AB 与A和B相关也与A及 B的相互 关系有关也可以写成 Y ∼ A + B+ A:B.
          Y ∼ (A + B + C )2 Y = β0 + β1 A + β2 B + β3 C + β4 AB + β5 AC + β6 BC 与多个一次变量相关, 同时也与它们的n个元素的组合有关, 这里的 n就是 ( )n 中的n.也可以写成Y ∼ A*B*C– A:B:C.
          1. 感谢博主的提醒以及limma包(?)的formula翻译。

            你的后两个公式的差别就在于是否与它们的相互作用有关。我个人觉着你使用第一个fit0就好了

            从博主的说法来看(fit1,fit0)的test应该是测试哪个模型更适合某个实验吧,当确定了合适的模型,再走标准流程?(我的统计学知识还是需要update)

            其实我之前的问题还想问interaction项的理解来着。在edgeR,如果要比较genotype导致对两种treatment不同差异表达的时候应该表示成(WT-Wt)-(wT-wt)这种形式?在WT-Wt和wT-wt这两组间寻找差异?
            如果我要找treatment对这两种基因型造成的一般性影响,就是说找WT-Wt和wT-wt这两组的共同部分,又该如何实现呢(难道就是DESeq里的fit0)。

            不知道我的理解有没有混乱。

          2. fit本身不计算结果,它只是对数据进行调整。contrast公式是用于最终计算差异结果的。如果要找共同的部分,以我的观点,我会分开计算,而后找共同的部分。

          3. 我看到grokbase上有个关于deseq的帖子

            http://grokbase.com/t/r/bioconductor/135ehek9kd/bioc-deseq-multifactor-design-paired-analysisz

            head(fit1)
            (Intercept) pair3 conditionTrt1 conditionTrt2 deviance converged
            gene0 10.188336 -0.30184662 0.1698471 0.188197344 0.0741656 TRUE
            gene1 9.992372 -0.32221027 0.1415406 0.001119304 0.1568869 TRUE
            gene10 6.893693 0.06795915 -0.4954365 -0.270402107 0.9072335 TRUE
            gene100 6.477847 0.04684995 0.1250743 -0.094116036 2.0611685 TRUE
            gene1002 4.463446 0.15038277 0.3186187 0.106512332 0.3701619 TRUE
            gene1003 4.090079 -0.05247162 -0.3560189 -0.054765642 1.9631768 TRUE
            pvalsGLM = nbinomGLMTest( fit1, fit0 )
            padjGLM = p.adjust(pvalsGLM, method="BH" )

            Yes, you are right. The 3rd and 4th columns of fit1 are the log2 fold changes for trt1 vs control and trt2 vs control.
            The p-values and adjusted p-values here are from a likelihood ratio test of
            a model including the condition and pair variables vs a model which only
            includes the pair variable.

            从作者的回复来看fit1还是会给出同时给出后面各组相对于control的fold change的值的?
            这些值之间应该是可以比较的吧?
            这些fold change值是和两两对比是一致的么

            两两比对的p值是显著性,那pvalsGLM,padjGLM是什么意义呢?

          4. 其它的问题我可能还真回答不上来。只能回答padjGLM是对p.value重新做一个计算,我们通常把BH算法出来的结果当成是FDR(其实不严格是,但大家都这么用)。

          5. 嗯,感谢博主的回复。
            还有

            如果要找共同的部分,以我的观点,我会分开计算,而后找共同的部分。

            这里我问过了,可以一步得到的,可能参考deseq2会更容易理解,作者之一现在推荐用deseq2。
            理解formula略伤身…

  9. 您好。我想问一下,我从公司那里获得了一个excle文档,里面有NM号对应的各组的counts数,软件是NTSEQ 6.0,我想转换成每个基因按照长度和测序量标准化后的变化的p值,请问这个标准化和p值计算应该用什么软件呢?

    1. 如果有count number了,我建议您使用edgeR或者DESeq来处理,参见http://pgfe.umassmed.edu/ou/archives/3777
      如果你一定要做rpkm,你可以使用biomaRt来注释,然后获得长度,然后按照cufflinks的文献说明自己一步一步地做。但是这样是非常困难的,除非你的数学功底很扎实。

  10. 博主你好。目前的问题是:我只有RPKM值,该如何找差异基因?之前你回答说用DESeq或者voom,我不是很理解,这两个软件和edgeR的什么区别造成了这个选择呢?

    1. 我想我的说法是edgeR/DESeq2/voom都可以,但是它们都是针对raw reads count来设计的。如果是RPKM值的话,可以做,但不如cufflink出的结果好。不过你的情况,只有使用edgeR/DESeq2/voom了。

  11. 欧老师您好,我问一个理论性的问题。
    请问tophat在做map的时候,(1)是只按照注释过的也就是genes.gtf去map还是全部基因组都map?理论上能给我解释一下么?
    因为我现在在做牛的。我发现用Ensembl的btau4.0(2007版本)和UMD3.1(2009版本)做的话,同样的参数,Cuffdiff的输出,大部分相同,但仍有不同,比如(2)后者比前者多了几个DE genes,(3)有些前者认为DE genes的但后者不认为是,(4)即使相同的基因(ensembl ID),但是基因名不同了(这个我可以理解),而基因的location也不同了(当然染色体也是相同),请问这个正常么?
    请您得闲时帮忙回答问题(1)-(4),谢谢!

    1. (1). 如果你提供了genes.gtf,tophat就不再去自己拼装transcriptome了,这会对拼装的结果中junction的部分影响巨大。如果你要分析alternative splicing,最好就是提供一下transcriptome。而且,不使用de novo会一定程度地提高map的效率和准确率。对于人,鼠等比较成熟的模式生物,提供transcriptome是比较好的选择。
      (2). 不同的注释版本是会有不同的。
      (3). 和前一问应该是一个问题。
      (4). 正常。

  12. 博主,记得你之前有过cufflinks和Deseq二者区别的简单比较,尤其是在重复样本的情况下。能否在重复下上述程序在估计差异表达方面最重要的区别! 非常感谢

    1. 两者区别挺大的。我只能简单的讲讲两者的应用环境。对于cufflinks,两两比较的进较的时候比较常用。对于edgeR/DESeq/voom,可以进行比较复杂的比较。cufflinks的假阳性比较低,但假阴性会高一些。

      1. 例如,我现在有三个处理组的二代测序数据。mapping之后的,DE等分析,你觉得用哪个程序比较好一点!

  13. 换句换说,我想问的问题或建议是,是基于cufflinks的两两比较好一点(这样cuffdiff就要做三次),还是要用DESeq之类的。从统计的角度看,好像是后者好一点。但是我个人对DESeq不熟,另外就是操作上Deseq用R语言写的,这种内存需要比较大的肯定不能再R自己的界面下运行,所以想请你粘贴个例子如何在Linux shell下如何调用DEseq中的函数的命令行?

    1. 我建议你使用cuffdiff做,cuffdiff可以把三组都一次性传上去,它会两两比较,你可以读一下它的文档。

  14. 您好,我最近做RNAseq,做的是植物材料,对照和虫子处理,做了分析后的能够获得2000个sig gene. 做volcano发现,很多p_adj_value都是一个值, 特别是比较小的值,所以做出了的volcano 比较难看, plot最上面有几行就是平行的,请问这样正常么?我用的是tophat2, cufflinks分析的,125bp, single end. 如果需要其他信息或者我没说明白,烦请联系我,非常感谢。帮我分析下。

  15. 博主,你好,我想请问一下用DEseq做差异组分析时,为什么用A_vs_B 跟B_vs_A 两组颠倒了下位置后筛选出来的差异个数就不同了,比较了下两次的结果,pvalue的值有差异,请问博主对这种问题有什么意见吗?

    1. 这是一个over和under的比较问题。我想问一下,你在比较的时候是从头从走了一遍还是只在最后一步重新比较了一下?你的design是怎么写的,最后是怎么比较的?

  16. Hi,

    你好,有个下游的问题请教下。请问如果拿到了DEGs, 如何画一个比较不错的调控网络。我做的是植物抗病。我现在有一个材料对两个胁迫的DEG dataset, 我想是否能够分别预测或者画出分别得网络,然后找到common pathways (当然我也可以找到common DEGs in both conditions) ,并用network表示出来。我有几个idea不知道是否可行,一个是我要下载所有其他人用这个材料抗病的表达数据,包括我的,构建一个大的这个种的specific抗病网络,然后提出我的部分。 但是好像比较复杂. 另一个是用KEGG的资源来构建,相对容易的多。但这些都还是想法,您是否知道有经验?或者 有啥建议,什么软件包之类的。谢谢!

    1. 一般来说,我们就的做法就是enrichment analysis。如果还是不能满足您的需要,可以考虑使用更多的数据,比如结合ChIP-seq, ATAC-seq, HiC-seq, histone modification epigenetics的数据,使用诸如GeneNetworkBuilder之类的包来建网。

  17. 您好,我做了细菌的RNA-SEQ,但是实验组和对照组只有三十几个差异基因,而且是在放宽筛选阈值的情况下(p<0.05),而且p值相对差异表达倍数较大,例如,差异为3倍时,p值就已经大于0.05,FDR已为1,请问这是这么回事呢

    1. 这就意味着组内差异太大了。两种方案,再做几组重复,2使用适当的normalization方法。您还没有讲你是用的什么差异分析软件。

  18. 博主好, 在用cuffdiff分析差异表达的时候,为什么用的是从tophat得到的accepted_hits.bam文件,而不是cufflinks分析得到的transcripts文件呢

    1. 那个transcripts文件只是整合了新发现的基因,并没有计数。要计数的话,还是需要从bam文件开始啊。

orange进行回复 取消回复

电子邮件地址不会被公开。 必填项已用*标注