使用SAGE-Seq对人类乳腺组织样本进行基因表达谱分析

用于标记映射和排序误差最小化的管道

为了分析单个基因的表达，我们使用SeqMap(Jiang和Wong 2008年)并提出一个映射管道（补充图S3），将标签与RefSeq基因对齐。这个映射管道允许我们将标签映射到RefSeq基因的线粒体、正反义转录本。如果一个转录物有多个CATG，那么最靠近poly（a）尾部（3′端）的一个被称为最佳标记(图1A). 如果一个标记映射到多个RefSeq位置，并且只有一个标记是最佳标记，则最佳标记被视为唯一映射的位置。否则，该标记是一个非唯一标记，其计数在映射的位置之间平均分配。意义标签被定义为映射到已知转录物外显子意义链的标签。反义标签是指不能映射为义标签，但能够映射到已知转录基因的反义链的标签(He等人，2008年).

在单独的窗口中打开

图1。

SAGE-Seq标签对齐和排序错误最小化。(A类)最佳标记是指紧邻poly（A）尾部的最3′端NlaIII位点（CATG）的标记。(B类)根据补充图S3中的标签对齐管道，样本N1的标签对齐统计。其他数据集的详细映射如补充电子表格1所示。(C类)标签映射期间显示的假定排序错误。所有列出的标签都是唯一映射到同一RefSeq基因的最佳标签“{“类型”：“entrez-notide”，“属性”：{“文本”：“NM_001010”，“term_id”：“1653962309”，“term_text”：“NM_001010“}}NM_001010号“（并非列出了在同一位置唯一映射到此基因的每个标记。）X（X）-axis是按标签计数的降序列出的标签序列。最有可能由测序错误引起的序列中的单基差异用红色标记。此特定示例的测序错误最小化步骤按以下方式进行：将所有这些标记的计数相加，并将其分配给标记“CATGGCGTCCGCCTGCTA”，并删除所有其他标记。

映射结果还可用于识别排序错误，如所示图1C。我们将唯一映射到相同位置相同基因的标签计数合并，以减少因测序错误（测序错误最小化）而产生的噪音和采样偏差，从而减少后续差异表达基因分析中的假阳性数量。参考基因组中的标签用作共识标签，以最小化测序错误。例如，假设标签“GCCGTGTCCGCCTGCTA”出现190793次，它正好映射到参考基因组。如果有3198个标签通过单个碱基对与此标签不同，在最小化测序错误后组合在一起，则总数为193961个；因此，单碱基对失配的比例为1.6%（3198/193961）。这相当于每基0.1%的测序错误率（17×0.001×[1−0.001]¹⁶= 1.7%). 这与Illumina高质量读数的估计一致(Shendure和Ji 2008). 以库N1为例，我们证明了大约76%的标签可以使用我们的管道进行唯一映射；其中6%为线粒体标签，46%为唯一RefSeq感观最佳标签，14%为唯一感观非最佳标签，10%为唯一反义标签(图1B; 其他库的映射结果的补充电子表格1）。所有后续分析都是在46%的独特感觉最佳标签上进行的。

正常和癌症转录组概述

细胞群体的基因表达模式在许多方面类似于生态系统中不同物种的物种群体，在我们的研究中，一个物种的个体就像一个转录物。在典型的生态系统中，有些物种很丰富，而大多数物种很稀少(Magurran 2003年). 同样，SAGE-Seq分析数据显示，大多数基因表达水平较低（罕见转录物），少数基因表达量较大（大量转录物）(图2A、B). 有趣的是，尽管很少表达的标签占大多数，但在考虑其群体（表达水平）时，高度表达的独特标签仍然占主导地位。通过绘制标签计数占总标签计数的累计比例作为唯一标签计数的函数，我们发现，尽管只有一个计数的唯一标签占总标签数的63%S公司（唯一标签的总数），它们只占N个（标签总数）(图2B). 问题是，这些低计数标签是由测序错误控制的假标签，还是以极低水平表达的真标签。如上所述，在将序列错误最小化后，可以根据映射信息识别和纠正由于序列错误导致的一个不匹配的标签。两个以上不匹配的标签仅占所有读取的0.01%（请参阅方法）。因此，这些低计数标签不能用测序错误来解释，因为它们比这些错误所能解释的要丰富得多。它们可能是低丰度转录物和核酸污染（可能是基因组DNA）的混合物，因为Sage-Seq制备协议不包括从RNA样品中消除基因组DNA的步骤。

在单独的窗口中打开

图2。

独特标签计数频率图和非参数经验贝叶斯方法。(A类)库N1（黑色）和N5（红色）中唯一标记计数的频率。X（X）-axis是观察到的标签计数年-axis是显示具有特定计数的唯一标记数的频率。(B类)描述库N1中唯一标签分布的饼图：62.6%的唯一标签的标签计数1、12.1%、计数2、5.5%、计数3和19.8%的计数大于3。外部图显示了唯一标记计数的累计分数。尽管62.6%的唯一标签具有计数1，但它们仅占总标签计数的3%。(C类)标签比例散点图。X（X）-axis是通过随机抽样库N1的10%获得的伪库1中标签的比例。Y（Y）-轴是通过随机抽样1%库N1获得的伪库2的平均比例。数据点通过以下方式获得。例如，以1×10的比例查找伪库1中的所有标记⁻⁶，然后计算这些标记在伪库2中的平均比例，例如给出1×10⁻⁵。这给出了（1×10）处的数据点⁻⁶, 1 × 10⁻⁵). 虚线为年=x个黑色符号表示使用最大似然估计值的比例，其中低表达和中等表达的标签（<100/百万）被高估，高表达的标签被低估（>1000/百万）。红色符号表示使用非参数经验贝叶斯方法计算的比例，该方法在低丰度和高丰度标签中具有不同排序深度的两个库之间具有改进的、更具可比性的校正比例。

