具有读映射不确定性的RNA-Seq基因表达估计

摘要

1简介

利用快速发展的测序技术，研究人员现在正在使用高通量测序仪，通过一种名为RNA-Seq（Nagalakshmi等。,2008; 王等。,2009). RNA-Seq是全基因组鸟枪测序的转录组类似物（Staden，1979)，关键区别在于RNA-Seq主要用于估计样本中转录本的拷贝数。该方案的主要步骤是：（i）从样品中分离RNA，（ii）通过逆转录和片段化将RNA转化为cDNA片段，（iii）高通量测序器[例如来自Illumina（加利福尼亚州圣地亚哥）、Applied Biosystems（加利福尼亚州福斯特城）和Roche 454 Life Science（美国康涅狄格州布兰福德）的测序器]用于从cDNA片段生成数百万个读取，（iv）使用比对工具将读取映射到参考基因组或转录集，（v）使用映射到每个基因的读取计数来估计表达水平。因为RNA-Seq的主要输出是读取计数，所以它们被称为“数字”基因表达测量，而不是来自微阵列的“模拟”荧光强度。

RNA-Seq是一种很有希望的微阵列替代品，因为初步研究表明，RNA-Seq-表达估计值具有高度可重复性（Marioni等。,2008)根据定量PCR（qPCR）和尖峰试验（Mortazavi等。,2008; 那加拉克什米等。,2008). 尽管这项技术还很年轻，但RNA-Seq已经足够成熟，可以用于酵母（Nagalakshmi等。,2008),拟南芥（李斯特等。,2008)，鼠标（Cloonan等。,2008; 莫塔扎维等。,2008)和人类（马里奥尼等。,2008; 莫林等。,2008). RNA-Seq的主要优点是其大的动态范围（跨越五个数量级）、低背景噪声、需要较少的样本RNA和检测新转录物的能力，即使在没有测序基因组的情况下（Wang等。,2009). 为了使该技术发挥其全部潜力，需要解决一些实验和计算挑战，包括处理读取映射不确定性、测序不一致性、估计潜在的新亚型（或者剪接转录物）表达水平和RNA-Seq读取的有效存储和校准。

在本文中，我们介绍了我们在解决读取映射不确定性的计算问题方面的工作。由于RNA-Seq读数并不跨越整个转录物，因此它们来源的转录物并不总是唯一确定的。同源基因家族、低复杂性序列和同一基因交替剪接异构体之间的高序列相似性是导致绘图不确定性的主要因素。此外，多态性、参考序列错误和测序错误要求在读取比对中允许不匹配和索引，这进一步降低了我们对映射的信心。由于这些因素，大量读取是多重读取：与参考基因组或转录集中的多个位置具有高得分比对的读取（Mortazavi等。,2008). 我们将把对多个基因的解读称为基因多重阅读并将该图解读为单个基因但多个亚型异构多重读取基因多重读取的可映射读取的比例不同，并取决于转录组和读取长度。对于我们分析的数据集，这一比例在17%（小鼠）和52%（玉米）之间，代表了RNA-Seq数据的很大比例。

以前有两种策略用于处理基因多重读取。第一个简单地丢弃它们，只保留用于表达式估计的唯一映射读取（Marioni等。,2008; 那加拉克什米等。,2008). 一种稍微复杂一些的方法，只使用唯一映射读取，通过外显子位置的分数来调整外显子覆盖率，从而产生唯一映射读取（Morin等。,2008). 第二种策略是通过独特的读数映射，将多个读数的一部分按比例分配给基因，从而“拯救”多个读数（福克纳等。,2008; 莫塔扎维等。,2008). 研究表明，与仅使用唯一映射读取（Mortazavi等。,2008). 最近的一种方法通过显式估计亚型表达水平来处理亚型多重读取，但不处理基因多重读取（Jiang和Wong，2009).

我们提出了一种在存在多重读取时估计表达式的方法，该方法以统计严格的方式处理映射不确定性。使用RNA-Seq数据的生成模型，我们通过代表真实映射的潜在随机变量统一了基因和亚型多重读取的概念。模型参数对应于异构体表达水平、转录物之间的读数分布和测序错误。我们使用期望最大化（EM）算法估计最大似然（ML）表达水平，并表明以前的救援方法大致相当于EM的一次迭代。我们的推理方法可以被认为是纠正SAGE测序错误方法的RNA-Seq模拟（Beissbarth等。,2004)和微阵列交叉杂交（Kapur等。,2008).

