用RNA-Seq和微阵列评估C57BL/6J和DBA/2J小鼠纹状体的基因表达

摘要

介绍

桑德伯格等人。[1]似乎是第一个使用微阵列检测两种近交系小鼠（C57BL/6[B6]和129SvEv）之间差异脑基因表达的小鼠。使用1.8倍的变化阈值，这些作者发现检测到的约1%的转录物存在差异表达。随后的研究检查了近交系小鼠或大鼠大样本中大脑基因的差异表达[2]，[3]，[4]，[5]以及隔离啮齿动物种群[6]，[7]，[8].重组自交系和F₂群体已被用来生成表达QTL（eQTL）图谱，而这些图谱又常常与行为数量性状位点（bQTL[4]，[5]，[9]微阵列技术和分析技术的改进使得可以非常准确地测量大脑基因表达的变化；重要的是，累积数据记录表明，自交系或来自自交系的群体之间的表达差异最大，例如选择性繁殖系[10]，实际上相当小（15%到30%）。在某种程度上，这些微小的变化反映了这样一个事实，即由于探针饱和，寡核苷酸阵列的杂交等温线往往不是线性的[11].由于选择性剪接的高度，这种效应在大脑中可能会更加明显[12]寡核苷酸阵列的第二个问题是SNP的影响，这是近交系或近交系衍生群体特有的问题[8]，[13]啮齿动物寡核苷酸阵列基于B6小鼠或Brown-Norway（BN）大鼠菌株的序列。即使是与B6或BN密切相关的近交系也可能存在数百万SNP的差异[14]，这反过来会导致显著的杂交伪影[13]。屏蔽SNP可以改善这种情况，但在某些情况下会导致从分析中删除探针（Illumina）或探针（Affymetrix）。寡核苷酸阵列的第三个问题是与预定义报告者/探针相关的注释和总结问题。最后，由于微阵列主要询问3′非翻译区（3′-UTR），因此它们提供的关于选择性剪接的信息相对较少。Affymetrix 1.0 Exon ST阵列收集关于选择性剪接的数据，但当用于检测差异选择性剪接时，由于每个探针的探针数量较少，因此对“SNP效应”特别敏感。

下一代测序（NGS）的出现和RNA-seq应用为检测差异基因表达提供了一个明确的替代寡核苷酸阵列，并有效地解决了上述问题。Mortazavi等人。[15]比较了短读Illumina/Solexa 1G平台产生的数据与Affymetrix 430 2.0阵列产生的数据之间的小鼠肝脏基因表达（仅B6菌株）。比较的短读指标是外显子模型每千碱基每百万映射读数（RPKM）。在基因水平上，这两个平台相当一致。然而，RNA-Seq也产生了新的数据。虽然90%的唯一定位阅读位于已知外显子内，但附加数据表明“新的和修订的基因模型，包括改变或附加的启动子、外显子和3′-UTR以及新的候选microRNA前体。”Morazavi等人。[15]他们还能够检测到超过100000个剪接位点，并注意到约3500个基因表达一个或多个交替的内部剪接。随后对小鼠和其他物种（如酵母）的研究证实，RNA-seq和微阵列在基因表达水平上产生了类似的数据[16]，[17]，[18]据我们所知，还没有研究使用RNA-seq来检测近交系之间的差异表达。然而，考虑到类似的应用，Bullard等人。[19]检查了最适合检测差异基因表达的各种标准化范式；他们得出结论，上四分位归一化引入的偏差最小。

目前的研究将RNA-Seq与两个微阵列平台（Affymetrix MOE 430 2.0和Illumina MouseRef-8 v2.0）进行了比较，以检测B6和DBA/2J（D2）近交系小鼠品系之间的差异纹状体基因表达。比较B6和D2菌株有几个原因。首先，B6是小鼠基因组的参考菌株，D2菌株是桑格小鼠基因组项目中正在进行完全测序的几个菌株之一。此外，还有大量的基因表达数据集对B6和D2菌株以及重组自交系和F菌株进行了比较₂来自这些菌株的杂交；这些数据集中包括来自纹状体的数据[8]，[20]，[21]最后，有广泛的表型数据比较B6和D2菌株[见22]包括行为、解剖和行为纹状体相关表型的数据[2]，[23]，[24]我们表明，通过使用严格的数据处理要求，RNA-Seq确定的差异表达与Affymetrix和Illumina平台在更多情况下一致，而不是仅与一个平台一致，褶皱变化方向冲突的情况很少。RNA-Seq的大动态范围检测到的基因比在微阵列上看到的多数千个。通过使用该技术获得的附加信息说明了RNA-Seq的价值。

