物理化学指纹图谱通过整合机器学习和结构分析揭示肠道病毒A71神经毒力的趋同进化决定因素

《Virus Research》：Physicochemical fingerprinting reveals convergent evolutionary determinants of enterovirus A71 neurovirulence through integrative machine learning and structural analysis

　　

Z-scale描述符（Z-scale descriptors）将肠道病毒A71（EV-A71）序列转换为连续的物理化学空间。五个位置显示出与神经毒力相关的趋同物理化学变化。逻辑回归（Logistic Regression）的受试者工作特征曲线下面积（AUC）为0.889，优于七种复杂机器学习算法（Machine Learning, ML）。SHapley Additive exPlanations（SHAP）揭示P2124疏水性是中枢神经系统（Central Nervous System, CNS）入侵的主要驱动因素。关键位置映射到三维（3D）蛋白质结构中的复制机器界面。

**研究背景与问题**
肠道病毒A71（Enterovirus A71, EV-A71）是无包膜的单股正链RNA病毒，属于小RNA病毒科（Picornaviridae），是手足口病（Hand, Foot and Mouth Disease, HFMD）的主要病原体。多数感染呈自限性皮肤黏膜病变，但部分毒株具有显著神经嗜性，可引发无菌性脑膜炎、脑干脑炎、急性弛缓性麻痹及致命性心肺功能衰竭。在亚太地区，EV-A71的反复大规模暴发给儿童健康带来沉重负担，使其成为全球卫生安全优先关注的病原体。传统研究主要依赖反向遗传学方法，针对个别突变（如VP1-Q145E、VP1-E164D、5′UTR-158C）在小鼠适应模型中进行验证。这类方法通量低，且孤立考察突变，无法捕捉全球异质性谱系中神经毒力进化的多维、组合特征。大规模基因组监测虽已产生前所未有的序列数据，但循环毒株的高维异质性使直接序列-表型映射变得复杂。机器学习（Machine Learning, ML）提供了一种从高维生物数据中提取模式的有力范式，但特征表示是关键瓶颈：传统独热编码或k-mer编码将氨基酸视为缺乏生化内容的类别标签，而Z-scale描述符（Z-scale descriptors）将每个残基编码为三维物理化学空间中的点，捕获疏水性（Z₁）、空间位阻（Z₂）和电子极性（Z₃），从而保留驱动功能后果的连续生物物理关系。本研究旨在回答以下问题：能否将全球循环EV-A71毒株的高维序列异质性缩减为一组紧凑的物理化学偏移信号，从而稳健区分神经毒力与非神经毒力分离株？

**研究内容与结论**
研究人员构建了一个整合的“物理化学指纹图谱”框架，包括：（i）将多聚蛋白序列嵌入Z-scale空间；（ii）通过严格的两阶段统计筛选识别与神经毒力相关的位置；（iii）在留一法交叉验证（Leave-One-Out Cross-Validation, LOOCV）下基准测试八个ML分类器；（iv）通过SHapley Additive exPlanations（SHAP）剖析模型归因；（v）将识别出的位置映射到三维蛋白质结构，以提供可实验验证的非因果生化假设的结构背景。主要结论：从267个全长基因组（7个CNS相关，260个非CNS，来自GenBank 1998-2024年全球分离株，涵盖中国、法国、美国等地）中，通过两阶段统计筛选（Mann-Whitney U检验和逻辑回归优势比）识别出五个核心位置（P2124、P997、P1743、P1711、P1246），其物理化学变化与CNS受累显著相关。L2惩罚逻辑回归（Logistic Regression）在LOOCV下获得最高AUC（0.889），敏感度0.857，特异度0.846。SHAP分析显示P2124疏水性（Z₁）是模型最主要贡献者，P1246极性（Z₃）与较低CNS几率相关。BrCr参考株成熟蛋白映射将P2124定位于3D RNA依赖性RNA聚合酶（RNA-dependent RNA polymerase, RdRp）的拇指结构域表面环（残基393），P1743位于3D RdRp的N端食指亚域（残基12），P1711位于3C蛋白酶（3C protease）的C端柔性环（约残基163），P1246位于2C ATP酶/解旋酶结构域（约残基135），P997位于2A蛋白（无实验结构）。系统发育分析显示，七个CNS相关毒株散布于至少四个遗传谱系，支持趋同进化假说。该研究发表在《Virus Research》。

**关键技术方法**（不超过250字）
研究人员从GenBank检索1998-2024年全球EV-A71相关核苷酸序列，经长度（≥7000 bp）和表型锚定（组织标本来源，映射至呼吸道、系统性、CNS层次）过滤，最终获得267个序列（7个CNS，260个非CNS，来源包括中国、法国、美国等）。使用MAFFT v7.490进行多序列比对，IQ-TREE v2.2.0构建最大似然系统发育树（1000次超快bootstrap）。将多聚蛋白序列通过Z-scale描述符（Z₁疏水性、Z₂体积、Z₃极性）嵌入连续物理化学空间。采用两阶段统计筛选：Mann-Whitney U检验保留20个候选特征（P<0.001），再通过逻辑回归优势比和生物学过滤（蛋白域功能背景、进化保守性、共线性排除）获得5个核心位置。基准测试8个ML算法（逻辑回归、随机森林、支持向量机、梯度提升机、极端随机树、k-近邻、高斯朴素贝叶斯、多层感知机），使用LOOCV和类权重调整，选取逻辑回归为最优模型。SHAP计算模型特征归因，并将五个位置映射至实验确定的EV-A71蛋白结构（PDB: 3N6L、8CNY、5GQ1等）。

