ACNE：一种估计Affymetrix SNP阵列等位基因特异性拷贝数的总结方法

摘要

动机：目前用于估计DNA拷贝数（CN）的算法借鉴了基因表达分析方法的概念。然而，单核苷酸多态性（SNP）阵列具有特殊的特征，如果将其考虑在内，可以提高整体性能。例如，等位基因之间的交叉杂交发生在SNP探针对中。此外，目前大多数CN方法都侧重于总CN，而研究表明，等位基因特异性CN对一些研究至关重要。因此，我们开发了一种估算高质量等位基因特异性CN的总结方法。

结果：该方法估计所有Affymetrix SNP阵列的等位基因特异性DNA CNs，直接处理SNP探针内探针之间的交叉杂交。该算法在计算DNA CN方面优于（或至少性能与）其他最先进的算法。它能更好地区分正常状态下的畸变，并能提供更精确的等位基因特异性CN。

可利用性：该方法在开源R软件包ACNE中可用，该软件包还包括对aroma.affmetrix框架的附加组件(网址：http://www.aroma-project.org/).

联系人： arubio@ceit.es

补充信息： 补充数据可在生物信息学在线。

1简介

基因组畸变参与不同疾病的发病机制，尤其是癌症（平克尔等。,1998; 波拉克等。,1999)。这些基因组异常可能导致细胞癌变，尤其是当抑癌基因或癌基因受到影响时。

DNA拷贝数畸变（CNA）是指基因组畸变，包括基因组某一部分（整个染色体、臂或片段）的扩增或缺失。CNA中的DNA拷贝数（CN）可能比正常状态（CN等于2）更大（扩增）或更小（缺失和纯合缺失）。杂合性丢失（LOH）是另一种基因组畸变，发生在缺少两条亲本染色体之一的基因组片段中。对于LOH区中的每个SNP，其中一个等位基因的CN等于零。LOH区域的CN可能是一条（失去一条亲本染色体），两条（复制中性LOH），甚至大于两条。在后两种情况下，具有LOH的区域的附加副本已由使用当前副本的单元重新生成。

单核苷酸多态性（SNP）阵列可用于研究CNA和LOH病例。虽然这些阵列的最初应用是基因分型，但它们也可以用于估计被询问SNP所在位点的拷贝数。还有一些基因组畸变并不影响总的CN或基因型，例如易位和逆转。使用SNP阵列无法识别这种畸变。

一些公司已经开发了这种阵列，但我们将重点关注Affymetrix（Affymetix Inc。，2009)平台。Affymetrix有一个从10K阵列开始的SNP阵列家族，它询问10000个SNP。随后，他们开发了100K和500K芯片组，均使用两种不同的限制性内切酶。较新的全基因组人类SNP阵列5.0（GWS5）和全基因组人SNP阵列6.0（GWS6）阵列也使用了两种酶，但它们在不同特征上与前人有所不同。GWS阵列除了SNP探针外，还有非多态性探针，用于研究DNA CN变异（CNV）并覆盖没有SNP的基因组部分。GWS5的SNP数量与500K和420000个非多态探针相同，而GWS6的SNP约为900 000个，非多态探针略少于950000个。每个芯片中的探针被分组为对应于单个SNP的探针。在同一SNP中与两个等位基因（A和B）互补的探针称为探针对。阵列中的探针与SNP的两条链互补。每个探针的长度为25 nt。10–500K阵列中的探针由完美匹配（PM）和错配探针（MM）组成，其中MM与中央核苷酸中的PM不同。为了给更多PM留出空间，GWS阵列没有MM。

提出了几种从SNP阵列中获取DNA CN信息的处理方法。使用的第一种算法直接继承自基因表达分析，如dChip（Li和Hung Wong，2001)和RMA（爱尔兰等。,2003)。还有一些特定于CN数据的最新数据，例如CNAG（Nannya等。,2005)、PLASQ（LaFramboise等。,2007)，斜体（Rigaill等。,2008)、CRMA（本特森等。,2008年b)，CN5（Affymetrix公司。，2008)和CRMA v2（本特森等。,2009).