方法

结果

RNA-Seq号

我们从10只B6和11只D2小鼠（三个Illumina GAIIx流细胞上的21条通道）中产生了单端RNA-Seq读数。每个样本的读取情况的高级摘要如所示表S1，包括总读取数和唯一映射读取数。为了便于总结与已知基因注释相关的RNA-Seq数据，我们首先根据Ensembl构建59个外显子注释组装一个基因级表示，如方法这产生了由Ensembl 59注释的36229个基因的汇总读取计数。其中12632个基因在所有车道上都没有读取，因此没有进一步考虑。另外7414个基因对至少一个B6和一个D2样本没有读取，被认为是不可靠的。如图所示图S1与其余基因相比，被移除的基因的读取计数往往很低。我们随后的分析侧重于剩余的16183个基因。

尽管读取计数的分布在泳道和流通池之间是可变的，但这些差异可以使用上四分位数缩放程序进行归一化[19]( 图1 )而不采用更极端的归一化程序，如分位数归一化[19]，[38]此外，上四分位缩放( 图1 )与仅基于总读取计数缩放（类似于RPKM度量）的分布相比，导致了同质分布[15](图S2).

在单独的窗口中打开

图1

每条车道/流单元每个基因的读取计数分布。

A）所示为每个基因的非标准化计数分布。然而，计数被一个常量除，该常量在此定义为以兆字节（MB）为单位的总唯一读取计数的平均值。这样做是为了便于比较箱线图。对于B），计数值通过上四分位校正总计数（MB）进行标准化。使用edgeR计算上四分位校正值[42]此外，这两个图都以log2比例显示。还要注意的是，X轴的标记形式为traint_lane_flowcell。

在我们手中，单因素泊松模型与之前提出的似然比检验相结合[19]，产生了倾向于反保守的结果。最简单的解释是，我们的生物复制品之间的变异性比固定分散泊松模型所能解释的要大，这种可能性以前已经提到过了[18]，[19]其他广义线性模型可用于解释过度分散，例如与拟似然拟合的泊松族之一[39]或负二项族之一[40]然而，由于我们的设计相对简单，我们决定使用精确的测试[41]用于edgeR Bioconductor包中实现的负二项分布[42].

在估计了常见的色散参数并应用精确测试后，错误发现率（FDR）控制[35]应用程序。我们发现1727个差异表达基因（DE）（q值<0.01），其中958个基因在B6中相对于D2更高表达，769个基因在B中相对于D22更低表达（分别为B6/D2>0或<0的对数倍变化）。在中演示数据集S1（附带说明文本S1)是基因、它们的折叠变化、平均读取计数的度量以及指示特别低表达（低读取计数）基因的标志。根据基因是否被确定为差异表达，对总结合外显子长度的分布进行了检查，结果表明，似乎没有与之前看到的情况相反的明显基因长度偏差[19]，[43]对每个四分位基因长度的DE基因计数的检查也没有显示任何基因长度偏差(图S3)。相反，如果我们检查为单因素泊松模型生成的相同图，其中包含偏移量中的上四分位缩放程序读取计数调整，我们可以看到证据表明较长的基因区域倾向于差异表达(图S4)这表明统计模型的选择可能是基因长度对差异表达影响的一个重要因素。

除了基因长度之外，我们还评估了基因组结构对结果的潜在影响。如Walter等人。[13]由于等位基因序列变异导致杂交效率的差异，微阵列分析可以检测假阴性或假阳性差异表达。RNA-Seq中可能会出现类似的现象，因为应用了启发式算法来加速重新校准算法。例如，SNP和碱基调用错误的组合可能导致读取映射错误或根本没有映射，从而可能导致差异表达基因中的菌株特异性偏差。为了解决这一问题，我们计算了16183个基因中每个结合外显子每千碱基的SNP数量，并评估了差异表达以证明方向性偏差。如中所示图S5就已知SNP的密度而言，D2>B6差异表达调用的分布相对于B6>D2调用向下移动，尽管分布在很大程度上重叠。为了进一步解决这个问题，我们合成了一个D2参考序列，并将D2读数重新对齐。我们发现，唯一映射读取数略有增加（每个D2样本的读取数增加不到0.2%），而多映射读取数增加最小（每个D2-样本的读取量增加不到0.1%）。此外，在总结、归一化和应用edgeR程序后，我们发现958个最初鉴定为DE且B6>D2的基因中有934个（97.5%）仍然是DE，但其余24个基因的状态发生了变化，这些基因不再被检测为DE，表明它们最初可能是假阳性。最初没有DE的36个基因的状态也发生了变化；根据重新校准的读数，这些被检测为D2>B6的DE，表明等位基因偏差可能导致假阴性结果。最后，最初检测到的5个基因的状态发生了变化，即DE D2>B6，未发现差异表达。总的来说，相对于RNA-Seq发现的1727个差异表达基因，假定的假阳性和假阴性差异表达调用的数量是最小的。改变状态的基因显示在数据集S1.