就像江和黄的统计方法(2009)，我们的模型可用于估计个体亚型表达水平。与他们的方法相比，我们的方法将基因多重阅读纳入统计模型，不需要知道哪些亚型共享外显子序列。此外，我们的模型可以适应跨转录本的非均匀阅读分布（例如3′偏差），这可以同时从数据中学习。

使用从实际数据中获得的参数进行模拟的结果表明，与仅使用唯一读取或救援策略的方法相比，我们的方法可以提供更准确的基因表达估计。我们表明，将基因表达水平估计为估计的亚型表达水平的总和可以提高基因水平的准确性。我们在实际数据集的读取分布中发现轻微的3′偏差，并确定考虑这种不均匀性可以提高准确性，具体取决于偏差的强度。最后，通过不同读取长度的模拟，我们表明当我们的方法用于处理多重读取时，RNA-Seq的最佳长度约为小鼠和玉米转录体的20–25个碱基。

1.1问题陈述

表达式分析的目标是估计转录组：在给定时间内细胞中所有表达的转录物及其频率的集合。RNA-Seq数据本身允许我们估算相对的样本中亚型的表达水平。结合准确的物理样本量或峰值测量，可以估计绝对表达（Mortazavi等。,2008)，但这是一个单独的问题，我们不会在这里讨论。

有两种相对表达的自然度量：成绩单分数和核苷酸的分数由特定基因或亚型组成的转录组。对于异构体我，我们将用τ表示这两个量_我和ν_我分别是。在亚型水平上，这些量由等式关联

(1)

(2)

哪里ℓ_我是亚型的长度，单位为核苷酸我在基因水平上，表达只是可能亚型表达的总和。为了便于记法，我们根据百万分之核苷酸（NPM）和百万成绩单（TPM），通过将ν和τ乘以10得到⁶分别是。

本文解决的问题是使用一组RNA-Seq数据来估计一组先前确定的亚型和基因的ν和τ值。我们假设RNA-Seq数据的形式为N个序列读取，每个长度我参考序列，但不一定是基因组坐标，被假设可用于所有M（M）亚型。如果需要，序列聚类或基因组坐标可用于将亚型分组为基因。

RNA-Seq的基本假设是，从异构体中获得的读取部分我是ν的函数_我假设读数均匀分布，所有读数都可以分配给亚型和poly（A）+mRNA，我们得到c（c）_我/N个方法ν_我作为N个→∞, 哪里c（c）_我是从异构体读取的次数我。即使读取值不是沿亚型长度均匀分布的，只要它们与ν成比例采样_我，这个结果仍然成立。

1.1.1与RPKM估计值的比较

我们将ν和τ值的估计问题与计算表达式的问题进行了比较每千字节读取数/百万映射读取数（RPKM；摩塔扎维等。,2008). 亚型的测量水平我，以RPKM为单位，定义为10⁹×c（c）_我/(N个_米ℓ_我)，其中c（c）_我是映射到异构体的读取次数我,N个_米是可映射读取的总数ℓ_我是亚型的长度我。在均匀分布读数的假设下，我们注意到RPKM测量值估计为10⁹× ν_我/ℓ_我，是τ的非标准化值_我。归一化因子为

它与转录组中转录本的平均长度成正比。当表达的平均转录长度为1 kb时，1 TPM相当于1 RPKM，大约相当于小鼠（Mortazavi）中每个细胞一个转录物等。,2008).

由于样本之间的平均表达转录长度可能不同，我们更喜欢使用ν和τ，而不是RPKM测量。例如，如果其他基因的表达发生变化，导致平均表达的转录物长度不同，则两个样本中转录物分数相同的亚型将具有不同的RPKM值。此外，即使其中一个物种的mRNA长度趋向于更大，两个物种之间的τ值也具有可比性。

2方法

2.1生成模型

我们使用RNA-Seq读取测序过程的生成模型估计基因和亚型表达水平。我们使用如中所示的有向图形模型（贝叶斯网络）图1。模型生成N个i.i.d.长度读数我。读取序列由R（右）_n个随机变量和是观测数据。每个读取与三个潜在的随机变量相关联，G公司_n个,S公司_n个和O（运行）_n个，分别表示读取的亚型、起始位置和方向。模型的主要参数为θ=[θ₀,…, θ_M（M）]，对应于表达式级别。此模型的完整数据可能性为

