利用共表达和RNA SEQ数据进行差异基因表达分析

 

MRFSEQ

As a fundamental tool for discovering genes involved in a disease or biological process, differential gene expression analysis plays an important role in genomics research.

High throughput sequencing technologies, e.g., RNA-Seq, are increasingly being used for differential gene expression analysis which was dominated by the microarray technology in the past decade. However, inferring differential gene expression based on the observed difference of RNA-Seq read counts has unique challenges that were not present in microarray-based analysis. The differential expression estimation may be biased against low read count values such that the differential expression of genes with high read counts is more easily detected. The estimation bias may further propagate in downstream analyses at the systems biology level if it is not corrected.

To obtain a better inference of differential gene expression, we have proposed a new efficient algorithm based on a Markov random field (MRF) model, called MRFSeq, that uses additional gene coexpression data to enhance the prediction power. Our main technical contribution is the careful selection of the clique potential functions in the MRF so its maximum a posteriori (MAP) estimation of can be reduced to the well-known maximum flow problem and thus solved in polynomial time. By using incorporating the loopy belief propagation technique, MRFSeq is able to not only call differentially expressed (DE) genes but also assign confidence scores to each inferred DE gene. Our extensive experiments on simulated and real RNA-Seq datasets demonstrate that MRFSeq is more accurate and less biased against genes with low read counts than the existing methods based on RNA-Seq data alone. For example, on the well-studied MAQC dataset, MRFSeq improved the sensitivity from 11.6% to 38.8% for genes with low read counts

抽象

动机：

 RNA-Seq越来越多地被用于差异基因表达分析，由于过去的十年里微阵列技术主导。然而, (基于观察到的RNA-Seq读数的差异)来推断差异基因表达,呈现出在基于微阵列的分析中所没有呈现的独特的挑战。

该差异表达估计是针对低读取计数值，高读取计数基因的差异表达更容易被检测到。估计偏差如果没有纠正在系统生物学级别的下游分析中可能会进一步传播。

结果：

为了获得更好的差异基因表达推论，我们提出了一种基于马尔科夫随机场的新型高效算法（MRF）模型，称为MRFSeq，它使用额外的基因共表达数据来提高预测能力。我们的主要技术贡献就是仔细选择MRF中的团聚势函数,所以它的最大后验估计可以降低到众所周知的最大流量问题,从而在多项式时间内解决。

我们在模拟和真实RNA-Seq数据集上进行大量的实验表明,与仅基于RNA-Seq数据的现有方法相比,  MRFSeq更准确，且对读取计数低的基因的偏倚更小。例如，在经过深入研究的MAQC数据集MRFSeq上低读数基因的敏感性从11.6％提高到38.8％。

 

可用性：MRFSeq以C语言实现，

可在http：// www.cs.ucr.edu/~yyang027/mrfseq.htm

联系方式：yyang027@ucr.edu或jiang@cs.ucr.edu

补充信息：补充数据可在

生物信息学在线。

2013年2月27日收到;于2013年6月6日修订;接受

2013年6月10日

1引言

在基因组学研究中已被广泛使用。RNA-Seq，下一代测序技术，最令人兴奋的应用之一，被用于揭示生物系统中转录组的复杂性（Wang等人，2009）。

许多前所未有的发现正在被RNA-Seq产生，例如新的亚型的推断，

反义调控模式的刻画和基因间表达模式的研究（Carninci等，2005; Graveley，

等人，2011; Nagalakshmi等人，2008; Trapnell等，2010）。

近年来，RNA-Seq在定量分析基因表达和转录变体发现方面发挥了重要作用。

在过去的十年中，这些应用大部分都是以基于微阵列的技术为主。

在这些定量分析中，RNA群体部分测序，获得的读取序列与参考基因组被对齐

。然后将对齐的读段分配给基因(基于它们共享的共同区域)。该分配给基因的阅读数被称为基因的读断数, 该基因已被证明几乎与基因的表达水平（是线性相关的Marioni等，2008）。
差异基因表达分析是鉴定基因在感兴趣的生物条件之间是否不同的基因表达。

鉴于RNA-Seq读取计数数据，可以通过检查观察到的读数差异是否显著,来检测差异表达 （DE）或同等表达（EE）基因来， 例如大于一些自然的随机变化。

为了测试RNA-Seq读取计数之间的差异的显著性, 首先将读取计数的分布被假定为泊松分布（Marioni等，2008; Srivastava和Chen，2010; Wang等， 2010）。

然而，泊松分布可能低估了读取计数的差异并在差异基因表达分析中,导致意外的假阳性（Nagalakshmi等，2008; 罗宾逊和史密斯，2007）。为了解决这个问题，对RNA-Seq读取数据应用于负二项式分布（Anders和Huber，2010; Robinson和Smyth，2007,2008; 罗宾逊等人，2010年），并已成为最先进的统计数据模型。

