生信分析流程转录组研究转录组分析

第2周：RNA-seq数据分析最佳实践

2019-01-15  本文已影响599人  TOP生物信息

原文来源：Conesa A, Madrigal P, Tarazona S, et al. A survey of best practices for RNA-seq data analysis[J]. Genome biology, 2016, 17(1): 13.

1. 摘要

  RNA测序（RNA-seq）具有广泛的应用，但没有统一的分析流程能适用于所有情况。我们回顾了RNA-seq数据分析的所有主要步骤，包括实验设计，质量控制，序列比对，基因和转录水平的定量，可视化，差异基因表达，可变性剪接，功能注释，基因融合检测和eQTL定位。

2. 背景

图1 RNA-seq信息分析的通用路线图。主要分析步骤列在预分析，核心分析和高级分析的上方。a预处理包括实验设计，测序设计和质量控制步骤。 b核心分析包括转录组分析，差异基因表达和功能分析。 c高级分析包括可视化，其他RNA-seq技术和数据整合。缩写：ChIP-seq染色质免疫沉淀测序，eQTL表达数量性状基因座，FPKM等于Fragments / (外显子长度*Mapped Reads)，GSEA基因富集分析，PCA主成分分析，RPKM等于total exon reads/ (mapped reads (Millions) * exon length(KB))，sQTL剪接数量性状基因座，TF转录因子，TPM每百万转录本。

3. 实验设计

较便宜的短SE reads通常足以用于研究注释良好的物种中的基因表达水平，而较长和PE reads优先用于研究缺乏注释的转录组。

饱和曲线（Saturation curves）：评估在给定测序深度下预期的转录组覆盖度

对测序实验进行适当规划以避免技术偏差与良好的实验设计同样重要，特别是当实验涉及需要分批处理的大量样品时。

4. RNA-seq数据分析

4.1 质量控制
4.1.1 Raw reads

涉及测序质量，GC含量，adaptor，不合适的k-mers和PCR重复的分析，进而检测是否存在测序错误，PCR人为误差或污染。
工具：FastQC、NGSQC、FASTX-Toolkit、Trimmomatic

4.1.2 Read alignment

比对率：the percentage of mapped reads，整体测序准确性和污染DNA存在的全局指标
原文此处举了一个例子，需要理解一下：

RNA-seq reads比对到人类基因组上，比对率期望是70％到90％，同时由于序列的相似性，一条reads可以同时比对到多个位置，因此存在多比对reads。
但是当reads比对到转录组时，比对率会小一些，因为落在未注释的转录本区域（存在但是暂时还没有发现的那些）的reads会被我们忽略。同时，由于存在一个基因有不止一个转录本，且不同的转录本共享外显子区域（常见）的情况，故一个reads可以比对到多个转录本上（常见），所以多比对reads明显增多。

覆盖度的均一性：不同的外显子上，正负链上

poly（A）针对性选择富集mRNA之后，如果reads主要集中在转录本的3'末端，可能表明起始材料中的RNA质量低。

比对上的reads的GC含量可以用来揭示PCR偏向性。我是这样理解的，下图就是一个例子，不同的GC含量对应的测序深度不一样，很可能就是由于GC含量的差异引起了PCR的偏向。

工具：Picard、RSeQC、Qualimap

4.1.3 Quantification

一旦计算出实际的转录本定量值，就应检查它们的GC含量和基因长度偏差，以便在必要时可以应用校正标准化方法。
如果参考转录组被很好地注释，还可以分析样品的RNA组成，以此来评估RNA纯化的质量。

工具：R包（NOISeq、EDASeq等）

4.1.4 Reproducibility（可再现性）

通过检查重复之间的可再现性和可能的批次效应来评估RNA-seq数据集的全局质量也是至关重要的。
如果基因表达差异存在于不同的实验条件之间，则应该可以预知相同条件的生物学重复将在主成分分析（PCA）中聚集在一起。

4.2 转录本鉴定

  当参考基因组可用时，RNA-seq分析通常是将reads比对到参考基因组或转录组以推断表达了哪些转录本。仅仅比对到已知物种的参考转录组不能发现新的、未注释的转录本，而是只关注量化分析。相反，如果生物体没有测序的基因组，则分析思路是首先将reads组装成较长的contig，然后将这些contig作为表达的转录组进行处理，将reads再次比对回去进行定量。在任何一种情况下，reads覆盖度可用于量化转录本表达水平（图1b）。

