RNA sequencing: the teenage years（中篇）

Posted 2023-03-30

参考技术A 前言

这一部分主要涉及常规RNA-seq的数据分析，提到了常见的数据分析流程。

RN A-seq数据分析

用于分析测序读长以确定差异表达的计算方法的数量在过去10年里大量增加，并且即使对于最简单的DGE分析来说，在分析实践中，每个步骤也存在着大量的差异。然而，每个步骤都可以使用不同方法，这些方法的不同组合会对从数据中得到的生物学结论产生重要的影响。这些工具的最佳组织取决于正在研究的特定生物学问题，以及可用的计算机资源。虽然有着尽可能多的排列组合，但是我们的重点在于研究，每个世界大在样本之间的差异表达的可能性的工具和技术。针对这个目标，我们可以将分析过程划分为4个阶段（FIG 2；TABLE 2）。

第1阶段是将一个测序平台产生的原始测序读长导入工具，并将这些读长回贴到转录组上。

第2阶段，是对每个基因或转录本相关的读长数目进行定量（表达矩阵） 。这一过程涉及一个或多个不同的比对(alignment)，组装(assembly)与定量(quantification)亚过程，或者是可以在单个步骤中从读长计数中，整体地生成表达矩阵。

第3阶段是通过过滤低表达特征来改变表达矩阵 ，这一步的 关键步骤是对原始读长计数进行归一化， 用于解释样本之间的技术差异。

第4阶段是样本组之间的统计建模与协变量(covariates)，以及计算与差异表达相关的置信统计量。

Figure2-差异基因表达的RNA-seq数据分析流程

Figure 2-差异基因表达的RNA-seq数据分析流程。差异基因表达(DGE)分析的第一步是原始RNA测序读长的FASTQ格式的数据，DGE的分析有多种方式。主流的分析流程有三种（用实线划的三个方框，分别用A，B和C表示），并且图上还列出了许多替代工具（用虚线表示）。

在A分析流程中，比对工具例如TopHat，STAR或HISAT2使用一个参考基因组来将读长回贴到基因组的位置上，然后使用一些定量工具，例如HTSeq和featureCounts，来将读长比对于基因的特征上。在归一化后（通常归一化的方式都内嵌到了一些分析工具，例如TMM），基因表达就通过一些计建模工具，例如edgeR，DESeq2和limma+voom进行计算， 计算结果是一些差异表达基因或转录本的列表， 这数据用于下一步的可视化和生物学解释。

在B分析流程中，使用一些较新的免比对工具，例如Kallisto与Salmon，这些工具会在一步操作中 组装转录组并对相应的转录本进行定量。 这些工具的输出结果通常是 转录本定量的一些估计值（例如tximport，TXI） ，然后通过与A分析流程中相同的归一化和统计建模，产生出差异基因或转录本列表。

在C分析流程中，第一步是比对读长（这一步的工具通常是TopHat，虽然有些分析方法也会用STAR与HISAT），接头使用CuffLinks来处理原始读长，再然后是使用CuffDiff2包来输出转录本丰度的估计值，以及一个差异表达基因或转录本的列表。

其它常用的工具还包括StringTie，这个工具使用TopHat（或类似工具）的输出结果来组装一个转录本模型，然后将结果输出到RSEM或MMSEQ中，用于估计转录本的丰度值，最后将转录本的丰度值输出给Ballgown来计算差异表达基因或转录本。而SOAPdenovo-trans这个工具则能同时对读长进行比和组装，其结果用于输入给RSEM或MMSEQ。

TABLE2-RNA-seq数据分析工具

第1阶段-测序读长的 比对(alignment)与组装(assembly)

测序完成后，分析的起点就是数据文件，这个数据文件包含了测序计数的碱基，这些数据文件通常是以FASTQ文件的格式存在。处理这些FASTQ文件最常见的第一步操作就是将测序读长 回贴到已知的转录组上（或已经注释的基因组上） ，将每个测序读长转换为一个或多个基因组坐标。这一过程可以使用多个不同的比对工具，例如TopHat，STAR或HISAT，它们都依赖于一个参考基因组。

由于测序的cDNA都源于RNA，而RNA有可能跨外显子边界，因此当与参考基因组（含有内含子与外显子）进行比对时，这些工具进行一个剪接比对后，测序读长之间会出现一些间隙。

