论文阅读(十六)：利用线性链条件随机场模型检测阵列比较基因组杂交数据的拷贝数变异

1.论文链接：Detection of Copy Number Variations from Array Comparative Genomic Hybridization Data Using Linear-chain Conditional Random Field Models

摘要：

拷贝数变异（CNV）约占人类基因组的12%。除了CNVs在癌症发展中的固有作用外，它们还被报道是复杂疾病易感性的基础。每个变异的范围可能从大约1000个核苷酸到不到5兆碱基不等。阵列比较基因组杂交（aCGH）允许跨基因组拷贝数的改变。利用aCGH数据分析cnv的关键计算挑战是检测拷贝数变化的段边界，并推断每个段的拷贝数状态。马尔可夫随机和更具体的条件随机为数据预处理、分割和拷贝数状态解码提供了一个框架。

关键词：拷贝数变化，aCGH，条件随机场

为了更多地了解可遗传拷贝数多态性（CNPs）在确定疾病表型中的作用，目前正在对CNPs进行系统的定位和编目。阵列比较基因组杂交（aCGH）允许跨基因组拷贝数的改变。利用aCGH数据分析拷贝数变异（CNVs）的关键计算挑战是检测拷贝数变化的段边界和推断每个段的拷贝数状态。本章提出了一种基于条件随机算法的新型无向图形模型，用于aCGH数据中的CNV检测。该模型将数据预处理、分割和拷贝数状态解码结合到一个统一的框架中。所描述的方法被称为CRF-CNV，在定义有意义的特征函数中提供了极大的存在性。因此，它可以将任意大小的局部空间信息生态型地整合到模型中。对于模型参数的估计，采用共轭梯度（CG）方法对似然值进行了优化，并在CG框架内开发了正向/反向算法。

16.1介绍

DNA序列中的结构变异，如可遗传的拷贝数改变，已被报道与许多疾病有关。还观察到体细胞染色体畸变（即，肿瘤样本中的基因扩增和缺失）在不同的癌症类型或亚型中显示出不同的临床或病理特征[3，6，14]。为了更好地了解遗传拷贝数多态性（CNPs）在确定疾病表型中的作用，需要对CNPs进行系统的定位和分类。在癌症样本中体细胞拷贝数畸变的鉴定可能导致重要癌基因或肿瘤抑制基因的发现。由于现有技术在评估拷贝数变异（CNV）方面的显著能力，最近研究界对研究可遗传的以及体细胞CNV产生了极大的兴趣[1，3，4，6，14，16，23]。

一般来说，CNV检测基本上有三种技术平台：基于阵列的技术（包括阵列比较基因组杂交（aCGH）[13，20]，以及许多其他变体），单核苷酸多态性（SNP）基因分型技术[1，14]和下一代测序技术[2]。基于阵列的技术测量疾病/测试样品相对于正常/参考样品的DNA拷贝数改变。原则上，对于每个基因（或基因组片段），每个个体从其父亲继承一个拷贝，从其母亲继承一个拷贝。副本总数为两份。然而，当个体中存在拷贝数突变时，拷贝总数可能是1（即，缺失），或三个或更多。对于每个DNA片段或“克隆”，aCGH可以间接测量通过荧光强度的拷贝数，通过荧光强度的拷贝数，它代表与片段杂交的DNA的丰度。测试和参考样品的强度比，以对数标度测量，预计与测试样品相对于参考样品（假设有两个拷贝）的拷贝数变化成比例，尽管在该过程中可能会从各种来源引入显著的噪声。基于阵列的技术主要用于大片段的扩增和缺失，尽管不同的实验平台和使用不同大小克隆的设计可能会产生非常不同的分辨率和基因组覆盖率[14]。其他两个平台，SNP基因分型和下一代测序，使用不同的技术来测量拷贝数的变化，并在近年来得到普及[2，17]。本章主要介绍aCGH数据分析。

不同的平台有不同的挑战。自然，人们应该通过利用不同数据集的特殊属性来为这些平台使用不同的方法。毫不奇怪，近年来已经提出了各种算法用于不同的数据。另一方面，所有这些研究的主要目标是鉴定共享相同拷贝数的连续片段集。因此，处理来自不同平台的数据的基本计算任务是相同的：将基因组分割成具有相同拷贝数的离散区域。完成上述分割步骤后，将为每个片段分配一个拷贝数。后一项任务被称为“分类”。对来自不同平台的数据进行分割的一个重要共性是克隆之间的空间相关性。许多现有的方法已经利用了这样的属性，通过利用相同的方法，隐马尔可夫模型（HMM），它可以方便地使用链结构建模空间依赖。结果表明，HESTERN的初步成功[1，4，16，23]。然而，所有这些障碍都有一个继承的限制，即，它们都是一阶Hynst，不能考虑长程相关性。

