基于效应方向和效应大小的稀有变异分类的聚类方法

《Algorithms》：A Clustering Approach for Rare Variant Classification by Effect Direction and Magnitude

　　

几项基于基因的检验（例如序列核关联检验（sequence kernel association test, SKAT））已被开发用于评估稀有单核苷酸变异（rare single nucleotide variants, SNVs）与疾病性状之间的关联。然而，这些汇总方法无法在关联区域内将潜在因果变异与无效变异区分开来。为解决这一局限性，研究人员提出了gvClust，这是一种聚类方法，通过将高斯混合模型（Gaussian mixture model, GMM）应用于多变异模型中的变异水平汇总统计量，将稀有变异分类为无效组和信号组。信号变异进一步根据其效应方向和大小被划分为风险亚组和保护亚组。研究人员在模拟研究中采用调整兰德指数（adjusted Rand index, ARI）、均方误差（mean squared error, MSE）和不同样本量、效应配置、结局类型及连锁不平衡（linkage disequilibrium, LD）结构下的簇数选择准确度对gvClust进行了评估。在模拟中，gvClust随样本量增加而性能提升，在大样本量下对连续性状的簇数确定实现了高准确度，并优于k-means聚类和仅初始化聚类。研究人员随后将gvClust应用于一项大型全基因组关联研究和meta分析中与血压性状相关的六个基因中的稀有变异。在实际数据应用中，gvClust识别出了明显的无效簇、风险簇和保护簇。这些结果表明，gvClust为分类关联区域内的稀有变异提供了一个实用框架，并可能有助于改善稀有变异信号的生物学解释。

## 论文解读：基于效应方向和效应大小的稀有变异分类的聚类方法

### 一、研究背景与目的
全基因组关联研究（genome-wide association studies, GWAS）已鉴定出大量与复杂疾病和数量性状相关的常见遗传变异，但低频率变异和稀有变异（minor allele frequency < 1%）作为未解释遗传力的重要来源，其作用机制尚不明确。现有基于基因的关联检验方法（如负担检验、序列核关联检验（SKAT）、SKAT-O、聚合柯西关联检验（ACAT-O））能够检测基因或区域内稀有变异的集合效应，但无法在显著关联区域内区分潜在因果变异与无效变异，也无法根据效应方向（风险或保护）和效应大小对稀有变异进行细分。为填补这一空白，研究人员提出了一种新的聚类方法（gvClust），利用高斯混合模型（Gaussian mixture model, GMM）对多变异模型获得的变异水平汇总统计量进行建模，将稀有变异分类为无效（null）、风险（risk）和保护（protective）组。该方法有助于深入理解稀有变异与疾病性状之间的生物学机制，为药物靶点发现和基因治疗设计提供支持。该论文发表在《Algorithms》上。

### 二、主要关键技术方法
研究人员采用的主要技术方法包括：（1）基因水平关联检验：使用负担检验、SKAT和ACAT-O识别与性状显著关联的基因区域（信号区域）；（2）多变异模型：对信号区域内所有稀有变异进行多变量回归，获得每个变异的效应估计（β系数）及其标准误（standard error, SE），并调整协变量；（3）高斯混合模型（GMM）：假设每个β系数来自一个无效簇（均值为0）和K个信号簇（各有不同均值），通过期望最大化（expectation–maximization, EM）算法迭代估计参数（簇均值、混合比例），初始值通过可变阈值法和K-means聚类确定，利用贝叶斯信息准则（Bayesian information criterion, BIC）在2至7个簇中选择最优簇数；（4）模拟研究：基于PhenotypeSimulator和sim1000G R包生成独立或具有连锁不平衡（linkage disequilibrium, LD）的稀有变异（取自1000 Genomes Project参考数据集，BTG3基因区域），设定六种不同效应方向和效应大小的场景，评估不同样本量（5000、15000、25000）和结局类型（连续/二分类）下的聚类性能，指标包括调整兰德指数（ARI）、均方误差（MSE）和簇数准确度，并与k-means和仅初始化聚类进行对比；（5）实际数据应用：数据来源于Surendran等人（2019）的大型GWAS和meta分析（包括超过800,000名个体，来自CHARGE、CHD Exome+、GoT2D:T2DGenes、ExomeBP四个联盟和英国生物银行UK Biobank），对与收缩压（SBP）、舒张压（DBP）、脉压（PP）和高血压（HTN）相关的六个基因（NPR1、DBH、COL21A1、NOX4、PLCB3、CEP120）中的稀有变异进行聚类，并利用综合注释依赖性耗竭（Combined Annotation Dependent Depletion, CADD）评分和后果类型进行注释。

### 三、研究结果
#### 3.1. Simulation Studies（模拟研究）
通过1000次重复模拟，研究人员评估了gvClust在不同条件下的性能。

##### 3.1.1. Simulations Without LD Structure（无LD结构模拟）
- 随着样本量增大，ARI值显著提高，例如在场景2中（1/3无效、1/3正向效应、1/3负向效应），连续性状下样本量为5000时平均ARI为0.60，增至25000时达0.96。
- 三种Z分数加权方案（等权重、逆标准误权重、组合权重）表现相当。
- 场景2和6（效应方向相反）的ARI值高于其他场景。
- 当样本量达25000时，连续性状在所有场景下的簇数确定准确度均高于0.9。
- MSE随样本量增大而减小。
- 方差分析（ANOVA）在所有重复中均显著（p<0.05），表明不同簇的均值存在显著差异；Mann-Whitney U检验显示无效簇的效应大小与0无显著差异（p>0.05）。
- 二分类性状的ARI值低于连续性状（例如场景1，无LD，样本量15000时连续性状ARI=0.66，二分类ARI=0.33）。