如果测序的物种没有一个可用的高质量基因组注释（含有已经知的外显子边界），或者说如果希望将测序读长与转录本（而不是基因）关联起来，那么可以使用比对的读长进行转录组的组装。一些组装工具， 例如StringTie，SOAPdenovo-Trans 使用 利用那些已经比对好的结果中的空隙(reads的gap) 来推测其外显子边界，以及可能的剪接位点。

当参考基因组注释没有或者是不完整时，或者是你感兴趣的组织（例如在肿瘤组织）中转录本异常的情况下，这些从头组装转录本的工具尤其好用。当使用的是双端测序和/或更长的测序技术时，这种转录组组装方法效果更好，因为这些测序技术有更大的可能性跨越了剪接位点(splice junctions)。但是，从RNA-seq数据中进行转录本的完整组装对于计算DGE来说，并不是一个必需的步骤。

最近，已经开发出了计算高效的“免比对”(alignment-free)工具，例如Sailfish，Kallisto与Salmon， 这些工具可以直接将测序读长与转录本进行关联 ，从而无需单独的定量步骤（参考后面的第2阶段部分）。这些工具在那些表征更高丰度（以及更长的）转录本方面表现得非常良好；然后它们在那些定量低丰度或短转录本方面表现不佳。

用于将测序读长回贴到转录本的不同的工具在 它们如何对测序的子集进行分配方面有着显著的差异，这会影响最终的表达估计值。 当有来自一个不同基因，伪基因或转录本的多个读长时，这种效果尤为明显。一项比对12个基因表达估计方法的比较显示， 一些比对方法低估了许多与临床相关的基因的表达，这主要是由于其并不精确的回贴读长(ambiguous reads) 所导致的。

在RNA-seq数据的计算分析中， 如何将多个回贴的读长(比对到多个位置的reads) 合理进行分析仍然是一个值得研究的方向。 通常的做法是将这些读长排除在下一步的分析之外（一种常见的做法是在定量前过滤掉这些reads，），但这可能会使结果产生偏（参考阶段2-转录本丰度的量化）。其它的估计包括生成“合并“表达特征， 这些特征包含了那些共同回贴后的读长的重叠区域， 以及要在随后的置信度计算中，对每个估计的基因的回贴的不确定性进行估计。（生成包含合并映射重叠区域的 “融合”表达特征 ，以及计算每个基因的映射不确定性估计，以用于 后续的置信度的计算。 ）

第2阶段-转录本丰度的量化

一旦读长被回贴到基因组的位置或转录组的位置，接下来的步骤就是 将这些读长分配到基因或转录本上，以确定它们的丰度。

不同的比较研究表明， 在量化步骤中采用的不同的方法对最终的结果影响最大，这种影响甚至超过了第1步中比对工具的选择。

对每个基因(即该基因所有转录本的亚型)测序读长丰度的量化依赖于 转录组注释来对那些重叠到已知基因上的测序读长进行量化。 但是，使用短读长对测序读长进行特定异构体分配来说还需要一个估计步骤， 因为许多读长并不能跨越剪接位点 ， 因此它们无法精确地分配给特定的异构体（isoforms）。（把短reads分配到特定isoforms则需要统计模型估计，尤其是很多reads不跨越剪接点，并且不能明确分配给特定isoform时。）

当一 个基因的主要表达形式 是在不同长度的转录本（isoform）之间进行转换时，那么即使在 仅研究基因层面的 差异表达分析的前提下，对这些转录本进行定量则会产生一个更加精确的结果。例如，在一个样本中， 一个主要的转录本 也许只有 另外一个样本中同样转录本的一半长度 ，但是前者的 表达量是后者的2倍 ，那么单纯地基因基因定量的工具无法区分这个转录本的差异表达。

常用的量化工具包括RSEM，CuffLinks，MMSeq与HTSeq以及前面提到的那些免比对工具。

一些基于读长计算的工具，例如HTSeq（或者是R equivalent，featureCounts） 通常会丢弃许多比对好的读长，包括那些回贴到多个位置的读长，或者是重叠多个表达特征的读长。 其结果就是，在随后的分析中 清除了那些同源和重叠的转录本。