RNA-seq和微阵列结果的比较

首先，我们比较了三个平台检测到的基因重叠。7211个基因对所有三个平台都是通用的( 图2 )。有202个基因是Affymetrix微阵列所特有的，529个基因是Illumina微阵列所独有的，4179个基因仅由RNA-Seq检测到( 图2 )。这些基因中的大多数是蛋白质编码的，尽管每个平台也检测到其他Ensemble基因类别中的一些基因，如数据集（S1，第2页，第3章).

在单独的窗口中打开

图2

纹状体表达的小鼠基因检测平台重叠。

对于RNA-Seq：当RNA-Seq数据与Ensembl中注释的所有基因对齐时，通过过滤器进行检测的基因。微阵列：筛选后检测到的基因至少有一个探针（Affymetrix）或探针（Illumina）。

接下来，我们比较了不同平台之间的差异表达。RNA-Seq共检测到1727个DE基因。RNA-Seq DE基因的最大子集（591个基因）仅在RNA-Seq-DE分析中检测到表达，因此其DE的证据也仅限于RNA-Seq数据。使用Affymetrix和Illumina微阵列检测到369和212个RNA-Seq差异表达的基因，但在这些微阵列分析中，这些基因中只有51%和31%是DE。 图3 总结了555个RNA-Seq DE基因的最终子集，这些基因也通过至少一个探针或探针在两个微阵列上检测到。其中222个基因是基于RNA-Seq的DE，而两个微阵列平台都检测到这些基因是“存在”的，但没有检测到DE。这表明差异表达可能不会延伸到3′UTR，这是由微阵列优先询问的。这些基因中有144个是使用所有三种平台的DE，表明可能在整个基因中存在高置信度差异表达。

在单独的窗口中打开

图3

RNA-Seq检测到差异表达（DE）平台重叠。

根据RNA-Seq数据，1524个DE基因中，555个也是检测Affymetrix MOE 430 2.0和Illumina MouseRef-8 v2.0。当在所有三个平台上进行查询时，144个在所有三种平台上显示DE（5个在所有这三个平台中显示DE，但在折叠变化方向上不一致），表明可能在所有外显子中都是DE的基因。125和54分别是Affymetrix和RNA-Seq上的DE或Illumina和RNA-Seq（其中3和2分别在折叠变化方向上不一致）。虽然在所有三个平台上都检测到，但222只在RNA-Seq上显示DE。

在所比较的555个( 图3 )，只有13个被标记为读取计数低。此外，与在至少一个阵列平台上发现差异表达的222个相比，仅在RNA-Seq上发现的平均对数转换计数略高（p=0.015；双尾吨-这表明许多仅在RNA-Seq分析中发现的DE基因可能确实为真阳性。

讨论

目前的研究使用了一种基因水平的总结方法，该方法已被证明提供了比单一碱基位置更可靠的值[44]或更小的间隔，因为序列位置之间观察到的不均匀性可能归因于随机六聚体启动[45]。与之前的报告一致，我们发现通过总读取计数调整读取计数差异似乎是不够的，但通过每条车道的上四分位进行额外的修正，可以提供足够的方法来规范化数据[19]。虽然使用泊松模型发现了相对较大的影响，但使用edgeR软件包，我们发现几乎没有证据表明存在基因长度偏差。我们表明，仅由于SNP的存在，在近交系小鼠品系的比较中有可能出现假阳性和假阴性，尽管对这种比较的影响并不严重。我们注意到，可能存在的任何偏差都是SNP密度的函数。对于更复杂和SNP密度更高的小鼠模型，例如异种库存[7]，[46]，SNP可能对某些基因的DE有更严重的影响。我们进一步注意到，这两个菌株之间的插入和删除（我们没有考虑）可能会导致类似的影响。在Sanger正在进行的普通近交系的测序基因组发布后，基因组结构对差异表达分析的影响可能不再是几个近交系小鼠的问题[14]，允许根据各自的参考序列和转录模型进行分析。