**研究结果**
3.1 系统发育离散性与趋同神经毒力获得相容
最大似然系统发育树显示，七个CNS相关毒株未形成单系群，而是散布于至少四个遗传谱系，Bootstrap值>80%，支持趋同进化候选物理化学谱的假设，但需更多CNS全长基因组验证。

3.2 基线特征与时间平衡
CNS组与非CNS组在采集年份中位数上无显著差异（CNS: 2013年，IQR 2012-2016；非CNS: 2014年，IQR 2010-2019；P=0.838），但地理分布有显著差异（P=0.014），反映CNS分离株的地理非随机性。

3.3 两阶段筛选识别五个核心物理化学位置
从3L维特征空间中，第一阶段Mann-Whitney U检验保留20个位置（P<0.001），第二阶段逻辑回归排序后经共线性与域背景过滤，得到五个核心位置（表3）。P2124_Z₁效应量最大（OR=4.28，95% CI 1.47-12.51），表明该位置疏水性每增加一单位，CNS受累几率增加。P1246_Z₃则与较低CNS几率相关（OR=0.50，95% CI 0.30-0.83）。

3.4 物理化学景观热图显示协同偏移
27个代表性样本（7个CNS，20个随机非CNS）的热图显示，CNS分离株在P2124和P997处趋向更高疏水性，P1246处趋向更低极性，形成协同物理化学模式，支持候选趋同信号。

3.5 逻辑回归在内部LOOCV中AUC最高
八种算法LOOCV基准测试显示，L2惩罚逻辑回归获得最高AUC（0.889），敏感度0.857，特异度0.846，F1分数0.632，马修斯相关系数0.601。线性分类器在该极端不平衡数据集（2.6%阳性）中优于集成和深度模型，与Z-scale值与结局的近似单调关系一致。

3.6 SHAP归因总结模型关联，而非致病机制
SHAP分析显示P2124_Z₁为模型最主要贡献者，高疏水性将预测推向CNS受累；P1246_Z₃相反，低极性将预测推向CNS受累；P997_Z₂、P1743_Z₂和P1711_Z₁贡献中等至适度。所有五个特征呈近似单调方向性，支持线性模型选择及CNS受累与多位置协同物理化学偏移的统计关联。

3.7 探索性列线图用于假设生成风险分层
列线图将逻辑回归模型可视化为基于点的评分系统，P2124_Z₁和P997_Z₂占据最宽评分区间，与优势比和SHAP排序一致。该列线图仅用于假设生成，而非验证的临床预测工具。

3.8 结构定位将关键位置映射至成熟病毒蛋白
三维结构映射（图10）显示：P2124位于3D RdRp拇指结构域表面环（残基393），P1743位于3D RdRp N端食指区域（残基12），P1711位于3C蛋白酶C端柔性环（约残基163），P1246位于2C ATP酶/解旋酶结构域（约残基135），P997位于2A蛋白（无实验结构）。这些结构背景为后续实验提供生化假设，但未证明物理化学变化改变病毒复制或组织嗜性。

**讨论与结论**
讨论部分强调，本研究识别的五个位置物理化学指纹是假设生成关联，而非确定的神经毒力决定因素。核心限制是仅7个CNS阳性全长基因组，导致阳性类极度稀缺（7:260），使LOOCV AUC估计可能乐观且不稳定，统计假阳性风险增加，分类器泛化能力受限。CNS全长基因组的匮乏源于EV-A71神经系统并发症罕见、临床测序常针对VP1而非全长、CSF病毒载量极低（Leon等报道0.06-0.55 reads per million）以及公共数据库缺乏神经表型注释。因此，所有结果（模型、五位置签名、结构解释、列线图）应视为假设生成，迫切需要更大、临床注释更完整的独立队列进行前瞻性验证。此外，二元CNS/非CNS表型将神经系统疾病谱压缩为单一标签，可能掩盖脑膜炎、脑炎、急性弛缓性麻痹等不同结局的决定因素；基于比对的坐标系统与BrCr标准编号存在偏移，成熟蛋白映射（尤其P1246至2C、P997至2A）在切割位点边界不确定时应视为近似；Z-scale仅捕获三个主要物理化学轴，未涵盖骨架柔性、局部上位效应、RNA二级结构或宿主免疫选择等特征。总体而言，结果支持EV-A71 CNS受累可能与非结构蛋白区域的物理化学偏移统计相关，而非证明复制机器适应或因果关系。在仅7个CNS阳性基因组的限制下，五位置签名提供了病毒基因组监测与结构假设生成之间的可解释桥梁，其当前价值在于假设生成和风险优先排序，长期价值取决于更大平衡队列的前瞻性验证及直接实验测试。

**结论**
本研究提示，EV-A71 CNS受累不仅与离散的氨基酸身份相关，还可能与非结构多聚蛋白中五个位置残基的连续物理化学性质相关。基于Z-scale描述符（疏水性、体积、极性）构建的简约逻辑回归模型在内部LOOCV中达到0.889的AUC；然而，该估计受限于仅7个CNS阳性基因组，不应视为可泛化预测的确凿证据。SHAP归因将P2124_Z₁（3D RdRp残基393处的疏水性）识别为模型主要贡献者，P1246_Z₃（2C解旋酶约残基135处的极性）为与较低CNS几率关联最强的特征，但这些是统计和模型归因发现，而非确认的机制。基于BrCr的成熟蛋白映射将所选位置置于结构可解释的背景中，并生成可测试的生化假设；未建立改变复制细胞器形成、神经嗜性或CNS选择性复制的证据。七个CNS相关毒株在不同谱系中的系统发育离散性与趋同进化相容，但该解释需要更多CNS全长基因组方能确认。因此，列线图和风险评分框架现阶段仅用于探索性优先排序，而识别出的位置需要前瞻性队列验证和直接实验测试。