非参数经验贝叶斯归一化

如果每个SAGE库的测序深度相同（即相同N个)，不同库中的标记数可以直接比较。然而尽管大多数样品经过Illumina测序的一个泳道，N个不同图书馆的数量从100万到1300万不等(表1). 因此，为了准确比较不同文库的基因表达模式，需要对标签计数进行标准化。一种直观的标准化方法是使用比例第页，定义为无，其中n个是唯一标记的计数。被称为人口频率的最大似然（ML）估计量(费希尔1922). 这种方法的缺点是第页测序数据中缺失的任何标记（未检测到的标记）被赋值为零，而它高估了第页低丰度和中丰度标记的数量，并低估了高丰度标记(图2C，黑色符号）。此外，Illumina基因组分析仪的高通量高质量读数可以很好地估计不同富集水平下标签的群体频率(图2A). 作为先验应用的标签频率异质群体表明，对标签富集的最佳估计应该小于观察到的数量，因为低丰度转录物的群体大于丰度转录体的群体(图2A). 因此，需要一种更复杂的SAGE-Seq数据规范化方法。

我们应用了非参数经验贝叶斯（NEB）方法(好1953；罗宾斯1956；Orlitsky等人，2003年)规范化具有不同排序深度的库（请参见方法）。与ML相比，NEB有两个优点。首先，ML简单地认为未检测到的标签为零，NEB估计未检测到标签的比例为对₀=n个₁/N个，其中n个₁是计数为1的唯一标记的频率。验证NEB估算值的准确性对₀，我们以1%的步长从1%到10%随机采样库N1，以生成10个具有不同测序深度的伪库。我们使用NEB估算对₀将每个伪库中未检测到的标签与它们在原始库中的各自比例进行比较，发现它们非常一致（补充图S2）。其次，NEB根据观察到的计数和唯一标签计数频率分布的性质来调整标签计数(图2A)，作为减少序列抽样偏差的经验先验(盖尔和桑普森1995). 它还通过检测到的标签的估计总比例将调整后的比例重新规范化为1–对₀（有关NEB和ML标准化之间的比较，请参阅补充材料部分“差异表达基因的算法比较”）。为了显示采样偏差的影响，我们从库N1中随机抽取10%（伪库1）和1%（伪库2）标签，生成两个序列深度相差10倍的伪库。当比较两个伪库中标签的比例时，我们发现NEB标准化后的比例更具可比性，而ML高估了这一比例第页用于低计数标签和低估第页用于测序深度较低的伪库2中的高计数标记(图2C).

正常和癌症转录组的多样性

基于序列的基因表达谱的一个优点是它同时测量许多基因的绝对表达水平。因此，我们可以获得细胞内以及库之间转录多样性的全局视图。我们使用了两种不同的测量方法来比较我们分析的图书馆中的转录本多样性。首先，我们使用了辛普森多样性指数（SID）(辛普森1949)以表征每个文库中的转录组多样性。SID捕获标记计数分布的方差，并且与排序深度无关（请参阅方法）。较高的值表示较高的多样性，这意味着标签数量在不同基因之间的分布更广。我们发现，在我们的数据集中，一般来说，癌症样本库的多样性高于正常样本库(图3A、B; Wilcoxon秩和检验，对= 0.07284; 这个对-由于样本数量有限，值处于显著性边界）。这种趋势可能是因为肿瘤表达更多的基因，要么是因为它们由更多不同的细胞群组成，要么是由于它们失去了维持组织和细胞类型特异性基因表达模式的正常表观遗传控制(图4D).

在单独的窗口中打开

图3。

正常和癌症转录组的多样性。(A类)辛普森多样性指数用于衡量文库内基因表达多样性。与正常组相比，癌症组的库内多样性较高。(B类)方框图描绘了正常和癌症样本的辛普森多样性指数。对=0.07284（威尔科森秩和检验）。(C类)定义为“1–Morisita-Horn相似性指数”的距离用于测量跨库的基因表达多样性。正常组的图书馆彼此更相似，而癌症图书馆则更为多样化。(D类)使用中定义的“距离”进行分层聚类C类分离正常和癌症库。

在单独的窗口中打开

图4。

SAGE-Seq标记映射和测序深度饱和度曲线。(A–C)三个选定基因家族表达谱的差异覆盖率：转录因子(A类)、GPCR(B类)和ABC运输车(C类).Y（Y）-axis列出了基因和x个-axis是平均基因表达指数（标准化标签计数的对数）。红色和蓝色分别表示传统SAGE和SAGE-Seq。SAGE-Seq在这些基因家族中检测到的基因比传统SAGE更多。(D类)唯一的最佳标记基因数(年-轴）与测序深度相关(x个-轴）。最佳标记基因的数量是由最佳标记映射的唯一基因的数量，如果多个标记映射到同一基因的最佳标记，则计算为一个。黑色和红色分别表示正常组和癌症组。符号”○” 以及““分别标记传统SAGE和SAGE-Seq。实体曲线（饱和度曲线）是通过对正常（或癌症）组中所有库的组合进行采样而得到的模拟结果，这些库描述了随着测序深度增加的趋势。传统的SAGE鉴定的最佳标签基因比SAGE序列少得多。SAGE-Seq显示，癌症样本（红色三角形）比正常样本（黑色三角形）具有更多独特的最佳标记基因。传统SAGE无法检测到这种差异（红圈与黑圈）。

第二种多样性是跨文库测量的，以研究不同个体的文库之间的基因表达多样性。为了询问A和B两个库的相似程度，我们使用了Morisita-Horn（MH）相似性指数C类_MH公司(A类,B类) (沃尔达1983)（参见方法），并将其距离计算为D类≡ 1 −C类_MH公司(A类,B类). 莫里西塔-霍恩指数比其他距离测量方法有几个优点。首先，MH指数受以下因素影响不大N个和S公司，这对于确定测量的差异不是由于测序深度的差异至关重要。其次，与基于皮尔逊互相关的距离相比，对于标准差接近零的数据，MH指数没有奇异性。

我们发现，癌症样本不仅在每个个体中差异更大（辛普森指数），而且在不同个体中差异也更大（MH指数）(图3C). 这并不完全令人惊讶，因为正常细胞的生理作用在不同个体中基本相同，而肿瘤基因多样，在体内没有功能作用；因此，没有选择压力来将其表型保持在一定范围内。基于MH索引的库的层次聚类表明，癌症库之间的差异更大（库之间的距离更大），并且它们与正常库也有非常明显的区别(图3D).

SAGE-Seq与传统SAGE的数据质量比较

在读取比对和测序错误最小化之后，我们进一步评估了SAGE-Seq分析全基因组基因表达的能力，并将其与传统SAGE进行了比较。随着测序覆盖范围的加深，SAGE-Seq在同一组内不同文库之间的数据相关性更高（补充图S4）。此外，传统SAGE只能检测比例为10的基因⁻⁵到10⁻³而SAGE-Seq显示了更大的动态范围（定义为检测到的基因富集范围），涵盖了10个数量级⁻⁷到10⁻²例如，编码转录因子的基因通常以中或低水平表达，SAGE-Seq在我们的样本中检测到大约1300个转录因子的表达（在人类基因组中总共1658个转录因子中）。传统SAGE检测到的大多数转录因子也由SAGE-Seq检测到，而384个转录因子仅由SAGE-Seq检测(图4A). 我们观察到GPCR和ABC转运体基因编码的类似现象(图4B、C)已知其在正常细胞和癌细胞之间有差异表达，且表达水平相对较低(Li等人，2005年；2009年院长).

为了确定测序的SAGE-Seq文库离饱和有多远，我们计算了检测到的唯一最佳标签基因（唯一映射的最佳标签）的数量与每个文库的测序深度的关系。唯一映射的最佳标签的数量是检测到的基因数量的一个良好指标。更深入的测序有望检测到更多的基因，直到检测到所有基因时达到平稳期。为了克服缺乏涵盖广泛测序深度的数据，我们合并了所有癌症（或正常）文库，并分析计算了在不同测序深度检测到的唯一标签基因的数量，以获得饱和度曲线（图4D; 见分析计算方法）。对于低于300万读数的测序深度，检测到的基因数量随着测序深度（快速生长区）的增加而急剧增加。该速率继续以低于500万次读取的较慢速度（缓慢增长区域）增长，直到正常和癌症样本的读取量都达到1000万次左右为止(图4D). 这表明SAGE-Seq的理想测序深度应在1000万读以上，每个库的最小期望序列深度为500万。Sage-Seq数据点（三角形）都靠近或位于生长缓慢的区域，在那里大多数转录组都被测序。传统的SAGE数据点（圆圈）仍处于快速增长区域，由于测序深度较低，超过一半的转录组未被检测到。图4D还表明癌症（红三角）样本中表达的基因比正常（黑三角）样本多，这与我们之前的发现一致，即癌症样本在每个库内和库间具有更高的转录多样性。

采样噪声和生物变异性

为了检测两种条件下差异表达的基因，了解基因表达变异的来源很重要。SAGE-Seq或任何使用测序技术的技术的变异性的主要来源是取样变异；大多数分析传统SAGE数据的算法使用各种方法解决了这一问题(Velculescu等人，1995年；Cai等人，2004年). 基于序列的转录组分析可以建模为二项式采样过程，替换近似泊松分布，因为使用当前技术(Kharchenko等人，2008年)，测序转录物是测序器上加载的cDNA总量的极小一部分。如果同一个文库被多次测序，泊松模型表明特定基因标签数的方差等于其丰度。

当检查正常文库中基因的经验方差与其各自的归一化计数（由NEB归一化并缩放到相同的测序深度N个=100万），观察到的差异确实取决于其基因表达水平（计数）(图5A; 带斜率的红色虚线α_正面≈2.0（对数图中）。然而，如果我们将随机二项抽样的方差-平均斜率表示为α_{苏格兰皇家银行}，根据泊松分布预计为1.0(苏格兰皇家银行代表第页安多姆b条诺米亚的秒采样，蓝色虚线输入图5A)，我们观察到过度分散(α_正面≈ 2α_{苏格兰皇家银行}> α_{苏格兰皇家银行}=1），这意味着观测数据的过度可变性明显大于随机参考模型（本例中为泊松模型）中的预期可变性。这表明基因表达变异性对平均表达水平具有非线性依赖性，这表明在我们的数据集中，过度分散可能是除采样变异外，生物个体变异的结果（见方法）。我们还观察到，在正常组和癌症组中，管家基因的子集和唯一映射的最佳标签的子集之间存在过度分散（数据未显示）。

在单独的窗口中打开

图5。

差异表达基因及其变异。(A类)去除噪声和归一化后，七个标准库的均值-方差图。红色虚线是对数图中最佳的线性拟合。斜率给出指数α_正面≈1.9蓝色虚线是通过采样引入的均值-方差线。(B类)鉴别差异表达基因的流水线：（1）测序误差最小化：标签对齐后，在相同位置映射到相同基因的标签被组合在一起；（2） NEB用于规范不同序列深度的不同库；（3） 过滤以删除少于两个库中计数≥3/百万的标记，然后是日志₂转型；（4） SAM用于检测差异表达基因。(C类)检测到差异表达基因(顶部)和激活的途径(底部)在SAGE-Seq和传统SAGE中。SAGE-Seq以1%的FDR识别约4000个差异基因，而传统SAGE以更宽松的截止值（10%的FDR）识别<200个差异基因。在对=0.001，SAGE-Seq确定了99条在乳腺癌中显著激活的通路，而传统SAGE仅显示32条。仅由SAGE-Seq识别而被传统SAGE遗漏的80条通路均为乳腺癌相关通路。(D类)重叠比率（定义为重叠基因数除以顶部传统SAGE中的基因数x个差异表达基因的百分比，其中x个在0和1之间更改。黑色符号表示实际数据（SAGE-Seq与传统SAGE）。这表明SAGE-Seq和传统SAGE之间的顶级差异表达基因列表几乎没有重叠。红色符号表示模拟（SAGE-Seq与向下采样的SAGE-Seq）。抽样SAGE-Seq意味着从每个SAGE-Seq库中对50 k个标签进行二项抽样；50000是传统SAGE的典型测序深度。模拟证实了与实际数据得出的相同结论：与传统SAGE相比，SAGE-Seq给出了不同的顶级差异表达基因列表。更深入的测序表明，传统SAGE与SAGE-Seq相比识别出不同组的差异表达基因，从而证实了我们的结论，即传统SAGE缺乏足够的测序深度。

差异表达基因分析

转录组分析的主要应用之一是识别不同样本之间差异表达的基因。在标签对齐和测序错误最小化后，我们的分析管道用于识别差异表达基因(图5B)首先应用非参数经验贝叶斯方法作为归一化步骤，以减少采样偏差，并将不同库引入相同的排序深度(N个=100万；归一化测序深度对差异基因没有影响，这与文库的测序深度不同。）。归一化后，所有库中计数≥3/百万in≤2的标签被丢弃。这有效地删除了大量非信息性标签，这些标签要么包含离群值，要么计数太低，无法检测具有统计意义的差异表达式，并在后续分析中节省了计算时间和存储空间。

然后应用对数变换获得表达指数，并解耦观察到的方差和基因平均表达水平之间的相关性(图5A). 从数量上讲，我们库中观察到的方差与表达式级别的平方成正比。根据统计学中的delta方法，对数变换是稳定方差的正确变换（参见方法）。另一种转换是电弧（电弧），这也是一个类似对数的变换，但具有在零处没有奇异性的优点(Huber等人，2002年). 补充图S5显示，在对归一化计数应用以2为基数的对数转换后，对于中等丰度和高丰度标记，表达指数的方差几乎与其平均值无关。最后，将SAM（微阵列显著性分析）算法应用于两组样本中的表达指数，以识别差异表达的基因(图5B；Tusher等人，2001年). 我们也尝试了标准吨-测试并发现许多误报，这些误报是由于低估了经验标准差而导致的极端情况吨值。SAM算法稳定方差以减少误报。此步骤中可以使用其他统计测试，而不是使用SAM，如Robinson和Smyth的调节吨-测试或Baggerly’s吨_w个测试(Baggerly等人，2003年,2004；Lu等人，2005年；Robinson和Smyth 2007). 分析差异表达基因的另一种选择是使用过度分散的模型，如过度分散的logistic回归或过度分散的对数线性模型(Baggerly等人，2004年；Lu等人，2005年). 然而，这些基于模型的方法是否可以扩大到SAGE-Seq数据的更深测序深度，需要通过更多数据的系统分析进行验证。

我们比较了SAGE-Seq和传统SAGE中正常和癌症的差异表达基因列表。SAGE-Seq和传统SAGE分别检测10052和4953个最佳标记基因的表达（即存在）（补充电子表格2），99%（4904）重叠。我们使用问-Storey和Tibshirani（2003）传统的SAGE序列不够深，不允许类似的对-值或FDR截止值，如SAGE-Seq。SAGE-Seq在1%FDR时识别出约4000个差异表达的最佳标记基因，而传统SAGE在10%FDR下检测到的差异表达基因不到200个(图5C). 更深入的测序可以提高SAGE-Seq的统计能力，以检测更多差异表达的基因。为了比较两个差异表达基因列表，我们根据它们的顺序检查了基因的秩吨-分数。前10%的基因吨-分数（传统SAGE的495个基因和SAGE-Seq的1005个基因）被用作差异表达基因列表，以比较这两种方法。根据先前的研究，SAGE-Seq检测到所有26个已知在正常和乳腺癌样本中差异表达的基因，而传统SAGE仅检测到四个（补充表S1）。

令人惊讶的是，当比较两种方法识别出差异表达的前10%的基因之间的重叠时，我们只识别出54个基因。进一步分析证实，传统SAGE和SAGE-Seq检测到的最高差异表达基因差异很大(图5D; 黑色符号）。许多因素都可能导致这种差异，例如文库制备方案和样本的差异。除了这些因素外，我们还观察到，SAGE-Seq检测到的最高差异表达基因通常以中等或低水平表达（～100/百万；参见补充图S7），传统的SAGE要么完全无法检测到，要么标签计数太低（两个或三个），无法显示具有统计能力的差异表达。SAGE-Seq中这些差异表达的标签不太可能来自基于观察到的标签计数的测序错误。这些数据表明，SAGE-Seq测序深度的增加导致检测到一组不同的差异表达基因。为了证明这一点，我们采用了模拟方法，因为使用从原始乳腺组织中分离的细胞数量有限的特定细胞群，不允许从同一样本中生成SAGE-seq和传统SAGE库。我们选取了14个SAGE-Seq库，并对其进行了抽样（二项式抽样），以达到传统SAGE的序列深度水平（～50000）。这些模拟文库的顶部差异表达基因也与原始SAGE-Seq文库几乎没有重叠(图5D，红色符号）。

排序错误最小化

在每个碱基0.1%的错误率下，没有错误的标签数量为（1−0.1%）¹⁷=0.9831353，一个错误的标签总数为17×0.001×（1−0.001）¹⁶= 0.01673003. 因此，在至少两个碱基中出现错误的标签群体为：1−0.9831353−0.01673003=0.0001346471。

辛普森多样性指数

辛普森多样性测量(辛普森1949)（SID）定义为：

哪里n个_我是的计数我第个标签和N个是标记计数的总数。SID=1表示一个标签支配系统的所有标签计数，这意味着没有多样性（最高支配）。SID值越大，多样性越高（优势度越低）。仿真证实，SID不受测序深度的强烈影响。

Morisita-Horn相似指数

Morisita-Horn相似指数，C类_MH公司(A类,B类)，在两个库之间，A类和B类，定义为：

哪里第页_我(A类)和第页_我(B类)是基因（或标签）的比例我用于库A类和B类分别是。MH相似性指数与测序深度无关N个.

抽样方差和生物方差

高通量测序技术被建模为具有替换的二项式采样过程。混合不同样本的生物可变性，由于采样，引入了另一层可变性。定义第页_n个、带计数的标签比例n个作为，第页=n个/N个，其中N个是标记计数的总数。从我们的层次模型来看，平均比例基因丰度是基因真实丰度的无偏估计量。然而，秒²(x个_我)，标度计数的观测方差x个_我基因的我是个体真实生物方差和抽样方差的相加（分析证明的补充材料），

使用非参数经验贝叶斯校正的归一化

我们使用简单的Good-Turing估计器（SGT）实现经验贝叶斯(盖尔和桑普森1995). 假设mRNA中所有唯一标签的总数为秒、和第页_我是标签的真实比例我，这是我们想要从数据中估计的。观察到的标签计数的经验贝叶斯估计第页是：

哪里n个_第页是带有count的标记数第页。因此，每个代表的所有标记的预期总机会第页次(第页≥1）为：(第页+ 1)n个_第页+ 1/N个，其中N个是测序深度(N个=n个₁+2个n个₂+ 3n个_三+…). 因此，样本中所有标签的预期总概率为：（2n个₂+ 3n个_三+…)/N个= 1 −n个₁/N个在SGT中，未检测到的标签的比例，对₀估计为

哪里n个₁表示计数为1的唯一标记的数量（频率）。SGT后修正的标签总数，N个*，是N个* = Σn个_第页第页*.计数基因比例的经验Bayes估计第页,第页_第页*，被重新规范化N个*作为

不可见标签的良品期估计方差

的方差对₀可以用以下公式计算：Var(对₀)=变量(n个₁)/N个².n个₁= Σ_我无_我(1-第页_我)^{N个− 1}在二项抽样近似的假设下。引入新的随机数x个_我:x个_我=1，如果我第个标签在排序深度仅使用一个标签进行排序N个和x个_我=0否则。然后：

和

因此，

它给出了：

这表明未发现标签的Good-Turing估计器是一个稳定的估计器，如补充图S2所示。

饱和度曲线计算

我们将mRNA中所有独特标签的总数定义为秒和第页_我作为标签的真实比例我.第页_我是根据数据估计的。因此，根据二项式抽样，未发现标签的平均值为：

其中，总和是可能的唯一最佳标记和N个是排序深度。因此，检测到的唯一最佳标记基因的数量是唯一最佳标记的基因总数减去n个₀(N个)如所示图4D.

方差稳定

对于我们数据中观察到的方差-均值关系，稳定方差的正确转换是基于delta方法的对数转换。我们提供以下证据。假设一个随机变量x个遵循平均值分布μ和方差σ².考虑转换克(x个). 泰勒展开式克(x个)大约μ第一个订单是克(x个) ≈克(μ) + (x个负极μ)克′(μ). 因此，转换后的变量克(x个)具有近似平均值克(μ)和近似方差Var[克(x个)] ≈σ²[克′(μ)]²。在我们的数据中，μ和σ²满足观察到的依赖性σ²∼μ²，产生Var[克(x个)] ≈μ²[克′(μ)]².假设转型克稳定方差，Var[克(x个)]是独立于μ，因此克′(μ) =c（c）/μ，其中c（c）是一个常量。在以下方面进行整合μ给出了稳定化变换的形式克应该是：克(x个)=对数x个.

使用的数据库

人类转录因子基因列表来自NCBI(网址：http://www.ncbi.nlm.nih.gov). 转到Entrez Gene并搜索人类转录因子。在筛选出非人类基因后，1658个编码转录因子的人类基因在这个列表中。七份正常和七份癌症原始SAGE-Seq数据已保存在GEO（登录号：。｛“type”：“entrez geo”，“attrs”：｛“text”：“GSE24491”，“term_id”：“24491”｝GSE24491标准).

我们感谢Andrea Richardson（Brigham and Women's Hospital）在获取乳腺肿瘤样本方面的帮助；黄海燕（Haiyan Huang）、李才（Li Cai）、王莫林（Molin Wang）和戴维·哈林顿（David Harrington）进行了有价值的讨论；喜欢尼克森的英语校对。这项工作得到了Dana-Farber女性癌症项目之友（X.s.L.）、NIH R01 1HG004069（X.s.L.）、NCI P50 CA89393 CA和R01 CA116235-04S1（K.P.）、AVON基金会（K.P.）、Terri Brodeur乳腺癌研究基金会（s.C.）和Susan G.Komen基金会PDF0707996（M.s.，R.M.）的支持。

作者贡献：Z.W.进行了计算分析。Z.W.、X.S.L.和K.P.设计了研究并撰写了手稿。Z.W.、X.S.L.、C.M.、A.S.和J.L.开发了分析方法。S.C.、M.S.、R.M.、M.B.和T.N.进行了实验并分析了数据。S.S.为研究提供了正常组织样本。