图2序列比对和转录本鉴定策略。常规RNA-seq分析的三种基本策略。 a可以使用带注释的基因组，并使用允许gap的比对工具将reads比对到基因组。下一个（新的）转录本发现和量化可以在有或没有注释文件的情况下进行。然后在功能上注释新的转录本。 b如果不需要新的转录本发现，可以使用无gap比对工具将reads比对到参考转录组。转录本鉴定和定量可以同时进行。 c当没有可用的基因组时，需要首先将reads组装成contig或转录组。为了定量，将reads比对回新的参考转录组，并且如（b）中进行进一步分析，然后如（a）中那样对新转录本进行功能注释。可在每个分析步骤中使用的代表性软件以粗体文本表示。缩写：GFF General Feature Format, GTF gene transfer format, RSEM RNA-Seq by Expectation Maximization

方框3.比对到参考序列
比对到参考基因组可以鉴定新基因或转录本，并且需要使用允许gap或剪接的比对工具，因为reads可以跨越剪接点。挑战在于正确识别剪接点，特别是当测序错误或与参考序列的差异（SNP, InDel等）存在时或者在寻找非典型连接和融合转录本（junctions and fusion transcripts）时。最受欢迎的RNA-seq比对工具之一TopHat遵循两步策略，其中未剪接的reads首先被比对以定位外显子，然后未比对的reads被分开（允许gap）并单独地比对以识别外显子剪接点。其他几个比对工具适用于识别SNP或indel（GSNAP，PALMapper，MapSplice），检测非典型剪接点（STAR，MapSplice），实现超快速比对（GEM）或比对长reads（STAR）。在比对过程中要考虑的重要参数是RNA-seq文库的链特异性，允许的错配数，reads的长度和类型（SE或PE）以及测序片段（fragments）的长度。此外，可以通过向一些reads比对工具提供注释文件来利用现有的基因模型，以便准确地定位外显子坐标并帮助识别可变剪接事件。基因模型的选择也会对量化和差异表达分析产生强烈影响。如果转录组注释是全面的（例如，在小鼠或人类中），研究人员可以选择直接比对到所有感兴趣基因的所有转录本序列的Fasta格式文件。在这种情况下，不需要gap比对，可以使用非剪接的比对工具，例如Bowtie（图2b）。比对到转录组通常更快但不能从头发现新的转录本。

4.2.1 Alignment

  当有参考序列可用时，有两种选择：比对到基因组或比对到带注释的转录组（图2a，b;方框3）。无论是使用基因组还是转录组作为参考，reads可以唯一地比对（它们可以仅被分配到参考序列中的一个位置）或者可以是多比对（multireads）。出现基因组multireads主要是由于旁系同源基因的重复序列或共有结构域。当比对到基因组上时，多比对reads通常占比对结果的较大一部分，不应该被丢弃。当参考序列是转录组时，多比对更频繁地产生，因为在基因组上唯一比对的reads将同样很好地比对到共享外显子的所有基因转录本类型上。基因组或转录组比对，转录本鉴定和定量成为研究基因不同转录本表达的重要挑战。

4.2.2 Transcript discovery

  使用Illumina测序技术提供的短reads鉴定新的转录本是RNA-seq中最具挑战性的任务之一。短reads很少跨越几个剪接点，因此难以直接推断所有全长转录本。此外，很难识别转录起始和终止位点。在任何情况下，PE reads和更高的覆盖度有助于重建低表达的转录本，并且重复（replicates）对于在信号检测的下游解决假阳性（即，比对错误或污染）问题是必不可少的。

4.2.3 De novo transcript reconstruction

  当参考基因组不可用或不完整时，RNA-seq reads可以使用诸如SOAPdenovo-Trans，Oases，Trans-ABySS或Trinity等软件重新组装成转录组（图2c）。一般而言，PE链特异性测序和长reads是优选的，因为它们提供更多信息。尽管组装低表达的转录本不可能，因为缺乏足够的覆盖度以进行可靠的组装，但是过多的reads也存在问题，因为它们会导致潜在的错误组装和增加的运行时间。对于样品间的比较分析，建议将来自多个样品的所有reads组合成单个输入，以获得巩固的contig（转录本），然后将短reads比对回去以进行表达量估计。

  无论是有参还是无参，使用Illumina测序得到的短reads完全重建转录组仍然是一个具有挑战性的问题，并且在许多情况下，从头组装会导致数十或数百个contig（只是完整转录本的一个片段）。新兴的长reads技术，如PacBio的SMRT，提供的reads足以对大多数基因的完整转录本进行测序，是一种很有前景的替代方法。（16年的综述文章，2年（多）以后再来看，发现实际确实也是这样啊！）