已经证明，在各种领域，条件随机ELD（CRF）始终优于HALTH，主要是因为CRF可以潜在地整合数据中的所有信息[21]。这一特性使得CRF对CNV模型数据特别有吸引力，因为人们可以使用来自一个区域的数据来定义特征函数，而不是使用单个或两个数据点来定义Hyndrome中的排放和转换。我们提出了一种基于CRF的新型无向图形模型，用于aCGH数据的CNV检测[25]。我们的模型有效地结合了数据预处理，分割和拷贝数状态解码到一个统一的艾德框架。我们的方法（称为CRF-CNV）在定义有意义的特征函数方面提供了很大的灵活性，因此它可以有效地将任意大小的局部空间信息集成到模型中。对于模型参数估计，我们采用了共轭梯度（CG）方法进行似然优化，并在CG框架内开发了有效的前向/后向算法。拷贝数解码的最后一步是基于Viterbi算法。该方法使用真实的数据与已知的拷贝数以及模拟数据与现实的假设进行评估，并与两个流行的公开可用的程序进行比较。实验结果表明，CRF-CNV优于基于贝叶斯隐马尔可夫模型的方法在两个数据集的拷贝数分配。当与非参数方法相比时，CRF-CNV在真实的数据上实现了更高的精确度，同时保持了相同的召回水平，并且当应用于模拟数据时，两种方法的性能相当。

本章的其余部分组织如下。在第16.2节中，我们简要概述了aCGH数据和从aCGH数据中检测CNV的现有方法。关于CNV的CRF模型开发和实施的详细信息见第16.3节。我们在两个数据集上的实验结果以及与其他程序的比较见16.4节。我们将在16.5节进行一些讨论，以此结束本章。

16.2 aCGH数据和分析

16.2.1 aCGH数据

虽然从理论上讲，我们的方法可以应用于不同实验平台的数据，但在本分析中，我们主要关注aCGH数据。在数学上，aCGH数据通常由一组克隆的log 2强度比的阵列以及每个克隆沿基因组沿着的物理位置组成。这里，每个克隆是阵列中的一个数据点，所有克隆根据它们在基因组上的物理位置排序。图16.1绘制了人癌细胞系和Snijders等人分析的正常参考DNA [18]。每个数据点代表一个基因组区段/克隆，y轴代表标准化的log 2强度比。基于aCGH的CNV检测的主要目标是将基因组分割成共享相同平均log 2比率模式的离散区域（即，具有相同的拷贝数）。在来自人类的正常样本中，所有常染色体的拷贝数通常为2。理想地，如果测试样品（例如，一个癌细胞系）也有两个DNA拷贝，如果它有一个单拷贝的增加（或单拷贝的丢失），那么这个值应该是0.585（或-1），等等。图16.1中的虚线从下往上分别表示一个、两个和三个拷贝的三个值。一种直接的方法是使用一些全局阈值来基于其log 2比率为每个克隆分配拷贝数。然而，如图16.1所示，aCGH数据可能非常嘈杂，不同片段之间的边界模糊。它还可能具有复杂的局部空间依赖结构。这些属性使得分割问题本质上很难。使用全球阈值的方法在实践中通常效果不佳。

16.2.2现有算法

通常，需要许多步骤来检测aCGH数据的拷贝数变化。首先，原始log2比率数据通常需要一些预处理，包括归一化和平滑。标准化是减少实验因素引起的系统误差的必要步骤。通常，输入数据通过利用测试和参考DNA中没有差异的一些对照数据集进行归一化。标准化的目的是使对照数据集中的中位数或平均log2比值为0。平滑用于减少由于随机误差或突变而产生的噪声。平滑方法通常使用滑动窗口来过滤数据，试图在处理突然变化和减少随机误差的同时对数据进行曲线处理。除了基于滑动窗口的方法之外，还提出了几种其他技术，包括分位数平滑，小波平滑等（参见[5，8]）。

分析aCGH数据的第二步称为分段，其目的是识别具有相同平均log 2比的分段。概括地说，有两个相关的估计问题。一是推断变异的数量和统计意义;二是准确定位变异的边界。因此，针对这两个估计问题，提出了几种不同的算法。Olshen等人[12]提出了一种基于递归循环二进制分割（CBS）算法的非参数方法。Hupe等人[9]提出了一种称为GLAD的方法，该方法基于中值绝对偏差模型，将离群值（数值上与其余数据相距甚远的观测值）与其周围的分段分离。Willenbrock和Fridlyand [24]使用真实的模拟模型比较了CBS（在DNACopy中实现）和GLAD的性能，得出的结论是，CBS总体上优于GLAD。在我们的实验研究中，我们采用了文献[24]中的模拟模型。在获得分割结果之后，需要后处理步骤来联合收割机具有相似平均水平的分割，并将其分类为单拷贝增益、单拷贝丢失、正常、多拷贝增益等。诸如GLADMerge [9]和MergeLevels [24]的方法可以对分割结果进行后处理，并相应地对它们进行标记。

