RNA-seq数据中的转录长度偏差混淆了系统生物学

摘要

背景

高通量测序很可能成为转录组分析的首选平台。测序平台查询整个转录环境的能力为了解生物系统中转录复杂性的水平提供了新的见解。转录测序也提供了新的机会，例如定量测量剪接变体[1]和单核苷酸多态性（SNPs）用于等位基因特异性表达，而无需任何先验知识。然而，这种新的细节水平需要仔细的统计建模，以提供新技术所承诺的好处。

不同的技术显示不同的数据特征，因此具有不同的优势和劣势。对数据的技术和统计属性的调查将揭示每种技术的优缺点。我们假设，使用统计方法检测样本之间的差异表达受到转录长度的影响，并且这种影响是标准RNA-seq过程固有的。

目前的RNA-seq协议在测序之前使用mRNA片段化方法，以获得整个转录物的序列覆盖率。简单来说，这意味着给定转录本的总读取次数与转录本的表达水平乘以转录本的长度成正比。换句话说，与表达相似的短基因相比，长转录本将具有更多的读映射。由于实验的威力与样本大小成正比，因此检测较长基因的差异表达的威力更大。这是数据的固有属性，不会被任何涉及对片段短于转录本的全长转录本进行排序的过程所改变。相比之下，微阵列的强度测量值仅与转录物的表达水平以及探针本身固有的任何特征（如GC含量）成正比[2,三]. RNA-seq数据中较长转录物的较高采样率的特征在识别样本之间差异表达基因的情况下变得重要。大多数用于检测差异表达的统计方法对读取次数较多的转录本具有更强的能力。因此，与同一样本中的长抄本相比，短抄本在统计上总是处于劣势。

在这里，我们探索了几个先前发布的高通量RNA测序数据集，并表明最广泛使用的RNA-seq协议与较短的转录物相比，可以检测到较长转录物中更多的差异表达。对于使用单个或一组诊断探针评估表达水平的微阵列平台，不存在这种偏见。我们还证明，在迄今为止所分析的所有实验中，无论是特定样本、平台还是统计分析，都存在此处提出的固有偏差。

结果

基于简化的计算，我们可以证明，使用非常简单的测试程序检测差异表达的能力取决于转录本的长度（参见方法）。为了实证研究RNA-seq和微阵列中转录长度偏差的影响，我们使用了三个不同的已发布数据集对全长转录进行排序。第一组数据比较了Illumina基因组分析仪和Illumiana微阵列的测序结果，并寻找人类胚胎肾脏和B细胞系之间的差异表达[4]. 第二种方法使用SOLiD测序比较小鼠胚胎干细胞和类胚体[5]第三个比较了Illumina测序和Affymetrix微阵列在人类肝脏和肾脏样本中的测序结果[6]. 为了确定哪些基因差异表达，每项研究都使用不同的统计方法计算显著性。在这里，我们使用这些统计数据来检查转录长度差异表达的行为。

对于每个平台，我们首先根据转录长度将所有基因装箱到相等的基因数箱中。接下来，我们根据相关出版物中定义的统计程序的截止值，将基因指定为差异表达（DE）。我们发现使用的特定统计截止值对结果没有实质性差异。然后，我们计算了每个箱子中DE基因的百分比​图11显示了每个实验中绘制为RNA-seq和微阵列转录长度函数的DE基因百分比。显然，无论平台、统计分析程序或差异表达的总体比例如何，检测DE的能力与RNA-seq的转录长度密切相关（图1a、c和​和1d）。1天). 正如预期的那样，来自两个不同平台的微阵列数据没有观察到这种趋势（图​（图1b1亿和​和1e第1页).

在单独的窗口中打开

图1

差异表达与转录长度的关系数据根据转录长度进行分类，并绘制使用统计截止值表示的差异转录百分比（点）。还绘制了线性回归曲线（直线）。一–e（电子）使用RNA-seq的所有数据和研究中的微阵列[4-6]分别是。（f）和克从微阵列数据中提取33%的高表达基因（蓝色交叉）和33%的低表达基因（红色三角形），用于分析出现在两个平台上的基因[6]. 回归显示了差异表达百分比与转录长度的显著趋势一,c（c）,d日和（f）和低表达基因克请注意，此图说明了不同实验之间的常见数据特征，而不是平台、方法或实验之间的比较。

