XGBoost-SHAP分析优先识别胰腺祖细胞中预测影响HNF4A结合基序的低频变异

《PLOS Genetics》：XGBoost-SHAP analysis prioritizes lower-frequency variants predicted to affect HNF4A binding motifs in pancreatic progenitors

　　

专业术语翻译 摘要 背景：非编码遗传变异是许多复杂性状的基础，然而预测其在特定细胞类型背景下的调控影响仍是一项重大挑战。转录因子（Transcription Factor, TF）结合是非编码变异影响基因调控的主要机制，但染色质可及性、变异频率和转录因子结合之间的关系复杂且非线性，限制了传统线性模型的性能。 方法：研究人员从七个诱导多能干细胞来源的胰腺祖细胞样本中生成了单核ATAC测序（single-nucleus Assay for Transposase-Accessible Chromatin using sequencing, snATAC-seq）数据，并鉴定了93,803个与可及染色质区域重叠的杂合单核苷酸多态性（Single Nucleotide Polymorphism, SNP）。通过二项检验对每个变异的等位基因特异性效应（Allele-Specific Effect, ASE）进行量化。研究人员训练了XGBoost、逻辑回归和正则化逻辑回归模型，利用染色质读段深度、峰值评分、变异次等位基因频率（Minor Allele Frequency, MAF）以及五种胰腺转录因子（HNF4A、HNF1A、NEUROD1、NKX6-1和PDX1）的结合预测来预测等位基因特异性效应。模型性能通过五折交叉验证的受试者工作特征曲线下面积（Area Under the Curve, AUC）以及在测试集和验证集SNP上进行评估。使用Shapley（SHAP）值分析来解释模型预测，并根据其调控特征贡献对变异进行聚类。 结果：XGBoost显著优于线性模型，尤其是在严格的等位基因特异性效应阈值下，这与它捕捉特征间非线性相互作用的能力一致。染色质读段深度、峰值评分和次等位基因频率是最强的预测因子。按读段深度和MAF十分位数分层显示，XGBoost在各十分位数中的整体预测性能均优于线性模型，表明其优势并非仅由这些变量解释。基于Shapley值谱对变异进行降维，揭示了七个具有不同染色质状态和变异等位基因频率的聚类。这些发现得到了七个聚类间等位基因频率分布高度显著的成对差异的支持。其中一个聚类富集了较低MAF的SNP，这些SNP显示出与HNF4A结合基序的强烈预测相互作用。由于MAF及其SHAP贡献是嵌入的输入，因此该聚类的低MAF特性是预期的；其高ASE富集度和独特的HNF4A相互作用结构是具有生物学信息价值的发现。联合注释依赖性耗竭（Combined Annotation Dependent Depletion, CADD）分析进一步支持在该聚类中HNF4A相关预测与预测有害性之间存在更强的关系。 结论：将单核染色质可及性数据与可解释的XGBoost-SHAP框架相结合，提供了一种在特定发育背景下优先考虑非编码变异并识别非线性调控关系的方法。在该框架内，低MAF聚类中HNF4A相关预测与CADD预测有害性之间的更强关系，与HNF4A结合的调控元件受到选择性约束一致，并提出一个假说，即低频非编码变异如何可能导致HNF4A相关青少年的成人发病型糖尿病（Maturity-Onset Diabetes of the Young, MODY1）家族内的表型变异。

论文解读文章

研究背景与目的

全基因组关联研究（GWAS）已鉴定出数千个与常见疾病相关的位点，但大多数关联变异位于非编码区，被认为通过细胞类型特异性的基因调控发挥作用。预测这些非编码变异的调控功能是当前基因组学的一项重大挑战。转录因子（TF）结合是变异影响基因调控的主要机制之一，但染色质可及性、变异频率和TF结合之间存在复杂的非线性关系，传统的线性模型难以有效捕捉。本研究旨在利用单核ATAC测序（snATAC-seq）产生的染色质可及性数据，结合可解释的机器学习框架，预测非编码变异的调控功能，特别是在胰腺祖细胞这一特定的发育背景下。研究人员重点关注五种关键的胰腺转录因子——HNF4A、HNF1A、NEUROD1、NKX6-1和PDX1，并探究它们与变异功能的关系。