在单独的窗口中打开

图1。

我们的方法使用的RNA-Seq数据的图形模型。

我们假设我们得到了所有M（M）转录组中可能存在的亚型。随机变量G公司_n个取[0中的值，M（M）]，值0表示“noise”异构体，它生成不映射到已知异构体的读取。我们让对(G公司_n个=我|θ)=θ_我，其中∑_我θ_我=1.

随机变量S公司_n个取值为[1，max_我ℓ_我]，其中Ş_我是亚型的长度我。我们打电话给对(S公司_n个=j个|G公司_n个=我)的读取起始位置分布（RSPD）。假设阅读在成绩单中统一生成，对(S公司_n个=j个|G公司_n个=我)=ℓ_我⁻¹这里，我们还假设mRNAs具有poly（A）尾巴，允许读取从亚型的最后位置开始并延伸到poly（A）尾巴。对于poly（A）-mRNA样本，RSPD变为对(S公司_n个=j个|G公司_n个=我)=(ℓ_我−我+1)⁻¹。我们在下一节中提出了一个非均匀RSPD模型。

随机变量O（运行）_n个是二进制的，带有O（运行）_n个=0表示读取顺序n个与其亲本同种型的方向相同，并且O（运行）_n个=1表示反向补码。这个随机变量允许我们对非特定于链的RNA-Seq协议进行建模，例如Mortazavi中使用的协议等。(2008). 对于此类协议，我们设置对(O（运行）_n个=0|G公司_n个≠0）=0.5。对于特定于股的协议，我们可以简单地设置对(O（运行）_n个=0|G公司_n个≠0)=1.

模型的观测数据是读取序列，我们用R（右）_n个随机变量。为了简化符号，我们总结了n个-用一组指示符随机变量读取Z轴_{日本[日本]}，其中Z轴_{日本[日本]}=1，如果(G公司_n个,S公司_n个,O（运行）_n个)=(我,j个,k个). 当协议特定于股时，我们使用变量Z轴_日本=Z轴_日本0.从亚型导出的读取序列的条件概率我>0由给定

其中γ^我是同种型的序列我,是它的反向补语w个_吨(一,b条)是位置特定的替换矩阵w个_吨(一,b条)是我们观察到性格的概率一在位置吨参考亚型序列中的相应字符为b条.由w个矩阵允许我们对读取位置和基相关替换过程进行建模。例如，在读取的最后位置，基本错误通常更频繁，因此w个_我可能给出比w个₁除了测序错误外，多态性和参考序列错误等因素可能会导致观察到的亚型序列和衍生阅读之间的替换。我们当前的模型没有区分这些过程，并通过w个矩阵。

从“噪音”亚型得出的读数(我=0）由位置无关的背景分布β生成

随机变量S公司_n个和O（运行）_n个在以下情况下不相关G公司_n个=0，因此我们设置对(O（运行）_n个=0|G公司_n个=0)=1,对(S公司_n个=1|G公司_n个=0）=1和ℓ₀=1。

2.1.1非统一RSPD

读取起始位置的分布通常是不均匀的（Wang等。,2009). 这种不一致性可能是由多种因素造成的，包括碎片协议和合成偏差（Dom等。,2008). 为了考虑这些不均匀性，我们允许在模型中使用经验RSPD。目前，我们使用的RSPD取决于给定起始位置沿异构体长度的分数。具体来说，我们使用

哪里电子海图_RSPD公司是[0,1]上的经验累积密度函数，表示为分段线性函数B类参数，=₁…ϕ_B类我们对所有亚型使用一个RSPD，因为通常没有足够的数据来估计单个亚型的这些分布。该方案允许我们模拟一般现象，如5′或3′偏差，但不考虑异构体特定效应（例如序列组成）。

2.2 EM算法推断

对于亚型表达水平估计，我们感兴趣的是推断模型参数θ=[θ₀, θ₁，…，θ_M（M）]. 在假设从转录组中均匀采样读取的情况下，这些参数对应于相对表达水平。根据这个假设，我们估计了ν_我通过并使用方程式(2)用于转换为τ_我.

给定RNA-Seq数据，我们通过找到最大化观测数据可能性的θ值来估计表达水平：

(3)

方程式(三)显示了我们的目标是找到混合模型的ML比例。我们使用EM算法（Dempster等。,1977)找到θ的ML值。我们假设模型的所有其他参数，即对(第页_n个|G公司_n个=我)要么是固定的，要么是提前估计的。对于固定对(第页_n个|G公司_n个=我)，观测数据的可能性[方程式(三)]是凹面的（参见补充材料)，从而保证EM算法能够找到ML值。然而，即使有无限的数据，参数也可能是不可识别的，这取决于亚型（Lacroix）之间的结构关系等。,2008). 这种情况对应于似然曲面中的平台。