为了进一步研究转录长度偏差，我们使用了Marioni等人（2008）的数据，并观察了RNA-seq和微阵列数据中出现的转录物。我们根据微阵列上测量的平均强度水平将基因分为三个相等的组。然后，我们计算了%DE作为高表达组和低表达组长度的函数（图​（图1f1英尺和​和1g）。1克). RNA-seq对低表达基因显示出更强的长度偏差，这在一定程度上有所改善，但对高表达基因来说仍然显著。我们认为高表达基因的斜率较低，因为几乎所有这些基因都有足够的能力在这个数据集中被称为差异表达，即使较短基因的p值较高。相比之下，使用微阵列平台未观察到高表达基因的显著趋势，而低表达基因的趋势很小。尽管对于RNA-seq，被称为DE基因的总数可能更大，但增加重复次数将增加微阵列检测到的DE基因的数量。对平台间DE绝对速率的仔细校准可能是进一步调查的主题，但此处确定的趋势与所使用的特定统计截止值相一致，并且对其稳健。

许多当前的统计方法都使用由基因长度归一化的表达水平度量。这提供了对表达水平的无偏测量，但也会以长度相关的方式影响数据的方差，从而导致对差异表达估计的相同偏差（参见玩具示例的方法）。为了证明这种效果，我们展示了从Marioni等人的数据的复制路线计算出的每个基因的样本平均值和方差。在图中​图2a2a个我们证明了样本均值和方差近似等于泊松随机变量的期望值。然而，当平均值除以转录本的长度时，关系变得更加复杂，数据不再是泊松。图​图2b2亿显示了相同的数据，每个转录本的标签计数除以转录本的长度。两条拟合曲线显示了三分之一具有最长转录本和最短转录本的数据的均值-方差关系。很容易看出，对于按长度归一化的基因，较短的转录物与较长的转录物相比，在相同表达水平下具有更大的差异。

在单独的窗口中打开

图2

均值-方差关系这里我们展示了Marioni等人的肝脏样本中跨车道的样本方差[6]数据绘制为每个基因平均值的函数（a）。接下来，我们有相同的数据，其中每个基因的标签数除以基因的长度（b）。红线拟合三分之一最短基因的均值和方差之间的线性关系，而蓝线拟合最长基因的线性关系。在绘图中一拟合非常接近平均值和方差之间的相等线（黑线），这是泊松过程的预期结果。在绘图中b条对于给定的表达水平，短基因比长基因具有更高的方差。

当在具有不同长度的基因或基因集之间进行比较时，RNA-seq数据中转录物长度偏差的后果变得最成问题。这最有可能发生在系统生物学中进行基因集测试时，与其他基因集相比，特定基因集具有长度偏差。如果一个集合包含的基因短于平均值，那么在差异表达中会出现低表达，而如果集合包含的是长于平均值的基因，那么在差别表达中该类别更有可能出现高表达。为了证明这种效应，我们使用Marioni等人的RNA-seq和微阵列数据寻找了过度代表的KEGG通路。在每次分析中，我们只使用两个平台上发现的基因，然后使用DAVID软件进行通路分析[7,8]. 我们发现肝脏和肾脏之间的差异表达有几个过度表达的途径。桌子​表11和​和22在p值0.1以下显示所有过度表达的通路（然而，经过多次测试校正后，只有前16个类别仍保留在微阵列数据中，其重要性与测序数据相同）。以粗体突出显示的类别不会在其他平台的列表中的任何位置出现过度显示。经过多次测试校正后，微阵列平台包含四条低于阈值0.1的通路，所有这些通路都可以在测序数据中找到。相比之下，RNA-seq数据包含九个类别，其中三个类别不包含在阵列数据的任何地方。图​图3三显示了与每个类别相关的基因的长度。第一个方框显示了基因在两个平台上显著的路径中的分布。第二个框给出了仅在测序平台中出现显著的基因（即，不出现在阵列平台的列表中的任何位置），第三个框是分析中所有转录物的长度分布。可以清楚地看到，仅在RNA-seq平台上过度表达的类别中的基因明显长于平均值。

在单独的窗口中打开

图3

KEGG通路中发现的基因长度显著高于差异表达基因图中的第一个方框表示在两个平台的四个重要类别中发现的基因长度。第二个方框是仅在测序数据中显著类别中发现的基因长度。第三个方框是这两种技术共有的所有基因的长度。可以看出，测序数据特有的类别往往具有更长的转录本。

讨论

RNA-seq数据中的转录长度偏差是采样过程的可预测结果，不能通过除以转录长度来校正（例如Cloonan等人（2008）或Sultan等人（2008）中的统计方法）。泊松随机变量的长度偏差是预期的，其中基因的预期读取计数与转录本的长度和表达水平成正比。换句话说，与短基因相比，长基因的采样率更高，因此无论使用何种特定测试，在给定的统计显著性下检测差异表达的能力更强。检测样本间差异表达的统计检验需要估计样本的平均值和方差。将平均值除以成绩单的长度可以消除该度量的偏差，但随后会将长度偏差引入方差中，问题仍然存在。类似地，与微阵列一样，对于高表达水平的基因，检测差异表达的能力更强。然而，我们在这里并没有关注这一现象，因为我们认为与转录长度偏差相比，表达水平是一个具有生物学意义的数量，而转录长度偏差是技术性的。

