胸衣：启用差异基因表达分析从头开始汇编成绩单

摘要

下一代测序技术使得在非模式生物中进行差异基因表达研究成为可能。对于这些研究，通过执行从头开始装配在RNA-seq数据上。然而，转录组组装产生大量的连接，在差异基因表达检测之前必须将其聚集成基因。在这里，我们介绍了Corset，这是一种使用共享读取和表达式对contigs进行分层聚类的方法，然后将读取计数汇总到集群中，以便进行统计测试。使用一系列指标，我们证明Corset的性能优于其他方法。胸衣可从购买https://code.google.com/p/corset-project网站/.

背景

下一代RNA测序，RNA-seq，是研究转录组各个方面的强大技术；它具有广泛的应用，包括基因发现、选择性剪接事件检测、差异表达分析、融合检测和识别SNP等变体以及转录后编辑[1],[2]. RNA-seq相对于较旧的技术（如微阵列）的优势之一是，它能够对非模型生物进行转录全分析，因为生成和分析数据不需要参考基因组和注释。当没有参考基因组可用时，转录组为从头开始直接从RNA-seq读取中组装[三]. 有几个程序用于从头开始转录组组装：绿洲[4]和跨深渊[5]延伸天鹅绒[6]和深渊[7]基因组组装器，以及特制转录组组装器，如Trinity[8]. 这些程序能够将数百万次的短读汇编成转录序列，称为连读。

RNA-seq的一个常见且生物学上重要的应用是识别在两种或多种条件下差异表达的基因[9]. 但是，对从头开始组装转录组具有挑战性，因为每个基因有多个连接。由于转录组组装子区分同一基因的不同亚型，并分别报告，因此出现了具有共享序列的多个连续序列。此外，他们经常报告的连续序列并非真正代表不同亚型，而是由诸如测序错误、重复、覆盖范围变化或二倍体个体或混合群体中的遗传变异等人为因素引起的。因此，转录组汇编程序经常报告转录本的片段版本，或仅因SNP或indel而不同的重复连续体。令人惊讶的是，模拟结果表明，即使在没有任何测序错误、SNP或选择性剪接的情况下组装数据，也可以为每个基因生成多个连续序列[10]. 因此从头开始组件通常较大；例如，具有8000万次读取的程序集可以生成数十万个连续[11].

不可避免的是从头开始转录组组装给差异表达分析带来了几个问题：i）读取不能与重复序列明确对齐，并且确定模糊对齐读取的来源容易出错；ii）差异表达测试的统计能力降低，因为读数必须分配给更多的连续数据，从而减少每个连续数据的平均计数；iii），多次测试的调整更为严重；以及iv），一旦发现差异表达的contigs，就很难解释，因为许多基因将多次出现在列表中。对基因而不是连续基因进行差异表达分析将克服这些困难。然而，根据一组从头开始组装的contigs并不简单，文献中也没有对其进行彻底的研究。

从一个从头开始组装转录组（图1)：首先组装RNA-seq读取，然后将读取映射回contigs，然后将contigs聚类为基因，然后总结每个基因簇的表达水平，并进行统计测试以检测差异表达。

图1

对非模式生物进行基于国家的基因级差异表达分析的管道。清除RNA-seq读取是第一步从头开始组装成连续序列。读取结果被映射回转录组，并且必须建立contigs和基因之间的关联（contigs的聚类）。然后估计每个基因的丰度。最后，对计数数据进行统计测试，以确定哪些基因是差异表达的。Corset在单个步骤中执行聚类和计数（虚线框）。

几项研究对该分析流程中的各个步骤进行了比较。例如，不同的从头开始装配前的装配工和步骤（如质量控制）已经过检查[12]–[15]. 同样，对差异表达的基于计数的统计测试方法的选择进行了评估[16],[17]. 然而，很少有研究比较甚至提出了从转录组集合中获取基因水平计数的途径[18],[19]到目前为止，只有一条自动化管道能够满足这一需求[20]; 由Trinity提供，用于运行RSEM[21]然后是edgeR[22]或DESeq[23]. 然而，对于汇编器的选择，此管道是不灵活的。

在本文中，我们介绍了Corset，一种从任何从头开始转录组组装。Corset将一组多映射的读取（每个读取报告多个对齐）从头开始组装转录组并根据共享读取和表达模式的比例对连续体进行分层聚类。如果群体之间的表达水平不同，表达模式允许在共享序列的基因之间进行区分，例如副录。使用映射读取，Corset然后输出基因级计数。然后，可以使用基于计数的框架（如edgeR和DESeq）轻松测试基因级计数的差异表达。我们证明，在一系列指标上，Corset与其他聚类方法相比始终表现良好。此外，由于它是一种独立于汇编程序的方法，它允许将来自不同来源的连接字和转录本组合在一起。它也更易于使用，在一次软件运行中包含了集群和计数步骤。

结果和讨论

胸衣簇连体并计数读数

对非模型生物体进行基因级差异表达分析的第一步是组装连接体，这可以使用多种软件进行。如前所述，这个过程产生每个基因的多个序列或连续序列。因此，下一步是将连续基因分组或聚类为基因，以便于下游差异表达分析。这个集群步骤是Corset的第一步。

Corset要求在转录组组装后，读取映射回contigs，允许读取映射到多个contigs（多映射）。然后，这些多映射读取被用作检测contigs之间序列相似性的代理，并提供有关contigs表达水平的信息。Corset还使用读取信息过滤出映射读取次数较少的contigs（默认情况下小于10次读取）。Corset的方法与其他用于聚类contigs的工具不同，因为大多数其他工具只使用来自程序集的序列信息。