4.3 转录本定量

  RNA-seq的最常见应用是估计基因和转录本表达量。该应用主要基于比对到每个转录本序列的reads数量，尽管有些算法如Sailfish依赖于reads中的k-mer计数而无需比对。最简单的量化方法是使用HTSeq-count或featureCounts等工具统计比对reads的原始计数。这种基因水平（而不是转录本水平）定量方法利用含有外显子和基因物理坐标的GTF文件，并且经常丢弃multireads。单独的原始reads计数不足以比较样品中的表达水平，因为这些值受诸如转录本长度，reads总数和测序偏差等因素的影响。RPKM（ total exon reads / (mapped reads (Millions) * exon length(KB)) ）是一种样本内标准化方法，将消除特征（feature）长度和文库大小（测序量）效应。RPKM及其随后的衍生指标FPKM（ fragments / (外显子长度*Mapped Reads)，类似于RPKM的表征转录本表达量的样品内标准化指标）和TPM（transcripts per million）是最常报道的RNA-seq基因表达定量指标。应该注意的是，RPKM和FPKM对于SE reads是等效的，并且FPKM可以使用简单的公式转换为TPM。在不同来源和组成的样品之间，TPM更具可比性。

  几种复杂的算法能够解决相关转录本共享大部分reads的问题，来估计转录物水平表达。Cufflinks使用TopHat（现在的HISAT更好）这类工具的比对结果来估计转录本表达量。Cufflinks可以利用PE reads，并且可以使用GTF信息来识别表达的转录本，或者可以仅从比对数据来推断转录本。量化转录组比对（有别于基因组比对）的表达量的算法包括RSEM（RNA-Seq by Expectation Maximization），eXpress，Sailfish和kallisto等。这些方法在转录本中分配多比对reads（我的理解是，如果一个reads比对到多个转录本，那么每个转录本都会计数），并输出针对测序偏差校正的样本内标准化值。此外，RSEM算法使用期望最大化方法返回TPM值。NURD提供了一种用低内存和计算成本估算SE reads的转录本表达量的有效方法。

4.4 差异基因表达分析

  差异表达分析（图1b）要求在样本之间比较基因表达值。RPKM，FPKM和TPM标准化了样品比较的最重要因素，即测序深度。这些标准化方法基于总计数或有效计数，并且当样本具有异质转录本分布时，即当高度差异表达的特征（features）可能使计数分布偏斜时，通常表现不佳。考虑到这一点的标准化方法是TMM，DESeq，PoissonSeq和UpperQuartile，它们忽略高度可变和/或高度表达的特征（features）。影响样品比较的其他因素包括转录本长度的变化和实验条件的变化，转录本序列上覆盖度的位置偏差（在Cufflinks中可校正），平均片段大小和基因的GC含量（在EDAseq包中校正）。 NOISeq R软件包包含各种诊断图，用于识别RNA-seq数据中的偏差来源，并在每种情况下应用适当的标准化程序。最后，尽管存在这些针对特定样本的标准化方法，但数据中仍可能存在批次效应。这些影响可以通过适当的实验设计最小化，或者通过批量校正方法如COMBAT或ARSyN去除。这些方法虽然最初是为微阵列（microarray）数据开发的，但已被证明可以很好地与标准化的RNA-seq数据一起使用。

原文接着介绍了几种概率分布和基于一些统计假设的工具

在没有生物重复的情况下，不可以进行总体推断，任何p值计算都是无效的。因此，在没有重复的情况下分析RNA-seq数据的方法仅具有探索价值。我们建议在允许的情况下，RNA-seq实验至少具备三个生物学重复。

方法的选择（甚至软件包的版本）可以显着影响分析的结果，我们建议详细记录所用程序的设置/参数和版本号，并考虑使用多个程序重复进行重要分析。

方框4.用于检测差异基因/转录表达的软件工具的比较
首先，当重复样本的数量非常小时，或者对于以非常低水平表达的基因，所有方法都需要特别谨慎。在这些工具中，limma在许多情况下表现良好，并且运行速度最快。DESeq和edgeR在排序基因方面表现相似，但在控制FDR（false discovery rate，伪发现率）时分别相对保守或过于自由。SAMseq在FDR方面表现良好，但在重复次数相对较高（至少为10）时，具有可接受的灵敏度（言外之意是说重复次数要高一些？）。 NOISeq和NOISeqBIO（NOISeq对生物重复的适应性）在以某种灵敏度为代价避免检测假阳性方面更有效，但在不同重复次数的情况下表现良好。 Cuffdiff和Cuffdiff2在比较中表现较差。这可能反映了这样的事实：检测转录水平的差异表达仍然具有挑战性，并且涉及将reads分配给可变转录本类型的不确定性。在最近的一项比较中，BitSeq与其他转录水平的软件包（例如Cuffdiff2）相比毫不逊色。最后，选择分析方法时的一个重要考虑因素是实验设计。虽然一些差异表达工具只能执行成对比较，但其他如edgeR，limma-voom，DESeq，DESeq2和maSigPro可以执行多重比较，接受不同的协变量或分析时间序列数据。