RSEM会使用期望最大化的方法来分配那些比对不明的读长，Kallisto这个无参比对工具会将比对不明的读长包括在它们相应的转录本计数中，从而导致结果偏倚。使用tximport包可以将转录本丰度估计转化为读长计数值(read count equivalents)，转录本丰度估计可以转换成等效的read计数。量化步骤中产生的结果通常会 合并为一个表达矩阵， 在这个矩阵中，每一行是表达特征（基因或转录本），每一列是样本名，表达矩阵中的值要和是实际的读长值，要么是一种估计丰度。

第3阶段-过滤和归一化

通常来说，定量后的基因或转录本计数结果 还需要过滤和归一化 ，从而用于 解释读长深度，表达模式以及技术偏倚。

过滤 去除在所有样本中都低丰度表达的基因 是很直接的方式，并且已经证明可以改善对真正差异表达基因的检测。

而对表达矩阵进行归一化则更加复杂。直接转换可以调整丰度值，以便能更加说明GC含量的差异以及读长深度。 早期用于归一化的方法就是RPKM，但这种方法现在已经淘汰， 它已经被那些能够校正样本之间更细微差异的方法所取代，例如 四分位数法或中位数归一法。

广泛的研究表明，归一化方法（normalization）的选择会对最终的结果以及生物学结论产生重要的影响。

大多数进行归一化的算法依赖于两个重要的假设：

第一，大多数基因的表达水平在重复的样本组之间保持不变；

第二，不同的样本组在总体的mRNA水平上不表现出显著差异。

当这些基础假设不成立时，那么就需要慎重考虑是否进行归一化，以及如何进行归一化。例如，如果一些基因 在一个样本中高表达 ，同时相同的基因，以及另外的一些基因在同一组中的 另外一个样本里正常表达 ，如果对 读长深度进行简单的归一化则明显不够， 因为相同数目的测序读长会 分配到第二个样本里更多的基因上面。

归一化过程，例如截断均值化M值(The Trimmed Mean of M-values,TMM)方法（它已经整合到的edgeR包中）就能解决这个问题。

选择一个合适的归一化方法或许很困难； 一种做法就是深度使用多种方法进行分析，然后比较它们结果的一致性。如果结果对于归一化方法高度敏度，则应该对数据进行进一步的探索，以确定差异来源。 但是，比较不同的归一化方法时，要谨慎确保这种归一化方法的比较不是为了选择与原始假设最兼容的归一化方法。

处理这些问题的一种方法的 spike-in control RNAs， 这种方法会引入一些外源已知的RNA序列，这些外源已知的 RNA序列已知，浓度已知， 在建库的过程中，将它们添加到样本中。

RNA-seq中的Spike-ins方法包括 外部RNA控制协会混合物 (External RNA Controls Consortium  mix, ERCCs)，s pike-in RNA突变物( spike-in RNA variants, SIRVs )与测序spike-ins(sequencing spike-ins, Sequins )。

由于预先知道spike-in的浓度，这些浓度直接与生成的读长数相关，因为就可以 对这些来自样本转录的表达水平进行校正。也有人指出，如果不进行spike-in控制时，那么就不能对那些有强烈倍数变化基因（变化大是基因）的实验进行分析。

然后在实践中，很难在预设水平上一致地整合spike-ins，并且它们 在基因水平上对测序读长数目进行归一化比转录本更加可靠 ，因为在一个样本中，每个 异构体（isoform）的表达水平不同 。

目前，尽管已发表的RNA-seq DGE实验中spike-in对照并未得到广泛使用，但随着单细胞实验的开展这一状况可能会改变，因为单细胞RNA-seq中spike-in应用广泛，当然前提是这个技术能进一步优化达到稳定的水平。

第4阶段-差异表达的统计建模

获得表达矩阵后，就可以 构建统计模型 评估哪些转录本发生了显著的表达改变。

有几个常用工具可以完成此任务； 一些基于基因水平的表达计数，其它的基于转录本水平的表达计数。

基因水平的工具通常依赖于比对的reads计数，并使用广义线性模型来进行复杂实验设计的评估。这些工具包括EdgeR， DESeq2 和 limma + voom等 工具，这些工具计算效率高并且彼此之间结果稳定性好。

评估差异isoforms表达的工具，例如CuffDiff，MMSEQ和Ballgown，往往需要更多的计算资源，并且结果的变化也更大。但是， 在差异表达工具应用之前的操作（即关于比对、定量、过滤和标准化）对最终结果的影响更大 。