胸衣的工作原理是基于共享读取对contigs进行聚类，但当观察到样本之间的不同表达模式时，会将contigs分开。这是使用聚合层次聚类算法实现的。任意两个contigs之间的距离是根据contigs之间共享的读取数定义的，因此共享读取的比例越低，距离越大（请参见材料和方法）。共享序列的基因，如旁白，可能距离很小，因为许多阅读是共享的。由于我们不希望这些连续体聚集在一起，Corset执行了一项测试，以检测在条件组（或实验组）中，这对连续体之间的相对表达水平是否恒定。如果两个contigs之间的相对表达式不恒定，则将两个Contig之间的距离设置为最大值。这被作为似然比检验纳入算法中，其中，零假设假设两个连续计数之间的比率在条件组中相等，而替代假设允许该比率随条件组而变化。此重叠比率测试的计数数据建模为泊松分布P（P）-阈值约<10^-5默认情况下应用（有关阈值的详细描述和理由，请参阅材料和方法）。

分离具有共享序列但不同表达率的重叠群的重叠群比率测试是Corset聚类算法的新特征之一。尽管这个功能可以被关闭，例如，为了确保差异剪接的亚型聚集在一起，我们发现它可以有效地从不同的基因中分离连接物（附加文件1：图S7）。例如，图2展示了人类ATP5J公司和GABPA公司位于相反链上但具有重叠UTR的基因。人类原发性肺成纤维细胞数据的集合产生了该区域的八个contigs（见材料和方法）。虽然每个基因都有单独的连接（连接1到3和8），但使用非传递方案会导致两个基因组装在一起的连接（连4到6）。当未进行连接性比率测试时，所有这些连接性被分配到同一个簇，并且在敲除和野生型条件（错误发现率（FDR））之间未检测到显著差异表达 = 0.053）。然而，检查成对连体细胞之间的连体细胞计数比率则说明了一个不同的情况（图2B） ●●●●。contig 3和contig 2在样本中的计数比率是恒定的，这意味着它们应该在同一个集群中。相比之下，重叠群3和重叠群4之间的重叠群比率在不同条件下显著不同，因此Corset将它们划分为不同的簇。当对所有成对组合进行测试时，这八个连续序列被分成四个不同的簇，差异表达的统计测试显示簇a和簇d在相反方向上有显著差异表达（FDR = 10^-11和10^-7）。

图2

紧身胸衣使用表达信息梳理不同基因的挫伤。（A）从包含这两个基因的人类基因组的一个区域收集的连体ATP5J公司和GABPA公司Trinity装配了8个连续体（底部轨道），如果未应用连续体比率测试，则将其分组为一组。包括这个测试可以让紧身胸衣将这个区域分成四个簇（框）。值得注意的是，连续4到6是假嵌合体，由重叠的UTR引起ATP5J公司和GABPA公司这些基因是差异表达的，如基础水平覆盖率所示，重复平均值（顶轨）。（B）在聚类时，Corset在计算连续线对之间的距离时检查条件之间的相等表达式比率：这里我们考虑连续线对2和3（顶部）以及连续线对3和4（底部）。绘制每个样本（点）与每个连续点对齐的读取次数的比率。可以看出，contig 2和contig 3在组间的表达率相同，因此聚集在一起，而contig 3和contig 4在不同条件下的表达率不同，因此分裂。此功能帮助Corset分离共享序列但来自不同基因的contigs。

一旦胸衣被应用于完整的数据集，代表基因的连续分组就会被报告，并在以后被称为聚类。Corset还报告与每个集群关联的读取计数。所有读取都被唯一地分配给一个集群（请参见材料和方法）；因此，每个读取只计算一次，即使读取最初是多映射到contigs的。读取计数表可以提供给基于计数的差异表达式程序进行统计测试。

在模型生物数据集上测试Corset

我们使用三个RNA-seq数据集测试了Corset相对于其他聚类和计数方法的性能：鸡雄性和雌性胚胎组织[24]，人类原发性肺成纤维细胞，具有和不具有小干扰RNA（siRNA）敲除HOXA1型[25]，以及在分批和恒化器条件下生长的酵母[26]. 我们选择了三种模式生物来比较我们的从头开始差异基因表达（DGE）结果与基于基因组的分析（此处称为真值数据集）相对应。在鸡的数据集中，我们测试了雄性和雌性之间的DGE。鸡基因之间的同源性，在性染色体上约为90%[24]为聚类算法提供了一个具有挑战性的测试。之所以选择人类数据集，是因为人类是注释最好的物种之一，而酵母用于评估聚类是否对剪接最少的生物体有益。每个数据集都是使用Trinity和Oases进行组装的，它们具有不同的底层组装策略，以确保结果一致。总的来说，六种不同的组件被用作评估紧身胸衣的起点。

胸衣聚类在精确性和召回率之间取得了良好的平衡

我们最初感兴趣的是将Corset产生的聚类与其他可用方法进行比较。Trinity和Oases的输出都提供了一些聚类信息，这是基于组装期间de-Bruijn图的划分（分别称为组件和轨迹）。基于序列相似性的独立工具也经常使用[27],[28]，CD-HIT-EST是一个流行的选择[29],[30].