4.5 可变剪接分析

  转录水平差异表达分析可以检测来自相同基因的不同转录本类型的表达变化，并且已经提出了使用RNA-seq进行可变剪接分析的特定算法。这些方法分为两大类。第一种方法将转录本表达估计与差异表达的检测相结合，以揭示总基因表达中每种转录本类型的比例的变化。BASIS，一种早期的方法，使用分层贝叶斯模型直接推断差异表达的转录本类型。CuffDiff2首先估计转录本表达量，然后比较它们的差异。通过整合这两个步骤，在进行统计分析以寻找差异表达的转录本类型时，将第一步的不确定性考虑在内。FDM（The flow difference
metric）使用来自比对上的外显子reads和连接处reads的比对累积转录本图来推断转录本类型，并使用JS散度来测量差异。最近，提出了一种新方法rSeqDiff，该方法使用分层似然比检验来检测差异基因表达，同时不进行剪接变化和差异转录本表达。如RNA-seq基因组注释评估项目（RNA-seq Genome Annotation Assessment Project）论文中所讨论的，所有这些方法通常受到短reads序列在转录本水平上准确鉴定的内在限制的阻碍。

  所谓的“基于外显子的”方法跳过不同转录本表达估计这一步，并通过比较样品之间基因的外显子和剪接点的reads分布来检测可变剪接的信号。该方法基于以下前提：可以在外显子及其连接处的信号中跟踪不同转录本类型表达的差异。DEXseq和DSGSeq采用了类似的想法，通过检测基因外显子（和剪接点）reads计数的显着差异来检测差异剪接基因。 rMATS通过比较连接处reads定义的外显子水平来检测外显子的差异使用。rDiff通过比较基因可变区域的reads计数来检测差异转录本表达，有或没有注释的可变转录本皆可。DiffSplice使用比对图来鉴定可变的剪接模块（ASM），并使用ASM的信号识别不同的剪接。这种外显子或连接方法的优点是它们在识别单个可变剪接事件时具有更高的准确性。如果研究的重点不在于所有转录本类型，而在于包含和排除特定外显子和这些外显子所包含的功能蛋白结构域（或非翻译区域外显子的调节性的特征区域（features）），则基于外显子的方法是合适的。

4.6 可视化
IGV

它提供了灵活的方式来显示外显子，转录本和连接点的reads丰度。可以隐藏内含子以更好地显示外显子上的信号，并且热图可以帮助多个样本信息的视觉比较。但是，RNAseqViewer比IGV慢。

RNAseqViewer

  用于差异基因表达分析的一些软件包（例如Bioconductor中的DESeq2或DEXseq）具有使结果可视化的功能，而其他专门用于可视化的工具的开发，例如CummeRbund（for CuffDiff）或Sashimi-plot，可用于可视化外显子的差异剪接。Sashimi-plot的优点在于，当样本数量较少时，它们的junction reads显示更直观，美观。可以使用SplicePlot获得用于剪接数量性状基因座（sQTL）的Sashimi，structure和hive图。可以使用SpliceSeq生成剪接图，SplicingViewer绘制剪接点和可变剪接事件。 TraV是一种整合数据分析的可视化工具，但其分析方法不适用于大型基因组。

Sashimi-plot
4.7 基因融合发现

融合基因（Fusion gene），是指将两个或多个基因的编码区首尾相连．置于同一套调控序列（包括启动子、增强子、核糖体结合序列、终止子等）控制之下，构成的嵌合基因。——百度词条

同源基因家族和高度多态性基因（例如HLA基因）产生的reads不能容易地唯一比对到它们在参考基因组中的最初位置。对于非常高表达的基因，RNA-seq的小但不可忽略的测序错误率将产生错误比对到同源基因座的reads。建议过滤高度多态性基因和同源基因对。还建议过滤不太可能参与基因融合的高表达基因，如核糖体RNA。