除了基于泊松或者基于负二项分布的方法之外，还提出了两个数据驱动的概率方法，baySeq（Hardcastle和Kelly，2010）和NOISeq （Tarazona等，2011）。此外，鉴于基因的带注释的或推断的mRNA转录物（或同种型），最近公布了用于检测转录物水平的差异表达的一些统计方法（格里菲斯 等人，2010;李和杜威，2011; Trapnell等，2013;郑和 陈，2009）。

 

由于可以通过简单地总结 其同种型的表达水平 来计算具有已知(或推断) 同种型基因的表达水平，所以这些转录水平的方法 可以用作检测同种型的差异表达替代方法（Trapnell等，2013），尽管这些方法的准确性显然取决于所提供的y亚型的质量。

 

尽管RNA-Seq数据的统计特性在上述统计方法中充分研究和考虑，但这些方法遭受以下问题。

首先，它观察到统计功效值随着读取计数值而增加（Oshlack和Wakefield，2009; Oshlack等，2010;Young等，2010）。

请注意，基因的读取计数是正比于基因表达水平乘以基因长度。

因此，与短和/或表现不佳的同行基因相比，长或高度表达的基因更可能被检测为DE基因。

DE基因检测中存在这种偏倚是不可避免的，甚至使使用于读取计数的数据, 标准化或重新标定（Oshlack andWakefield，2009; Young等人，2010）。

据了解，如果，DE基因的选择偏倚未校正的，可能导致偏向下游分析（Oshlack和Wakefield，2009; Oshlack等，2010; Young等，2010）。

其次，基因表达之间的依赖关系不是在这些方法中使用。

基于微阵列数据的基因表达分析，基因共表达模式的先验知识, 已被用于改善算法的性能

用于检测表型相关的通路（Rahnenfuhrer等，2004），

寻找重要的途径监管机构（Sivachenko等人，2005），

鉴定差异基因表达模式（Jacob等，2012）和

微阵列数据的分类（Rapaport等，2007）。

 尤其要获得更准确魏和李（Wei and Li，2007）提出了一个DE基因的推断

马尔可夫随机场（MRF）模型，整合基于微阵列数据的gammagamma模型