正如Willenbrock和Fridlyand [24]所指出的，同时执行分割和分类更可取。合并这两个步骤的一个简单方法是使用线性链HMM。潜在的隐藏状态是真实的拷贝数。给定一个状态，log 2比率可以使用高斯分布来建模。从一种状态到另一种状态的转变揭示了相邻克隆之间拷贝数变化的可能性。给定观察到的数据，标准算法（前向/后向，Baum-Welch和Viterbi）可以用来估计参数和解码隐藏状态。近年来，针对aCGH数据提出了几种HSPs变体[7，16]。Lai等人。[11]已经表明，在给定足够的信噪比的情况下，Hynomial对小像差的表现最好。Guha等人。[7]提出了一种贝叶斯HMM，可以对参数施加生物学意义的先验。Shah等人。[16]通过引入离群值和位置特异性先验，扩展了贝叶斯HMM，可用于对遗传拷贝数多态性进行建模。

请注意，所有这些模型都是一阶Hynomial，无法捕获长期依赖性。直观地说，考虑高阶HHT来捕获信息局部相关性是有意义的，这是从aCGH数据观察到的重要属性。然而，考虑更高的阶数将使HISTORY更加复杂和计算密集。

16.3 aCGH数据的线性链CRF模型

为了克服障碍的局限性，我们提出了一个基于线性链条件随机场（CRF）理论的新模型[10，21]。CRF是无向图模型，设计用于计算给定输入变量X的输出随机变量Y的条件分布。术语“随机场”，相当于图论中的“马尔可夫网络”，已广泛应用于统计物理和计算机视觉。CRF也被称为条件马尔可夫网络[22]。因为CRF是条件模型，变量X之间的依赖关系不需要明确指定艾德。因此，可以通过使用输入变量X来定义有意义的特征函数，这些特征函数可以有效地捕获局部空间依赖性。CRFs已被广泛应用于语言处理、计算机视觉和生物信息学，与包括HPLF在内的有向图模型相比，具有显著的性能。

一般来说，线性链CRF（见图16.2）被定义为条件分布。

其中：

16.3.1特征函数

16.3.2参数估计

边缘分布可以通过组合前向和后向变量获得：

16.3.3评价方法

16.4实验结果

16.4.1真实样本

16.4.2模拟数据

尽管真实数据的结果显示，当匹配范围指数放宽到2时，CRF-CNV的性能优于CBS，并且与CNA-HMMer相当，但由于数据集规模非常小，实验是有限的。为了进一步评估CRF-CNV的性能，我们使用从文献[24]获得的模拟数据集测试了这三种算法。数据集由500个样本组成，每个样本有20条染色体。每条染色体包含100个克隆。每个克隆可能受到六种可能的拷贝数状态的影响。作者通过从原发性乳腺肿瘤组织的数据库中采样片段，并使用几种机制（例如，每个样本中癌细胞的比例，给定拷贝数状态下强度值的变化）来控制噪声水平，从而生成这些样本。通过使用文献中的模拟数据，我们进一步评估了CRF-CNV的性能。

为了训练CRF-CNV，我们像往常一样将500个样本分为三组。这次，训练集第1组包含样本1-50，验证集第2组包含样本51-100，测试集第3组包含样本101-500。我们使用了之前讨论的相同网格搜索方法来获取超参数{aj}、u和σ^2。对于每组固定的超参数，我们使用CG方法来获取参数θ。最后，我们使用Viterbi算法来解码第3组样本中可能的隐藏拷贝数状态标签，并与CNA-HMMer和CBS的结果进行比较。此外，我们还比较了CRF-CNV对第2组和第3组的预测，以查看在新的测试数据下，我们的模型可能会因为从小样本中推断出的次优参数而遭受多大的性能下降。为了便于比较，CBS和CNA-HMMer的结果分别针对这两组进行了展示。我们还使用第1组作为训练数据，为CNA-HMMer分配信息性先验。