我们根据CD-HIT-EST和组装者自己的聚类来评估Corset的聚类。对于鸡，组装了300000多个contigs，而对于人类，组装了100000多个Contig（表1). Trinity和CD-HIT-EST报告了大量集群，例如，鸡数据集上有超过200000个集群。默认情况下，Corset会删除支持它们的读取数非常少的contigs，以便在所有情况下都提供较少的集群（附加文件1：表S1）。这使得聚类列表更易于管理，而不会影响检测差异表达的灵敏度。Oases产生的簇也比CD-HIT-EST少，但它将许多不相关的重叠群组合在一起，最大的簇包含数千个重叠群（表1).

表1与Corset相比，各种集群选项的集群数量统计

使用精确度评估聚类（真阳性/（真阳性 + 假阳性））和回忆（真阳性/（真阳性 + 假阴性）从头开始程序集。阳性和阴性的计算方法是，取所有成对的contigs组合，并评估contigs是否正确地放置在同一个簇中（真阳性），是否正确地分为不同的簇（真阴性），是否错误地放置在相同的簇中（假阳性）或是否错误地分开（假阴性）[31]. 真相信息是使用适当的参考基因组注释得出的（参见材料和方法）。由于映射读取次数较少，Corset筛选出的轮廓也被筛选出来，用于评估竞争方法。

我们发现，CD-HIT-EST的精确度普遍较高，但召回率较差。相比之下，Oases的聚类在召回方面表现良好，但在所有情况下的精确度都在零左右。从概念上讲，这再次表明Oases将许多不相关的contigs分组到同一个集群中（过度聚类）。Trinity的聚类结果显示了精确性和召回率之间的更好平衡。在所有情况下，紧身胸衣的召回率和汇编程序的聚类精度都优于CD-HIT-EST（图三)这表明它在精确性和召回率之间提供了良好的平衡。此外，在六个集合中的两个集合中，Corset是最精确的（鸡胚-酵母菌和酵母菌-绿洲）。

图3

不同聚类方法的性能比较。对于汇编程序自己的集群（Trinity或Oases）、CD-HIT-EST和Corset，我们显示了针对召回的精确度。精确度是真阳性与真阳性加假阳性的比率，召回率是真阳性和真阳性加伪阴性的比率。我们显示了六个不同组件的结果：（A）与Trinity组合的鸡数据；（B）用绿洲收集鸡数据；（C）与Trinity组合的人类数据；（D）用Oases汇编人类数据；（E）用Trinity汇编酵母数据；和（F）用Oases收集酵母数据。X表示完美聚集。

聚类性能不仅受聚类算法的选择影响，还受汇编程序的选择和RNA-seq数据集的质量的影响。我们简要研究了装配质量的某些方面如何影响聚类召回和精度。我们发现，随着集合中基因碎片的增加，召回率降低。来自一个没有序列共享的普通基因的轮廓不太可能被任何算法聚集在一起，而几乎完全冗余的轮廓应该总是聚集在一起。大多数基因属于这两个极端之一（附加文件1：图S2B）。汇编程序的聚类和Corset的行为与预期一致，对具有完全重叠的连接性序列的基因提供了接近完美的召回，而当连接性断开时则为零召回。然而，CD-HIT-EST即使对于没有片段的基因也无法实现良好的召回（附加文件1：图S4）。完全脱节的连接部分似乎决定了任何聚类算法都可以实现的最佳召回的上限。

我们发现，当基因共享序列时，如旁系、公共结构域、重叠的UTR或重复序列，就会发生较差的聚类精度，即来自不同基因的重叠群被分组在一起。在某些情况下，这也可能导致嵌合体连接被错误组装（例如，图2). 以前已经证明，在从头开始转录组组件[32]我们还观察到，在我们的组合中，假嵌合体的比例很高，鸡和人的嵌合体约占5%至15%，酵母的嵌合体占40%（附加文件1：图S3）。Oases对共享序列或具有嵌合连接的基因的聚类精度始终低于Corset和CD-HIT-EST。Trinity比Corset稍微差一点。对于没有共享序列的基因，Corset和CD-HIT-EST（附加文件1：图S5）。

这些结果表明，集群性能受底层组装质量的影响（这反过来取决于数据集），但Corset集群在一系列组装质量上是稳健的。

聚类对差异基因表达结果的影响

低精度类似于过聚类，使得一些差异表达的基因无法检测，因为具有不同相对表达的基因组合在一起。此外，集群的功能注释变得模糊。然而，糟糕的回忆类似于欠聚类。它延长了聚类列表的总长度，这有几个后果：不方便进行后续研究（如基因本体论），导致更大的多重测试更正，并增加统计不确定性。为了评估这些对差异表达结果的影响程度，我们使用每个聚类选项进行了基因级差异表达分析。管道中的其余步骤，包括基于计数的丰度估计，在每种情况下都是相同的（参见材料和方法），并使用edgeR进行DGE测试。使用基于基因组的绘图方法，将显著差异表达的簇与差异表达的基因进行比较。如果聚类与基于基因组的分析中的差异表达基因相匹配，则该聚类被视为真正阳性。不管用于从基于基因组的分析中生成真正差异表达基因的统计测试如何，Cuffdiff 2[25]（图4)或edgeR（附加文件1：图S10），我们在Corset与其他contig聚类选项的比较中发现了类似的结果。

图4

聚类对差异基因表达排名的影响。唯一真正差异表达聚类的累积数量与从头开始分析结果如图所示。唯一的真阳性指的是只计算排名表中出现多次的基因的第一个实例。在所有情况下，Corset的性能都与CD-HIT-EST和汇编程序自己的集群相同或更好：（A）与Trinity组合的鸡数据；（B）用绿洲收集鸡数据；（C）与Trinity组合的人类数据；（D）用Oases汇编人类数据；（E）用Trinity汇编酵母数据；和（F）用Oases收集酵母数据。为了进行比较，我们还显示了无聚类的结果，其中分析是在连续的水平上进行的，而不是在聚类的水平上。