与Illumina或Affymetrix微阵列平台相比，RNA-Seq提供了更多“检测到”的基因，这些基因具有可靠的信号，不受SNP或注释问题的影响。我们认识到在分析描述RNA-Seq的背景水平方面存在困难，并建议今后进行其他实验，精确测量所选模型系统中计数范围相对较广的基因的mRNA丰度，类似于早期研究中所做的工作[15]。尽管平台存在差异，但我们发现，总体而言，基于折叠变化方向和显著性值，微阵列和RNA-Seq的一致性相对较好，58%的RNA-Seq-DE基因在所有平台上被检测，确定在相同方向上差异表达，错误发现率为0.01。通过RNA-Seq和至少一个微阵列平台发现这些基因中最大比例的DE在所有三个平台上都是DE。有趣的是，Affymetrix MOE430 2.0比Illumina更符合RNA-Seq数据。Affymetrix和测序数据之间的高度相关性在SAGE样数字基因表达（DGE）的背景下已经观察到[47].

所有三个平台询问的RNA-Seq-DE基因中比例最大的是那些仅在RNA-Seq中差异表达的基因。我们发现，相对于那些与至少一个微阵列平台一致的基因，这些基因往往具有更大的平均读取计数，这进一步表明RNA-Seq在这两个微阵列平台上的效用。与两种微阵列技术共同差异表达的144个基因可能代表注释基因中相对同质表达的实例，因为来自这两个微阵列平台的探针应倾向于偏向基因的3′UTR端。在这方面，仅在差异表达方面与一个平台一致或仅在RNA-Seq中发现的基因可能代表转录亚型丰度的差异。这些类别中每一种的差异都说明了RNA-Seq相对于微阵列的优势，在这种情况下，RNA-Seg计算整个基因的DE，而不仅仅是在基因内的单个探针（集合）位置。据估计，每个小鼠基因的平均替代转录本约为2.5个[48]图中有选择性剪接，探针位置明显影响微阵列数据的解释。

另一个潜在的跨平台敏感性生物学评估可以通过跨平台检测Y染色体基因来完成。然而，对本研究中用于分析的Ensembl注释的审查表明，只有20个基因注释位于Y染色体上。其中，如果我们检查通过SNP过滤器的探针、背景上的检测以及我们分析中使用的独特Ensemble映射，那么Illumina阵列上只剩下7个探针，询问Y染色体上的5个基因，Affymetrix阵列上没有探针。我们确实注意到，在Illumina阵列检测到的5个基因中，其中4个也被RNA-seq检测到（唯一的例外是Illuminia阵列上最低可接受信号水平）。然而，由于带注释的基因如此之少，而且阵列上的覆盖率如此之低，因此目前使用这种类型的分析进行评估并不能提供信息。

此外，目前存在许多针对RNA-Seq的实验方案选择，每种方案都有其自身的优势和对下游分析的影响。例如，为了去除高度丰富的rRNA分子，要么富集含有poly-A的序列，要么耗尽rRNA[49]，[50]在我们的RNA-Seq实验中，我们使用了poly-a富集程序，因为目前的rRNA耗竭策略效率较低[50]然而，使用不同的选择程序可能会影响与微阵列的一致性，这是未来研究的一个有趣的主题。

尽管目前已有许多分析RNA-Seq数据的方法，而且还有更多的方法在继续开发，但可以提出的最基本和最基本的问题是，一个基因/转录物是否在两组之间有差异表达。基因水平的差异表达为进一步的实验奠定了基础，这些实验旨在识别和量化样本之间的转录亚型[44]，[51]和从头开始未标记表达区域的识别[15]，[52].方法存在以推断转录亚型的选择性剪接和相对定量[44]，[51]，并且配对-end序列的生成应该允许对观察到的任何新的选择性剪接事件有更大的信心[53]。未来应继续使用该技术来评估小鼠菌株间选择性剪接的细节。

支持信息

致谢

脚注

工具书类