4.8 小RNA

  sRNA的长度通常为18-34个核苷酸，它们包括miRNA，小干扰RNA（siRNA），PIWI互作RNA（piRNA）和其他类型的调节分子。由于缺乏复杂性，sRNA-seq文库很少以常规RNA-seq文库那样的深度进行测序，通常测序量为2-10百万reads。 sRNA-seq数据的生物信息学分析不同于标准RNA-seq方案（图1c）。首先去除adaptor序列，并计算得到的reads长度分布。在动物中，通常存在22和23个核苷酸的峰，而在植物中存在21和24个核苷酸冗余reads的峰。例如，用于预测和分析sRNA种类的工具miRTools 2.0默认使用长度为18-30个碱基的reads。阈值取决于实际应用，并且在研究miRNA的情况下通常在19-25个核苷酸的范围内。

  与标准RNA-seq一样，sRNA reads必须使用标准工具与参考基因组或转录组序列比对，如Bowtie2，STAR或Burrows-Wheeler Aligner（BWA）。此外，有一些比对工具（例如PatMaN和MicroRazerS）被设计用于比对具有预设参数值范围的短序列，该预设参数值范围适合于短reads的最佳比对。比对本身可以在有或没有错配的情况下执行，后者更常用。另外，可以去除超出预先设定位置个数值的reads。在研究miRNA时，通常允许在每个基因组上有5-20个不同的比对。然后简单地计数sRNA reads以获得表达值。然而，研究者还应该验证他们的sRNA reads没有被降解的mRNA显着污染，例如，通过检查miRNA文库是否显示出高表达基因（如GAPDH或ACTB）的基因组区域有较高覆盖度。

  进一步的分析步骤包括与已知的sRNA的比较和sRNA的从头鉴定。有用于此目的的特定工具，例如分别用于动物和植物miRNA的miRDeep和miRDeep-P，或用于反式作用siRNA（phasi RNA？）预测的UEA sRNA Workbench工具。miRTools 2.0，ShortStack和iMir等工具也可用于全面注释sRNA集和鉴定不同类别的sRNA。

总结一下：

4.9 使用RNA-seq结果进行功能分析

  标准转录组学研究的最后一步（图1b）通常是分子功能或通路的描述，其中涉及差异表达的基因（DEGs）。起初为微阵列技术开发的功能描述的两种主要方法是（a）比较DEGs与剩下的基因组基因，探讨DEGs中哪些功能更突出，以及（b）基因富集分析（GSEA），它的基础是根据差异表达的测量值给转录组排序。 RNA-seq偏向性如基因长度使这些方法直接应用于计数数据变得复杂，因此提出了RNA-seq特异性工具。例如，GOseq估计差异表达结果的偏差效应（如基因长度），并调整功能富集测试中使用的传统超几何统计量来解释这种偏差。相似的，Gene Set Variation Analysis (GSVA)或SeqGSEA包也结合了类似于GSEA的剪接和实施富集分析。

  功能分析需要为研究中的转录组提供足够的功能注释数据。诸如Gene Ontology，Bioconductor，DAVID或Babelomics等资源包含大多数模式物种的注释数据。然而，在从头组装或重建转录组期间发现的新转录本将缺少一些功能信息，因此为了更好地对那些结果进行功能分析，注释是必需的。通过在蛋白质数据库（如SwissProt）和含有保守蛋白结构域的数据库（如Pfam和InterPro）中搜索相似序列，可以利用直系同源序列对编码蛋白质的转录本进行功能注释。诸如Gene Ontology（GO）之类的标准术语词汇表的使用允许直系同源序列的功能信息在一定程度上具有可交换性。诸如Blast2GO等热门工具允许针对各种数据库和特定的专业术语表对完整的转录组数据集进行大量注释。通常，从RNA-seq数据重建的50-80％的转录本可以用这种方式用功能性术语注释。然而，RNA-seq数据还显示转录组的一些重要部分缺乏蛋白质编码能力。这些长非编码RNA的功能注释更具挑战性，因为它们的保守性通常不及蛋白质编码基因。Rfam数据库包含注释最全面的RNA家族，如核糖体或转运RNA，而mirBase或Miranda专门研究miRNA。这些资源可用于短非编码RNA的基于相似性的注释，但是对于其他RNA类型（例如长非编码RNA），还没有标准的功能注释程序可用。

5. 与其他数据类型整合

5.1 DNA测序
5.2 DNA甲基化
5.3 染色质特征
5.4 MicroRNAs
5.5 蛋白质组学和代谢组学
5.6 多种数据的整合和可视化

6. 展望

6.1 单细胞转录组测序
6.2 全长转录组

最近在学习转录组和小RNA，所以就学习了这一篇综述。初次接触转录组，有些术语不太理解，可能翻译有错误。如果您发现了错误，请帮我指出来，谢谢！！！

上一篇下一篇

猜你喜欢

热点阅读