我们模型的EM算法的实现很简单，详细信息在补充材料算法的一个关键方面是E步：计算Z轴_尼日（或Z轴_{日本[日本]})随机变量，给定当前参数值θ^(吨)。对于统一RSPD和特定于绞线的协议，此计算为

(4)

2.2.1通过读取对齐进行近似

由于读取、异构体和读取起始位置的数量很大，因此计算所有的预期值是不现实的Z轴_日本因此，我们通过只允许少量Z轴_日本变量不为零。为了确定要考虑哪些变量，我们将读取与异构体序列对齐，最多保持所有对齐x个不匹配。出租π^x个_n个表示所有的集合(我,j个)读取对齐对n个最多x个失配和噪声对（0，1），我们通过以下公式近似似然函数

方程式(4)类似地近似为

(5)

注意θ^(吨)_我/ℓ_我与亚型的RPKM成正比我，我们观察到(5)大致相当于Mortazavi的救援计划等。(2008)当一个read被认为与π中的所有候选位置一致时_n个并且没有考虑噪声亚型。因此，救援方案可以看作是EM算法的一次迭代，初始参数为θ⁽⁰⁾，仅从唯一映射读取派生。

为了加快计算速度，我们还过滤了与大量位置对齐的读取，并调整θ估计值以解释这些丢弃的对齐（请参阅补充材料).

3结果

3.1模拟实验

在缺乏我们知道所有真实mRNA数量的样本的RNA-Seq数据的情况下，我们进行了模拟以验证我们的方法，并评估其相对于其他RNA-Seq-表达水平估计器的性能。我们使用我们的生成模型来模拟两个参考转录集的读取，一个来自小鼠，另一个来自玉米。我们的方法以及其他两种方法的估计值与每个模拟中的样本表达值进行了比较。

3.1.1比较方法

我们将我们的方法的基因表达水平估计值与之前使用的两种处理多头的方法的基因表达水平估计值进行了比较。这些先前的方法都不能估计亚型表达水平，因此我们只比较了基因表达估计。第一种也是最简单的方法，独特的，仅根据唯一映射读取来估计表达式。如果一个read仅与单个基因的亚型对齐，则认为它是唯一映射的。此方法将表达式级别计算为

(6)

哪里c（c）_我^联合国是将该映射读取到基因的唯一映射数我和c（c）^联合国是唯一映射读取的总数。然后通过以下公式计算成绩单分数

(7)

哪里是基因的有效长度我对于具有单一亚型的基因，就是那个亚型的长度。对于选择性剪接基因，我们取是与基因亚型外显子相对应的所有基因组区间的联合长度我如ERANGE软件包（Mortazavi）中所述等。,2008).

第二种方法，我们称之为营救，按照τ的比例将多重读取分配给基因_我^联合国值。这种方法是在Mortazavi中引入的等。(2008)目前用于ERANGE软件包以及RSAT软件包（Jiang和Wong，2009). 这个营救该方法计算基因计数我作为

(8)

其中π_n个是读取的一组基因指数n个地图。仅映射到τ为^联合国=0在这些基因中平均分配。ν的值_我^营救和τ_我^营救然后进行与等式类似的计算(6)和(7).

3.1.2模拟程序

我们首先从Mortazavi中描述的小鼠肝脏RNA-Seq数据推导出模拟参数等。(2008). 核苷酸分数（ν）由该数据计算得出营救方法（以减少对我们方法的偏见）和小鼠UCSC基因（Hsu等。,2006)作为参考基因集。对于具有多个亚型的基因，这些部分在亚型之间随机划分。我们估计了位置特定替换矩阵(w个_吨)来自唯一映射读取和来自不可映射读取的“噪声”读取模型（β），它们占数据的46.2%。我们通过运行Bowtie对齐器（Langmead等。,2009)最多允许两个失配。使用这些参数、统一的RSPD和非转移特异性模型，我们用我们的生殖模型模拟了小鼠参考基因集中长度为25的3000万次读取。

