MADAP，用于解释一维基因组注释数据的灵活聚类工具

摘要

在分析来自高通量技术的大规模基因组注释数据时，一个反复出现的任务是识别噪声信号谱中的峰值或簇。这些应用的例子包括基于转录起始位点谱的启动子定义、基于ChIP-ChIP数据的转录因子结合位点定位以及从全基因组SNP谱中识别数量性状位点（QTL）。这种分析的输入是一组与计数或强度相关的基因组坐标。输出由离散数量的峰值组成，峰值具有各自的体积、延伸和中心位置。为此，我们开发了一个灵活的一维集群工具，称为MADAP，我们将其作为web服务器和独立程序提供。一组参数使用户能够针对特定问题自定义程序。web服务器以文本和图形形式返回结果，对于中小型应用程序以及需要本地安装的大型应用程序的评估和参数调整非常有用。用C编写的程序++可以从免费下载ftp://ftp.epd.unil.ch/pub/software/unix/madap。MADAP web服务器可以访问http://www.isrec.isb-sib.ch/madap/.

简介

各种实验性基因组注释技术为散布在完整或部分基因组上的染色体位置提供计数、概率或强度值。这些技术包括cDNA和标记测序协议，以绘制mRNA的5′端和3′端(1,2)ChIP芯片分析揭示转录因子结合位点和表观遗传标记，以及用于各种基因型-表型关联研究的高密度SNP分析平台。处理此类数据时经常出现的一个问题是通过某种信号检测和噪声滤波算法识别单个簇（也称为峰值或孤岛）。基因组注释数据的分析在许多方面对任何聚类软件提出了挑战。基因组坐标与可用的定量测量值通常在分析的染色体范围内分布不均匀。此外，水平轴本质上是离散的，在大多数情况下，目标物体的大小和形状基本上是事先未知的。此外，实验方案中的位置特异性读数可能是几个重叠的生物和技术过程的函数。例如，在启动子定位的情况下，cDNA 5′末端的数量被认为主要反映转录起始事件，但技术伪影或cDNA合成提前终止也可能导致该信号。

峰值识别算法有多种变体(三,4)并举例说明了它们在MS数据和色谱剖面或时间序列分析中的应用(5). 然而，对这些方法在基因组注释数据中的应用的描述很少。现有方法经常缺乏灵活性，以R中的mclust包为例(6)，实施附加约束是困难的。这些方法通常基于关于信号生成过程和峰值形状的或多或少明确的物理假设，因此不能直接移植到新的应用程序。此外，标准聚类算法容易受到非典型分布的干扰。通常，没有实验金标准可用于评估结果。在这种情况下，最终的参考仍然是应用于代表性示例的人类直觉。

我们最初开发了MADAP程序，用于从全长cDNA到基因组序列的mRNA 5′末端图谱推断启动子(7). 最近，我们发现MADAP在解释ChIP芯片数据方面很有用。考虑到可能会有更广泛的用途，这里介绍的web服务器主要用于快速评估MADAP对新类型数据的适用性。找到合适的参数设置可能是将MADAP应用于新基因组注释问题的最大障碍。请注意，以前建立的参数可能不适合具有不同特征的新数据集，例如背景噪声。该服务器还可用于小型生产应用程序。然而，对于大规模基因组注释任务，我们建议使用本地安装的程序，该程序可以从我们的FTP站点免费下载。

用户输入和算法描述

MADAP网络服务器将上传的文件作为输入，该文件包含一组制表符分隔的数字，例如表示染色体上的位置。描述数据点的数字出现的频率与此给定位置的测量强度（或强度）相对应，例如全长转录本5′端的次数(8)在该位置观察到。或者，可以在gff格式的文件中提供数据(www.sequenceontology.org/gff3.shtml网站). 在这种情况下，每个特征的频率分数必须在gff文件的第6列（“核心”）中指示。

MADAP的功能是确定描述输入数据集的最可能模型。使用具有中心、标准偏差和相对频率的正态（高斯）分布的混合来建模数据点，这种方法也称为混合建模。虽然隐藏在输入数据中的簇的分布形状通常不为人所知，也不一定类似于正态分布，但我们观察到，如后面进一步描述的MADAP算法可以很好地处理大多数类型的分布。使用正态分布的优点是，在给定正态分布集的中心位置和标准偏差的情况下，可以很容易地计算在特定位置观察到单位事件（例如一个cDNA 5′端）的概率。正态分布的数目、中心位置和相对频率最初是从数据中推导出来的。使用标准期望最大化（EM）方法(4)MADAP优化了零到多个簇的中心位置及其标准偏差。这种方法的一个已知缺点是，远离聚类中心的孤立点（杠杆点）可能会对选择最能描述点簇的模型产生不良的影响。为了控制对孤立点的干扰，我们在混合模型中添加了额外的非高斯均匀“背景”成分，从而减少了对模型选择的负面影响。因此，MADAP不同于标准算法(9)在添加背景分布和指定一组附加约束的可能性中，将在下面对模型的优化步骤进行总结。