研究方法

研究人员使用了来自7个诱导多能干细胞（iPSC）来源的胰腺祖细胞（iPSC-PPC）样本的snATAC-seq数据。主要关键技术方法包括：1）通过snATAC-seq生成染色质可及性图谱，并鉴定与可及峰重叠的杂合SNP；2）利用二项检验量化每个SNP的等位基因特异性效应（ASE）作为模型训练的表型；3）采用XGBoost梯度提升树模型，并与广义线性模型（GLM）和弹性网络正则化逻辑回归模型（GLMNET）进行比较，预测ASE状态；4）运用Shapley（SHAP）值分析解释XGBoost模型的预测，并对变异进行聚类；5）使用FIMO工具计算五种TF的位置权重矩阵（PWM）得分，评估变异对TF结合序列的影响。样本来源于iPSCORE队列，所有分析的变异在iPSCORE群体中MAF≥5%。

研究结果

XGBoost模型在预测ASE方面优于线性模型：通过比较XGBoost、GLM和GLMNET三种模型在不同ASE p值阈值下的表现，研究发现XGBoost在所有阈值下均表现出更优的性能，尤其在最严格的ASE阈值（p=1e-5）下优势最为显著。XGBoost的优势在于其能够捕捉读段深度、峰值评分、MAF和TF PWM得分等特征间的非线性交互作用。按读段深度和MAF十分位数分层分析证实，XGBoost的整体性能优势并非由这些单一变量驱动。

snATAC-seq读段深度、峰值评分和SNP MAF对预测结果至关重要：SHAP分析表明，snATAC-seq读段深度是模型中最重要的特征，其次是SNP MAF、峰值评分和各TF的PWM得分。这表明染色质可及性和变异频率是比单纯的TF结合预测更强的ASE预测因子。

UMAP聚类揭示SNP MAF与HNF4A PWM得分的非线性相互作用：基于SHAP值的UMAP聚类将SNP分为7个簇。其中，簇1、2和6主要包含低读段深度、低ASE比例且富集于抑制/静息染色质的SNP。簇7则富集于高读段深度的活性启动子区域，具有最高的ASE比例。簇3和簇4在读段深度和染色质状态上相似，但簇3的ASE比例是簇4的7.48倍。簇3的独特之处在于：其SNP具有相对较低的MAF，并且展现出涉及HNF4A基序得分的强SHAP交互模式，特别是MAF与HNF4A参考等位基因PWM得分之间的交互。进一步的逻辑回归分析确认，MAF与簇3成员身份的强关联独立于读段深度和峰值评分。此外，在簇3中，HNF4A结合破坏的SHAP得分与CADD预测有害性之间的正相关性显著强于其他簇，表明该簇中影响HNF4A结合的变异可能具有更重要的功能后果。

讨论与结论

本研究证明了将小样本量的snATAC-seq数据与可解释的XGBoost-SHAP框架相结合，在特定细胞类型中优先排序功能性非编码变异的可行性和生物学价值。XGBoost因其捕捉非线性关系的能力而优于线性模型。SHAP分析不仅揭示了各特征的贡献，还通过聚类发现了具有不同调控环境的变异群。特别地，一个由较低MAF变异组成的簇（簇3）表现出与HNF4A结合基序的独特预测交互模式和高ASE富集度，这为理解HNF4A相关MODY1的家族内表型变异提供了潜在机制假说。该框架为实验验证提供了可测试的候选变异，并为将预测方法扩展到整个非编码基因组奠定了基础。研究结论明确指出，整合单核染色质可及性数据与可解释的XGBoost-SHAP框架，为在特定发育背景下优先处理非编码变异并识别非线性调控关系提供了一种有效方法。