另一个具有相同大小和读取长度的模拟集是从玉米基因集生成的。我们选择玉米来评估基因组重复性的变化如何影响表达估计。使用的玉米基因集是从网址：http://maizesequence.org，删除重复的mRNA。每个玉米基因的表达水平都是从小鼠肝脏基因表达估计值中取样并替换的。

我们使用Bowtie映射模拟读取，每个对齐最多允许两个不匹配。表1总结了真实数据集和模拟数据集上的读取映射结果。基因多重读取分别占小鼠和玉米数据集中可映射读取的17%和52%。阅读映射到的基因和亚型的数量分布如所示补充图6.

表1。

不可匹配、唯一映射、映射到多个基因或在三个RNA-Seq数据集中过滤的读取片段

数据集	%未映射	%独一无二的	%多个	%已过滤
鼠标实数	46.2	44.4	9.2	0.2
鼠标模拟	47.6	43.2	8.7	0.6
玉米模拟	47.5	25	27.1	0.4

在单独的窗口中打开

3.1.3仿真结果

这个独特的,营救和相对长度单位方法以与输入相同的比对集运行，以估计基因表达水平。将估计值与样本值进行比较，即根据每个亚型的样本读取真实计数得出的表达式估计值[通过方程式(6)和(2)具有真实读取计数]。我们用样本值而不是模型参数进行比较，因为如果观察到所有潜在变量，样本值是可以做出的最佳估计。与模型参数的比较受到采样误差的影响，采样误差同样影响所有方法，并掩盖了它们之间的差异。

我们使用了两种度量表达式估计误差的方法。首先，作为对异常值稳健的精度的一般度量，我们计算了样本值的中值百分比误差（MPE）。其次，我们计算了估计值与样本值显著不同（误差百分比>5%）的基因比例。我们将第二个统计数据称为误差分数（EF）。对于这两种度量，我们只考虑样本ν或τ分别不小于1 NPM或1 TPM的基因。我们还计算了假阳性率，这是样本τ低于1 TPM的基因的分数，估计τ至少为1 TPM。表2给出了小鼠和玉米模拟数据估计值的MPE和EF统计数据。估计值与样本值的散点图见补充图1和2.

表2。

的错误独特的,营救和相对长度单位根据小鼠和玉米RNA-Seq数据的模拟，相对于样品表达值的估计基因表达水平

			NPM（ν）或TPM（τ）中的样本基因表达
			[1, 10)	[10, 102)	[102，10三)	[10三，104)	[104，105)	全部
小鼠RNA-Seq数据的模拟
		N个	5577	5240	1028	114	9	11968
	MPE公司	独特的	18.9	18.7	19.1	19.9	20.7	18.8
		营救	2.8	1.1	0.8	0.7	1.2	1.6
ν		相对长度单位	2.3	0.8	0.4	0.2	0.3	1.1
	EF公司	独特的	93.9	96.2	96.5	100	100	95.2
		营救	26.9	6.1	6.4	7.9	33.3	15.9
		相对长度单位	18.8	2	0.8	0	0	9.7
		N个	6279	4025	886	111	15	11316
	MPE公司	独特的	29.6	29.2	30.9	32.8	32.1	29.6
		营救	12.6	6.8	6.1	5.9	5.8	8.2
τ		相对长度单位	2.6	1	0.4	0.3	0.4	1.5
	EF公司	独特的	93.7	93.9	95.6	99.1	100	94
		营救	79.5	73.2	72.2	69.4	66.7	76.6
		相对长度单位	27.8	6.2	1.1	0	0	17.7
玉米RNA-Seq数据的模拟
		N个	8934	4737	988	119	14	14792
	最大允许功率	独特的	86.8	87.8	88.7	88.1	85.9	87.3
		营救	11.3	3.3	0.9	0.6	0.7	6.6
ν		相对长度单位	3.7	1.2	0.6	0.5	0.4	2.3
	EF公司	独特的	97.3	97.3	97.5	93.3	100	97.3
		营救	65.8	42.6	22.7	11.8	7.1	55
		相对长度单位	40.5	16.5	6.4	2.5	21.4	30.2
		N个	9210	4931	1040	113	12	15306
	MPE公司	独特的	86.1	84.2	85.2	80.5	96.3	85.5
		营救	21.3	11.8	8.9	8.5	7.7	16