高通量测序可以使用不同的处理方法来确定转录水平，例如大规模并行签名测序（MPSS）、基因表达序列分析（SAGE）、基因表达式帽分析（CAGE）。这些方法只计算每个转录本的一个序列，因此不会受到转录长度的影响。然而，在许多情况下，研究人员希望通过对整个RNA序列进行测序来检查转录组的全部复杂性，因此我们推测这些无偏见的方法将构成RNA-seq实验的一小部分。由于片段大小是转录长度偏差的基础，而不是读取的碱基数，因此当前平台上读取长度的改进不会改变转录长度偏差。

使用外显子水平分析作为减少转录长度范围的一种方法可能不会显著减少偏差。虽然在查看人类基因组中基因和外显子的长度时，基因的长度明显长于外显子，但对数中的四分位间距（IQR）₂长度相似（基因IQR=1.45，外显子IQR=1.23）。这意味着长度加倍的基因数量与长度加倍的外显子数量相似，这意味着偏见也同样强烈。扩大测序深度将提高检测差异表达的能力。然而，由于转录长度偏差是基因之间的相对度量，因此不会影响呈现的结果。

目前，我们建议的解释基因间长度偏差的唯一方法是使用固定长度窗口方法，窗口大小小于最小的基因。在这种方法中，可以计算每个窗口的聚集标记数，并评估差异表达。进一步的扩展可以将每个转录本的多个窗口组合成一个单一的度量，方法类似于在RMA微阵列算法中组合多个探针[9]但这一建议需要进一步探索和分析。尽管如此，由于需要在窗口级别进行分析，这将需要丢弃部分数据或为每个基因引入可变数量的窗口。此外，窗口的较小尺寸将需要每个样本的较大总读取次数，以实现统计显著性，从而在以基因为中心的分析中降低到小基因的等效水平。

重要的是要理解，使用统计截止值从聚集的标签计数生成DE基因列表，不可避免地会对具有较长转录本的基因更加敏感。因此，包含具有不同长度分布的基因的基因集和本体类在基因集测试中可能看起来被错误地低估或高估。由于基因集测试是许多系统生物学实验和大型生物医学项目（如国际癌症基因组联合会）的一个组成部分，因此长度偏差将显著影响目前正在使用RNA-seq的许多应用。将需要复杂的统计方法来开发一个新的分析框架，为不同的转录长度提供类似的错误发现率。

结论

目前可用的RNA-seq和微阵列技术的不同优势使得这些平台对转录组的综合分析具有互补性。使用高通量测序来询问全长转录物的一个技术特征是，相对于类似表达的短转录物，较长的转录物产生更多的读数。这种更高的抽样意味着，与短转录本相比，长转录本检测差异表达的统计能力更强。使用三个已发布的RNA-seq数据集，我们已经证明，正如我们假设的那样，更长的转录物更有能力检测差异表达。这些数据集使用不同的样本、平台和统计方法。这种偏差在微阵列数据中不存在。如果无法解释，检测差异表达能力中的转录长度偏差可能会混淆系统生物学和基因集测试方法。了解新技术中的技术问题将导致开发更复杂的分析方法。

方法

在这里，我们表明，在一些非常简单的假设下，当测试给定表达水平的两个样本之间的差异时，较长的基因比较短的基因更重要。让X（X）是图书馆中与特定成绩单对应的读取次数。的预期值X（X）与转录本总数成比例N个乘以基因的长度L（左）

其中c是比例常数。假设数据分布为泊松随机变量，则方差等于平均值。

例如，我们可以使用t检验来测试相同库大小的两个样本之间特定基因计数的差异是否与零显著不同

(1)

哪里D类是两个样本的观测平均值之差瑞典（D）是的标准误差D类.

t检验的功效取决于E（D）/S.E.（D）=δ本质上是t分布的非中心性参数

可以看出，这与长度的平方根成正比。因此，对于给定的表达水平，测试对于更长的转录长度更为重要。

按基因长度除法

考虑到上面的简单设置，我们可以看到表达水平除以基因长度的效果。在这里

和

鉴于我们假设X（X）作为泊松随机变量分布，我们可以看到，一旦我们除以长度，分布就不再是泊松分布μ' ≠变量(μ').

使用上述t检验的相同格式，我们发现我们恢复了等式1中相同的长度相关性。

哪里

t检验的威力再次取决于E（D）/S.E.（D）=δ哪里

我们得到了与等式1相同的测试结果，等式1仍然有平方根L（左）依赖。

缩写

DE：差异表达；IQR：四分位范围。

竞争性利益

作者声明，他们没有相互竞争的利益。

作者的贡献

AO和MJW构思了这个想法并撰写了论文。AO分析了数据。

审核人意见

鸣谢

工具书类