SNP阵列的整个处理流程由几种低层方法组成，即背景去除、归一化（均衡多个阵列的信号水平）、汇总（提取与每个等位基因的CN成比例的信号）、，基因组后处理（处理与探针所在序列的大小和其他属性相关的信号偏差）和分割（识别具有相同畸变的基因组相邻区域）。在这里，我们重点关注SNP阵列低水平分析的探针汇总步骤，并提出等位基因特异性CN估计（ACNE）。ACNE是一种基于非负矩阵分解（NMF）的多样本摘要方法，它直接处理交叉杂交问题，从而得到等位基因特异性CNs（ASCN）。ACNE使用CRMA v2在背景和偏移消除后给出的数据，并可用于替代CRMA v2中已经实现的问题汇总模型。

使用SNP阵列可以获得ASCN，即任何SNP的每个等位基因的拷贝数。ASCN本身在生物学上与某些疾病相关（达菲等。,2008)。ASCN对于识别存在LOH的区域也很重要。LOH被证明是重要的，因为它有助于识别基因组中使肿瘤抑制因子失活的异常区域。如图所示，当正常组织受到一些污染，即肿瘤不是100%纯时，ASCN也有助于解释CN图。

上述一些方法提供了ASCN（dChip和CRMA v2的摘要方法的简要描述包含在补充材料)。然而，结果不如总CN精确。性能较差的主要原因是每个等位基因信号之间的串扰。反过来，这是因为交叉杂交。尽管交叉杂交（由于与不同于其靶基因序列的杂交而添加到一个探针的信号）发生在所有芯片类型上，无论潜在靶基因来自何处，它都因阵列类型而异。在SNP阵列中，等位基因之间存在交叉杂交，因为探针对内的探针几乎相同（它们只在与SNP对应的核苷酸上不同）。两个等位基因信号之间的特殊约束可用于计算和控制这种特殊串扰。例如，表达式数组不可能有类似的约束。CRMA v2执行全局串扰校准。在这里，我们建议在每个SNP的基础上纠正交叉杂交。这种方法，如所示第3节，提供了对总CN和等位基因特异CN的更可靠估计。

2方法

如前所述，MMs不包括在最新的Affymetrix SNP阵列中，此处不考虑它们。因此，我们在处理问题时仅提及PM。

2.1数据

这里使用的数据（CEL文件）来自国际HapMap项目和两个不同的公开数据集。第一个数据集由来自国际HapMap项目（Altshuler）的31个样本组成等。,2005; 国际HapMap联盟，2003)在Mapping250K_Nsp上杂交（在补充材料).

第二个数据集由一项乳腺癌研究的24个肿瘤样本组成的一个子组（哈弗蒂等。,2008)（数据保存在NCBI-GEO中，注册号为GSM182833-44、GSM182848-59和GSE7545）。这些阵列在Mapping250K_Nsp上杂交。

最后，第三个数据集来自前列腺癌研究，它有74个肿瘤样本在GWS6上杂交（刘等。,2009)（GSM37452-44、GSM37452-609和GSE14996）。

初始低级步骤（背景和偏移量删除）是使用CRMA v2中的方法（使用默认参数）执行的。它们包括等位基因之间偏移和全局串扰的校准，以及探针序列效应的标准化。

2.2总结模型

如前所述，SNP探针对的探针几乎相同，只是SNP位点的核苷酸不同。因此，如果探针第页和第页对应于相同的探针对_个人账户将接近_{英国皇家银行}.图1说明了属于同一探测对的探测如何具有相似的相似性。

2.2.1初始化Φ和C类

初始ASCN。我们的算法从初始化C类基于基因型的天真估计的矩阵。使用强度矩阵作为参数计算CN矩阵，并假设C类_艾岛+C类_Bi公司=2，即假设总CN为2。然后，对于每个SNP，如果大多数探针对的A等位基因探针的强度高于B等位基因的强度的2倍，则SNP将被分配C类_艾岛=2和C类_Bi公司=0（基因型AA），基因型BB反之亦然。否则，分配SNPC类_艾岛=C类_Bi公司=1（AB型）。的参数其他等位基因信号的2倍是根据现有的HapMap基因型数据的经验结果选择的。在比较预测基因型和验证的HapMap基因型时，此设置提供了最小的基因型错误调用百分比。

探针亲和力初始化。初始化后C类，我们推导了一个稳健的初始化Φ步骤如下：选择具有不同等位基因特异性CNs的两个不同样本，我们求解方程组Φ。使用不同的样本多次重复相同的过程。最后，计算的Φs已计算。我们发现，使用50个这样的随机对，在计算时间和准确性之间是一个很好的权衡，同时提供了对Φ使用中位数有助于抵御数据中存在的异常值（发出错误信号的探针）和初始假设不具有中性CN的样本。