我们研究了三种不同的方法来评估聚类对DGE结果的影响。首先，我们检查了唯一真阳性簇的累积数量作为簇总数的函数（图4). 唯一的真阳性是指当一个基因有多个聚类时，只计算排名靠前的聚类。在四个案例中，Corset的表现优于其他替代品（鸡-三位一体、鸡-酶类、人-绿洲、酵母-绿洲），在其余两个案例中表现同样好（人-三位一、酵母-三位居）。该指标因报告给定基因的多个簇而受到惩罚（即回忆不良）。检查这个指标的另一个版本也很有启发性，它不会以这种方式惩罚：唯一真阳性的数量是唯一假阳性数量的函数（图5). 在这种情况下，精度更高的聚类算法比汇编程序的聚类更好，后者在召回方面表现更好。作为聚类的最终评估，我们观察了真相分析中差异表达基因与从头开始组件（表2). Corset与基于基因组的真相分析一致性最高。

图5

聚类对差异基因表达受体操作特性（ROC）曲线的影响。唯一的真阳性差异表达簇与唯一的假阳性簇从头开始分析结果如图所示。唯一的阳性指的是只计算排名表中出现多次的基因的第一个实例。在所有情况下，紧身胸衣的性能都类似于或优于CD-HIT-EST和汇编程序自己的集群：（A）与Trinity组合的鸡数据；（B）用绿洲收集鸡数据；（C）与Trinity组合的人类数据；（D）用Oases汇编人类数据；（E）用Trinity汇编酵母数据；和（F）用Oases收集酵母数据。为了进行比较，我们还显示了无聚类的结果，其中分析是在连续的水平上进行的，而不是在聚类的水平上。

表2基因水平日志中的皮尔逊相关 ₂ 折叠更改

DGE结果也说明了将连续基因聚类为基因的一般重要性；与聚类分析相比，对无聚类的contigs进行差异表达分析后，对于相同数量的唯一真阳性结果，其列表要长得多（图4). 这甚至是来自酵母的燕麦组合的情况，酵母是一种几乎没有选择性剪接的有机体（图4F） 强调了从组装中消除冗余的重要性，即使对于预期存在最小选择性剪接的基因组也是如此。从所有指标来看，Corset是最佳或接近最佳的可用方法。这表明Corset实现的精确性和召回之间的平衡转化为更准确的DGE结果。

Corset允许组合多个转录组

理想的聚类工具将允许组合来自不同来源的转录组，因为来自同一基因的多个转录本将聚集在一起，而不管其来源如何。然而，这仅适用于独立于从头开始汇编程序。虽然一些出版物使用CD-HIT-EST组合多个转录组组装[12],[14],[29]，我们已经表明CD-HIT-EST不是最有效的contig聚类工具。然而，胸衣提供了一种方便的方法来聚类来自不同来源的contigs。首先将读取分别多重映射到每个转录组，然后在一次Corset运行中一起处理所有bam文件。

不同的汇编程序各有优缺点，将几个从头开始汇编程序。为了演示胸衣的实用性，我们将人类数据集中的Trinity和Oases集合聚集在一起。这一组合数据集有效地将contigs的数量增加了一倍。Corset能够处理这种冗余度，得到一个簇数为37741的组合转录组，分别少于Trinity或Oases集合的43664和38477。此外，该联合转录组包含注释了约200个单独使用任一组分转录组未检测到的额外基因的连接。相比之下，CD-HIT-EST在相同的组合数据集上生成了115980个集群。

组合转录组的另一个应用是当部分组装的基因组或注释可用时。补充从头开始使用基于基因组的数据组装数据有几个优点：i）它增加了已知转录序列的数量，例如，因为注释中的基因或基因区域很少或没有读取覆盖范围，而这些基因或区域在组装中不存在；ii）从头开始如果组装的contigs与已知基因聚集在一起，则可以很容易地进行注释；以及iii）如果引用注释中的转录本与两者重叠，则允许将程序集中的断开片段聚集在一起。我们通过使用Corset将人类RNA-seq数据集的Trinity转录组与人类Ensembl版本73注释相结合，证明了这一最终益处。我们随机抽取了50%、25%、12%和6%（约10万、50万、25万和12.5万份成绩单）的完整合奏转录组，以模拟部分注释。Trinity contigs在聚类召回方面有显著改进，且不影响聚类精度（图6A） ●●●●。

图6

补充a 从头开始 用额外的转录组组装。补充a从头开始用任意一个连续的组件（A）部分注释或（B）相关物种提高了从头开始组装连体。我们显示了针对Trinity contigs计算的召回率和精确度。（A）我们从Ensembl的人类注释中随机抽取100%、50%、25%、12.5%和6%的转录本，以模拟部分注释，将人类RNA-seq读取映射到每个集合，并使用Corset将读取与映射到Trinity集合的读取聚集在一起。（B）我们将人类RNA-seq读码映射到黑猩猩、猩猩、猕猴、绒猴和丛林宝宝的Ensemble注释上，然后使用Corset将这些读码与映射到Trinity汇编的读码进行聚类。”两个图中的“无”表示三一议会本身。