在第一个模型初始化步骤中，程序为每个可能数量的集群生成几个初始模型。簇的最小和最大数目可以由用户指定的参数定义；簇的初始数量还受到输入数据中不同数据点的数量的限制。簇的中心位置最初归因于频率最高的数据点。相邻“集成范围”内的数据点被包括在每个初始集群中。后台控制的第二个参数从模型初始化的所有位置减去用户定义的常量。使用完整数据再次计算后续步骤。

在接下来的步骤中，每个初始模型都使用EM算法进行迭代进化。簇的初始中心将被优化，直到数据可能性稳定为止，如果达到最大迭代次数，则丢弃这些中心。在第三步中，EM步骤生成的模型需要符合用户定义的模型约束。这些约束包括归属于一个簇的最小数据点数量和相邻簇的峰值之间的最小距离。移除模型中不符合的簇，并使用簇数减少的模型重复前面描述的EM步骤。如果没有剩余的集群，则拒绝该模型。记录满足所有约束的模型。

在对所有初始模型进行优化后，最终模型被选为数据似然最高的模型。为似然计算提供了两种变体：混合模型下的通常似然(4)，以及将每个数据点归因于该位置密度最高的集群后计算的可能性（请参阅web服务器上的解释文档）。

由于计算资源的限制，MADAP的web服务器版本具有一些限制，特别是在初始集群的数量方面。建议用户要么将其数据集拆分为假定包含少于50个集群的范围，要么在本地计算机上安装MADAP。要分析我们基础设施上的较大数据集，请联系作者。可以在MADAP web服务器的站点上找到对算法及其参数的进一步描述。

用户定义的示例和参数

在下文中，我们将通过两个示例描述如何使用MADAP从全长cDNA数据或ChIP芯片数据定义启动子。请注意，这是一个部分探索性（无监督）数据分析问题。给定基因组范围内转录起始位点的谱，没有客观的方法来回答它们代表多少启动子的问题。共识的答案可能来自与新方法和数据的互动学习过程。因此，我们为MADAP算法配备了各种参数，使用户能够在所需的方向上指导分区过程。web服务器上的默认参数对应于转录起始站点（TSS）任务的优化值。特别是，我们尝试将本应用中选择的参数与生物信号的假设联系起来，如下所示，按估计的重要性顺序排列。为了提高TSS映射的鲁棒性和精确度，我们要求每个簇至少有10个计数（此处为cDNA 5′端）（参数n个= 10). 我们假设平均约70%的全长转录物在其“主要”TSS的20 bp内启动，因此最好用初始高斯分量的固定标准偏差来描述(d日= 20). 替代启动子被定义为最小距离为50 bp的相邻TSS(第页= 50).

图1显示了第二个演示文件上MADAP web服务器的输出，与ENSEMBL基因组查看器显示的相应基因组注释对齐(11).

在MADAP web服务器上，结果以图形方式显示在输入数据直方图的顶部。在总图中，确定簇的分布显示为在相应簇的近似位置绘制的数字。每个簇的详细视图可以确定中心和估计的标准偏差。除了推断集群位置的可视化表示外，MADAP服务器的输出还包括文本格式的可解析输出文件。MADAP的主要结果在一个文件“输出”中报告，该文件包含所用参数的概要以及在上述模型下发现的簇的描述。补充文件旨在帮助跟踪MADAP的迭代行为，包括带有迭代步骤大纲的“组件”文件，以及恢复优化模型属性的“摘要”文件。程序执行期间遇到的最终问题会在“错误”文件中报告。

结论

这里我们为集群程序MADAP提供了一个web服务器，该程序最初是为确定TSS而开发的(7). 尽管MADAP使用内部正态分布，但其设计用于模拟任何形状的非接触分布，并证明在这方面具有显著的鲁棒性。

其他人利用cDNA全长测序数据或5′标记（5′SAGE，CAGE）进行启动子定位，并可能开发出TSS聚类问题的替代解决方案。然而，据我们所知，这些方法都没有通过web服务器明确描述或公开。

MADAP原则上是通用的，足以解释任何来源的数据，根据以中心位置、体积和延伸为特征的有限数量的簇。在初级基因组注释数据的范围内，我们设想在分析转录物的3′端和推断polyA信号时，以及在从ChIP芯片或平铺阵列获得的数据中，扩展MADAP的使用。

补充数据

补充数据可在NAR Online上获取。

确认

支付这篇文章的开放存取出版费用的资金由瑞士政府提供。

利益冲突声明。未声明。

参考文献