τ		相对长度单位	4.6	1.5	0.6	0.5	0.3	2.8
	EF公司	独特的	97.2	96.7	97.1	98.2	100	97
		营救	89.4	88.3	85.8	82.3	91.7	88.8
		相对长度单位	47.5	18.8	6.1	4.4	16.7	35.1

在单独的窗口中打开

给出了不同表达水平的基因以及表达至少为1 NPM（ν）或1 TPM（τ）的所有基因的误差度量。粗体值表示估计值显著(对<0.05）更准确，通过配对Wilcoxon符号秩检验进行评估。除了一个类别，相对长度单位比其他方法准确得多。对于高度表达的类别（10⁴−10⁵NPM）玉米，相对长度单位实际上在νEF方面的表现比营救我们将这种奇怪的现象归因于该类别中少数基因（14个）中的几个重复基因。

对于鼠标营救和相对长度单位两种方法都将大多数基因的σ估计在百分之几以内相对长度单位估计值通常更准确（成对Wilcoxon符号秩检验，对=7.2 × 10⁻²⁷³). 这个独特的该方法给出的估计值要差得多，总体MPE为18.8%，EF为95.2%。的准确性营救和相对长度单位随着更多数据可用于为这些基因分配多重读取，对于高表达的基因而言，该值最高。EF统计表明，这些方法的准确性差异较大相对长度单位产生最小数量的异常值估计值（总νEF为9.7%，相比之下为15.9%营救). 在τ估计方面，两种方法之间的差异甚至更大。这种差异主要是由于以下事实营救不估计单个亚型表达水平，而是使用单个“有效”长度来计算τ。最后，假阳性率独特的,营救和相对长度单位方法小鼠数据分别为1.8%、0.97%和0.89%。

对于玉米数据集，广泛的近期基因paralogy使得RNA-Seq分析具有挑战性，因为唯一的映射读取在所有可映射读取中占少数。尽管如此相对长度单位这组数据的估计值为2.3%，相比之下营救估计。然而，μ估计值的EF远高于小鼠，分别为30.2%和55.0%相对长度单位和营救分别是。玉米数据集的假阳性率也要高得多独特的,营救和相对长度单位方法假阳性率分别为4.2%、7.7%和4.4%。

3.3最佳读取长度的确定

排序技术可以产生不同的读取次数和长度。假设一个经济模型中RNA-Seq的成本与测序的碱基数量成正比，我们希望确定在表达估计中允许最高精度的读取长度。给定固定的吞吐量，可以产生大量的短读或少量的长读。读取时间越长，由于映射不确定性导致的估计误差越小，而读取次数越多，由于采样导致的误差越小。最佳读取长度可以在这两个错误源之间实现平衡。