最后，我们将这个想法扩展到使用相关物种的参考注释[33]. 人类RNA-seq读取独立映射到黑猩猩、猩猩、猕猴、绒猴和丛林宝宝的Trinity汇编和Ensembl版本73转录序列，相关物种的映射严格性较低（见材料和方法）。我们的结果表明，即使使用相关物种的注释，聚类召回也会得到改善，而且不会损失精确度（图6B） ●●●●。这种策略的有效性将取决于物种之间的差异，在这个例子中，物种之间的差异从大约600万年前（我的；黑猩猩）、15我的（猩猩）、32我的（猕猴）、46我的（狨猴）到68我的（丛林婴儿）[34]. 此外，这种改进还将取决于注释的完整性（这些物种的转录本数量是人类Ensembl注释的一半到三分之一）。

Corset总结了对基因级计数的读取

如前所述，Corset不仅可以对连续体进行聚类，还可以以可直接由DGE分析软件使用的读取计数的形式为每个聚类生成表达水平。此功能是为了方便起见而提供的，因为它用一个步骤取代了聚类contigs和估算丰度的两步过程。

我们将软件计数方面的性能与其他三条用于基因级计数估计的管道进行了比较：（1）RSEM[21]; （2） 从每个簇中选择最长的contig作为代表序列并映射到该序列；（3）将每次读取映射到所有contigs，只允许进行一次对齐，然后在集群内聚合计数（参见材料和方法）。在每种情况下，都使用了Corset聚类。我们发现，总的来说，所有计数方法都产生了与Corset相似的结果；Corset和其他方法产生的计数之间的Pearson相关性一直很高（附加文件1：表S2）。尽管RSEM和紧身胸衣报告了多达95%的集群的相同计数（附加文件1：表S3），我们发现少数集群的计数存在显著差异。在这些集群中，RSEM报告的计数较少（附加文件1：图S11A）。此外，对于这些集群，胸衣计数在生物复制之间的变异性较小（附加文件1：表S4和图S11B），表明它们更加精确。作为对差异的最终评估，我们将RSEM聚类级别计数与通过对映射到基因注释的读取（附加文件1：图S12A）。这证实了RSEM低估了一小部分集群的计数，导致下游分析中缺少真正的差异表达（附加文件1：图S12B）。

因此，我们没有发现证据表明，使用Corset生成的计数比其他程序（如RSEM）进行基因级分析存在劣势。相反，我们发现了一些微妙的线索，即简单的基于计数的丰度估计方法在检测差异表达方面可能更稳健从头开始组装转录组，而不是试图首先在转录水平推断丰度的方法。然而，我们重申，所有计数方法都给出了类似的结果。

结论

最近，人们开始关注使用RNA-seq分析非模式生物转录体的潜力，一些研究探索了分析管道中的各个步骤，例如清洁读取的影响[13]，数字标准化[35]，不同的装配商[12]以及装配后的改进[36]. 然而，在RNA-seq的许多应用中，感兴趣的结果不是转录组本身，而是检测样本之间差异表达的基因。为此，关于如何从组装的转录组进展到差异表达基因列表，目前几乎没有工具甚至指南。本研究提出了一种新的聚类算法Corset从头开始收集contigs并生成基因级计数。本研究也是第一次比较不同管道的DGE分析从头开始组装的转录体，证明进行基因水平而非转录水平的分析是有利的，即使是对酵母等物种，也要进行最小的选择性剪接。

然而，执行基因级分析需要从头开始将连接聚集成基因。在这里提出的算法之前，聚类选项仅限于汇编程序自己的分组或基于序列的聚类工具。此外，尚不清楚这些方法在从头开始组装转录组；具体来说，它们处理诸如基因片段化、错误嵌合体和重复等问题的能力，所有这些都可能导致contigs的聚类错误。

我们的方法Corset利用了序列相似性和可用于聚类contigs的表达数据。表达数据的使用意味着Corset有能力将带有错误嵌合体的Paralogue和contigs分为不同的簇。基于相对表达分离连续序列的一个可能结果是，基因的差异剪接异构体可能被分割成单独的簇。对于基因级差异表达分析，我们认为这没有任何缺点；其中一个或两个亚型簇应被检测为差异表达。但是，如果用户希望确保异构体聚集在一起，可以关闭此选项。总的来说，我们发现Corset提供的聚类在我们检查的几个指标中表现得比其他方法更好。类似地，Corset以计数形式提供的表达式数据得出的结果与所有替代估计值相等，有时略为准确。

因此，Corset在单个软件工具中提供了新的方法，有效地取代了从从头开始将转录组组装到基因水平计数。胸衣很容易运行，因为不需要对bam文件进行索引或排序，并且可以处理单端、成对读取或混合读取。最后，Corset提供了一种方便的方法来合并来自不同从头开始装配、参考注释或基因组引导装配。我们相信这些特征将对非模式生物的RNA-seq分析非常有益。

材料和方法

数据集集合

我们使用来自鸡、人和酵母三种模式生物的公开可用RNA-seq数据进行差异基因表达分析。所有数据集均由Illumina HiSeq 2000的100-bp配对读取数据组成。对于每个数据集，我们修剪了读取[37]然后进行了三项分析：从头开始使用Oases和Trinity组装的转录组，以及一种基于基因组的分析——“真相”——用于比较。来自艾尔斯的鸡数据集等. [24]，Short Read Archive（SRA）登录号SRA055442，包含约12亿次读取。对于从头开始分析我们只使用了一条数据通道（大约3.2亿次读取），因为整个数据集计算量太大，无法组合。然而，所有的数据都被用于基于基因组的“真相”分析。该数据集由八个样本组成——雄性和雌性胚层，以及雄性和雌性第4.5天性腺组织，一式两份。Trapnell发布的数据集等. [25]，基因表达综合登录GSE37704，来自人类原发性肺成纤维细胞，siRNA敲除HOXA1型。该数据集包含三个敲除副本和三个总计超过2.31亿次读取的控件。最后，我们包含了一个酵母数据集，SRA登录号为SRR453566到SRR4536571，发表于Nookaew等. [26]. 该数据集包含大约3600万次读取。三个复制品在分批条件下生长，三个在恒化器条件下生长。