(Kendziorski等，2003; Newton等人，2001）和提取来自KEGG路径的基因共表达网络（Kanehisa and Goto，2000）等 ，DE基因可以通过MRF模型的最大后验（MAP）估计来确定。

他们的实验结果表明额外的基因共表达信息可以帮助检测基因表达的更微妙的变化（例如，局部的已知途径中的干扰）并显著改善DE基因的整体预测准确性（Wei和Li，2007）。

但是由于连续微阵列的强度值和不连续的RNA-Seq读数差异，MRF（Wei and Li，2007）的模型不能立即应用于RNA-Seq数据。

而且，由于MAP估计问题对于MRF模型通常是NP-Hard（Boykov，2001）

MRF模型（Wei and Li，2007）通过启发式算法解决方法，迭代条件模式，它提供了一个没有置信度分数的近似最佳预测。

在这项工作中，我们提出了一种新型的MRF模型MRFSeq，将RNA-Seq读数与先前的基因共表达网络知识结合起来,来推断DE基因。

 

不同于（Wei and Li，2007）中的MRF模型，我们仔细选择了MRF模型中这个集团聚函数,, DE基因的MAP估计可以降低到众所周知的流量网络方面最大流量问题水平。根据Kolmogorov和Zabih的工作（Kolmogorov和Zabih，2004）。

由于最大流量问题是多项式时间可解的，我们的MRF模型可以用多项式时间精确求解。此外，我们介绍一个loopy置信传播方法（Mooij，2007; Weiss，2000）来计算推断DE或EE基因的置信度。

 

我们广泛的模拟实验和真实的RNA-Seq数据表明,,,,,MRFSeq通过在不损失精度的情况下获得相当高的灵敏度，实现了大幅改善的总体估计性能。
一个对预测结果的详细分析表明，通过MRFSeq预测的DE基因在不同的读取计数值上的分布更均匀,相比于的现有方法所获得的计数。

因此，MRFSeq可以帮助缓解DE基因对低读数计数基因的选择偏倚。

我们的分析进一步表明，大多数DE或EE基因可以单独由RNA-Seq数据正确预测，也可以通过MRFSeq正确预测，这意味着在差异分析结果中,使用基因共表达的先验知识不会引入新的偏差分析结果。

此外，我们比较MRFSeq与最近出版的转录水平方法，Cuffdiff 2（Trapnell et al。，2013）使用来自UCSC hg19的注释转录组的   真实的RNA-Seq数据（Meyer

等人，2013）。比较显示MRFSeq比Cuffdiff 2敏感得多。

本文的其余部分安排如下。

第2.1节定义了我们算法中使用的术语和符号，同时2.2节提供了MRF模型的构成和其团势函数。2.3节中显示了最大流量问题的减少。

实验结果在第3节中描述，其中还包含MRFSeq与现有差异表达分析方法（包括edgeR）DESeq（Anders和Huber，2010）， baySeq（Hardcastle和Kelly，2010）NOISeq（Tarazona等，2011）和     Cuffdiff 2（Trapnell等，2013）之间的比较（Robinson等，2010）。

 特别是第3.4节比较了低读数和低基因的方法的性能显示MRFSeq

不仅实现了显着的整体更高的准确性，但也提供了较少的偏见预测。

第4节给出了一些结论性意见。由于页面限制, 在正文中, 用于每个预测可信度的计算, Loopy置信传播方法以及一些数字和表格都被省略了，并在提供补充材料。

2 材料和方法

2.1术语和符号

假设G = {g1，g2，...，gn}是待测差异表达的基因

和X = {x1，x2，...，xn}二元随机变量，使每个Xi属于{0,1}表示基因gi的DE状态。

随机变量如果xi = 0，基因gi是一个DE基因; xi = 1，基因gi是一个EE基因。

假定两个随机变量xi和xj是相关的;两个基因gi和gj形成一对共表达基因。

配置x是随机变量X的0-1赋值。

假设在在感兴趣的两个条件A和B中,分别有p和q重复。

让读数计数aij和bij是读数的数量，与之对齐基因gi分别在条件A和B的第j个重复。

对于每个基因gi，两组读取计数    RA,i = {ai1，ai2，:::，aip}和

RB,i = {bi1，bi2，:::，biq}

所有读数映射到参考基因组后,从两个条件A和B的所有重复中重新总结出来。

对于观察到的读取差异计数, 流行的统计测量为错误的发现率[FDR，例如p值被修正, 多重测试（Benjamini and Hochberg，1995）]和先验概率。

目前的统计学方法通过独立检验来推断DE基因

对于每个基因，如果其读取计数的差异测量值超过一定的阈值（Oshlack等，2010）。在我们的方法中，DE基因是观察到读数与基因共表达的差异的最大化似然函数的配置决定，而不需要事先知道阈值。

MRFSeq使用，但不是仅限于, DESeq的FDR qi（Anders and Huber，2010）

读取计数RA,i和RB,i的差值测量，其中qi &[0,1]。为了提高我们算法的计算效率，通过将区间[0,1]划分为相同长度的间隔为0.05的20区间来进一步离散FDR qi(被进一步离散)。

让yi 2&{1，2……20}表示观察到的差异qi所属的间隔，并且

Y = {y1，y2，...，yn}是所有离散化的FDR的集合。

该给定其观察值Y隐藏变量X的联合概率,由MRF模型表示，一种能够捕获随机变量的统计依赖性的图形模型（Kindermann和Snell，1980）

在下一小节中有描述。

给定X对Y的联合概率，估计基因的DE状态实际上涉及两个推论问题。

首先是MAP估计问题，即搜索一个配置x使得Pr (x |Y）被最大化。

算法为MAP估计问题将在本节后面讨论。

其次是边际概率问题，即计算概率Pr(xi|Y）作为每个基因gi上配置的致信水平。该边际概率问题的loopy置信传播方法在补充材料中给出。

2.2 MRF模型

令H = (Vx，E）是一个无向图,表示G的共表达式网络，

使每个节点vxi  Vx与随机变量xi ￠X相关联以及节点vxi 共享的每个边(i,j)

并且编码两个相关的随机变量xi和xj的依赖性。

假定两个变量xi和xj被假定为相关的, 如果基因gi和gj是共表达的。

为了确定哪对基因是共表达基因，使用COXPREdb（Obayashi和Kinoshita，2011）中定义的相关系数ci,j作为两个基因gi和gj之间的基因共表达的测量。

如果ci,j大于阈值p，则认为两个基因是一对共表达基因。

我们在整个这项工作中使用p=0.5，因为它文献被广泛使用（Patil等，2011; Watson，2006）。在我们的模型中，我们认为每个基因的DE状态应该

依赖于其观察到的计数差异和其共表达基因的DE状态。

换句话说，我们可以假设每个随机变量在条件上独立于H中 由非相邻顶点索引的变量。因此，满足以下性质得到：

其中N(vxj)表示H中vxi的邻居。

由Hammersley-克利福德定理（Besag，1974）

给定Y的随机变量的X的联合分布可以被分解为一个集团潜在函数的形式

T， 在给定的集合中配置的正函数在图H使得

为了模拟共表达基因之间的成对依赖性，我们可以使用一个MRF模型，它只含潜在的函数大小最多2个.这种类型的MRF被称为成对MRF（Besag,1986），并将用于我们的工作。有两种潜在函数在我们的MRF模型中被使用。一个是一元函数    评分随机变量xi与其观察证据yi的兼容程度。另一个是成对的潜在函数   测量每对相关变量xi和xj之间的统计相关性。通过定义潜在函数，

给定Y的X的联合分布可写为

其中Z是标准化项以确保联合概率PrjXjY）总和为1.

令一元函数定义如下：

为了计算一元函数，在我们的MRF模型中, 应该给出两个先验概率之间的比例

应该作为已知参数。为了估计参数，首先合成10 000个DE基因和10000个EE基因的4个重复（每个条件2个）的读取计数。

我们对DE和EE基因读取计数的模拟, 遵循在DESeq模拟研究中使用的相同的步骤（Anders和Huber，2010）。

对于DE基因，在两个条件之间观察到的读取计数的log2倍数变化率是从均值为0和方差为0.7的正态分布随机抽取。

对于EE基因，平均值设为0，方差为0.2。

读取计数的仿真后，先前引入的离散化FDR被计算为观察到的合成读数的差异。。

假设有myi DE基因和nyi EE基因，在这个模拟中离散化的FDR是yi。

我们进一步假设xi成立的两个背景概率相等，即

Pr(xi= 1）= Pr(xi = 0）。

通过贝叶斯准则，先验概率的比率得到如下：

对称地，我们有

对于每对共表达基因gi和gj的成对函数 ，

在COXPREdb中定义的相关系数ci,j（Obayashi和Kinoshita，2011）被用作xi和xj之间统计相关性的度量。因此，成对的潜在函数定义如下：

这完成了X的联合分配的规范。为了方便起见, X的联合分布可以在方程（2）的两边取负对数改写为：

这里是一个常数，

          

每个  是一个单项，

每个  是一个两项

E(X Y)     伪能量函数

最大化联合概率Pr(X Y)的配置实际上是使伪能量函数E(X Y)配置最小化（Besag，1986

2.3 MAP估计

与启发式方法不同，迭代的条件模式用于近似Wei和Li（2007）的MRF模型的MAP，我们在本小节中表明，通过仔细设计MRFSeq中的潜在函数MRFSeq的MAP估计问题不再是NPHard问题，因为它可以减少到流量网络上的最大流量问题并且在多项式时间内得到最优解决。

如果一个随机变量xi被一个配置x所反转，

给xi状态赋值违反了其先验概率，

即

对于一个倒置的随机变量xi，

我们定义  代替 

 反演的代价。

如果xi和xj的分配状态是不同的,,即两个相关变量xi和xj被称为由配置x分开。

对于一对分离的变量xi和xj 而不是ci,j。

分离的成本是ci,j。 Kolmogorov和Zabih（2004）

证明当伪能量函数的E(X Y)的成对项是子模块，

即满足以下属性：

寻找使伪能量函数最小化的配置可以通过寻找最小化反转和分离总成本的配置来完成。也就是说，MRF模型上的MAP估计问题可以降低到在流程新网络H0上的最大流量（或最小切分问题）

这样, H0的最小切割对应于MRF模型的MAP估计，并且切割的饱和容量恰好是反转和分离的总成本。 很容易验证我们的配对项是子模块

 从我们的MRF模型（其图表示为H = Vx，E）到流动网络的减少

可按如下方式完成。

H0的节点是H节点和两个附加节点的联合，源节点s和宿节点t。

H的每个无向边(i,j)被转换成两个有向边(i,j)和伴随容量(i,j)。

对于每一个节点xi，两个有向边(s,i) (i,t)被加到E’

边(s,i)的容量是

 否则，边(s,i)容量是0。

 对称地

边(i,t)的容量是

否则，边(i,t)的容量为0。

 运行标准最大流量算法后，

例如Edmond和Karp算法（Edmonds and Karp，1972）

在流网络H’上获得最小割

其中, Vs是与s相邻的节点，Vt是与t相邻的节点。

它表示0-1分配，使得对应于Vs的节点的所有随机变量被分配1，

并且对应于Vt的节点的所有随机变量被分配0。

然后，如果xi为1，则将基因gi推断为DE基因，

xi为0，则将基因gi推断为EE基因

2.4 RNA-Seq数据集

两个公开可用的人类RNA-Seq数据集，即

MAQC数据集(Bullard等，2010; Shi等，2006）和

Griffth数据集（Griffith等，2010）将被用作基准数据集评估我们选择的差异基因表达分析方法的表现。

每个数据集都与一个额外的qRT-PCR数据集相关联，以验证基因的DE状态。

 MAQC数据集由两个样本组成，

 Ambion的人类大脑参考RNA（脑）和Stratagene的人类通用参考RNA。

 每个样品提供7个重复和总共4500万个长度为35bp的单端RNA-Seq读数。

MAQC数据集的读数从7100万个中唯一获得通过ReCounts（Frazee等人，2011）校准的映射读数。

格里菲斯的数据集是由DE的qRT-PCR验证或由ALEXA-Seq（Griffith等，2010）突出显示的表达基因制成的。

它包含96个和1.98亿个跨越两种人类结肠直肠癌细胞系的双末端读数

只有氟尿嘧啶的抗药性不同。

为了平衡两个样品中的测序深度，如（Tarazona等，2011），如读取的库大小设置为每个条件大约1亿次读取。

从SRA数据库           下载了MAQC数据集的原始RNA-Seq读断

(Leinonen等，2011），

而从ALEXASeq的FTP站点下载了Griffith(格里菲斯)的数据集的RNA-Seq读取

跨平台的基因关联是用BioMart（Zhang等人，2011）。

 在下游分析步骤中丢弃不匹配的基因。

为了获得Griffith数据集的读数，使用Tophat（Trapnell等，2009）将原始RNA-Seq读数与人类基因组UCSC hg19（Meyer等，2013）的高覆盖率组装进行比对，其中两个不匹配被允许和映射到多个位置的读取被删除。

最后，使用R软件包对格里菲斯数据集中每个基因的读数进行总结

Bioconductor的基因组特征和RSamtools以及来自Ensembl（版本60）的基因组注释信息（Flicek等，2011），只有外显子读数。

 为了公平比较，一个伪读数，1，应用于读数为零的所有基因，以避免一些统计计算中的零分问题。

2.5评估标准

根据Bullard等人的评估方法。 （2010年），

我们所有的实验结果都根据

精度（PRE）   PRE = TP /(TP + FP) *100％

灵敏度（SEN）SEN = TP /(TP+FP)  *100％，

TP是真阳性的数

FP是假阳性的数

FN是假阴性的数

为了结合这两个评估指标，F评分（FS）（van Rijsbergen，1979），定义为

 

FS =2 *(PRE x SEN)/(PRE t+SEN)*100%

被用作我们测试中预测方法总体性能的度量。

3实验结果

3.1差异基因表达分析方法的选择

将我们的方法与现有的基因差异分析方法进行比较[王闯1] , Tarazona等人提出了相同的选择标准。 （2011年）。

但是，Fisher的精确测试（Fisher，1922年）（Tarazona等，2011）进行了比较，在这里被排除在外，因为它的性能显示远低于其他方法的表现。

 最后，有四种方法包括edgeR（Robinson等，2010），

DESeq（Anders和Huber，2010），

baySeq（Hardcastle和Kelly，2010）和

选择NOISeq（Tarazona等，2011）在我们的测试中进行比较。

 请注意，NOISeq有两个版本，NOISeq_real和NOISeq_sim，

版本NOISeq_real用于我们的实验,  因为在我们的仿真和真实数据集的重复次数总是大于1个。

 在这些方法中需要一些合理的截断值（MRFSeq除外）来决定统计差异测量的显著性。

为了获得可比较的性能分析方案, 我们的实验采用了文献中截断值。

 更具体地说，DESeq中选择的FDR 0.1用于DESeq和edgeR。

对于NOISeq和baySeq, 我们分别选择0.8和0.999的概率，正如在Tarazona等人的工作中所做的那样。 （2011）。

 在两个难度等级下进行实验以比较我们的方法MRFSeq与其他选择的方法。

在第一级，基准数据集的所有读取计数都是根据DESeq仿真研究中假设的相同分布生成的。

在第二级，所有读取计数的基因是从两个真实数据集中 MAQC和格里菲斯数据集累积的，并可能包含低读取计数。

除了与基因水平方法的比较之外，MRFSeq也与最近公布的转录水平方法  Cuffdiff 2- 在两个RNA-Seq数据集上进行了比较。

3.2模拟研究

3.2.1模拟实验

 我们的仿真实验遵循（Wei和Li，2007）的框架。

下载了与MSigDB中186个KEGG途径相关的所有基因集（Subramanianetal。，2005）。

 然后使用COXPREdb（Obayashi。）定义186个基因组的共表达网络 (Kinoshita，2011），并且他们形成了186个无向图。

如果其共表达网络中的边数少于节点数，则丢弃基因集。

过滤后，保存由2194个不同基因组成的37个基因组。

37个共表达网络（参见补充表格S1）合并为由2194个节点和8512个边缘组成的全局网络。

 所有的方法在真正的DE基因的五种不同丰度水平下进行测试。

 绩效评估分为五类，其中每类代表10％的丰度水平间隔，使得五类覆盖DE基因的丰度水平，范围从0至50％（Wei和Li，2007）。

 在五个级别中的每一级中，我们随机选择10个路径的组合来形成这些集合真正的DE基因，同时将其余基因保持为真正的EE基因，使得真正的DE基因的百分比在该水平的范围内。

 通过遵循相同的步骤来模拟DESeq中使用的读取计数，随机获得10个基准数据集和读取计数的10个不同组合。

模拟读取计数范围从25到401。

 所有的方法都适用于50个基准数据集。

由于页面限制，总结第一区间[0,10）的绩效评估在表1中,所有五个时间间隔的完整评估在补充表S2中给出。

为了方便读者，精度灵敏度曲线也提供在补充图S1。

3.2.2比较模拟数据的方法

MRFSeq显然具有最好的F分数(即总体性能)和显着提高了其他方法的灵敏度。

其    F-分数是      14.2,   6.8,  6.8,  3.7和五个区间中第二高的 4.8％

而它对灵敏度的改进是13.6，10.6，9.1， 1.7和                    2％

 尽管baySeq在区间[30,40]和[40,50]中提供了相近的灵敏度分数，但它不能获得可比较的精确分数，因此总体表现较差。

虽然在区间[0,10）中, 灵敏度实现了相当大的提高 ，MRFSeq将精度提高至少6.9％。在另外四个区间中，MRFSeq的精度略低于其他方法。

MRFSeq的精确度与这些区间内的最佳精度之间的差异分别为2.3,1.4,2.2和1％，这实际上小于标准偏差。

MRFSeq的灵敏度和F分数的标准偏差大于其他方法的标准偏差。

这是因为MRFSeq的表现与共表达网络上真正的DE基因的拓扑分布有些相关。

由MRFSeq实现改进的数量可能取决于拓扑分布

 尽管如此，仿真结果表明共表达数据通过增加真实DE基因的覆盖面可以显著地帮助改善差异基因表达分析

3.3真实RNA-Seq数据的性能

3.3.1 MAQC数据集上的实验

除了之前的模拟研究之外，推断DE基因的性能在MAQC数据集上被评估。

 以前，塔拉索纳等人。 （2011年） 在不同的重复次数（或泳道）上测试选定的方法，从每个条件2次重复至7次重复，在MAQC数据集中查看排序深度

如何影响方法的性能。

 结果表明, 增加测序深度会降低(除NOISeq外的所有选定方法)的测序精度。

 为比较性能并理解MRFSeq的精确度会随着测序深度的增加而改变，在我们的工作中使用了由Tarazona等人设计的实验。

 考虑不同数量的重复，使得读取的文库大小在两个样本的每一个中从14变化到4,500万读数。

通过qRTPCR的归一化阈值循环值(CT)测量MAQC数据集中基因的表达水平。

为了验证MAQC数据集的真正DE基因，如果其CT值的log2倍变化率（LR）大于某个阈值b，则基因被定义为真正的DE基因，

 例如0.5或2，

 而如果一个基因的LR小于阈值a，那么它就是一个真正的EE基因，

例如0.2（Bullard等，2010）。

任何基因的LR在两个阈值a和b之间的,  被认为是边界基因。

在以前的研究中，所有的边缘基因都被丢弃了（Bullard等，2010; Tarazona等，2011）。

由于qRT-PCR的检测限制,低表达的基因在一些qRT-PCR测定中可能不存在。

如果被检测到的基因至少不能在qRT-PCR分析的四分之三中出现

那么它将在在至少一种qRT-PCR测定中也将被除去（Bullard等，2010）。

与以前的研究不同，我们不会抛弃那些边缘基因。

为了进一步测试低读数计数基因的推断能力,低表达基因也保存在我们的实验中。这给了我们总共836个基因。

如果LR大于阈值b，我们将基因定义为真正的DE基因。

否则 该基因是一个真正的EE基因。

 当阈值b设定为0.5时，有669个真正的DE基因，

当   b是     2.0时，有373个真正的DE基因。

836个基因的共同表达网络 形成836个节点和2426条边的图。

在DE基因的这两种不同的丰度水平（或LR值）上, 所有的方法都经过测试

预测结果再次综合评估精度，灵敏度和F分数在图Figure1中。

3.3.2比较MAQC数据集的性能

与仿真研究中的结果类似，在真实DE基因的丰度水平上, MRFSeq显着提高了灵敏度数和F-分数。

 对于b = 0.5和2时分别考虑的所有测序深度,灵敏度的提高至少9.2和8.8％。

 在获得最佳灵敏度分数的同时，除了表现出极高精度的NOISeq外，MRFSeq还可与其他方法的精度进行比较。

请注意，虽然NOISeq在所有方法中的精确度最高，但它的敏感性远低于其他方法的分数，其整体表现（以F分数衡量）受此影响。

随着序列深度的增加，NOISeq的精度仍然保持稳定，而所有其他方法失去一些精度。

对于b的两个值(0.5和2)，当重复次数从两次增加到七次时, DESeq的精度

分别下降4.0和4.7％。

当BaySeq的精度下降5.4％和6％，而edgeR失去2.0和2.6％同时，

MRFSeq精确度仅下跌1.6和2.8％。

以下事实可以解释MRFSeq的相对较小的精度损失

如果不以较强的致信度预测，许多假阳性可以通过MRFSeq方法使用共表达信息来消除。

因此，MAQC数据集上的这些结果表明，共表达信息不仅有助于获得更多真实DE基因的覆盖范围, 而且随着测序深度的增加也有助于相对保持精度稳定。

 此外，在差异基因表达分析方面,它可以帮助减少我们对深度覆盖的RNA-Seq数据的依赖。

3.3.3考虑置信度分数

就像DESeq和edgeR的FDR或baySeq和NOISeq的先验概率一样，

MRFSeq估计每个预测的DE基因的置信度（即边际概率）而且 可以应用置信度阈值来选择输出的DE基因（而不是遵循MAP估计算法）。

当致信度被应用不同的阈值时,我们对MAQC数据集上MRFSeq的性能感兴趣。

为了计算致信度分数，loopy置信传播算法在MAQC数据集中的所有七个重复中运行。

为了比较MRFSeq与其他方法的性能，如（Tarazona等人，2011）中所做的那样，每个选定的方法描绘了每个点代表一定阈值下的精度和灵敏度的精密度灵敏度曲线。

为了描述DESeq和edgeR的精度-灵敏度曲线，选择了从10 6到1的FDR阈值范围。

请注意，此FDR截止范围涵盖了在实际中使用所有阈值，这些FDR阈值之间产生灵敏度值在45%和100％之间。

对于其他不使用FDR的方法, 在相同的范围内（即45-100％）导致灵敏度的等效阈值被用于绘制精确灵敏度曲线。

精度灵敏度曲线在补充图S2表明，一般来说，MRFSeq比其他方法提供更准确的置信度分数。

注意，与MAP估计算法不同，使用由loopy置信传播算法获得的边际概率来推断DE基因需要额外的知识来选择适当的置信度（边际概率）阈值。

此外，loopy置信传播算法是一种启发式算法，因此不能保证正确的边际概率。

因此，MRFSeq使用MAP估计来选择DE基因和loopy置信传播算法仅用于估计每个选定的DE基因的置信度分数。

3.3.4比较格里菲斯数据集上的表现

 格里菲斯的qRT-PCR数据数据集由94个蛋白质编码基因的193个外显子测定组成。

 与MAQC数据集的LR不同，应用双尾t检验从Griffith数据集的qRT-PCR数据中鉴定真正的DE基因。

如果t检验小于0.05，则认为t检验的P值是显着的（Griffith等，2010）。

在此标准下，83个真正的DE基因和11个真正的EE基因被确定和被使用在测试选定的方法。

从COXPREdb中提取94个基因的共表达网络形成94个节点和25条边的图。

格里菲斯数据集上方法的性能见附表S3。

MAQC数据上, MRFSeq仍然有最好的整体表现，尽管其相对其他方法的改进不那么显著。

请参见补充材料进行详细讨论。

3.4低读数计数基因的表现

3.4.1读数低的基因

为了理解这些方法对具有不同读取计数水平的基因如何执行，对真实数据集的预测进一步分析。

 数据集中的基因分为两类，读数低的基因和读数正常的基因。

在（Bullard等，2010）中，如果基因在两个条件的每个重复中具有少于10个读数，则称该基因具有低读数。

否则，该基因据说有一个较好的读数。

由于格里菲斯的数据集只包含具有较好读数的基因，因此我们下面仅仅考虑MAQC数据集。

在MAQC数据集中的836个基因中，有453个基因具有低读断数，383个基因具有较好的读断数。

baySeq和NOISeq方法为标准化基因长度提供了额外的选项。

这两种具有归一化基因长度的方法分别表示为baySeqlen和NOISeqlen。

为了进一步研究低读断数基因的标准化效果，baySeqlen和NOISeqlen也适用于MAQC数据集。

通过为LR值选择阈值b = 0.5,表2比较了不同方法低读数的基因的预测结果。

补充表S4列出了(低读数或较好读数的)真实的和预测的DE基因的详细数量。

3.4.2低读断数的基因有显著的改善

 对于读数低的基因，MRFSeq的灵敏度为38.8％，而次优灵敏度仅为11.6％。

同样，MRFSeq的F-分数达到53.3％，而第二个最好的F-分数只有20.4％。

除了这些重要的改进，

 MRFSeq的预测表明读数低的基因与较好读数的基因之间具有更平衡的模式。

MRFSeq的RTP1|h是43.1％，而其RPPl|h为42.7％。

（由edgeR获得）第二好的RTP1| h和RPP1|h仅分别为13.0和13.3％。

这一结果表明，所有其他方法都针对低读数的基因。

他们预测的大部分DE基因来自具有较好读断数的基因。

对读数低的基因的基因长度应用规范化之后，baySeq和NOISeq的性能略有改善。

然而，长度标准化并不能真正改善低读数基因的整体表现，或纠正选择偏倚。

3.5与Cuffdiff2进行比较

与使用原始读取计数的基因级方法不同，Cuffdiff 2需要将读数映射到给定的基因转录物作为输入来调用差异基因表达（Trapnell等人，2013）。

 为了评估Cuffdiff 2在MAQC和Griffith数据集上的表现,(Trapnell等，2013）中使用Tophat

将两个真实数据集的RNA-Seq读数映射到注释的转录组UCSC hg19。

FDR值的阈值为0.1被用于称为Cufflink 2的DE基因。

两个数据集的MRFSeq和Cuffdiff 2的预测精度总结在补充表S5中，

分别用MAQC和Griffith数据集的LR阈值b = 2和截断P值为0.05。

精度-灵敏度曲线也在补充图S3中说明。

该表显示MRFSeq通过实现更高的灵敏度而显著地具有更好的F分数(并且因此整体性能),而当Cuffdiff 2获得一个更好的精度。

当我们考虑完整的FDR频谱或将FDR值限制在0.1以下时，

精确灵敏度曲线也表明MRFSeq比Cuffdiff 2具有更好的整体性能。

详细分析表明，Cuffdiff 2预测的真实DE基因的LR值比MRFSeq小。

在MAQC数据集中，有290个真正的DE基因，LR值从0.5到2。

MRFSeq的预测涵盖了290个基因中的171个， 而Cuffdiff 2只能检测到140个真正的DE基因。

 在格里菲斯的数据集中，9个真正的DE基因与P值相关，P值可以测量LR值从0.005到0.001之间差异的显著性。

尽管通过MRFSeq所有9个真正的DE基因都是可以预测的，但Cuffdiff 2只能识别4个DE基因。这个结果与（Trapnell等，2013）的讨论一致。

 一般来说，Cuffdiff 2可能会报告较低的LR率的DE基因较少

 因为它由于片段数量的不确定性，在表达中控制了变量 。

3.6 DESeq和MRFSeq预测的一致性

 如果一个基因仅通过使用RNA-Seq数据正确地预测其DE状态，但通过MRFSeq错误预测，那么该基因被定义为不正确的倒转。

虽然我们上述结果表明，使用基因共表达的先验知识显著提高了差异基因表达分析的整体准确性，并有助于缓解对低读数的基因的偏倚，

 但它提出了这个问题，如果先前的知识可能会引入一些新的预测偏差。

 在本小节中，我们估计了MRFSeq预测中错误倒转基因的数量与一种流行的RNA-Seq根本方法DESeq的预测中的错误倒转基因的数量并

分析了在共表达网络中更可能错误反转的基因类型。

比较DESeq和MRFSeq在上述仿真和实际数据集上的详细预测结果。

在40个模拟基准数据集中，通过DESeq正确预测的73 619个基因中，只有3092个（4.2％）被MRFSeq错误地反转。

在MAQC和Griffith的数据集中，被DESeq正确预测，被MRFSeq错误地反转的基因分别只有16个（3.5％）和0（0％）。

一般来说，大部分被DESeq正确预测的基因，在MRFSeq预测中保持正确。

此外，我们观察到错误倒转的基因，相比其他基因，在基因共表达网络中具有更高的边缘度。

在补充图FigureS4中显示了在基因共表达网络中，所有基因平均边缘的比较和不正确反转基因的平均边缘的比较。

 边缘度差异的意义在于通过使用单尾t检验证实（Goulden，1956）在仿真和MAQC数据集上的t检验的P值分别为5.1 x10 和5.1 x10 。

 然而，由于基因共表达网络通常具有众所周知的无标度特性，因此只有少量基因具有较高的边缘度（Carlson等，2006; Stuart等，2003）。

该性质应该限制不正确倒转的基因的数目，因此大多数仅基于RNA-seq数据正确预测的DE或EE基因(通过例如DESeq)在MRFSeq的结果中得到良好保留。 

4结论和讨论

在这项工作中，我们提出了一种新的统计方法MRFSeq，它将RNA-Seq数据和共表达信息结合起来，并有效地获得DE / EE基因的MAP估计。

在补充材料1.3节中提供了我们的图形模型的所带来的改进好处的讨论。

 使用仿真和真实数据 的大量的实验，我们已经表明MRFSeq能够利用共表达信息和这个额外的信息可以帮助提供更准确的和更小偏差的差异基因表达分析。

显然，我们改进的性能（特别是对读取次数低的基因）的关键取决于是否存在一个高质量的基因共表达网络（见补充材料1.4节的讨论）。

最后，MRFSeq使用DESeq的DE分析结果。

去研究当其他工具被用于DE分析结果时，MRFSeq如何执行，这将是有趣的。

补充材料第1.5节，比较了MRISeq的性能与NOISeq合并的变体的性能。

我们计划使MRFSeq具有灵活性

所以在不久的将来它可以与任何DE分析工具结合使用。

我们的实验使用了COXPREdb（Obayashi和Kinoshita，2011），它由七种模型生物的共表达数据组成。

人们也可以考虑使用共表达数据的其他来源

如ACT（Manfield等，2006），         ATTED-II（Obayashi等，2007），

CSB.DB（Steinhauser等，2004），    CoP（Ogata等，2010）等。

 或从公开可用的表达数据构建定制基因共表达网络等

作为GEO（http://www.ncbi.nlm.nih.gov/geo/），

ENCODE（http://encodeproject.org/ENCODE/），

modENCODE（http:// www.modencode.org/）等，

 特别是对于COXPREdb未覆盖的生物体（或组织）。

此外，阈值p用于从给定的基因共表达网络中 提取共表达基因对，也可能会影响我们算法的性能。

根据文献（Patil等，2011; Watson，2006）和MAQC数据的一些初步测试，我们在我们的实验中经验地设定p = 0.5。

显然，较高的p可能会降低MRFSeq的灵敏度，

而  较低的p可能会降低MRFSeq的精度。

我们计划探索不同的共表达网络（包括p的选择）对MRFSeq性能的影响，并研究在未来工作中选择最优p的自动方法。

致谢

作者非常感谢匿名的审稿人建设性的意见。

 该研究得到了由加利福尼亚大学河滨分校，总统讲座教授基金， 国家科学基金会的资助的部分支持