处理总CN等于2的例外情况在上述初始化步骤中，假设总CN等于2。使用与探针亲和力初始化类似的过程，并改变不同探针的不同样本的作用，我们改进了CN矩阵的初始化。

这些稳健初始化的结果如下Φ⁽⁰⁾和C类⁽⁰⁾。此初始化步骤是随机的，因此每次运行可能会产生略有不同的结果。然而，在实现中，随机种子是固定的，因此结果在数值上是确定的。补充图S3结果表明，如果种子改变，解的精度不会改变。使用统一正随机数的标准初始化过程在这个特定的应用程序中不能很好地工作。补充图S4显示了本节中获得的亲和力和CN的初始值的质量。

2.2.2剔除异常值

有时，强度矩阵中会出现异常值，NMF算法通常对其敏感。为了识别和控制异常值，计算误差矩阵，作为初始强度矩阵和估计的初始强度矩阵之间的差值，等于Ŷ⁽⁰⁾=Φ⁽⁰⁾C类⁽⁰⁾然后，使用中值绝对偏差估计器稳健估计每个探针的样本间SD。任何大于跨样本SD的10倍（对于正态分布，大约为6个SD）的探针强度误差都被视为异常值。这些异常值被赋予相应的值Ŷ⁽⁰⁾.

2.2.3 NMF优化

我们认为当两个连续矩阵之间的绝对差最大时，估计值收敛|C类^(n个+1)−C类^(n个)|，小于ϵ（默认为0.01）。在极少数情况下（～1%），当算法不收敛时，我们在n个_最大值（默认为10）次迭代。上述NMF优化的伪代码为

1. n个← −1

2. 重复

3. n个←n个+ 1

4. Φ ^(n个+1)← 最大(YC公司^(n个)+, 0)

5. C类 ^(n个+1)← 最大值(Φ^(n个+1)+Y（Y）, 0)

6. 规格化（Normalize）Φ^(n个+1)使其列的总和相同

7. 重新计算C类^(n个+1)

8. 比例C类^(n个+1)使得列和的中值为2

9. 重新计算Φ^(n个+1)

10. 直到（最大值(|C类^(n个+1)−C类^(n个)|) < ϵ或n个≥n个_最大值)

算法每个步骤的复杂性在探测次数上是线性的(K（K）)和样品数量(我)。有关实际基准测试结果，请参阅痤疮包裹。

2.2.4 GWS阵列的应用

即使在技术探针复制的情况下，也可以认为探针的亲和力受周围探针的影响（兰登等。,2009)。因此，可以假设每个关联都是独立的，并使用算法进行计算。此外，具有完美复制的知识可用于验证优化结果。出于这些原因，我们选择独立处理完全复制的探针，并像前几代芯片类型一样处理GWS阵列。

方程式(8)与PLASQ（LaFramboise）模型极为相似等。,2007)。尽管如此，PLASQ假设属于某些组的探针（正向链、反向链、不同等位基因和不同数量的错配）的亲和力是相同的，而ACNE则不是。此外，ACNE和PLASQ之间的主要区别在于，在ACNE中，计算每个数据集的亲和力和CN，而在PLASQ中，需要控制样本来计算亲和力，然后这些样本用于计算其他数据集中的CN。

2.3稳健扩展

ACNE提供缩放步骤（伪代码中的步骤8），以便获得接近两个的总CN。然而，通过使用以下更复杂的方法完成算法，可以改善结果。为了减少计算时间，此步骤不包含在主迭代循环中。

假设大多数样本的总CN之和等于2，则继主算法之后的建议调整将比例因子与每个等位基因相匹配。通常，其他方法（如CRMA v2或dChip）使用参考（或所有样本）总CN的中值或修剪平均值来实现模拟缩放。

由于可能存在总CN不等于2的位点样本（因为在某些样本中SNP所在的区域中存在CNA），因此方程中所示的线性系统(13)必须稳健地求解，以承受异常值的存在。在这种情况下，离群值是CN不同于2的样本。我们使用迭代加权最小二乘法（IWLS）求解方程组。如果实验中的一组样本已知为法线，则可以将其用作参考和C类^T型应该只包括与这些样本对应的行。