基于基因组的“真相”分析

为了衡量不同聚类和丰度估计算法的性能，我们使用基于基因组的分析导出了一个“真值”集。

确定以下各项之间的对应关系从头开始组装的连接和参考注释基因，我们使用BLAT将组装的连接与注释对齐[38]（最小长度为200个底座，最小一致性为98%）。嵌合体contigs被视为来源不明。我们确定嵌合连接是指那些匹配两个或多个真基因（如上所述）且基因之间重叠少于100个碱基的连接。在其他情况下，如果一个contig与多个基因对齐，则将其分配给具有最长对齐长度的基因。当将差异表达的“真相”基因与从头开始聚类，我们将一个聚类分配给与其大多数邻接物相同的“真”基因。在“真值”集合中找不到的任何连续体或簇都被排除在显示的结果之外。由于读取映射次数较少而被Corset移除的轮廓也被排除在外。

为了计算“真”微分表达式，首先使用TopHat 2.0.6版映射读取[39]人类、鸡和酵母基因组的hg19、galGal3或sacCer3版本。在所有情况下，我们都向TopHat提供了基因注释（人类的RefSeq、鸡的Ensembl（v.70）和酵母的酵母基因组数据库），以支持剪接位点检测。这些相同的基因注释通过“gffread-merge”进行处理，以给出位点级注释。运行Cuffdiff 2.1.1以检测差异基因表达（-单位选项）。我们使用“gene_exp.diff”中的“significant”位点作为真阳性。作为cuffdiff 2的替代方案，我们还使用基于基因组的边缘R分析来定义真相（结果显示在附加文件中1：图S10）。EdgeR的运行方式与从头开始组装（请参阅下面的“统计测试”）。

从头开始装配

Oases 0.2.06（Velvet版本1.2.07）用于组装人类和酵母数据，其中kmer长度为19、23、27和31。对于鸡肉数据集，我们使用了31、41、51、61和71的kmer长度。使用Trinity-r2012-10-05创建鸡Trinity组件，使用Trinity-r2013-02-25创建人和酵母组件。所有情况下都使用了默认参数，最小连接长度为200个碱基。其他文件1：图S1、S2和S3显示了装配质量。

映射

读取已映射到从头开始使用领结将部件作为成对的发动机对齐[40]. 对于只允许一次对齐的单个映射，我们使用领结选项--最好的。对于多映射对齐，我们使用了选项--全部。当映射到相关物种时，我们使用领结设置--所有-m 6-n 3-e 1000-X 1000，以允许大量不匹配。对于人类数据集，这导致了30%（丛林宝宝）到70%（黑猩猩）的读对映射，相比之下，Trinity程序集的读对对应约为75%。

群集

我们使用CD-HIT-EST和默认参数对转录组进行聚类。对于汇编程序集群，我们从汇编fasta文件中的contig名称中提取集群。例如，对于Trinity，contig“comp1_c2_seq3”属于集群“comp1_c2”。对于绿洲，“Locus_1_Transcript_3/10_Confidence_0.000_Length_268”属于集群“Locus_1”。为了获得Corset聚类，我们将读取多次映射到转录组，并以实验组作为参数执行Corset(-克选项）。对于图中所示的微分表达式结果4和5和表2，我们使用下面描述的“单映射然后求和”方法估计计数。

丰度估算分析

对以下四种方法进行比较，以评估哪种方法的DGE结果最佳。在所有情况下，聚类都是相同的，使用Corset生成，实验组通过-克选项和使用-0米选项（以便报告所有contigs）。使用edgeR进行统计测试。

RSEM（RSEM）

使用命令“convert-sam-for-RSEM”将多映射bam文件转换为RSEM所需的格式。转录组使用“rsem-prepare-reference--no-polyA--no-bowtie--transcript-to-gene-map”进行准备，Corset聚类作为参数传递。使用“rsem计算表达--bam-配对端”估计基因丰度，并从“.genes.results”文件中提取“expected_counts”。

代表性连续法

选择最长的contig代表每个簇。读数被单一映射回这些重叠群。使用samtools-idxstats命令计算映射到每个代表contig的读取数。因为这些计数数据是每次读取的，所以我们将其除以二，以获得每个片段的计数。

单映射然后求和

我们使用samtools-idxstats将读取映射到所有contig，并计算每个contig的重叠数[41]. 为了获得基因水平的计数，我们对一个集群内所有contigs的计数进行了汇总。因为这些计数数据是每次读取的，所以我们将其除以二，以获得每个片段的计数。

紧身胸衣

我们将读取映射到转录组，并使用上述选项执行Corset。

统计测试

使用edgeR处理集群级计数数据。对于鸡的数据，我们用四个条件组（两种性别和两个时间点）对数据进行建模，如Ayers在al. [24]，但只测试了稍后时间点的雄性和雌性之间的差异。其他数据集各有两个条件组（共六个样本有三个重复），并对这些组之间的差异进行统计测试。在所有情况下，我们都使用edgeR GLM框架进行分段离散估计[42]. 以相同的方式对所有人进行统计测试从头开始程序集。使用Cuffdiff 2（图4和5)和edgeR（附加文件1：图S10）。虽然这些结果给出了一个稍微不同的重要真值基因列表，但将Corset与其他聚类方法进行比较的结果是相似的。