我们使用固定吞吐量和不同读取长度进行了模拟，以分析这种权衡。为了评估物种和组织之间的差异，我们从小鼠肝脏、小鼠大脑和玉米转录组进行了模拟。至于小鼠肝脏模拟，小鼠大脑模拟的参数是根据Mortazavi的数据确定的等。(2008). 我们将吞吐量固定为7.5亿个碱基（大致相当于小鼠肝脏数据），并生成长度为20到40个碱基的读取集。考虑到这种固定的吞吐量，每个集合中的读取次数从1875万次不等(我=40）至3750万(我=20). 对于小鼠肝脏，我们进行了额外的实验，其吞吐量是其他三个实验的两倍（15亿），一半（3.75亿）。根据相应组织和物种的表达水平参数模拟读数，位置无关错误率为1%。路线中允许的最大不匹配数为1我≤22，较长读数为2。对每个阅读长度进行了五次模拟，并使用我们的方法估计了表达水平。图2显示了这些集合上估计值的平均τMPE。ν和EF测量的趋势相同（未显示数据）。我们的结果表明，20到25个碱基之间的读取长度对所考虑的转录组和吞吐量具有最高的准确性。

在单独的窗口中打开

图2。

给定固定的基础吞吐量，基因表达估计精度随读取长度而变化(T型). 曲线为（1）小鼠肝脏，T型=375 × 10⁶，（2）小鼠肝脏，T型=750 × 10⁶，（3）小鼠肝脏，T型=1.5 × 10⁷，（4）小鼠大脑，T型=750 × 10⁶（5）玉米，T型=750 × 10⁶τMPE是根据真实水平至少为1 TPM的所有基因的真实表达值计算的。

4讨论

在本文中，我们证明了RNA-Seq的统计模型可以用于解决读取映射的不确定性，因此，与简单的方法相比，可以产生更准确的基因表达估计。通过模拟真实数据，我们证实了我们的方法提高了小鼠和玉米实验的准确性。精确度的提高对于重复基因组来说最为显著，例如玉米，它会产生大量的多重读取。此外，我们关于τ估计准确性的结果表明，估计单个亚型表达水平对于计算转录分数的基因水平估计至关重要。

我们选择关注基因水平的准确性，因为我们的初步工作表明，异构体水平的估计仅适用于具有当前读吞吐量的最高表达基因。当基因多重读数可以忽略不计时，我们预计我们的亚型估计值将与Jiang和Wong方法得出的结果相比较(2009)，因为它基于类似的模型。当基因多重读数显著时，我们已经证明我们的方法比营救方法，在Jiang和Wong中使用(2009)用于分配基因多重读取。

我们的最佳读取长度结果可能对未来的RNA-Seq分析具有重要意义。我们表明，给定固定的测序器吞吐量，在已知基因组的生物体上进行RNA-Seq实验以估计基因表达的最佳读取长度出奇地短。我们确定玉米和两个小鼠组织的转录体的长度在20到25个碱基之间。这一结果表明，当目标应用程序是基因表达时，测序技术开发人员应该争取更多的读取次数，而不是更长的读取时间。较长的读取当然有助于减少映射的不确定性，但使用能够处理这种不确定性的方法，可以通过大量的短读取获得更高的准确性。进行RNA-Seq实验的研究人员通常应使用测序器试剂盒产生的完整读取长度，但如果在以读取长度换取读取数量的试剂盒之间进行选择，则应选择产生更多读取的试剂盒。

我们的最佳阅读长度结果取决于以下假设：（i）转录组中所有可能的异构体都是已知的，（ii）只需要基因水平的表达估计，（iii）所有基因的准确性是主要目标，而不是特定子集的准确性。当这些假设不成立时，较长时间的阅读可能是有益的。此外，最佳读取长度可能取决于目标种类和吞吐量，因此应针对其他条件进行额外的实验，以确定合适的读取长度。

虽然我们的准确度结果是基于模拟的，但我们认为这是对表达估计方法的公平测试，因为模拟的是遵循标准假设的RNA-Seq数据。模拟的主要假设是，读取在转录组中均匀生成，或至少与ν成比例_我对于每个亚型（在转录本长度不一致的情况下）。所有当前的方法，包括我们自己的方法，都依赖于这个假设。模拟的一个替代方案是与微阵列估计值进行比较，这是为了验证RNA-Seq技术的总体效果（Marioni等。,2008; 莫塔扎维等。,2008). 然而，微阵列评估不能被视为基本事实，事实上，其精度低于RNA-Seq（Wang等。,2009). 理想情况下，我们将使用RNA-Seq数据集评估准确性，该数据集的亚型水平已通过qPCR确定，但由于qPCR的相对低吞吐量性质，很难获得这样的数据集。

由于我们的主要动机是读取映射的不确定性，因此我们方法的有效性取决于RNA-Seq数据集中的多重读取数。随着定序器每年产生更长的读取长度，人们可能预计多重读取的数量会显著下降。然而，有几个原因可以解释为什么在阅读跨越整个RNA分子之前，多重阅读仍然是相关的。首先，较长的读操作和成对的读操作不会像预期的那样减少多重读操作的数量。我们对小鼠转录组上75个碱基读取和配对的75个碱位读取（200个碱基插入）的模拟分别产生了10%和8%的基因多重读取。与25个基本读取的17%多重读取相比，这些仍然是重要的分数。其次，我们的结果表明，RNA Seq的最佳测序仪应该产生大量的短读数（20-25个碱基），在这种情况下，多读数是高度相关的。第三，亚型多重阅读甚至对于较长的阅读长度也是显著的，因为替代亚型通常共享其序列的很大一部分。由于准确的表达估计取决于识别单个亚型的表达，因此正确处理亚型多重读取仍然是一个重要问题。

我们未来的工作将包括一些改进，以更紧密地模拟RNA-Seq读取生成过程。例如，我们将扩展该模型以合并测序器质量分数、成对读取、可变读取长度（例如在454数据中），并允许在读取中使用索引。RNA-Seq数据的其他不一致性，例如对某些碱基组成的偏差（Dohm等。,2008)，也将得到解决。

补充材料

致谢

参考文献