##### 3.1.2. Comparison of Simulations with and Without LD（有/无LD模拟比较）
- 使用sim1000G包模拟的60个稀有变异显示中度或低LD：约29.9%的变异对相关系数在0.01-0.05之间，最大相关系数为0.26。
- 在场景2、4和6（效应方向相反且均值差异较大）中，有LD与无LD的ARI、MSE和簇数准确度相当（例如场景2，样本量15000，连续性状，组合权重：ARI 0.91 vs 0.9；MSE 0.0023 vs 0.0022）。
- 在场景1、3和5（信号簇均值差异较小）中，有LD时ARI和准确度较低，MSE较高（例如场景1，样本量15000，连续性状：无LD时ARI=0.66，有LD时ARI=0.51；MSE分别为0.0033和0.0043；簇数准确度分别为0.911和0.715）。这表明LD对性能的影响在簇均值差异较小时更显著。
- 对比gvClust、仅初始化聚类和k-means：gvClust在所有场景中均表现最优，中位聚类准确度最高且变异最小；k-means在多数场景中性能较差、变异大；仅初始化聚类性能次于完整gvClust，证实了EM迭代步骤的重要性。

#### 3.2. Application to a GWAS on BP Traits with Rare Variants（应用于血压性状稀有变异GWAS）

##### 3.2.1. Identification of Blood Pressure Trait-Associated Genes（血压性状关联基因鉴定）
- 使用SKAT检验，在MAF<0.01条件下，NPR1、DBH、COL21A1和NOX4四个基因与至少一个血压性状显著关联（p<2.5×10^-6）；PLCB3和CEP120在MAF<0.05下关联。
- 各基因包含的稀有变异数：NPR1（13个）、DBH（29个）、COL21A1（26个）、NOX4（9个）。
- SBP关联NPR1、DBH、PLCB3；DBP关联DBH、PLCB3；PP关联COL21A1、NOX4、CEP120；HTN关联DBH、NPR1、PLCB3（基于单变异检验p<1×10^-4）。

##### 3.2.2. Annotation of Rare Variants Within Blood Pressure Trait-Associated Genes（血压关联基因内稀有变异注释）
大多数变异为错义变异，CADD评分≥20。

##### 3.2.3. Clustering of Rare Variants（稀有变异聚类）
- 第一种策略（每个基因-性状单独聚类）：对11个基因-性状组合进行聚类，每个组合识别出2-3个簇。约70%的稀有变异归入无效簇。例如，NPR1基因与HTN关联：12个变异中8个为无效簇，3个变异（rs140425746、rs61757359、rs61758562）形成负效应簇（平均标准化β=-2.05），1个变异（rs116245325）形成单独正效应簇（效应大小4.33），类似模拟场景2。
- 第二种策略（同一性状的合并基因聚类）：对与PP相关的三个基因（COL21A1、NOX4、CEP120）共55个稀有变异进行聚类，识别出5个簇（ANOVA p=1.39×10^-19）。无效簇包含43个变异，其余变异分为两个负效应簇（平均效应-0.097和-0.0438）和两个正效应簇（最强效应0.334和中等效应0.173），类似模拟场景4。

### 四、讨论与结论
研究人员提出的gvClust方法基于变异-性状关联的汇总统计量对信号基因区域内的稀有变异进行聚类。模拟研究表明，样本量增大能提升聚类性能，LD对性能的影响在簇均值差异较小时更为明显，且gvClust在连续性状上优于二分类性状。与k-means和仅初始化聚类相比，gvClust因采用基于模型的EM框架而表现更优。实际数据应用显示，大多数与血压性状相关的稀有变异被归为无效簇，这提示自然选择可能是塑造人类基因组稀有变异的主要力量；同时，更稀有的变异倾向于归入负效应簇。

本方法的一个关键优势是能够将稀有变异分类为信号变异和无效变异，区分相反效应方向或不同效应大小的变异，从而更准确地识别真实遗传信号。此外，该方法可整合多种注释信息对信号变异进行优先排序，且计算效率高，适用于大规模测序数据分析。研究人员指出，尽管多变异模型已考虑了LD，但LD仍会影响聚类性能，当样本量足够大（如N=25000）时算法表现良好。未来可扩展至基于大样本GWAS汇总数据的应用，以增强鲁棒性并减少假阳性。

**研究结论**：简而言之，所提出的聚类算法能够根据SNP-性状关联的汇总统计量，识别不同效应大小的风险性和/或保护性稀有变异。该方法可便捷地应用于新兴的大规模稀有变异GWAS汇总统计量，将性状关联的稀有变异识别并归类为无效组和具有离散效应大小的信号组。因此，该方法可能有助于识别基因组区域中潜在因果的稀有变异簇，最终帮助理解人类复杂性状的遗传结构，为药物靶点发现和基因治疗设计提供支持。