表16.3显示了CRF-CNV、CBS和CNA-HMMer分别预测的第2组和第3组的段数总和，并将其与已知的段数总和进行了比较。有趣的是，对于这个模拟数据，CBS和CNA-HMMer预测的段数都较少。CRF-CNV在第2组预测的段数较少，而在第3组预测的段数较多。然而，段数并不能提供完整的情况。因此，我们使用F-measure来检查每种方法在第2组和第3组的边界预测准确性。表16.4显示了不同方法、不同组别和不同匹配范围的F-measure。正如预期的那样，随着D从0增加到4，所有方法和两个数据组的F-measure都有所增加。CBS在第2组和第3组的结果一致，CNA-HMMer在第2组和第3组的结果也一致，这并不奇怪。有趣的是，CRF-CNV在第3组的性能也非常接近其在第2组的性能。这个特性是理想的，因为它说明了CRF-CNV的稳健性。新测试数据的性能几乎与用于选择最佳模型参数的验证数据的性能相同。请注意，训练数据和验证数据的大小也非常小。可以预期，在小规模训练数据集的情况下，我们的方法可以可靠地预测在相同实验条件下生成的其他数据。就三种方法的性能而言，CNA-HMMer比CRF-CNV更准确，而CBS在完全匹配的情况下表现最差。然而，当我们通过增加D的值来放宽匹配标准时，CBS和CRF-CNV都比CNA-HMMer表现得更好。CNA-HMMer和CRF-CNV的结果与真实数据的结果一致。CBS的性能明显高于真实数据观察到的性能。然而，这可能归因于模拟过程，因为CBS被用来分割来自原发性乳腺肿瘤数据集的145个样本[24]。

16.5结论

近年来，检测拷贝数变异的问题引起了广泛关注，并提出了许多计算方法来解决这一问题。在这些计算发展中，CBS获得了极大的流行，并且已被证明在模拟数据上通常比其他算法表现更好[24]。然而，正如原始论文（以及我们的实验重新发现的）所示，CBS在标准的Coriell数据集中通过光谱核型分析识别的拷贝数变化上报告了更多的假阳性[18]。另一种常用的分割技术是隐马尔可夫模型（HMM）。HMM方法的优势在于能够在一个框架内执行参数（即均值和方差）估计和拷贝数解码，并且随着观察次数的增加，其性能有望提高。然而，几乎所有专门用于aCGH的HMM都是一阶马尔可夫模型，因此无法在数据中纳入长距离空间相关性。

我们提出了一种基于条件随机场理论的新计算模型。我们的模型在定义任意邻域上的有意义特征函数方面提供了极大的灵活性。因此，可以纳入更多的局部信息，并期望得到更稳健的结果。在本章中，我们使用能够将平滑有效地纳入我们统一框架的稳健充分统计量定义了特征函数。我们还开发了在共轭梯度方法内的有效前向-后向算法，用于高效计算模型参数。我们使用真实数据以及模拟数据评估了我们的方法，结果表明，当允许小的偏移时，我们的方法在两个数据集上的表现都优于贝叶斯HMM。与CBS相比，我们的方法在真实数据集上的假阳性更少。在模拟数据集上，我们方法的性能与CBS相当，CBS已被证明是三种流行的分割方法中最好的。

与其他条件随机场（CRF）类似，为了训练我们的模型，必须依赖一些训练数据。为了实用，像CNA-HMMer这样的贝叶斯隐马尔可夫模型（HMM）也需要训练数据来正确分配信息性先验。我们认为这个问题并不像最初看起来那么严重，主要有两个原因。首先，正如我们的实验所示，我们的算法确实非常稳健，即使找不到模型参数的最佳估计，也能持续表现良好。例如，我们在模拟数据集的分析中使用了一种简化的程序，通过随机选择一个子集进行训练。从理论上讲，从这种程序估计的参数可能严重依赖于这个特定子集，并且不一定全局最优。然而，表16.4中显示的结果表明，新测试数据的性能几乎与用于调整参数的验证数据相同。目前正在进行10折交叉验证以进一步验证这一观察结果。此外，我们的算法所需的训练规模非常小，这从真实数据和模拟数据中都可以看出。我们的算法在Coriell数据上报告的段数没有假阳性的主要原因是我们通过训练学习了数据中的结构。CNA-HMMer有无训练数据时性能的显著差异也支持这一观察结果。其次，确实存在一些经过实验验证的数据（如Coriell数据集）。随着时间的推移，预计会有越来越多的此类数据，因为最终，对于任何预测方法，都必须通过实验验证其中的一些预测。此外，我们方法的稳健性表明，我们可能只需要为每个特定平台训练我们的算法。然后可以将这些参数用于同一平台上生成的未来数据。

在计算成本方面，CRF-CNV有两个独立的部分：训练时间和预测时间。训练需要密集的计算来优化对数似然并确定超参数。此外，还可以执行k折交叉验证，这将需要更多的计算时间。相比之下，一旦参数被估计，预测阶段就相当高效。幸运的是，我们算法的训练阶段只需要一个小数据集，这使得算法仍然具有实用性。我们的算法的可能扩展和应用可以应用于其他高通量技术以检测拷贝数变化。