Corset算法

我们的软件接受一组bam格式的多映射读取对齐作为输入（每个样本一个或多个文件）。然后，算法按以下方式进行：

• 1 对每个读取对齐进行解析，并提取读取ID和contig ID。对于每次读取，我们都存储它映射到的一组contigs。

• 2 具有10个或更少读数的轮廓将被过滤掉。这一步骤对算法来说并不重要，但具有减少最终聚类总数以及每个聚类的重叠群平均数量的效果，这可以简化分析中的后续步骤。

• 三。 解析读取的数据并形成超级集群。每个超级集群包含与同一超级集群中的另一个contig共享一个或多个读取的所有contig。

• 4 然后，对于每个超级聚类，我们执行类似于中的算法的聚集层次聚类[43]，但具有如下所述的距离和连杆。使用层次聚类而不是其他聚类方法，因为它易于计算。

4.1我们使用指标创建距离矩阵：

$$\mathit{<trans\; data-src="distance">距离</trans>}<trans\; data-src="=">=</trans>\left\{\begin{array}{c}\hfill \mathrm{<trans\; data-src="1">1</trans>}<trans\; data-src="-">-</trans>\frac{{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathit{<trans\; data-src="ab">ab公司</trans>}}}{<trans\; data-src="min">最小值</trans>\left({\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathrm{<trans\; data-src="a">一</trans>}}<trans\; data-src=",">,</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathrm{<trans\; data-src="b">b条</trans>}}\right)}<trans\; data-src=",">,</trans>\mathit{<trans\; data-src="contig">康蒂格</trans>}\mathit{<trans\; data-src="ratio">比率</trans>}\mathit{<trans\; data-src="the">这个</trans>}\mathit{<trans\; data-src="same">相同的</trans>}\hfill \\ \hfill \mathrm{<trans\; data-src="1">1</trans>}<trans\; data-src=",">,</trans>\mathit{<trans\; data-src="contig">康蒂格</trans>}\mathit{<trans\; data-src="ratio">比率</trans>}\mathit{<trans\; data-src="different">不同的</trans>}\hfill \end{array}\right.$$

(1)

哪里，R（右） _一 是该映射到contig的读取总数一所有样本，以及R（右） _ab公司 是映射到两个重叠群的读取总数一和contigb条，跨越所有样本。因此，距离在零和一之间有界，零表示一对冗余的连续线，一表示没有相似性。”Contig ratio是指contigs的表达一和b条根据条件组的测量结果彼此成比例。我们假设，当两个连续序列起源于同一基因，并且没有选择性剪接时，这是真的。或者，如果连续基因不是来自同一个基因，或者存在选择性剪接，那么它们的表达不一定是成比例的，如果一个连续基因差异表达就会发生这种情况。我们使用“连续比率测试”来测试这些场景，该测试按以下方式进行。让第页 _{美国国际期刊} 是指向contig的地图读取次数一在条件下我，用于j个第个复制。然后我们将该映射的读取次数近似为contig一，在条件下我作为：

$${\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="ai">人工智能</trans>}}<trans\; data-src="=">=</trans>\mathrm{<trans\; data-src="1">1</trans>}<trans\; data-src="+">+</trans>{\displaystyle \underset{\mathrm{<trans\; data-src="j">j个</trans>}}{<trans\; data-src="\sum ">\sum </trans>}}\left({\mathrm{<trans\; data-src="r">第页</trans>}}_{\mathit{<trans\; data-src="aij">美国国际期刊</trans>}}<trans\; data-src="+">+</trans>\mathrm{<trans\; data-src="0.5">0.5</trans>}{\mathrm{<trans\; data-src="r">第页</trans>}}_{\mathit{<trans\; data-src="abij">阿比吉</trans>}}\right)$$

(2)

共享读取术语，第页 _阿比吉 ，用于避免重复计算读取次数。添加一个作为偏移，以确保X（X） > 0

然后将连续计数建模为泊松分布。请注意，我们使用泊松模型计算速度：

$${\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="ai">人工智能</trans>}}<trans\; data-src="~">~</trans>\mathrm{<trans\; data-src="Pois">波阿斯</trans>}\left({\mathrm{<trans\; data-src="\mu ">\mu </trans>}}_{\mathit{<trans\; data-src="ai">人工智能</trans>}}<trans\; data-src="=">=</trans>{\mathrm{<trans\; data-src="f">（f）</trans>}}_{\mathrm{<trans\; data-src="i">我</trans>}}{\mathrm{<trans\; data-src="\mu ">\mu </trans>}}_{\mathit{<trans\; data-src="bi">双</trans>}}\right)$$

(3)

$${\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="bi">双</trans>}}<trans\; data-src="~">~</trans>\mathrm{<trans\; data-src="Pois">波阿斯</trans>}\left({\mathrm{<trans\; data-src="\mu ">\mu </trans>}}_{\mathit{<trans\; data-src="bi">双</trans>}}\right)<trans\; data-src=",">,</trans>$$

(4)

哪里μ _人工智能 是连续计数的平均值一在条件下我和（f）是一个比例常数。定义（f） _我= μ _{人工智能/} μ _双 作为contig之间比例表达的真正度量一和b条在条件下i、。我们想检验零假设，H（H） ₀ ：f _我  = （f） _我'  = f、，比例常数与条件无关，H（H） ₁ ：f _我  ≠ （f） _我' .

