背景
微阵列是一种同时测量许多基因表达水平的有效技术，但这种方法存在一些局限性。低表达基因的表达估计值通常不可靠，因为真正的信号被交叉杂交效应掩盖了[1,2]. 此外，阵列的设计依赖于基因结构的注释，因此该方法对于发现新的剪接事件并不理想。最近开发的一种替代方法，称为RNA-Seq，有潜力克服这些困难[三]. RNA-Seq使用超高通量测序[4]确定大量cDNA片段的序列。生成的序列（读取）可以是长序列（>100核苷酸）或短序列，具体取决于平台[4]. 目前流行的两个短读平台是Illumina的Solexa[5-11]和应用生物系统公司（ABI）的SOLiD[12]. 每一个都可以在一次运行中生成数千万次的短读[5-12]. 在本文中，我们只考虑短读RNA-Seq。
RNA-Seq产生的读数首先通过计算机程序映射到基因组和/或参考转录本。然后，RNA-Seq的输出可以通过一系列“计数”进行汇总。也就是说，对于基因组中或假定转录本上的每个位置，它给出了一个计数，表示从该位置开始映射的读取数。例如（为了简化起见，我们缩短了基因并读取），如果一个具有单一亚型的基因具有序列ACGTCCCC，并且我们有12个ACGTC读取，8个CGTCC读取，9个GTCCC读取和5个TCCCC读取，那么这个基因可以通过计数12、8、9、5的序列进行总结。
RNA-Seq数据的定量推断，例如计算基因表达水平[7]和亚型表达水平[13]，基于这些计数。为了有效地利用数据，为这些计数建立一个适当的统计模型是至关重要的。当前的分析方法显式或隐式地假设一个朴素的恒速泊松模型，在该模型中，来自同一亚型的所有计数都是从泊松分布中独立采样的，单个速率与亚型的表达水平成比例[7,13,14]. 不幸的是，我们发现这个模型不能很好地拟合真实数据（见结果），需要一个更精细的模型。
为了更好地建模计数，考虑可变速率的泊松模型是很自然的；也就是说，来自亚型的计数仍然被建模为泊松随机变量，但每个泊松随机变数具有不同的速率（平均值）。通过检查不同组织计数之间的相似性（见结果），可以看出泊松比不仅取决于基因表达水平，还取决于读数的位置。因此，我们将速率建模为基因表达水平和从该位置开始读取的“排序偏好”的乘积。这种排序偏好是一个显示在此位置生成读取的可能性的因素。
多姆等。[15]发现GC丰富的区域往往比AT丰富的区域有更多的读取，但我们发现纯粹基于GC内容的模型工作得很差（附加文件1)。通过回顾如何在微阵列中处理相关问题，可以获得一些关于如何建模测序偏好的线索。微阵列中的每个基因都有一组探针，每个探针都可以连续测量基因表达水平。同一组测量值采用不同平均值的高斯分布建模，每个平均值都是基因表达水平和探针对cDNA序列亲和力的乘积。Naef和Magnasco[16]提出了仅依赖于探针序列的探针亲和力模型：
哪里ω我是探针的亲和力我,K（K）是探针的长度，I(b条伊克=小时))当k个第个基本对是字母小时，否则为0，α和β千赫是我们想要估计的参数，以及ε是高斯噪声，因此可以通过正则线性最小二乘法估计参数。该模型的主要特点是，它考虑每个位置出现的字母，而不仅仅是每个字母出现的总次数。这个简单的线性模型可以解释Affymetrix寡核苷酸阵列数据集中44%的亲和力差异。为其他阵列或数据集开发了类似的模型[17-20].
在RNA-Seq实验中，cDNA合成通常由随机启动启动。根据其序列，mRNA片段可能形成阻碍引物结合的二级结构。此外，引物通常由非随机侧翼序列标记，该序列可能根据mRNA序列优先与mRNA相互作用。由于这些影响，结合的概率取决于核苷酸序列和方案。合成后，cDNA与连接子连接，扩增，然后测序。在这些步骤中，cDNA的二级结构和协议的细节会再次影响效率。因此，协议和局部序列上下文可能对读取mRNA片段的可能性有很大影响。因此，在特定方案下，我们可以根据局部核苷酸序列预测测序偏好，至少可以部分预测。
结果和讨论
数据集和过度分散
本文使用了三个全基因组RNA-Seq数据集。前两个是由Illumina的Solexa平台生成的，第三个是由ABI的SOLiD平台生成的。第一个数据集[7]由来自三种小鼠组织的7900万、7600万和7000万个读数组成：大脑、肝脏和骨骼肌。每个读数的长度为25。第二个数据集[11]由来自10种不同人类组织和5种乳腺上皮或乳腺癌细胞系的1200万到2900万个读数组成。每个读数的长度为32。我们使用其中九种组织或细胞系的数据，并将它们合并为三组（第一组为脂肪、大脑和乳房，第二组为结肠、心脏和肝脏，第三组为淋巴结、骨骼肌和睾丸）。每个组包含6100万到7700万个读取。第三个数据集[12]由两种细胞系（类胚体（EB）和未分化小鼠胚胎干细胞（ES））的1600万个高质量读数组成。每个原始读数是35个核苷酸，但有些被截断为30或25个核苷酸以确保高质量。我们将这三个数据集分别称为Wold数据、Burge数据和Grimmond数据，这与最初生成数据的研究小组一致。正如我们刚才所描述的，这三个数据集中的每一个都包含几个代表不同组织、组或细胞系的子数据集，总共有八个子数据集：三个（组织）代表Wold数据，三个（组）代表Burge数据，两个（细胞系）代表Grimmond数据。在我们的所有处理和计算中，上述子数据集是单独考虑的；也就是说，一次只分析一个子数据集。
首先，从原始数据集中提取计数数据。材料和方法中描述了详细的程序。简言之，我们将读取映射到所有RefSeq基因的所有亚型，然后为了避免歧义，我们只计算唯一映射到RefSeq中只注释了一个亚型且与其他基因不重叠的基因的读取，我们称之为“非重叠单转录基因”。此外，我们只使用表达水平最高的前100个基因的计数来拟合我们的模型，因为它们具有最高的信噪比（参见附加文件1详细信息）。
两项证据清楚地表明，计数违反了恒速泊松模型。首先，数据严重过度分散。泊松分布的一个基本性质是均值和方差相等。如果方差大于平均值，则数据被称为过度分散，泊松假设不合适。表列出了每个子数据集中前100个基因的方差-均值比（也称为“Fano因子”）的最大值、中位数和最小值。所有比率都远大于1。其次，基因计数的“模式”（相对值）在同一数据集的不同子数据集中惊人地保守。图显示基因中的计数阿波（载脂蛋白E）的所有三个组织的沃尔德数据。尽管计数的绝对值在不同组织中相差100倍，但不同组织的变化模式高度一致。沃尔德数据的其他基因以及伯格和格里蒙德数据的基因也是如此。这有力地证明，来自同一基因不同位置的计数并不是从同一分布中取样的。相反，计数的分布似乎取决于其序列在转录本中的位置。这迫使我们考虑更复杂的模型。Hansen也描述了阅读率偏差强烈依赖于局部序列的观察结果等。[21]这是我们在审查论文时注意到的一项独立工作。
基因读取数阿波在Wold数据的不同组织中. (一)大脑(b条)肝脏(c（c）)骨骼肌。每条垂直线代表从该位置开始的读取计数。灰色线是UTR区域的计数，再加上100 bp。在这里，内含子被删除，外显子连接成一个整体。只显示了一条基因链上的计数；另一条链上的计数在不同组织中显示出相似性。Nt：核苷酸。
表1
| | 方差与平均值之比 |
|---|
| |
|
|---|
| 数据集 | 子数据集 | 最大值 | 中间带 | 最小值 |
|---|
| 沃尔德 | 大脑 | 248 | 36 | 21 |
| 肝脏 | 1,503 | 48 | 19 |
| 肌肉 | 2,088 | 34 | 18 |
| | | | |
| 伯格 | 第1组 | 835 | 78 | 14 |
| 第2组 | 1,187 | 102 | 28 |
| 第3组 | 1,593 | 112 | 20 |
| | | | |
| 格里蒙德 | 电子束 | 24,385 | 806 | 47 |
| 锿 | 9,162 | 345 | 22 |
泊松线性模型及其性能
对于核苷酸j基因的我，我们想对从这个核苷酸（表示为n个ij公司)取决于该基因的表达水平（表示为μ2)以及围绕该核苷酸的核苷酸序列（带有长度的序列K（K）表示为b条ij公司1,b条ij公司2, ⋯,b条ijK公司,). 我们假设n个ij公司~泊松(μij公司)，其中μij公司是泊松分布的速率，以及μij公司=ωij公司μij公司，其中ωij公司是排序偏好，可能取决于周围的序列。作为一种简单的方法，我们对偏好使用线性模型，从而得出泊松率：
哪里ν我=对数(μ我),α是一个常数项，I(b条ijk公司=小时)等于1，如果k个第个周围序列的核苷酸是小时，否则为0，以及β千赫是字母的影响系数小时发生在k个第个位置。该模型使用大约3K个参数来模拟排序偏好。为了拟合上述模型，我们迭代优化了基因表达水平和泊松回归系数（材料和方法）。
我们将我们的模型应用于八个子数据集中的每一个子数据集。作为局部序列上下文，我们在读取的第一个核苷酸之前使用40个核苷酸，在它们之后使用40个核酸（即读取的前40个核苷酸；请参阅附加文件1出于选择该区域的原因）。因此，我们的模型使用3×80=240个参数来建模排序偏好。与每个子数据集中的样本大小（约100000个计数）相比，这是一个相对较小的数字。
在线性回归中，可以用回归来解释的方差百分比，表示为R（右）2，用于测量菲特的质量。在泊松回归中，我们可以用偏差代替方差，并定义：
哪里天是拟合模型的偏差，以及天0是空模型的偏差[22]. 在我们的例子中，空模型是假设相同排序偏好的朴素模型。决赛R（右）2表中列出了我们获得的值粗略地说，这个简单的线性模型可以解释大约40%到50%的方差。
表2
| |
R（右）
2
|
|---|
| |
|
|---|
| | 泊松线性 | 市场 |
|---|
| |
|
|
|---|
| 数据集 | 子数据集 | 80个核苷酸一，非交叉验证 | 80个核苷酸一，交叉验证 | 40个核苷酸一，交叉验证 | 40个核苷酸一，交叉验证 |
|---|
| 沃尔德 | 大脑 | 0.52 | 0.51 | 0.51 | 0.70 |
| 肝脏 | 0.51 | 0.50 | 0.50 | 0.70 |
| 肌肉 | 0.48 | 0.46 | 0.46 | 0.59 |
| | | | | |
| 伯格 | 第1组 | 0.43 | 0.42 | 0.42 | 0.52 |
| 第2组 | 0.37 | 0.35 | 0.35 | 0.46 |
| 第3组 | 0.45 | 0.42 | 0.42 | 0.54 |
| | | | | |
| 格里蒙德 | 电子束 | 0.47 | 0.40 | 0.40 | 0.58 |
| 锿 | 0.45 | 0.39 | 0.37 | 0.54 |
图显示了线性模型中的所有系数。每个系数的渐近标准误差约为0.002，因此几乎所有系数都具有统计意义。这并不奇怪，因为我们的样本量远远大于参数的数量。在这种情况下，更重要的是系数的大小。通常，图中央部分的系数的绝对值大于两侧的系数，后者接近于零。这表明，read第一个位置周围的核苷酸对测序偏好有更大的影响。这是合理的，因为这些核苷酸往往形成一个读取局部二级结构的头部，它只涉及几个核苷酸，因此很容易预测。虽然更远的核苷酸可能会形成带有读头的非局部二级结构，但由于它涉及的核苷酸太多，并且可能因情况不同而有显著差异，因此很难预测其结构。
不同数据集中泊松线性模型的系数当我们将周围序列视为读取第一个核苷酸之前的40个核苷酸和之后的40个核酸时，八个子数据集中泊松线性模型的系数。位置-1、0、1分别表示read的第一个核苷酸之前的核苷酸、read的首个核苷酸和read的第二个核苷酸。核苷酸的颜色编码：红色，T；绿色，A；蓝色，C；黑色，G。核苷酸T（红色）的系数是基本水平，所以它们总是零。(一)沃尔德数据中的系数。子数据集的形状编码：矩形，大脑；三角形，肝脏；圆形，骨骼肌。(b条)伯格数据中的系数。子数据集的形状编码：矩形，组1；三角形，第2组；圆，第3组。(c（c）)Grimmond数据中的系数。子数据集的形状编码：矩形，EB；下面是如何读取这些系数的示例。在沃尔德大脑数据中，读取的第一个核苷酸（（a）中位置0处的蓝色矩形）中的C系数为0.82。这意味着如果核苷酸T被C取代，那么测序偏好将增加到e（电子）0.82=2.27倍。Nt：核苷酸。
同一数据集的每个子数据集中的系数都非常相似，尽管它们在不同的数据集中有显著差异。这有力地证明了这些系数是有意义的，而不仅仅是随机的。
虽然很难从生物学角度解释每个系数的大小，但我们可以通过使用的协议解释数据集之间系数的主要差异。Wold和Burge数据都是使用Illumina平台生成的，因此它们的曲线看起来很相似，尤其是在中部。然而，在Wold数据中的cDNA合成之前，mRNA被切割成大约200个核苷酸片段，但在Burge数据中没有。较短的mRNA片段不太可能形成非局部二级结构。因此，沃尔德数据的系数曲线应该有较轻的尾部。Grimmond的实验使用了ABI的平台进行测序，并在测序之前向合成的cDNA添加了不同的连接子，因此整个曲线看起来与Wold和Burge的数据完全不同。
我们的泊松线性模型表明，至少37%到52%的不均匀性可以由序列差异来解释。然而，这个百分比可能低估了局部序列上下文可以解释的偏差部分，因为简单的线性模型无法捕捉到许多其他影响。在线性模型中添加更多的预测因子是可能的，特别是添加二核苷酸组成可以显著改善拟合（附加文件1)，但我们倾向于考虑非线性模型，以便更好地理解计数的不均匀性在多大程度上是系统偏差而不是随机噪声。
MART模型及其性能
尝试过支持向量机和神经网络等方法（附加文件1)，我们选择了MART（多元加性回归树）作为非线性模型的最终选择。MART是Friedman提出的一种梯度树增强算法[23,24]. “gbm”包中提供了MART的一个版本[25]第页，共页[26]. 此外，为了避免非线性模型通常出现的过度拟合，我们使用交叉验证和R（右）2在测试数据中。
关于使用MART和估计交叉验证的详细信息R（右）2在材料和方法中给出。在此分析中，我们使用较短的周围序列。对于Wold和Burge数据，我们在reads的第一个核苷酸之前使用了25个核苷酸，在reads之后使用了15个核苷酸，对于Grimmond数据，我们使用了之前的15个核苷酸和之后的25个核苷酸。这些是泊松回归模型中系数较大的区域（附加文件1)。使用较短的周围序列可以降低输入数据的维数，从而缩短训练时间并减少过度拟合的机会。
最终交叉验证R（右）2表中列出了我们获得的值八分之七R（右）2值大于0.50，其中两个高达0.70。与线性模型相比，R（右）2增加0.10至0.20，显示MART模型的威力。图给出了两个方法如何执行的示例。图显示基因计数阿波在原始数据中，分别用泊松线性模型拟合计数和用MART拟合计数。很容易看出，MART更符合计数。因此，我们建议在根据数据进行任何统计推断时，应使用MART模型，而泊松线性模型仅用于为MART选择合理的周围序列区域。我们还注意到，与使用泊松线性模型确定的拟合数相比，使用MART确定的拟合计数沿基因的变化更快，但在这两种情况下，变化都不如原始数据中的剧烈。实际上，这两种方法拟合的计数的方差与平均值之比分别为55和91，均小于原始计数的127。这表明我们的两个模型仍然给出保守拟合。
配件数量阿波基因黑色垂直线表示沿阿波基因（UTR和另外100个核苷酸被截断）。(一)沃尔德大脑数据中的读取次数（真值）。这与图1a的中心部分（黑色垂直线）相同。(b条)使用泊松线性模型的拟合读数计数。我们使用前100个基因中的其他99个基因来训练线性模型，然后用它来预测阿波。此预测具有（交叉验证）R（右）2= 0.54. (c（c）)使用MART的拟合读数计数。我们使用前100个基因中的其他99个基因来训练MART，然后用它来预测阿波。此预测具有（交叉验证）R（右）2= 0.69.
我们的高潮R（右）2结果表明，从局部序列中可以预测测序偏好中至少50%到70%的不一致性。
我们使用最高度表达的基因训练的模型可以用于预测其他基因的测序偏好。例如，我们使用仅使用前100个基因训练的MART模型预测了Wold数据的大脑样本对所有独特基因的偏好，结果总结如下R（右）2（图)。正如预期，R（右）2对于表达水平较低的基因来说，其变异性较小，因为不可预测的随机性在平均值较小的泊松分布中占较大比例。平均值R（右）2对于高表达或中等表达的基因，高于0.5（每百万个映射序列外显子的每千碱基读数（RPKM）>30）R（右）2RPKM>1的基因为阴性，表明我们的模型始终优于统一模型。请注意，在这些数据中，1 RPKM代表平均每核苷酸仅0.034次读取。
的箱线图R（右）2世界大脑数据中的独特基因根据RPKMs，我们将至少有一个读码的基因分为六组：<1、1-5、5-15、15-30、30-100和>100；每组分别包含4205、3320、2807、1330、1094和383个基因。注意，在这些数据中，1 RPKM代表平均每核苷酸0.034个读取，RPKM>30的基因被认为相对丰富，RPKM<1的基因即使用于转录检测也不可靠[7].
我们模型的应用
我们的结果可能有助于从RNA-Seq数据进行定量推断。为了减少由于读取速率的不一致性而导致的基因表达估计偏差，我们建议在我们的MART模型下，通过沿着基因的读取总数除以测序偏好总和（SSP）来估计单个亚型基因的表达。相比之下，标准估计值将读取数除以基因长度，这相当于在所有测序偏好设置为1的统一模型下除以SSP。
为了测试新方法，我们首先将使用Wold RNA-Seq数据的小鼠肝脏子数据集估计的基因表达水平与使用Kapur使用的相同组织的Affymetrix微阵列数据估计的基因表示水平进行了比较等。[27]. 对于RNA-Seq数据，我们在统一模型和MART模型下估计基因表达水平，对于微阵列数据，我们使用稳健多芯片平均值[28]. 所有非重叠的单转录基因都包含在比较中，结果由Spearman的秩相关系数汇总。对于所考虑的所有基因，与统一模型相比，使用我们的MART模型将秩相关从0.771增加到0.773，这表示一个非常小的改进。
我们对测序偏好的高度预测模型未能导致基因表达评估的更显著改进，原因是什么？我们认为，答案是，当一个基因较大时，当将测序偏好中的显著局部差异在多个位置上求和，以生成整个基因的SSP时，这些差异将被消除。在这种情况下，MART模型下的单一共享点与统一模型下的单个共享点相差不大，并且新的估计值与通常的估计值几乎相同。为了查看新的估计值在与标准估计值不同的情况下是否会导致改进，我们首先通过折合变化量化两个估计值之间的差异，定义如下：
沃尔德数据中跨基因的平均折叠变化仅为1.02；因此，新估计值的性能如此接近标准估计值也就不足为奇了。一致地，当我们检查100个折叠变化最大的基因时（平均而言，这100个基因中的折叠变化为1.10），等级相关性显示出更大的改善，从0.095到0.198，即108%的相对变化。
表给出了不同数据集的1号染色体基因、外显子和连接的平均折叠变化。我们发现，根据测序偏好、测序平台和生成数据的实验室的平均区域大小，折叠变化可以大大大于1。例如，Grimmond数据显示，跨基因的平均折叠变化为1.25。因此，我们预计新的估计将显示该数据有更大的改善。为了看看情况是否如此，我们注意到卡普尔等。[27]计算来自小鼠胚胎样本的Affymetrix微阵列数据的基因表达水平，我们可以使用这些数据评估Grimmond EB数据的新估计值和标准估计值。对于所有考虑的基因，秩相关系数从标准估计的0.439增加到新估计的0.469，相对变化为6.9%。根据SSP的折叠变化，我们进一步将这些基因分为五类，每一类包含约20%的所有基因。表显示了每个箱子中基因表达水平的秩相关系数。很明显，在具有较大倍变化的基因中会出现较大的改进。对于20%的折叠变化最小的基因，改善幅度仅为0.1%左右，而对于20%折叠变化最大的基因，提高幅度约为26%。最显著的是，在折叠变化最大的100个基因中，秩相关从0.323变为0.526，相对提高了62.8%。这些结果表明，我们基于建模测序偏好的新估计可以显著改善基因表达估计。
表3
| 平均排序偏好的平均倍数变化 |
|---|
|
|
|---|
| 用于训练模型的数据集 | 基因 | 外显子 | 交叉点（读取长度=35） | 交叉点（读取长度=100） |
|---|
| 沃尔德 | 1.02 | 1.12 | 1.13 | 1.07 |
| 伯格 | 1.18 | 1.32 | 1.37 | 1.28 |
| 格里蒙德 | 1.25 | 2.17 | 2.34 | 1.73 |
表4
| 折叠更换箱 | 均匀模型SCC | 我们的MART模型的SCC | 相对改进 |
|---|
| (1.00, 1.09) | 0.465 | 0.466 | 0.1% |
| (1.09, 1.19) | 0.437 | 0.444 | 1.4% |
| (1.19, 1.33) | 0.413 | 0.434 | 5.1% |
| (1.33, 1.53) | 0.481 | 0.520 | 8.2% |
| (1.53, 4.82) | 0.389 | 0.490 | 26.0% |
接下来，我们研究了测序偏好的结合是否可以改善对异构体特异表达水平的推断。我们修改了Jiang的异构体特异性表达估计等。[13]假设每个外显子的平均计数与外显子SSP成比例，而不是外显子长度。图显示了RefSeq基因的四种亚型Clta公司在鼠标中。在均匀模型下[13]Grimmond EB数据的亚型表达分别为21.6%、53.4%、8.95%和16.0%（总和为100%）。当考虑到排序首选项时[13]分别为15.5%、52.9%、10.8%和20.7%。基于新表达式级别和序列首选项的新计数更适合数据（未显示数据）。
RefSeq基因的四种亚型Clta公司在鼠标中。此图是使用CisGenome浏览器生成的[36]. 顶部显示了小鼠4号染色体中的碱基位置和外显子的灰色块。底部显示四种亚型，外显子放大。第一种亚型的外显子1的尾部比其他三种亚型少6 bp。第二个亚型有7个外显子，而第三个亚型同时缺失了外显子5（54 bp）和外显子6（36 bp），第四个亚型缺失了外隐子6。
回到沃尔德数据，我们从表中注意到外显子SSP的倍数变化为1.12，这表明MART模型和统一模型在外显子水平估计方面可能存在足够的差异。为了评估这两个模型在外显子水平估计方面的性能，我们将我们对亚型表达水平的估计与Pan中给出的估计进行了比较等。[29]他使用定制微阵列研究了10个小鼠组织中的3126个“盒式”选择性剪接（AS）事件。每个组织中的每个AS事件都由七个探针靶向，然后计算选择性剪接外显子排除值百分比（%ASex）作为汇总统计。在蒋的论文中等。[13]介绍了他们估算亚型表达水平的方法，他们比较了Pan的ASex百分比等。[29]根据三种小鼠组织（肝脏、肌肉和大脑）的统一模型计算%ASex。特别是，他们根据两个标准选择了AS事件的子集：一个标准要求基因的适度表达水平和相对狭窄的%ASex置信区间；另一种则需要适度比例的外显子排除亚型。我们使用相同的基因子集，考虑到MART预测的测序偏好，并使用他们的方法计算ASex百分比。结果总结见表对于几乎每个基因子集，当我们考虑测序偏好时，皮尔逊相关系数较高，平均相对改善约为7.2%。这表明我们的MART模型为亚型表达水平估计提供了有意义的改进，即使是对于不一致性最小的Wold数据也是如此。
表5
| 选择标准 | 组织 | 所选AS事件数 | 采用统一模型的PCC | 通过我们的MART模型进行PCC | 相对改进 |
|---|
| 1 | 肝脏 | 472 | 0.48 | 0.50 | 4.2% |
| 肌肉 | 451 | 0.40 | 0.45 | 12.5% |
| 大脑 | 699 | 0.36 | 0.40 | 11.1% |
| | | | | |
| 2 | 肝脏 | 228 | 0.60 | 0.60 | 0% |
| 肌肉 | 194 | 0.48 | 0.51 | 6.3% |
| 大脑 | 298 | 0.44 | 0.50 | 13.6% |
综上所述，我们发现决定我们的模型能带来多大改进的主要因素是褶皱变化的幅度。因此，我们期望我们的方法可以应用于涉及短序列元素的许多其他问题。在新的亚型发现中，一个当前非常感兴趣的问题，关键是要考虑沿该区域读取的相对计数。例如，一个碱基读数比其周围区域多的区域表明存在一个新的外显子。然而，如果仅仅因为该区域比其周围区域具有更大的测序偏好，所以该区域具有更多读取数，则可能会产生误导。需要进一步努力将我们的方法融入到当前的异构体发现算法中。
虽然MART模型可以更好地估计排序偏好，因此可用于统计推断，但泊松线性模型的主要目的是选择合适的K（K）用于MART模型。然而，我们仍有可能从中获得更多信息，尤其是从系数图中（如图)。例如，如果曲线中间部分的系数具有较大的绝对值，这可能表明测序偏好的差异在实验中反复扩大，很可能是通过多轮PCR，我们可能需要使用更多的mRNA样本，而不是进行过多轮PCR。另一个例子是，如果系数曲线有重尾，这应该表明信使核糖核酸/cDNA倾向于形成复杂的非局部二级结构，这也是不利的，我们可能需要将信使核糖核酸片段成更小的片段和/或选择具有适当长度的更好的连接体。经验丰富的技术人员知道实验的所有细节，可能会对偏差的主要原因提供更多的解释，甚至是精确定位。这可能有助于改进RNA-Seq的方案。
材料和方法
从原始读取数据中提取计数数据
首先，我们从UCSC基因组浏览器网站下载[30]RefSeq基因序列[31,32]（2007年7月，Wold和Grimmond数据为鼠标mm9，Burge数据为人类hg19）。然后，我们将读取映射到RefSeq基因的所有亚型。对于Illumina数据，我们使用SeqMap直接绘制了25或32个核苷酸读取[33]，允许两个不匹配。对于ABI数据，我们使用了与补充图中描述的相同的策略第页，共页[12]其中，分别对35、30和25个核苷酸合格读进行三轮定位。在每一轮中，我们都使用SOCS[34]作为映射工具。映射后，我们选择了RefSeq中只注释了一个亚型并且与其他基因不重叠的基因，并将其称为“非重叠的单一亚型基因”。为了避免歧义，我们只保留了指向唯一位点的读码，而该位点位于唯一基因内。然后，我们计算了从这些独特基因的每个位置开始映射的读取数，从而得出计数数据。由于读取长度较短，一些位置与其他位置具有相同的局部序列（读取长度），因此我们的计数方法总是为它们分配一个零计数。这可能会影响我们的分析结果。然而，即使读取长度仅为25，这些位置也不到所有位置的2%，因此它们不应显著改变我们的分析。
之后还要执行几个步骤。避免注释中的UTR歧义和排序中的边界偏差[三]，我们截断了所有UTR和两端的100个核苷酸。然后，我们在截短后丢弃了太短（少于100个核苷酸）的基因。最后，在计算RPKM测量的基因表达水平后[7]，我们丢弃了所有基因，除了表达水平最高的前100个基因。这些顶级基因的计数是我们用于拟合模型的唯一计数。这些顶级基因的读取在所有明确定位的读取中占相当大的比例，因此为测序偏好提供了足够的信息。相比之下，低表达的基因没有或只有很少的读取，而适度表达的基因在相当大比例的位点上通常为零；因此，关于它们的测序偏好的信息是有限的。
每个子数据集中前100个基因的计数数据可在R包“mseq”中获得[35]，可在CRAN（综合R档案网络）中公开获取。
拟合泊松线性模型
我们使用以下策略来拟合我们的泊松回归模型：
1.初始化
，其中L（左）我是基因的长度我.
2.查看
作为已知偏移量，拟合泊松回归模型以获得
和
。这是一个标准算法，R的“glm（）”[26]实现它。
3.更新
，其中W公司我是基因所有核苷酸的测序偏好之和我也就是说，
.
4.跳到步骤2，除非偏差减小小于1%。
在上面，步骤2给出了α和β千赫鉴于，很容易证明步骤3给出了ν我鉴于α=和β千赫=因此，上述过程通过反复优化偏好参数和基因表达水平来最大化可能性。
R程序包“mseq”中提供了执行此程序的R代码[35].
使用MART和估计交叉验证的策略R（右）2
使用MART和估计交叉验证的策略R（右）2包括以下步骤：（1）将100个基因随机分为5组。在交叉验证的每一次折叠中，使用其中一个作为测试集，使用其他四个作为训练集。（2） 在每个折叠中，对于训练数据集中的每个基因，将每个计数除以该基因中计数的平均值。结果数字被认为是该位置的排序偏好。为了避免第3步中麻烦的零偏好，我们将零计数替换为一个小数字（计算中为0.5）。（3） 获取这些首选项的对数。（4） 使用周围的序列作为输入，使用这些日志首选项作为输出来训练MART。我们用于MART的参数是：交互深度=10，收缩率=0.06，树数=2000（该方法对参数的选择是稳健的；附加文件1)。此外，由于高表达基因的方差较小，我们对其对数偏好赋予了更大的权重。基因对数偏好的权重我设置为N个我/L（左）我，其中N个我是该基因的总读取次数，以及L（左）我是这个基因的长度。（5） 使用经过训练的MART预测测试数据的日志偏好。（6） 获得基因表达水平的最大似然估计。也就是说，假设一个基因的长度是L（左），日志首选项为一1, ...,一L（左），计数为n个1, ...,n个L（左），则基因表达水平为
（7） 根据步骤5中的对数偏好和步骤6中的基因表达水平计算偏差。同时计算零偏差。（8） 对所有五个折叠重复步骤2至7。（9） 计算最终交叉验证R（右）2，这是五个折叠中的偏差之和，大于零偏差之和。
R程序包“mseq”中提供了执行此程序的R代码[35].