如果对于任何特定的SNP，大多数样本都有中性CN，尽管所有样本都可能是肿瘤样本，那么这种归一化步骤都很有效。如果大百分比的样本出现像差，则鲁棒归一化方法无法识别正常样本。对于这些情况，必须在实验中添加额外的正常样本，否则CN估计值将在受影响的区域中产生偏差。我们已经在模拟中测试（未显示数据），击穿点位于～30%，即至少70%的样品必须正常工作在特定SNP（不是全基因组）。

2.4下游台阶

一旦执行了原始CN估计，就可以应用许多后处理方法。例如，可以执行片段长度补偿，以根据DNA片段的长度（Nannya等。,2005)。此外，总CN也可以根据GC含量效应进行校正（Nannya等。,2005)。最后，使用分割方法将样本划分为具有常见CNs的区域，例如循环二进制分割（CBS）（Olshen等。,2004)。所有这些过程都可以在汇总步骤后应用，但它们不是ACNE的一部分，也不用于第3节.

2.5实施

建议的ACNE方法可在痤疮中实现的包R（右）（R开发核心团队，2010)。除了提供一个低水平的估计之外痤疮包还包括高级aroma.afmetrix公司框架（Bengtsson等。,2008年a)这使得ACNE可以应用于非常大的Affymetrix SNP数据集。有关更多信息，请访问http://www.aroma-project.org/.

3结果

为了显示结果的准确性，我们比较了ACNE和CRMA v2（Bengtsson等。,2009)Affymetrix公司的CN5。，2008)和dChip（Li和Hung Wong，2001)。我们将评估分为两部分。首先评估总CN，然后评估ASCN。对于ASCN，我们将ACNE与CRMA v2和dChip的等位基因特异性版本进行比较。CN5未包含在该比较中，因为它需要配对样本来估计ASCN。

3.1总DNA CN结果

通过在CN发生变化的肿瘤样本中选择基因组的一部分，并比较位于变化点两侧的SNP的总CN，对总CN结果进行评估。我们选择了三个CN发生变化的不同地区。图2a显示了乳腺癌研究样本GSM182834中2号染色体上靠近150 Mb的一个这样的区域。此区域的副本从一个更改为两个。其他两个区域显示在补充材料，其中一个是基于Liu的GWS6数据等。(2009).

位于变化点之前的SNP被认为具有一定的总CN（删除区域），而在跳跃后的SNP具有更大的总CN.对于任何阈值，我们都有真阳性、TP、（正常区域中的SNP高于阈值）、假阳性、FP、（删除区域中高于阈值的SNP）、真阴性、TN、，（删除区域的SNP低于阈值）和假阴性（FN，正常区域的SNPs低于阈值）。FP速率FPR=FP/（FP+FN）和TP速率TPR=TP/（TP+TN）在构成ROC曲线的不同阈值下进行评估。本特森也使用了这种评估方法等。(2009)并在相应的补充材料.

对于此特定示例，所选区域位于140到160 Mb之间。我们还设置了一个148到152 Mb之间的“安全”区，因为很难知道跳跃的确切位置，所以分析中没有包括该区域。图3表明ACNE优于其他方法（CRMA v2、dChip和CN5）。CRMA v2和dChip具有类似的性能。在本研究中，CN5结果落后于其他方法。

3.2等位基因特异性CN结果

我们将我们的方法与通过使用dChip（'AS-dChip'）和CRMA v2（'AS-CRMA v2'）的等位基因特异性版本获得的数据进行了比较。

ASCN的验证比较困难，因为没有“基本事实”可供比较。我们采取了两种方法来验证我们的方法：分析B等位基因信号（β）的分数曲线图和使用HapMap的基因分型数据分析杂合和纯合调用。

3.2.1更明显的等位基因B组分