条件比例常数的估计我从连续计数中获得，即：

$${\stackrel{<trans\; data-src="^">^</trans>}{\mathrm{<trans\; data-src="f">（f）</trans>}}}_{\mathrm{<trans\; data-src="i">我</trans>}}<trans\; data-src="=">=</trans>\frac{{\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="ai">人工智能</trans>}}}{{\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="bi">双</trans>}}}$$

(5)

公共比例常数估计为：

$$\stackrel{<trans\; data-src="^">^</trans>}{\mathrm{<trans\; data-src="f">（f）</trans>}}<trans\; data-src="=">=</trans>\frac{{\displaystyle {<trans\; data-src="\sum ">\sum </trans>}_{\mathrm{<trans\; data-src="i">我</trans>}}}{\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="ai">人工智能</trans>}}}{{\displaystyle {<trans\; data-src="\sum ">\sum </trans>}_{\mathrm{<trans\; data-src="i">我</trans>}}}{\mathrm{<trans\; data-src="X">X（X）</trans>}}_{\mathit{<trans\; data-src="bi">双</trans>}}}$$

(6)

我们可以使用检验统计量的似然比检验来检验无效假设：

$$\mathrm{<trans\; data-src="D">D类</trans>}<trans\; data-src="=">=</trans><trans\; data-src="-">-</trans>\mathrm{<trans\; data-src="2">2</trans>}\left(<trans\; data-src="ln">自然对数</trans>{\mathrm{<trans\; data-src="l">我</trans>}}_{\mathrm{<trans\; data-src="0">0</trans>}}<trans\; data-src="-">-</trans><trans\; data-src="ln">自然对数</trans>{\mathrm{<trans\; data-src="l">我</trans>}}_{\mathrm{<trans\; data-src="1">1</trans>}}\right)<trans\; data-src=",">,</trans>$$

(7)

其近似为分布在n个 _条件 -零假设下的1个自由度。在这里n个 _条件 是条件总数，我 ₀ 是零假设下的可能性我 ₁ 是替代假设下的可能性。

拒绝零假设的任何一对连续体被定义为具有“连续体比率差异”，其距离将增加到最大值1。我们发现在D上设置一个阈值（相当于P（P）-阈值为10^-5。阈值和条件数之间的关系参数化为D类 _门槛  = 15 + 2.5 × n个 _条件 。此关系仅为近似关系，在以下情况下有效n个 _条件  < 10.这种近似不应影响聚类，因为我们发现DGE结果在很大范围内是稳健的P（P）-值（附加文件1：图S7）。

4.2层次聚类通过将距离最小的两个contigs合并在一起进行。然后使用下面的链接标准更新与此新簇对齐的读取数，并重新计算距离矩阵（如步骤3所示）。请注意，Corset使用的链接不同于标准链接方法，例如单个链接，因为它依赖于距离矩阵之外的信息：

$${\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathrm{<trans\; data-src="a">一</trans>}<trans\; data-src="’">’</trans>}<trans\; data-src="=">=</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathrm{<trans\; data-src="a">一</trans>}}<trans\; data-src="+">+</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathrm{<trans\; data-src="b">b条</trans>}}<trans\; data-src="–">–</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathit{<trans\; data-src="ab">ab公司</trans>}}$$

(8)

$${\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathrm{<trans\; data-src="a">一</trans>}<trans\; data-src="’">’</trans>\mathrm{<trans\; data-src="c">c（c）</trans>}}<trans\; data-src="=">=</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathit{<trans\; data-src="ac">交流电</trans>}}<trans\; data-src="+">+</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathit{<trans\; data-src="bc">公元前</trans>}}<trans\; data-src="-">-</trans>{\mathrm{<trans\; data-src="R">R（右）</trans>}}_{\mathit{<trans\; data-src="abc">美国广播公司</trans>}}$$

(9)

其中contigs一和b条那些是合并到集群中的吗a’。R（右） _{美国广播公司} 是映射到所有contigs的读取数一,b条、和c（c）.

4.3反复重复步骤4.1和4.2，直到将所有contigs分组为单个簇，或者当前最小距离增加到距离阈值以上。然后输出集群和每个集群的读取数。对齐到多个簇的读取被随机分配到它们对齐到的其中一个组。这只占我们测试中100-bp配对读取的1-5%。

我们的结果相对于距离阈值的选择是稳健的。默认值0.3是根据经验选择的，因为它对DGE结果来说稍微好一些（附加文件1：图S9），但与0.1和0.9之间的任何阈值相比，结果没有显著差异（附加文件1：图S8和S9）。与阈值相关的稳健性可以用距离接近0或1的大多数连续线对来解释（例如，附加文件1：图S2B）。

默认值P（P）-似然比测试的阈值，10^-5，被选为高水平的多重测试。该值是围绕典型注释中预期的基因数量设计的。我们再次发现，我们的结果在很大范围内对该参数的选择是稳健的，10^-3到10^-8个（附加文件1：图S7）。

我们的软件是开源的，可以作为C++源代码tar ball从[44]. 它已经在Linux x86和Mac OS X 10.7操作系统上编译和测试。使用Intel Xeon E7-8837的一个内核完成代码所需的时间从5分钟到5小时不等，通常比其他管道更快。在最坏的情况下，内存消耗小于60GB，程序会解析超过200GB的bam文件。内存需求高于其他聚类和丰度估计工具，但远低于从头开始我们测试的数据集的集合。