托比·迪伦，霍金;吉勒姆·里盖尔;保罗·费恩海;纪尧姆·布尔克 基因组数据峰值检测的约束动态规划和监督惩罚学习算法。 （英语） Zbl 1502.92008年 J.马赫。学习。物件。 21，第87号论文，40页（2020年）. 摘要：基因组数据中的峰值检测涉及到将DNA序列读取的计数分段到染色体的不同位置。目标是检测计数较高的峰值，并过滤计数较低的背景噪声。针对该问题的大多数现有算法都是针对特定数据类型中的模式定制的无监督启发式算法。我们针对这个问题提出了一个有监督的框架，使用带有学习惩罚函数的最优变化点检测模型。我们提出了第一个动态规划算法，该算法保证计算相邻段平均参数之间有约束的变点检测问题的最优解。实现该算法需要选择惩罚参数，该参数决定估计的分段数。我们展示了如何G.里格尔等[“使用最大裕度区间回归学习变化点检测的稀疏惩罚”，in：Proc.30th ICML，172-180（2013）]可用于选择该惩罚。我们将我们算法的最终实现与具有两种不同模式（宽H3K36me3数据和尖锐H3K4me3数据）的标记ChIP-seq数据集基准中的几个基线进行了比较。虽然基线非监督方法仅为单个模式提供准确的峰值检测，但我们的监督方法在所有数据集中都达到了最先进的准确性。我们提出的动态规划算法的对数线性定时使其可扩展到目前常见的大型基因组数据集。我们的实现可在CRAN上的PeakSegOptimal R软件包中获得。 引用于4文件 MSC公司： 92D10型 遗传学和表观遗传学 92D15型 与进化有关的问题 68T05型 人工智能中的学习和自适应系统 49升20 最优控制与微分对策中的动态规划 关键词：非凸的;受约束的;优化;变点;分段 软件：R（右）;CRAN（起重机）;峰值分段最佳;峰值SegDisk;PeakSeg公司;绿色荧光粉;同位素;科赛格 PDF格式BibTeX公司 XML格式引用 \textit{T.D.Hocking}等人，J.Mach。学习。第21号决议，第87号论文，40页（2020年；Zbl 1502.92008年） 全文： 链接 参考文献： [1] H Ashoor、A H´erault、A Kamoun、F Radvanyi、VB Bajic、E Barillot和V Boeva。HMCan：一种使用ChIP-seq数据检测癌症样本中染色质修饰的方法。生物信息学，29（23）：2979-29862013。 [2] IE Auger和CE Lawrence。段邻域的最佳识别算法。《公牛数学生物学》，51:39-541989年·Zbl 0658.92010号 [3] T Bailey、P Krajewski、I Ladunga、C Lefebvre、Q Li、T Liu、P Madrigal、C Taslim和J Zhang。ChIP-seq数据综合分析的实用指南。《公共科学图书馆·计算生物学》，9（11）：e10033262013年。 [4] A Barski、S Cuddapah、K Cui、T-Y Roh、DE Schones、Z Wang、G Wei、I Chepelev和K Zhao。人类基因组组蛋白甲基化的高分辨率分析。细胞，129（4）：823-8372007。 [5] JD Buenrostro、B Wu、HY Chang和WJ Greenleaf。ATAC-seq：一种全基因组染色质可及性检测方法。2015年分子生物学现行协议。 [6] ENCODE项目联盟。DNA元素百科全书（ENCODE）的用户指南。《公共科学图书馆·生物学》，9（4），2011年。 [7] C Gao、F Han和C-H Zhang。分段常数模型的极大极小风险界。预打印arXiv:1705.063862017年。 [8] EL Greer和Y Shi。组蛋白甲基化：健康、疾病和遗传的动态标志。Nat Rev Genet，13（5）：343-3572012年5月。 [9] J Hardwick和QF Stout。最优约化等渗回归。发布前arXiv:1412.28442014年。 [10] TD Hocking和G Bourque。多样本和细胞类型中同时监督检测峰值的机器学习算法。InProc.公司。太平洋生物计算研讨会，第25卷，第367-378页，2020年。 [11] TD Hocking、G Rigaill和G Bourque。PeakSeg：计数数据中峰值检测的约束最优分割和监督惩罚学习。InProc.公司。第32届ICML，第324-332页，2015年。 [12] TD Hocking、P Goerner-Potvin、A Morin、X Shao、T Pastinen和G Bourque。使用可视标签和监督机器学习优化ChIP-seq峰值检测器。生物信息学，2016年。 [13] TD Hocking、G Rigaill、P Fearnhead和G Bourque。约束变化点检测的对数线性时间算法。预打印arXiv:1703.033522017年。 [14] TD Hocking、G Rigaill、P Fearnhead和G Bourque。用于基因组数据中约束变化点检测的广义函数修剪最优分割（GFPOP）。预印arXiv:1810.001172018年。 [15] H霍夫林。融合套索信号近似器的路径算法。《计算与图形统计杂志》，19（4）：984-1006，2010年。编号10618600。统一资源定位地址http://www.jstor.org/stable/25765384。 [16] B Jackson、JD Scargle、D Barnes、S Arabhi、A Alt、P Gioumousis、E Gwin、P Sangtrakulcharoen、L Tan和TT Tsai。一种在区间上对数据进行优化分区的算法。IEEE信号处理快报，12:105-1082005。 [17] 肖恩·朱厄尔（Sean W Jewell）、托比·迪伦·霍金（Toby Dylan Hocking）、保罗·费恩黑德（Paul Fearnhead）和丹妮拉·维滕（Daniela M Witten）。钙成像数据的快速非凸反褶积。生物统计学，2019年2月。ISSN 1465-4644。doi:10.1093/biostatistics/kxy083。统一资源定位地址https://doi.org/10.1093/生物统计学/kxy083。kxy083。 [18] NA约翰逊。基因组学问题的有效模型和算法。斯坦福大学博士论文，2011年。https://purl.stanford.edu/jq411pj0455。 [19] R Killick、P Fearnhead和IA Eckley。具有线性计算成本的最佳变化点检测。《美国统计协会杂志》，107（500）：1590-15982012·Zbl 1258.62091号 [20] 勒巴比尔。通过模型选择检测高斯过程均值中的多个变化点。信号处理，85（4）：717-7362005·Zbl 1148.94403号 [21] R Maidstone、TD Hocking、G Rigaill和P Fearnhead。关于大数据的最优多变点算法。《统计与计算》，27:519-5332017年。ISSN 1573-1375·Zbl 1505.62269号 [22] P Mair、K Hornik和J de Leeuw。R中的同位素优化：池邻违规算法（PAVA）和活动集方法。统计软件杂志，32（5）：1-242009。 [23] A Cleynen和E Lebarbier。泊松和负二项速率模型的分段：惩罚估计量。ESAIM:PS，18:750-7691014年·Zbl 1310.62041号 [24] G里格尔。优化多变化点检测的删减动态规划。预印arXiv:1004.08872010。 [25] G Rigaill公司。一种修剪动态规划算法，用于恢复具有1到kmax变化点的最佳分段。《法国社会统计杂志》，156（4），2015年·Zbl 1381.90094号 [26] G Rigaill、TD Hocking、J-P Vert和F Bach。使用最大裕度区间回归学习变点检测的稀疏惩罚。InProc.公司。第30届ICML，第172-180页，2013年。 [27] V Runge、TD Hocking、G Romano、F Afghah、P Fearnhead和G Rigaill.gfpop：一个用于单变量图约束变化点检测的R包。预发布arXiv:2002.036462020。 [28] MB Rye、P Strom和F DrablöS。一个人工管理的ChIP-seq基准表明了当前寻峰程序的改进空间。核酸研究，39:gkq11872010。 [29] MJ Schell和B Singh。简化单调回归方法。《美国统计协会杂志》，92（437）：128-1351997年·Zbl 0890.62035号 [30] AM Szalkowski和CD Schmid。芯片序列峰值调用算法的快速创新正在超越基准测试工作。生物信息学简报，12（6）：626-6332011年。 [31] EG Wilibanks和MT Facciotti。ChIP-Seq峰值检测算法性能评估。《公共科学图书馆·综合》，2010年第5卷第7期。 [32] Y-C姚明。通过Schwarz标准估计转换点的数量。《统计与概率快报》，6（3）：181-1892988年2月·Zbl 0642.62016号 [33] NR Zhang和DO Siegmund。一种改进的贝叶斯信息准则及其在比较基因组杂交数据分析中的应用。生物统计学，63:22-322007年·Zbl 1206.62174号 [34] Y Zhang、T Liu、CA Meyer、J Eeckhoute、DS Johnson、BE Bernstein、C Nusbaum、RM Myers、M Brown、W Li等。基于模型的ChIP-Seq分析（MACS）。《基因组生物学》，9（9）：R1372008。 [35] 补充图1：背景水平随样本McGill0028而变化●●●●●●●●●●●● ●●●●●●●●●●●● ●●●●●●●●●●●● 麦吉尔0009●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0027●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0011●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0024●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0107●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0029●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●● 麦吉尔0322●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0102●●● ●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0007●●●●●●●●● ●●●● ●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0005●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 麦吉尔0004●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● [36] 而不是每个背景段的新段平均值（如我们提出的模型）。我们 [37] 在同一ChIP的不同背景/噪声区域中观察到不同的平均水平 [38] seq样本，这表明统一背景/噪声平均模型不会 [39] 在这些数据中是合适的。左：我们计算了每个样本和基因组的平均覆盖率 [40] 无峰值标签区域（这意味着生物学家观察到没有峰值 [41] 在该区域，因此它必须仅包含背景噪声），并观察到平均值为 [42] 基因组各区域之间差异很大。例如，在示例McGill0028中 [43] 背景/噪声区域的平均覆盖范围为1.78至6.76，这表明 [44] 统一/恒定平均值不是这些数据的良好模型。右：对于一个样本， [45] 我们计算了整个基因组的峰值，然后计算了每个峰值的平均覆盖率 [46] 附近背景的平均覆盖率。我们观察到附近背景的平均值 [47] 随着峰值平均值的增加而增加，这表明均匀/恒定平均值 [48] 不是这些数据的好模型。29 [49] 补充图2：CDPA失败示例的详细信息 [50] 考虑计算6个数据点集的5段上下约束模型 [51] y=[3,9,18,15,20,2]，使用泊松损失'（yi，m）=m−yilogm无约束PDPA计算模型=[3,9,16.5,16.5,20,2]，其中 [52] 如下所示。GPDPA计算的第一个成本函数是第一个数据的泊松损失 [53] 点：C1,1（u1）='（3，u1）=u1−3 logu1（26） [54] C1的最小值为3，因此它的minless算子对于u2≤3是凸的，对于 [55] u2≥3：( [56] 第二个成本函数是前两个数据点的总泊松损失：C1,2（u1）='（9，u1）+C1,1（u1 [57] C1，2的最小值是6，所以它的minless算子对于u2≤6是凸的，对于 [58] u2≥6：( [59] 截至数据点2的两个部分中的最佳成本是：( [60] 注意，在上述等式中，数据点1和2之间的非递减变化为 [61] 由min-less运算符C1≤，1强制执行。截至数据点3的2个分段的最佳成本定义为： [62] atu2≈15.41，因此在数据点3之前的两个段中的最佳成本简化为：( [63] 说明：上图中的成本函数如下图中的黑色曲线所示。该部分 [64] 左边（u2≤15.41）的蓝色是第一个数据后非递减变化的成本 [65] 点（t1=1）。右侧紫色部分（u2≥15.41）是非递减成本 [66] 在第二个数据点（t2=2）之后进行更改-43.2 [67] 平均参数u3（如上图中的灰色函数所示）： [68] 因为最佳成本是针对可能的变化点和所有可能的平均值计算的 [69] 值，则GPDPA能够计算最优解（出现在atu3=15，u2=31 [70] 18，u1=6，t1=2，但在算法计算所有 [71] 数据点C5,6）。相反，CDPA计算非递减的标量最小成本 [72] 第一个数据点后发生变化（上图中的红色圆圈，u2=13.5，u1=3，t1=1）， [73] 并放弃在第二个数据点（其中 [74] 最终达到最佳）。因此，CDPA是一种贪婪算法。32 [75] 补充图3：峰值检测精度的ROC曲线 [76] 我们想将我们提出的算法与其他算法的峰值检测精度进行比较 [77] 来自生物信息学文献，该文献通常使用 [78] 默认参数设置。因此在典型分析中，要控制假阳性 [79] 根据 [80] 某种可能性或显著性阈值。例如，MACS算法使用q值 [81] 阈值参数，HMCanBroad使用finalThreshold参数（更高的阈值 [82] 导致更多误报）。描述：为了解释评估中的这个修剪步骤，我们使用了Receiver Operat [83] ing特征（ROC）曲线分析。对于每个阈值参数，我们计算了 [84] 使用标签的假阳性率和真阳性率会导致 [85] ROC曲线。曲线下面积（AUC）通过改变阈值参数来计算 [86] eter在其整个范围内（从具有许多假阳性和真阳性的完整峰值列表， [87] FPR=TPR=0）的完全修剪/空峰列表。请注意，即使使用 [88] 每个模型中预测的峰值数量最多，并非所有标签都达到了最大值 [89] TP和FP，因为最大峰值列表不一定预测所有标签中的峰值。 [90] 因此，我们将每个ROC曲线线性外推至TPR=FPR=1，以便计算 [91] AUC（下图中的虚线）。1.0 0.8 [92] 不一定是等级制度。例如，GPDPA计算从0到 [93] 每个问题有9个峰值；具有1个峰值的最佳模型中的峰值可能不是 [94] 存在于具有2个或更多峰值的最优模型中。 [95] 补充图4：惩罚之间的测试AUC比较 [96] 功能CDPA（以前最好的）GPDPA（建议的）PDPA（无约束的） [97] 可视化算法之间是否存在任何显著差异。每个点 [98] 表示单个列车/测试分段的精度。显示平均差值 [99] 以及配对t检验的p值。我们建议使用两种方法来选择变更点模型中的段数 [100] 模型复杂度函数C。oracle模型的复杂性相对复杂 [101] 统计论据驱动的表达；线性模型复杂度的简单度量 [102] 段数。在本节中，我们从测试AUC的角度比较了这些方法。我们考虑学习由常数惩罚组成的单参数惩罚函数 [103] λ、 给定线性或oracle模型复杂度函数。我们期待神谕 [104] 由于其统计动机，模型复杂性导致测试AUC更高。这个 [105] 补充图4的第二行绘制了一个测试AUC值的分布 [106] 模型复杂性函数。与我们的预期相反，对于这三个人 [107] 变化点算法（CDPA、GPDPA、PDPA）、线性和预言惩罚显示 [108] 测试AUC无显著差异（0.00004至0.0006的平均差异，p值>0.05 [109] 成对试验27）。这些数据表明，尽管统计数据表明 [110] oracle模型的复杂性，并不比简单线性模型的复杂性更精确 [111] 用于实际基因组数据中的峰值检测。我们还比较了使用一个或多个参数学习惩罚函数， [112] 给出了线性模型的复杂度函数。我们预计罚款的测试AUC会更高 [113] 具有多个参数的函数。补充图4的第一行绘制了34 [114] 一个惩罚函数的测试AUC值相对于另一个的分布。同意 [115] 按照我们的预期，对于所有三种转换点算法，多参数惩罚 [116] 函数的测试AUC稍大（平均差值为0.002到0.006，p值<0.07）。 [117] 这些数据表明，最好学习多参数惩罚函数 [118] 用于基因组数据中的精确峰值检测。35 [119] 补充图5：峰值模型惩罚/阈值训练方法的测试准确性|无监督|监督 [120] 每个面板显示七个ChIP-seq数据集中的一个，按模式/实验（Broad [121] H3K36me3）、标签（AM）、细胞类型（免疫）和标签数量。很明显 [122] 监督模型通常比无监督和多参数模型更准确 [123] 模型通常比单参数模型更准确。 [124] 8.实施细节 [125] 8.1. GPDPA伪码 [126] 在本节中，我们为我们提出的广义剪枝给出了更详细的伪码 [127] 动态规划算法（GPDPA）。我们提出以下数据结构 [128] 和计算子程序：•FunctionPiece：表示一段aCk，t（u）成本函数的数据结构 [129] 上述数据结构和子例程用于以下伪代码 [130] 描述了用K解决上下约束转换点问题的GPDPA [131] 部分。 [132] 算法2广义剪枝动态规划算法（GPDPA）1：输入：数据集∈Rn，最大分段数K∈{2，…，n}。2：输出：最优段均值矩阵U∈RK×K和ends∈{1，…，n}K×K3：计算最小值和最大值。4:C1,1←OnePiece（y1，y，y）5：用于2吨起的数据点：6:C1，t←OnePiece（yt，y，y）+C1，t−1 7：用于2到K的分段：用于数据点起：//动态编程8:MinLessOrMore←MinLess if kis even，else MinMore 9:min prev←MinLessOrMore（t−1，Ck−1，t−1） [133] 算法是Ck，t，应该初始化为空FunctionPieceList的aK×narray [134] 物体。计算第4-6行的alltoccurs的C1，t。动态编程更新发生在第7-11行的for循环中。第9行使用 [135] 计算临时FunctionPieceList的MinLess（或MinMore）子程序min prev≤≥ [136] （表示函数Ck−1、t−1或Ck−）。第10行设置临时功能 [137] PieceList最小新到唯一可能变更点的成本ift=k；否则，它使用 [138] MinOfTwo子例程，用于计算每个可能的最佳变更点的成本 [139] 平均值。第11行添加了数据点的成本，并存储生成的FunctionPieceList [140] inCk，t.最佳段的解码意味着U（aK×实数Karray）和 [141] endT（aK×Karray of integers）变量出现在第12-18行的for循环中。对于 [142] 给定型号sizek，解码从第13行开始，使用ArgMin子例程求解 [143] u∗=arg minuCk，n（u）（前一段endt0和meanu0的最佳值为37 [144] 也已返回）。现在我们知道u*是最后（k）段的最佳平均值 [145] 从数据点0+1吨发生。这些值存储在Uk、kandTk、k（第14行）中。而且 [146] 我们已经知道分段k−1的最佳平均值为u0。请注意，0=∞标志 [147] 表示相等约束处于活动状态（第16行）。其他段的解码 [148] s<k使用FindMean子例程继续（第17行）。这需要最好的成本 [149] 模型插入到数据点0，找到存储u*成本的FunctionPiece， [150] 并返回上一段endt0和meanu0的新最佳值。平均值 [151] 段-1的末尾存储在Tk，s（第18行）。算法2的时间复杂性为O（KnI），其中MinLess的复杂性为 [152] 和MinOfTwo子例程，它在间隔数上是线性的（FunctionPiece [153] 对象），用于表示成本函数。有病理合成数据 [154] 间隔数I=O（n）的集合，这意味着时间复杂度为O（Kn2）。 [155] 然而，实际数据集中的平均间隔数在经验上是I=O（logn），因此 [156] 算法2的平均时间复杂度为O（Knlogn）。 [157] 8.2. Min-less算法 [158] 以下子程序用于实现MinLess子程序。•GetCost（p，u）：采用FunctionPiece对象p和平均值的算法 [159] 类似。38 [160] 算法3MinLess算法。1：输入：前一段endtprev（整数）和fin（FunctionPieceList）。2：输出：FunctionPieceListfout，初始化为空列表。3：前期成本←∞4：新下限←下限（fin[0]）。5:i←0；//从左侧的FunctionPiece开始6:whilei<长度（fin）：//继续，直到右侧的FunctionalPiece 7:FunctionPiecep←fin[i]8:if prev cost=∞：//在此间隔中查找min。9：候选平均值←OptimalMean（p） [161] 操作员。该算法初始化prev cost（第3行），这是一个使用的状态变量 [162] 在第8行上决定算法应该查找局部最小值还是交集 [163] 以有限的成本。由于prev成本最初设置为∞，因此算法首先如下 [164] 凸函数从左到右分段，直到找到局部极小值。如果没有最小值 [165] 在给定的凸函数片（第15行）中，只需将其推到 [166] 新的FunctionPieceList（第16行）。如果最小值出现在间隔内（第10行） [167] 存储成本和平均值（第11-12行），并使用 [168] 上限以该平均值结束（第16行）。然后算法开始寻找 [169] 另一个具有相同成本的FunctionPiece，通过计算凸的较小根 [170] 损失函数（第20行）。如果发现FunctionPiece在区间（第21行）中具有根， [171] 按下一个新的常量FunctionPiece（第22行），算法继续搜索 [172] 最低要求。在算法的末尾，如有必要，将推送常量FunctionPiece 39 [173] （第28行）。该算法的复杂度为O（I），其中是FunctionPiece的数量 [174] 对象内嵌。实现min-more运算符的算法是类似的。而不是 [175] 从左到右搜索，它从右到左搜索。而不是使用小根 [176] （第21行），它使用大根。 [177] 8.3. 实施细节 [178] 我们发现一些重要的实施细节： [179] 对于包含重复的数据序列的权重，使用数据的游程长度编码和相应的损失函数在计算上是有利的。例如，如果数据序列5,1,1,0,0,5,5被编码为n=4 countsyt5,1,0.5，并具有相应的权重swt1,3,2,2，则平均µ的泊松损失函数为'（yt，wt，µ）=wt（µ−ytlogµ）。 [180] 平均成本文本将Ck，t函数定义为总成本。然而，对于非常大的数据集，成本值将非常大，从而导致数值不稳定。为了克服这个问题，我们使用平均成本实现了更新规则。对于权重Wt=Pti=1wi，计算平均成本的更新规则为 [181] 对数（平均）空间中的区间对于非负计数数据yt∈{0,1,2，…}的泊松模型，没有可能的平均µ值小于0。因此，我们使用log（µ）值在FunctionPiece对象中实现间隔。例如，我们存储的不是Φ∈[0,1]，而是对数Φ∈[-∞，0]。 [182] 寻根对于泊松损失的ComputeRoots子例程，我们使用牛顿寻根。对于较大的根，我们求解logµ+bµ+c=0（线性为µ→∞），对于较小的根，则求解ax+bex+c=O（x=log¦Μ，线性为asx→−∞和µ→0）。当成本接近零（绝对成本值小于10−12）时，我们将停止查找根。 此参考列表基于出版商或数字数学图书馆提供的信息。其项与zbMATH标识符进行启发式匹配，可能包含数据转换错误。在某些情况下，zbMATH Open的数据对这些数据进行了补充/增强。这试图尽可能准确地反映原始论文中列出的参考文献，而不要求完整或完全匹配。