图2显示了难以分析的肿瘤样本的CN和β图。CN图(图2a） 显示三个不同的区域：从染色体开始到150 Mb（CN小于2），从位置150到180 Mb（中性CN），以及从180 Mb到染色体末端（CN再次低于2）。存在部分CN的一个可能原因是，该样本不是纯肿瘤，而是正常组织和肿瘤组织的混合。正常组织的存在导致β图中出现新的点云。在一次缺失中，肿瘤样本仅显示一个等位基因a（B），如果在正常组织中，同一SNP同时具有两个AB等位基因，那么等位基因B（a）将由于两个组织的混合而发出一些信号。在第二个区域，存在LOH，中性CN从150 Mb到180 Mb。在这种情况下，肿瘤组织将有AA或BB呼叫。在这种情况下，正常组织中带有杂合子叫声的SNP的β信号也将移向1/2，但这种影响比之前的情况弱，因为肿瘤样本的信号大两倍。

图2b、 c和d显示了使用ACNE、AS-dChip和AS-CRMA v2的三个区域中β的不同云。我们的结论是，在这个特定的样本中，染色体的第一和第三段有一个缺失，中间区域有中性CN的LOH。这个例子说明了ASCN的估计如何澄清对特定样本的理解。

ACNE比AS-dChip或AS-CRMA v2更好地分离不同的等位基因特定状态，β图显示位于预期位置的点云更紧密。

3.2.2改进杂合子呼叫的检测

除了研究β图外，我们还使用了HapMap提供的基因分型信息来评估ACNE。在正常样本中，每个等位基因的CN可以是（2,0）、（1,1）或（0,2）。我们分析了HapMapping250K_Nsp数据集的一些样本，这些样本具有可用SNP的基因型。该基因型在正常样本中产生ASCN。

ACNE、AS-dChip和AS-CRMA v2的估计ASCN如所示图4此图直观地显示ACNE比其他两种方法性能更好。结合ASCN和基因型信息，可以量化每种方法的性能。

使用ACNE，LH在[0，1]范围内。对于AS-CRMA v2，此值可能略超出此范围，因为负估计值为C类_一个和C类_B类是允许的。如果SNP是杂合的，其LH接近1，如果是纯合的，LH接近0。因此，LH可以用来辨别SNP是否是纯合的。设置不同的LH阈值，并使用HapMap中的基因型信息作为金标准，可以为每个阈值创建一个列联表，从而创建ROC曲线。该曲线告诉我们计算LH的不同算法的质量，以及ASCN的质量。图5显示了使用ACNE、AS-dChip和AS-CRMA v2方法的ROC曲线。除此曲线外，我们还包括了正常样本的等位基因B组分密度图(补充图S9)。我们截断了上面ACNE的曲线图，因为有大量的SNPβ非常接近于零和一。这些密度图表明，与AS-dChip和AS-CRMA v2相比，使用ACNE的点云更紧密，其质心更接近理论值。

4讨论

与以前的最先进的摘要方法相比，ACNE摘要方法最大的改进是它能够更准确地估计每个等位基因的CN值。我们还观察到，与CRMA v2、dChip和CN5相比，总CN有所改善。ACNE优于其他探针的原因之一是它能够估计SNP中每个探针的交叉杂交。

对于GWS6阵列，如中所述第2.2.4节，一个SNP探针对只有三个或四个重复。用其他算法代替ACNE可以很好地解决这个问题。如所示补充材料，尽管使用ACNE比CRMA v2的ASCN结果更好，但这些阵列的总CN改进不如以前的Affymetrix芯片类型中的明显(补充图S10和S11).

ACNE的复杂性与探针数量呈线性关系(K（K）)和样品(我)。在多达5000个样本的模拟中计算时间几乎完美地证实了线性关系。换句话说，ACNE可以很好地根据样本数进行缩放。

最近Staaf等。(2008)提出了新的分割算法，不仅使用总CN作为输入数据，还使用等位基因B（β）的分数作为输入数据。这些算法主要应用于Illumina数据，因为β图的噪声较小。准确估计β将有助于这些算法和其他算法使用Affymetrix阵列识别基因组的不同片段。

5结论

本文描述了一种从任何Affymetrix基因分型阵列估计ASCN的新算法。初始预处理步骤借用自CRMA v2。具体地说，该算法主要集中在摘要方法中。使用NMF，它可以估计每个等位基因的CN。ROC分析表明，ACNE不仅在ASCN方面，而且在总CN估计方面都优于其他最先进的方法。这些改进，尤其是在ASCN中，使得增加分割算法的信息以发现基因组中的不同畸变成为可能。

基金：纳瓦拉大学；研究基金会（Fundación para la Investigación Médica Aplicada）；NCI（授予U24 CA126551）。

利益冲突：未声明。

参考文献

作者注释