通过机器学习技术识别多转录组预测生物标志物，以探索肺癌的自然疗法

《Digital Discovery》：Identification of multi-transcriptomic prognostic biomarkers to explore natural therapeutics for lung cancer integrating machine learning

　　肺癌仍然是全球癌症相关死亡的主要原因，这凸显了开发新型治疗策略的迫切需求。细胞周期依赖性激酶1（CDK1）作为细胞周期的核心调节因子，已被发现是肺癌腺癌的致癌驱动因素和潜在靶点。本研究旨在整合转录组学、机器学习（ML）和先进的计算机模拟方法，以识别针对CDK1的天然产物衍生物抑制剂。为了识别差异表达基因，我们首先分析了四个不同的数据集（GSE19804、GSE10072、GSE18842和GSE10799）。蛋白质-蛋白质相互作用网络和拓扑分析表明CDK1是一个主要的枢纽基因（pKHG），在细胞周期和p53通路中富集。目标验证确认了CDK1的过表达、预后意义、免疫浸润联系及其突变关联。此外，通过基于ML的生物活性（pIC50）预测，从9667种天然植物化合物库中筛选出了潜在的先导分子。然后，选定的先导分子通过分子对接、分子动力学模拟、ADMET分析和结合自由能计算（MM-GBSA）进行了进一步评估。在选定的植物化学物质中，CID_14218027（?7.69 kcal mol?1）、CID_487089（?6.80 kcal mol?1）和CID_174880（?6.70 kcal mol?1）显示出最高的结合亲和力得分（GLIDE_XP得分）和稳定的分子相互作用。此外，MD模拟证实了配体-蛋白质复合物的构象稳定性，支持它们作为CDK1抑制剂的潜力。这一综合性的从组学到计算机模拟的流程将CDK1确定为一个稳健的治疗靶点，并突出了具有有利药理和物理化学特性的天然产物衍生物抑制剂。因此，这些发现为加速精准药物发现提供了一个可行的框架，后续正在进行实验验证。

1. 引言

肺癌是全球最致命的癌症，其每年导致的死亡人数超过了乳腺癌、结肠癌和前列腺癌的总和。2022年，肺癌导致了约182万例死亡，占所有癌症死亡人数的18.7%，远超过结直肠癌（9.3%）和乳腺癌（6.9%）。同年，新增了约248万例肺癌诊断，年龄标准化死亡率约为每10万人16.8例，这一趋势在人类发展指数较高的国家更为明显。如果目前的趋势保持不变，随着人口老龄化和增长，到2050年，新病例数可能增至约462万例，死亡人数可能增至约355万例，表现为总病例数和死亡人数的增加。肺癌之所以如此致命，是因为诊断往往较晚、分子异质性较高以及现有治疗方法的敏感性较低。尽管已经开发出许多FDA批准的合成药物，包括细胞毒性化疗药物、分子靶向药物和免疫检查点抑制剂，但肺癌患者的整体预后仍然具有挑战性。然而，药物间的相互作用、脱靶毒性、效力不足和药代动力学问题限制了这些发现转化为有效疗法。尽管标准合成药物在癌症治疗中发挥了重要作用，但这些药物通常受到副作用、毒性、耐药性和生产成本高的限制。特别是传统化疗虽然有效，但经常会引起临床上显著的副作用（如骨髓抑制和胃肠道毒性），从而降低生活质量并限制剂量强度。最初对治疗有反应的患者最终会发展出获得性耐药性（例如，在EGFR-TKI(osimertinib治疗后的耐药性），这推动了治疗需求的升级。因此，探索如细胞周期控制等替代弱点变得至关重要：细胞周期依赖性激酶是一种核心的有丝分裂激酶/酶，在肺癌中报告呈过度表达并且与预后相关。异常的致癌信号通路、DNA修复机制、细胞死亡和不受控制的细胞周期进展共同促进了肿瘤的起始和发展。这些分子异常表明，细胞周期调节因子、有丝分裂蛋白和其他重要的调节因子是治疗干预的有希望的靶点。转录组学和多组学研究有助于识别具有临床相关性的关键调节基因和信号通路。基于网络药理学的方法，包括蛋白质-蛋白质相互作用（PPI）网络分析、枢纽基因识别、转录组学因素和通路富集分析，是一些有效的发现关键分子调节因子的方法。

另一方面，天然化合物——更具体地说，来自药用植物的植物化学物质——已成为合成药物分子的一个有吸引力的替代品，因为它们的化学多样性和进化出的生物活性可以提供与现有合成库不同的新骨架和机制。先前的研究表明，超过60%的抗癌药物是天然的，包括著名的药物如紫杉醇、喜树碱或长春新碱。IMPATT数据库是一个大型数据库，包含4010种药用植物和17,967种植物化学物质的信息及其性质，这些信息在其他数据库中往往不完整。因此，这个数据库是肺癌药物的有用来源。如今，为了减少时间和成本，科学家们在通过实验（体内和体外）验证之前，考虑使用计算机模拟方法进行筛选。本研究开发了一个全面的计算工作流程，整合了转录组学、基于机器学习（ML）的先导筛选、分子对接和动态模拟研究，以筛选大量化合物库中的化合物，开发出毒性较低且疗效增强的治疗药物分子。先前的研究主要使用转录组学和机器学习（ML）方法在LUAD中识别潜在的生物标志物和治疗靶点或预测对已批准和实验性药物的响应。类似地，也有研究使用LUAD mRNA和突变谱开发了ML模型，预测对现有靶向和化疗药物的敏感性，在数十种药物上取得了良好的预测性能。与这些主要关注生物标志物/药物靶点识别、验证和现有药物响应预测的研究不同，本研究结合了基于LUAD转录组学的药物靶点识别和传统的基于ML的QSAR建模，预测植物衍生化合物的pIC50值，随后进行对接和分子动力学（MD）模拟，以优先选择有效的天然候选物。因此，本研究基于各种公开可用的GEO数据集中的转录组信息构建了PPI网络，使用传统的ML回归模型预测生物活性（pIC50），通过计算对接后的MM-GBSA评估结合强度，进行大规模MD模拟和药代动力学分析，创建了一种全面的基于天然产物的肺癌治疗药物发现方法。与传统的以对接为中心的研究相比，这种基于生物学的信息驱动的工作流程允许对疾病相关的靶点进行优先排序，并从大型天然化合物库中更系统地选择先导化合物。因此，本研究为优先选择针对CDK1的植物化学物质提供了一个全面的、以转化为导向的计算策略，并支持未来对所识别候选物的实验验证。本研究的完整指南见图1。

2. 材料与方法

2.1. 靶点识别

在本节中，我们对肺癌基因表达谱进行了系统分析，以揭示差异表达基因（DEGs），并突出主要的枢纽基因（pKHG），这对于理解疾病的潜在机制至关重要。

2.1.1. 微阵列表达数据集获取

我们使用了来自美国国家生物技术信息中心（NCBI）基因表达组学（GEO）数据库的四个公开可用的微阵列数据集（https://www.ncbi.nlm.nih.gov/geo/）。选定的数据集包括GSE19804、GSE10072、GSE18842和GSE10799，其中包含了肺癌和相应正常组织的基因表达谱。这些数据集的详细信息，包括平台和样本分布，总结在表1中。

表1 本研究中使用的肺癌GEO数据集摘要

GEO数据集
样本数量
癌症
对照组
平台
参考编号 GSE19804
120
60
60
GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
41,42

GSE10072
107
58
49
GPL96[HG-U133A] Affymetrix Human Genome U133A Array
43

GSE18842
91
46
45
GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
44

GSE10799
19
16
3
GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
45

2.1.2. 基于微阵列基因表达数据集的DEG分析

为了识别肺癌中的差异表达基因（DEGs），我们使用了NCBI提供的在线工具GEO2R（https://www.ncbi.nlm.nih.gov/geo/geo2r/）分析了四个选定的数据集（GSE19804、GSE10072、GSE18842和GSE10799）。每个数据集中的样本被分为“对照组”和“癌症组”，数据通过log2转换进行了标准化。使用limma包检测DEGs，同时使用Benjamini–Hochberg假发现率（FDR）方法进行调整以控制多重检验。调整后的p值<0.05且logFC>+1.0的基因被认为是显著上调的，而logFC<-1.0的基因被认为是显著下调的。

2.1.3. 共同DEGs（cDEGs）的识别

最后，我们通过交叉所有四个数据集中的DEG列表来识别共同的DEGs（cDEGs），这些基因是与肺癌相关的潜在候选基因。cDEGs使用R语言中的“dplyr”包提取，并通过R包“ggplot2”和“ggvenn”构建的Venn图可视化它们的重叠部分。

2.1.4. 通过PPI网络分析识别关键枢纽基因

确定了共同的DEGs（cDEGs）后，我们研究了相关蛋白质之间的相互作用，以揭示肺癌中的主要分子靶点。蛋白质-蛋白质相互作用（PPI）网络的筛选是通过STRING数据库进行的，只考虑了经过实验验证的相互作用。为了集中分析cDEGs之间的直接相互作用，我们将最小相互作用得分定义为低置信度（0.150），并且没有在初始壳层中添加任何额外的相互作用者，以减少分析结果中的间接相互作用。然后在Cytoscape（v3.10.4）中进行了可视化，每个节点代表一个蛋白质，每条边代表两个蛋白质之间的相互作用。为了揭示最显著的蛋白质，我们应用了CytoHubba应用程序，并在八种不同的拓扑方法中分析了枢纽基因：度数、最大邻域组件（MNC）、最大团中心性（MCC）、最大邻域组件密度（DMNC）、边缘渗透组件（EPC）、瓶颈、EcCentricity和接近度。通过这些拓扑分析，选择了相互作用百分比最高的蛋白质作为最重要的枢纽基因（KHGs），这些基因在识别肺癌相关网络中的关键中心靶点方面具有较高的价值。

2.1.5. KHGs的转录和转录后调控分析

为了从上游水平理解所识别KHGs的调控过程，我们进行了联合分析，以识别转录和转录后调控因子。转录因子的预测是在JASPAR数据库中使用TF-靶标相互作用进行的，miRNA-靶标相互作用则在miRTarBase中进行。调控网络是借助网络分析平台NetworkAnalyst（https://www.networkanalyst.ca/）构建的，该平台结合了这些相互作用并能够进行拓扑分析，以发现必要的调控因子。然后使用Cytoscape（v3.10.3）进一步可视化和分析这些网络，其中节点是调控因子或靶标基因，边是调节基因的相互作用。根据网络的拓扑特征，还对重要的调控因子进行了优先排序，从而揭示了可能调控肺癌相关KHGs的转录因子和miRNAs。

2.1.6. KHGs的GO和通路富集分析

我们首先使用Gene Ontology（GO）和KEGG通路富集分析，通过DAVID数据库来理解所识别DEGs的生物学功能。GO根据分子功能、生物过程和细胞组分对基因进行分类，而KEGG则识别相关的信号通路。对于GO富集分析，最小基因计数为2且p值小于0.05的术语被认为是显著的，KEGG通路则使用DAVID中的默认EASE得分进行选择。为了保证结果的可靠性，DAVID富集的GO术语和KEGG通路还借助Enrichr和GeneCloudOmics进行了进一步验证。在DAVID中找到的术语以及在DAVID和其他验证平台中也经常出现的术语被保留下来，使我们能够获得一套全面可靠的功能类别，据此将肺癌KHGs联系起来。

2.2. 目标验证

为了验证在肺癌中识别的主要枢纽基因（pKHG），我们专注于可用数据库中的肺癌腺癌（LUAD）数据集。pKHG的转录组学和蛋白质组学表达分析

pKHG是从候选枢纽基因中选出的，这是基于拓扑评分和富集分析的结合结果，确保了生物学上最相关和在网络中居中心位置的目标被优先用于下游验证。为了增强我们发现的可信度，我们接下来使用多个公开可用的资源在LUAD中验证了pKHG的表达模式。首先使用Tumor Immune Estimation Resource (TIMER 2.0)（https://timer.comp-genomics.org/）59根据TCGA数据检查了肺腺癌细胞组织和正常对照组之间pKHG的表达情况。为了进一步验证，我们还采用了GEPIA2平台（https://gepia2.cancer-pku.cn）60，整合了TCGA和GTEx数据集。箱线图参数设置如下：对数2倍变化系数（log2 FC）截止值=1，p值截止值=0.01，抖动大小=0.4，并且值经过对数2（TPM + 1）转换以便于可视化。GEPIA2的“Stage Plot”模块也被用来确定基因表达与LUAD病理分期（I–IV）之间的相关性。此外，我们还应用了UALCAN数据库（https://ualcan.path.uab.edu）61：（i）该研究使用TCGA RNA-seq数据验证了正常组织和LUAD组织之间mRNA表达的差异；（ii）该研究使用Clinical Proteomic Tumor Analysis Consortium (CPTAC)数据测试了pKHG的蛋白质表达水平。这种综合方法在转录组和蛋白质组学水平上都是一致的。

2.2.2. pKHG的生存分析

为了研究pKHG在LUAD中的预后作用，在GEPIA2阶段进行了Kaplan–Meier生存分析。中位截止值50%被用来将患者分为高表达组和低表达组。测量了总体生存期（OS）和无病生存期（DFS），并使用对数秩检验（log-rank test）来测试两者之间的差异。生存图在X轴上显示时间（以月为单位），在Y轴上显示生存概率（%），虚线表示95%的置信区间。这一讨论提供了关于pKHG在LUAD中的预后价值的宝贵信息。

2.2.3. 与pKHG相关的免疫浸润

我们使用TIMER2.0中的“Immune-Gene”模块单独探索了pKHG与LUAD上免疫细胞浸润之间的相关性。我们重点关注CD8+ T细胞、巨噬细胞和CD4+ T细胞，以确定最强的正相关和负相关。应用了解卷积算法来估算免疫细胞浸润，并使用经过肿瘤纯度调整的Spearman相关性来提供相关系数和显著性值。所得到的关联结果以热图的形式呈现，显示了这种肺癌枢纽基因与哪些类型的免疫细胞有最紧密的关联。

2.2.4. pKHG的突变和变异研究

为了继续研究与癌症相关的基因组变化问题，我们使用了cBioPortal平台（https://www.cbioportal.org/）62。利用这个资源系统地分析了已识别的主要关键枢纽基因（pKHG）的基因变异模式。在TCGA PanCancer Atlas研究中的数据分析（包括32种肿瘤类型和10,957个患者样本）是通过“Query by Gene”模块完成的。得到的“Cancer Type Summary”总结了pKHG在各种癌症中的突变和拷贝数变化类型。此外，还应用了“Mutations”模块来生成基因中确切突变位置的示意图。

2.3. 基于pKHG的计算机辅助药物发现

为了探索治疗机会，我们对经过验证的pKHG进行了基于结构的药物发现。这包括检索和准备目标蛋白的3D晶体结构、配体收集、虚拟筛选配体、分子对接（ADMET）、分子动力学（MD）模拟、MMGBSA、PCA（主成分分析）和3D-FEL（自由能景观）。

2.3.1. 目标蛋白的检索与准备

pKHG的晶体结构数据来源于RCSB蛋白质数据库63。为了便于后续的蛋白分析，使用BIOVIA Discovery Studio 202164移除了所有的异质原子、配体和水分子。然后使用Swiss-PDB Viewer 4.1.065对蛋白结构进行了能量最小化处理。

2.3.2. 化合物库的构建

为了检测可能的抑制剂，我们借助IMPPAT 3.0数据库31构建了一个通用的植物化学库。我们选择了33种具有传统药理意义的药用植物，并提取了这些植物根、叶和种子等部位的植物化合物。从这些植物中获得了9667种植物化合物，将它们收集起来形成一个全面的化合物库，用于未来的虚拟筛选和分子对接研究。

2.3.3. 基于物理化学性质的配体筛选

虚拟筛选（VS）是当代药物发现的重要组成部分，它提供了一种计算方法，可以在庞大的化学库中找到可能的生物活性化合物66。借助基于网络的ADMETlab 2.0平台67，药代动力学和毒性相关参数可以帮助这一筛选过程。为了缩小筛选范围，Lipinski提出的五规则（RO5）68常被用作判断化合物药物拟合性的标准。根据这一规则，化合物具有更高的口服生物利用度可能性：即它不超过五个氢键供体，不超过十个氢键受体，分子量小于500 Da，且logP值小于五。满足这些要求的植物化合物被认为是进一步计算和实验研究的良好候选物。重复的化合物被排除在外，剩下来的化合物具有可用的3D结构，将用于后续研究。

2.3.4. 基于机器学习的化合物生物活性预测（pIC50）

有效的抑制浓度（pIC50）值是一个重要的指标，表示将生物靶点的活性降低一半所需的浓度69。随着时间的推移，许多研究致力于改进这些预测策略，使得更容易聚焦于最有希望的药物候选物70–72。在这一背景下，本研究应用了机器学习（ML）技术来预测选定化合物对pKHG的生物活性。

2.3.4.1. 数据集的整理与准备

为了构建数据集，我们从ChEMBL数据库中获取了生物活性记录和化学结构，这是一个用于QSAR研究的知名平台。然后通过仅保留报告了对pKHG具有IC50值的条目来整理化合物。在选择最终的ChEMBL数据库时，我们主要关注了蛋白质类型、生物体、化合物活性的可用性，包括标准类型（IC50）和标准单位（nm）。这种精细的筛选产生了一个可靠且生物学上相关的分子数据集，作为回归建模和后续分析的基础。

2.3.4.2. 分子描述符的计算

分子描述符在QSAR建模中起着核心作用，因为它们将化合物的结构和物理化学特性转化为可测量的值，从而能够识别有意义的模式73。在本研究中，使用RDKit为每个化合物计算了Morgan指纹和ECFP4描述符。这些描述符提供了化合物结构特征的详细数值概况，这对于构建预测模型以估计QSAR建模中的生物活性至关重要。

2.3.4.3. 数据集划分为训练集和测试集

最终包含987个化合物的数据集通过scikit-learn的train_test_split函数进行了划分，根据类别标签进行分层，以保持两个子集中原始的类别分布。数据集被分为789个分子在训练集中和198个分子在测试集中，这在化学信息学和其他机器学习任务中是常见的做法，通常使用20–30%的数据来测试算法。随机或分层随机划分在化学信息学建模中作为基线被广泛使用，尽管最近的研究强调更严格的分割策略（例如基于支架或聚类的分割）可能为化学数据集提供更现实的外部泛化评估74,75。

2.3.4.4. 机器学习模型的开发与验证

为了进行预测建模，我们选择了最佳的回归模型，该模型因其速度快、准确性高、可解释性强以及测试和训练数据集的确定系数（R2）值而受到重视。为了提高模型的效果，应用了递归特征消除（RFE）来去除信息量较少的特征，减少噪声并使模型更易于解释。在这里，我们使用optuna来调整超参数。模型的预测性能使用广泛采用的指标进行了评估，包括平均绝对误差（MAE）、均方误差（MSE）、均方根误差（RMSE）和确定系数（R2）。这些指标是根据先前发表的公式计算的。

2.3.4.5. 生物活性（pIC50）预测的模型实现

在模型开发之后，我们应用了经过训练和微调的模型来预测收集到的化合物对pKHG的pIC50值。那些预测的pIC50值最高的化合物（>6.5），显示出更大的潜力，被考虑进行进一步的研究。然后通过分子对接研究分析了这些选定的化合物，以评估它们与pKHG的结合亲和力和相互作用。

2.3.5. 通过AutoDock Vina进行分子对接

分子对接在发现和优化潜在药物候选物方面非常重要。它可以用来估计两个分子（如蛋白质和配体）之间相互作用的可能性，并提供关于它们结合行为和潜在效力的信息76。因此，在本研究中，我们使用PyRx77来进行所选化合物与目标蛋白质的分子对接。PyRx提供了一个易于使用的图形界面，与AutoDock Vina结合使用，使得对接过程更加简便。它还结合了AutoDock Vina和Open Babel，提供了一个完整的分子对接和分析套件78。最后，我们通过PyMOL和BIOVIA Discovery Studio 2021分析了配体与目标受体之间的相互作用，使我们能够更好地理解它们的结合模式。

2.3.6. 通过Schr?dinger软件进行分子对接验证

在AutoDock Vina中的首次对接之后，根据结合亲和力将候选化合物进一步筛选，并选择了一些化合物进行GLIDE模块的再评分。首先，这些潜在的命中化合物使用Schr?dinger套件的LigPrep模块生成并进行了优化79。蛋白质结构进行了预处理，网格分别通过蛋白质准备向导和受体网格生成模块生成。这包括使用Prime校正缺失的氢原子，使用Epik在生理pH 7.0 ± 2.0下关闭侧链间隙和环，以及创建合适的质子化状态。最后，使用Maestro（版本11.8.012）软件的GLIDE模块在Extra Precision（XP）模式下将验证过的配体重新对接到目标蛋白质的活性位点中，确保只有可靠的目标在结合相互作用方面得到更详细的研究。

2.3.7. 嵌接后的MMGBSA

为了更好地理解配体-蛋白质相互作用，使用了对接后的精细化方法，特别是分子力学广义Born表面面积（MM-GBSA）。虽然对接分数提供了结合亲和力的初始值，但它们没有包括所有的溶剂化效应或熵贡献，这些可能会影响配体结合的稳定性和强度80。MM-GBSA通过提供更准确的结合自由能（ΔG）估计来克服这些缺点，从而更准确地反映了配体与受体之间的亲和力81。在这项研究中，MM-GBSA分析仅应用于在extra precision（XP）模式下获得良好结果的配体，使用了Schr?dinger套件中的OPLS4力场（PRIME模块）。这种选择性方法是必要的，因为并非所有配体在XP对接过程中都会产生可靠或有意义的分数，这由于结合模式和结构灵活性的差异82。每个配体的结合自由能使用以下公式计算：
ΔGbind = ΔGcomplex ? (ΔGreceptor + ΔGligand)

在这里，ΔGcomplex代表配体-蛋白质复合物的自由能，ΔGreceptor是受体的自由能，ΔGligand是配体的自由能。更负的ΔG值对应于更强的结合相互作用83。分子动力学模拟与模拟后的MM-GBSA计算

分子动力学（MD）模拟使用Schr?dinger的Desmond模块进行，以研究蛋白质-配体复合物的稳定性和构象动态。系统准备使用了System Builder向导完成。每个复合物都溶解在正交模拟盒内的简单点电荷（SPC）水模型中，确保蛋白质表面与盒子边缘之间的最小距离为10 ?。为了模拟生理离子条件，系统被中和，并补充了0.15 M的Na+和Cl?离子。进行了100 ps的能量最小化处理，以消除不利的接触，然后在300 K和1 atm的条件下，通过NVT和NPT集合进行平衡。生产MD模拟持续了100 ns，时间步长为2 fs。轨迹以20 ps的间隔保存，总共获得了5000帧数据。我们使用了OPLS3e力场，该力场为生物分子和类似药物的化合物提供了改进的参数化。模拟轨迹分析了RMSD、RMSF、SASA、回转半径、氢键以及主成分动态，以评估构象稳定性。然后使用gmx_MMPBSA包计算了选定复合物的基于MM-GBSA的结合自由能（ΔG_bind）。这种方法整合了范德华力、静电能和溶剂化能成分，以及溶剂可及表面积（SASA）的贡献，提供了配体亲和力的热力学估计。Desmond轨迹文件使用Schr?dinger的工具转换为GROMACS兼容格式，而相应的拓扑文件通过InterMol将*.cms文件转换为*.gro和*.top格式。结合相关的自由能（ΔGBind）根据以下关系确定：
ΔGbind = 〈GPL〉 ? 〈GP〉 ? 〈GL〉

在这个方程中，〈GPL〉、〈GP〉和〈GL〉分别对应于复合物、未结合蛋白质和自由配体的平均自由能。因此，总结合能量由以下公式给出：
ΔGbind = ΔEMM + ΔGSOLV ? TΔS

2.3.9. 药代动力学（ADME & 毒性）评估

在完成分子对接和模拟研究后，有必要检查表现最佳的化合物是否具有适合实际药物开发的特性。药代动力学分析基于SwissADME网络工具（https://www.swissadme.ch/）进行，该工具可以预测重要的ADME参数（吸收、分布、代谢和排泄）以及物理化学性质，包括溶解度、亲脂性和分子柔韧性。在这项分析中，筛选出在效力和药代动力学可行性方面具有最佳平衡的化合物。为了补充这些结果，使用了ProTox-III服务器（https://tox-new.charite.de/protox_III）来预测不同类型的毒性，如肝毒性、肾毒性、心毒性和神经毒性。它使用深度机器学习算法训练，以便对毒性估计更有信心。

3. 结果

3.1. 目标识别

3.1.1. 差异基因表达分析及DEG选择

为了研究肺癌中的转录变化，分析了四个独立的GEO微阵列数据集（GSE19804、GSE10072、GSE18842和GSE10799）。我们发现肺癌与正常组织之间存在广泛的转录变化，因为在所有四个数据集中许多基因都显著上调或下调。表S1列出了每个数据集中的这些基因。DEG可视化采用火山图形式呈现（图2A），显示出肺癌状态和正常组织中显著上调和下调基因的明显差异。为了寻找强烈的分子特征，我们进一步识别了四个数据集之间的共同DEGs。交集分析显示存在共同的DEGs（cDEGs），这些基因被整合到一个韦恩图中（图2B）。表S2列出了这些cDEGs的详细列表。这些共享基因可能是肺癌发病机制中的关键候选者，并被优先考虑进一步分析其功能。图2

(A) 火山图显示每个数据集中显著上调和下调基因的分布。(B) 韦恩图显示四个数据集之间DEGs的重叠部分以及共同的DEGs（cDEGs）。

3.1.2. 通过PPI网络分析识别关键枢纽蛋白

利用共同DEGs（cDEGs）构建了蛋白质-蛋白质相互作用（PPI）网络，该网络包含375个节点和349条边，平均节点度为1.86，平均局部聚类系数为0.299。预期的边数为222，PPI富集的p值为2.55 × 10?15，这意味着观察到的相互作用远超过偶然发生的概率。所有这些PPI网络都显示在图3A中，其中节点代表单个蛋白质，边代表实验验证的相互作用。图3

肺癌中共同DEGs（cDEGs）的蛋白质-蛋白质相互作用（PPI）网络。(A) 整体PPI网络，节点表示蛋白质，边表示实验验证的相互作用。(B) 使用八种不同的拓扑方法（Degree、MNC、MCC、DMNC、EPC、Bottleneck、EcCentricity和Closeness）识别的前10个枢纽基因。为了识别网络中最具影响力的蛋白质，我们应用了CytoHubba插件进行了八种拓扑分析。每种方法都确定了前10个枢纽基因。这些分析的结果总结在图3B中，突出了在多个指标中始终作为中心枢纽出现的蛋白质。这些关键枢纽基因（kHGs）可能在肺癌生物学中发挥关键作用，并代表进一步研究的潜在目标。在所有识别的kHGs中，CDK1表现出最多的蛋白质-蛋白质连接，因此其相互作用网络在图3A中特别突出。

3.1.3. kHGs的转录和转录后调控分析

在使用多种拓扑算法（八种不同方法）识别出枢纽基因后，我们观察到尽管许多枢纽基因在大多数方法中都是共有的，但某些基因仅由个别算法识别出来。每种中心性方法强调网络的不同属性，因此捕捉到了网络生物学的互补方面。为了减少任何单一算法的偏倚，并确保不遗漏潜在相关的候选者，我们考虑了从所有方法获得的枢纽基因的并集（n = 26）。为了更好地理解识别出的kHGs是如何被调控的，我们进行了综合分析，重点研究了它们与转录因子（TFs）和微小RNA（miRNAs）的相互作用。如图4A和图4B所示，网络分别突出了kHG–TF和kHG–miRNA的相互作用。从网络的拓扑评估中，我们确定了五个关键的TFs（FOXC1、GATA2、YY1、E2F1和HINFP）以及五个主要的miRNAs（hsa-miR-192-5p、hsa-miR-92a-3p、hsa-miR-193b-3p、hsa-miR-215-5p和hsa-miR-155-5p）作为中心调控因子。这些分子似乎在转录和转录后水平上控制基因表达中起着关键作用。图4

3.1.4. KHGs的GO和通路富集分析

为了阐明肺癌中识别出的关键枢纽基因（KHGs）的生物学意义，在DAVID、EnrichR和GeneCloudOmics中进行了基因本体学（GO）和KEGG通路的富集分析。在DAVID中检测到的所有术语以及那些在所有验证平台中都存在的术语都被保留下来，从而得到了肺癌KHGs的完整功能类别。富集的类别被分类为生物过程（BP）、细胞组分（CC）、分子功能（MF）和信号通路，相应的详细列表在SI（表S3）中呈现。在GO类别中，观察到与细胞周期相关的功能富集显著。有丝分裂细胞周期进展、纺锤体组织、染色体分离和检查点信号通路在生物过程（BP）中高度富集，AURKA、CDK1、CCNB1、NDC80和CDC20等基因被反复提及。这些结果表明，有丝分裂保真的破坏是肺癌进展的关键标志之一。细胞组分（CC）组中的KHGs富集到已知的有丝分裂结构，例如着丝粒、中心体和有丝分裂纺锤体，而NDC80、BUB1B和AURKA等KHGs被识别为这些富集的主要驱动因素。这表明肺癌细胞中的染色体分离机制存在结构失调。在分子功能（MF）中，激酶活性（蛋白激酶和丝氨酸/苏氨酸激酶）和结合ATP的功能最为富集，AURKA、CDK1、NEK2和TTK占据了中心位置。还出现了多种组蛋白激酶活性，这些枢纽基因可能在有丝分裂相关的表观遗传调控中发挥作用。与GO结果一致，KEGG通路分析显示了14条高度富集的通路，其中细胞周期通路（p = 3.83 × 10?7）与多个KHGs（CDK1、CCNB1、CCNB2、BUB1B、CDC20和TTK）的关联最强。与肿瘤抑制功能和细胞应激反应相关的通路也显著富集，包括p53信号通路和FoxO信号通路，这些通路也包括调节凋亡、衰老和基因组稳定性的基因。

3.2. 目标验证

基于PPI网络中枢纽基因排名的综合证据及其在重要通路中的强烈参与度，CDK1被优先考虑作为下游验证的主要关键枢纽基因（pKHG）。

3.2.1. 在肺癌中评估CDK1的转录和蛋白质组表达水平

为了验证CDK1作为主要关键枢纽基因，我们首先使用TIMER2.0评估了其在泛癌症数据集中的表达谱。如图5A所示，CDK1在多种肿瘤类型中的表达显著升高，包括肺腺癌（LUAD），p < 0.001，表明具有很强的统计显著性。使用GEPIA2在LUAD中的进一步验证确认了CDK1的上调。箱形图分析（图5B-i）显示LUAD组织中的CDK1转录水平显著增加，而阶段图（图5B-ii）显示CDK1的表达在病理阶段（I–IV）中仍然显著高（p = 0.000579）。这些发现表明CDK1的过表达在整个疾病进展中得以维持。此外，基于UALCAN的分析在转录组和蛋白质组水平上都提供了一致的证据。TCGA RNA-seq数据显示LUAD肿瘤中的CDK1表达明显高于正常样本（图5C-i），CPTAC蛋白质组数据同样证实了肿瘤组织中CDK1的蛋白质表达高于对照组（图5C-ii）。总体而言，这些结果表明CDK1在LUAD中在mRNA和蛋白质水平上都显著上调，增强了其作为生物学相关的主要关键枢纽基因的候选资格。图5

肺腺癌（LUAD）中CDK1的转录和蛋白质组表达谱。(A) 使用TIMER2.0对TCGA数据集中的CDK1表达水平进行的泛癌症分析，显示LUAD中的表达显著升高（***p < 0.001）。(B) 使用GEPIA2验证LUAD中的CDK1表达：(i) 肿瘤与正常组织的箱形图分析；(ii) 显示病理阶段（p = 0.000579）中的表达阶段图。(C) 基于UALCAN的LUAD中CDK1的验证：(i) TCGA RNA-seq数据确认肿瘤中的转录水平较高；(ii) CPTAC蛋白质组数据显示原发肿瘤中的蛋白质水平升高。

3.2.2. CDK1在肺癌中的生存分析

为了探索CDK1表达在LUAD中的临床相关性，我们使用GEPIA2平台进行了Kaplan–Meier生存分析。根据中位切割值将患者分为高表达组和低表达组。在整体生存（OS）分析中，CDK1表达高的患者与表达低的患者相比预后明显较差（图6A）。差异具有统计学意义（log-rank p = 2.6 × 10?5），风险比（HR）为1.9，表明CDK1的表达升高几乎使LUAD中的死亡风险翻倍。与此一致，无病生存（DFS）分析也表明CDK1水平高的患者复发更早，DFS更短（图6B）。这种关联具有显著性（log-rank p = 0.027；HR = 1.4），进一步证实了CDK1过表达在疾病进展中的不利作用。这些发现共同强调了CDK1作为LUAD（左下肺腺癌）中潜在的预后标志物的作用，其过表达与生存率降低和复发风险增加有关。图6。

根据CDK1表达分层对LUAD患者进行Kaplan–Meier生存分析。（A）总生存期（OS）；（B）无病生存期（DFS）。

3.2.3. 评估与CDK1相关的免疫浸润

利用TCGA（肿瘤遗传组学数据库）的信息，我们研究了CDK1表达与LUAD中免疫细胞浸润之间的联系。CD8+ T细胞浸润分析数据显示，CDK1表达水平与T细胞CD8+浸润之间存在显著的正相关关系，这一关系通过多种去卷积算法得到证实。热图（图S1A）中有一条明显的带状区域显示了这种关联，该带状区域显示出显著的统计学意义。其中，T细胞CD8+（naive_XCELL）浸润与CDK1表达之间的正相关性最强，Rho值为0.235，p值为1.35 × 10^-7。另一方面，T细胞CD8+（EPIC）与CDK1表达之间的负相关性最强，Rho值为-0.207，p值为3.46 × 10^-6。在LUAD中的CD4+ T细胞浸润方面，热图（图S1B）同样显示了一条橙色带状区域，代表了LUAD特异性的关联。CD4+ T细胞（Th2XCELL）与CDK1表达之间的相关性最强，Rho值为0.813，p值为2.96 × 10^-117。而T细胞CD4+（central memory_XCELL）与CDK1表达之间的负相关性最强，Rho值为-0.322，p值为2.27 × 10^-13。在LUAD中的巨噬细胞浸润方面，热图（图S1C）同样显示了一条橙色带状区域，表明了相关的关联。其中，巨噬细胞（M1_CIBERSORT）浸润与CDK1表达之间的正相关性最强，Rho值为0.367，p值为3.39 × 10^-17。巨噬细胞（M2_TIDE）与CDK1表达之间的负相关性最强，Rho值为-0.346，p值为2.45 × 10^-15。综合这些分析，我们发现CDK1表达与LUAD中的免疫细胞浸润之间存在强烈的关联。特别是，CD8+ T细胞和CD4+ T细胞亚群之间存在正相关和负相关，而CD4+ Th2细胞（一种体液免疫相关亚型）与CDK1表达之间的负相关具有条件依赖性。此外，CDK1表达与M1巨噬细胞浸润（促炎/抗肿瘤表型）呈正相关，而与M2巨噬细胞浸润（免疫抑制/促肿瘤表型）呈负相关。这种M1和M2巨噬细胞之间的转换，也被称为巨噬细胞极化，是肿瘤微环境巨噬细胞的功能重编程。总的来说，这些结果表明CDK1可以影响LUAD的免疫微环境，并可能对肿瘤进展和治疗敏感性产生影响。

3.2.4. 研究CDK1的突变和改变

为了研究人类癌症中CDK1的基因组特征，我们使用TCGA Pan Cancer Atlas的数据查询了cBioPortal数据库。研究发现，子宫内膜癌（UCEC）中的CDK1改变频率最高，接近8%，主要由扩增引起。其次是皮肤黑色素瘤，约为3%，胆管癌和子宫内膜癌为2.5–2.8%，乳腺浸润癌中的改变频率较低（图S2A）。总体而言，突变是各组中最常见的改变类型，结构变异（紫色）、扩增（红色）和深度缺失（蓝色）较少见，但在LUAD和BRCA等癌症中较为显著。这些改变的存在，尤其是扩增，使得CDK1成为一部分LUAD患者的重要致癌驱动因素。这表明针对CDK1的靶向治疗（如选择性CDK抑制剂）可能是对携带这些特定基因改变的患者的有效治疗策略。通过“Mutations”模块的进一步分析，生成了CDK1突变的示意图（图S2B），将改变映射到297个氨基酸的蛋白质序列上。错义突变（绿色）占主导地位，而截短（黑色）、剪接（棕色）和融合（紫色）变异较为少见。详细突变分布见表S4。这些数据显示，CDK1在LUAD中可以发生突变，尤其是预计具有功能损伤的错义突变。这项分析在来自癌症基因组图谱（TCGA）的肺腺癌（LUAD）患者队列中识别并描述了三个非同义的体细胞突变（S178L、I136N和E57V）。尽管这些突变目前尚未在主要临床数据库（OncoKB和CIViC）中被标注为已知驱动因素，但计算机模拟功能预测工具表明其中两个突变（I136N和E57V）可能是致病的。这些突变可能成为LUAD患者对CDK抑制剂敏感性的候选生物标志物。

3.3. 基于结构的药物发现，针对CDK1

3.3.1. 目标蛋白的结构制备和优化

CDK1蛋白的三维晶体结构（PDB ID：6GU7）是从RCSB蛋白质数据库获得的。在将其用于计算实验之前，按照方法部分所述进行了结构优化处理，包括去除杂原子和水分子，然后进行能量最小化，以确保获得稳定且适合后续分子建模分析的构象。

3.3.2. 化合物筛选和选择结果

在这项研究中，我们从IMPPAT数据库中选择了33种天然药用植物作为化合物来源，以识别潜在的CDK1抑制剂（表S5）。最初从数据库中检索了近9,667种化合物及其IMPPAT ID和SMILES。Lipinski的五规则（RO5）分析显示，在5,786种植物化学物质中，有4,236种的违规项为零（0），表明其具有良好的药物相似性特征。去除重复项后，从所有选定的药用植物中保留了2,113种独特的化合物（表S6）。然后选择这些化合物进行进一步的计算分析。

3.3.3. 使用机器学习方法进行生物活性预测（pIC50）

3.3.3.1. 模型开发、优化和验证

为了计算我们选定的化合物对CDK1的生物活性，我们选择了包含3,348种实验验证的抑制剂和4,536种CDK1活性数据的ChEMBL数据集（CHEMBL308），这些数据于2025年11月20日获取。目标是开发一个经典的回归模型，以便我们可以根据训练好的模型计算化合物库的IC50和pIC50值。根据我们的研究，随机森林（RF）、LightGBM、XGBoost和Stacking模型表现最佳。评估指标显示，RF、LightGBM、XGBoost和Stacking的R2值分别为0.177、0.117、0.241和0.214（图S3A）。因此，基于这些发现，我们观察到LightGBM在测试数据和训练数据之间的R2值差异最小，表明其预测结果最准确。LightGBM回归模型由于出色的预测准确性、计算效率以及处理大量数据的能力而受到了广泛关注。作为决策树的集成模型，LightGBM非常适合捕捉非线性结构-活性关系，同时能够在大型分子描述符和指纹集上快速训练。LightGBM模型使用优化的超参数在训练集上进行训练，包括n_estimators、num_leaves、max_depth、learning_rate、subsample_bytree、min_child_samples、reg_alpha和reg_lambda等参数。此外，为了提高模型的泛化能力和减少过拟合，还对模型进行了测试集的评估，使用了标准的QSAR回归指标：平均绝对误差（MAE）、均方误差（MSE）、均方根误差（RMSE）和决定系数（R2）。这些参数全面评估了模型的预测结果与实际实验数据的吻合程度。测试模型的表现令人印象深刻，MAE、MSE、RMSE和R2值分别为0.520、0.544、0.737和0.758（见图S3A）。总体而言，LightGBM的准确性相对较高，因此我们最终选择了它来预测选定天然化合物的生物活性（pIC50）值。

3.3.3.2. 应用于生物活性计算和化合物筛选

一旦验证了模型的预测性能，我们就使用优化后的机器学习模型对一个包含2,113种植物化合物的库进行了虚拟筛选，预测它们的pIC50值。共有380种化合物的预测pIC50 > 6.5（≈ IC50 < 1 μM），表明它们可能对目标具有抑制作用，并将它们列为后续验证的优先对象（见表S7）。由于实验值和预测值在散点图中非常接近（见图S3B），这表明QSAR模型具有很强的预测准确性，大多数数据点都紧密地沿着回归线分布。这些结果表明该模型在识别有潜力进一步研究的分子方面非常有用，这些分子可以通过分子对接或实验研究来确认其作为潜在药物靶点的可行性。

3.3.4. 通过AutoDock Vina进行分子对接分析

在先前选定的380种植物化学物质中，共有358种的3D SDF结构从PubChem数据库中提取，并使用PyRx中的AutoDock Vina工具进行了分子对接，以CDK1蛋白（6GU7）为目标。对接过程使用以蛋白质活性位点为中心的网格框进行，坐标为X = 22.0351, Y = 12.6323, Z = 4.7655。此外，还将AZD5438这种已知的CDK1抑制剂与AZD5438的对接亲和力进行比较，并将其对接亲和力作为我们的截止值（?8.8 kcal mol?1）。最后，我们选择了29种化合物进行进一步的分析与评估，以确定它们的结合强度（表S8）。

3.3.5. 使用Schr?dinger软件进行分子对接验证

为了进一步提高准确性，使用Schr?dinger的Maestro环境中的GLIDE模块对选定的29种化合物进行了额外的精度（XP）对接。在最初对接的化合物中，Top 5配体（CID_115196、CID_5317764、CID_14218027、CID_487089和CID_174880）在AutoDock Vina中的结合亲和力低于?9.2 kcal mol?1，被选为这一验证步骤的候选对象。它们的XP对接结果显示，结合亲和力范围在?6.80至?7.85 kcal mol?1之间（表2）。特别是，CID_115196、CID_5317764和CID_14218027的Maestro XP对接得分分别为?7.85、?7.70和?7.69 kcal mol?1，表明它们与目标蛋白的相互作用最强。此外，CID_487089（?6.80 kcal mol?1）和CID_174880（?6.70 kcal mol?1）的结合亲和力也优于参考化合物AZD5438（?5.15 kcal mol?1）。因此，这些Top 5化合物被认为是进一步计算的候选药物分子。

3.3.6. 嵌接后的MMGBSA分析

为了评估基于对接分数的筛选结果，我们进一步使用MM-GBSA方法估计了选定蛋白-配体复合物的结合自由能（ΔG_bind）。分析显示，CID_174880展现了最有利的结合能（?50.44 kcal mol?1），表明在所有评估的化合物中预测的结合稳定性最高。同样，CID_487089也显示出较强的结合亲和力，ΔG_bind值为?46.78 kcal mol?1。CID_14218027（?36.28 kcal mol?1）和CID_5317764（?36.15 kcal mol?1）的结合能适中但稳定。相比之下，CID_115196在测试的植物化学物质中显示出最弱的结合能（?27.81 kcal mol?1）。值得注意的是，参考抑制剂AZD5438显示出较强的结合自由能，为?44.58 kcal mol?1，与其报道的抑制活性一致。特别是，几个候选化合物，尤其是CID_174880和CID_487089，显示出比对照组更有利的MM-GBSA能量，而其他化合物如CID_14218027和CID_5317764的结合性能相当。总体而言，这些发现表明排名靠前的化合物具有等于或优于参考抑制剂的结合亲和力，支持它们作为进一步稳定性和动态分析的潜在候选者（表3）。最后，基于上述发现，我们选择了前三种化合物进行进一步分析，以验证它们的结合强度和稳定性。表3显示了与对照组（AZD5438）相比的top-docked化合物的MM-GBSA结合自由能。

3.3.7 蛋白质与配体之间的分子相互作用评估

所有对接的配体都通过氢键和疏水相互作用的结合，在CDK1活性位点表现出稳定的结合，这与对照化合物AZD5438的相互作用行为非常相似。对照化合物形成了四个方向良好的氢键（THR14、LYS130、ASN133和LEU83），以及多样的疏水接触，建立了一个强力的基准模型。在植物化学物质中，6GU7_487089显示出最具竞争性的相互作用模式，与关键残基如LYS130、ASP86、LEU83和ASN133形成了多个氢键，并且与ILE10、ALA31和LEU135有密集的疏水相互作用，表明其结合模式高度稳定且紧密。同样，6GU7_174880通过四个氢键（LEU83、LYS89、ASP86和GLU81）以及与ALA31、ILE10和LEU135的π-烷基和烷基相互作用表现出强烈的锚定作用，尽管其疏水网络相对不如6GU7_487089广泛。相比之下，6GU7_14218027显示出平衡但稍弱的结合模式，氢键主要集中于LYS33、LYS89和GLN132，疏水相互作用涉及ALA31、VAL18和PHE80。总体而言，虽然AZD5438保持了最多样化的相互作用模式，但6GU7_487089在结合稳定性和相互作用密度方面是最接近的竞争者，其次是6GU7_174880，而6GU7_14218027显示出相对中等的相互作用强度。这些发现表明某些植物化学物质，特别是6GU7_487089，有潜力模仿对照配体在CDK1活性位点的结合效率。完整的相互作用细节在表S9中提供，相应的3D和2D相互作用图在图7中展示。

3.3.8 通过重新对接验证对接协议

对接协议的验证通常是通过将共结晶的配体重新对接回目标蛋白质的活性位点来进行的，以检查对接方法的可靠性。这个过程的准确性是通过计算实验确定的配体位置与重新对接后的姿势之间的均方根偏差（RMSD）来衡量的。在这项研究中，使用Maestro平台将共结晶配体AZD5438重新对接到CDK1的结合口袋中（PDB ID: 6GU7）。最佳对接姿势在与原始配体构象叠加时显示出1.360 ?的RMSD值（图8）。由于这个RMSD值远低于通常接受的≤2.0 ?的阈值，它表明对接协议能够准确再现实验观察到的结合模式。总体而言，这一验证确认了对接设置的可靠性，并适合进一步分析结合亲和力和分子相互作用。图8展示了共结晶配体（绿色）和重新对接的共结晶配体（蓝色）之间的重新对接过程和叠加情况，使用的是GLIDE对接方法。

3.3.9 分子动力学模拟

为了理解所有四个CDK1复合物——6GU7_487089、6GU7_174880、6GU7_14218027和6GU7_AZD5438（对照）在类似于人体条件下的行为，对每个复合物进行了100纳秒的分子动力学（MD）模拟。通过这种模拟，我们可以跟踪结构随时间的演变，以及配体如何保持与蛋白质的结合。为了清楚地了解它们的稳定性，我们分析了一些重要指标，如RMSD、RMSF、回转半径（Rg）和溶剂可及表面积（SASA）。此外，还进行了主成分分析（PCA）和自由能景观（FEL）映射，以描述每个复合物在100纳秒轨迹中访问的最稳定构象状态。

3.3.9.1 均方根偏差（RMSD）

均方根偏差（RMSD）评估了对接复合物的结构完整性和构象变化，较高的值表明稳定性降低，较低的值反映系统性能稳定。RMSD分析显示，无配体的原型蛋白质（6GU7）显示出最高的平均结构偏差（3.245 ?），这是在没有配体情况下最高的构象灵活性。相比之下，所有结合配体的复合物都表现出明显较低的平均RMSD值，表明配体结合后结构稳定性增强。在筛选的化合物中，CID_487089的平均RMSD最低（2.446 ?），其次是CID_14218027（2.470 ?），表明这些化合物在结合口袋中有强烈的稳定作用。CID_174880显示出与CID_14218027相当的稳定性，平均RMSD为2.476 ?。对照配体的平均RMSD为3.002 ?，与原型形式相当，意味着稳定的作用较弱。总体而言，平均RMSD结果确认CID_487089在MD模拟期间与6GU7形成了最稳定的复合物（图9A）。

3.3.9.2 均方根波动（RMSF）

通过RMSF分析在模拟的平衡阶段检查了系统的残基级动态。正如预期的那样，无配体的原型蛋白质显示出最高的平均波动（1.571 ?），反映了其在没有结合配体时的更大结构自由度。相比之下，所有结合配体的复合物都表现出明显较低的平均RMSF值，表明配体相互作用有助于形成更稳定和有序的残基环境。在筛选的化合物中，CID_174880显示出最强的稳定效应，其平均RMSF值最低（1.181 ?），表明结合紧密且有效地限制了残基的移动性。对照化合物显示出适度的灵活性降低（1.390 ?），表明有合理的但不太明显的稳定作用。同时，CID_487089和CID_14218027显示出中等程度的稳定作用，平均RMSF值分别为1.253 ?和1.375 ?，表明结合相互作用平衡，减少了灵活性但不会过度限制蛋白质动态（图9B）。

3.3.9.3 溶剂可及表面积（SASA）

SASA反映了蛋白质的溶剂暴露表面，较低的值通常表示配体埋藏得更深，可能意味着结合更强。对6GU7蛋白质复合物的分析揭示了配体之间的溶剂可及表面积（SASA）存在显著差异。6GU7_487089复合物显示出最高的平均SASA值（148.17 ?2），表明溶剂暴露表面相对较大，暗示结合构象较为宽松。相比之下，6GU7_174880复合物显示出最低的平均SASA值（126.15 ?2），表明配体在结合口袋中埋藏得更深，形成了更为紧凑和稳定的蛋白质-配体界面。6GU7_14218027（130.14 ?2）和对照复合物（131.81 ?2）显示出中等的SASA值，反映了结构排列较为紧凑且溶剂暴露平衡。总体而言，6GU7_174880在结构紧凑性方面显示出最有利的SASA特征，突出了其作为最有效稳定配体的潜力（图9C）。

3.3.9.4 回转半径（Rg）

回转半径（Rg）是分子动力学模拟中蛋白质结构紧凑性和稳定性的关键指标。在这项研究中，对照复合物（6GU7/control）显示出相对较高的平均Rg值（4.60 ?），表明构象较为松散。同样，6GU7_487089复合物在结合配体的系统中显示出最高的Rg值（4.69 ?），表明紧凑性较低且结构稳定性较弱。相比之下，6GU7_14218027（4.40 ?）和6GU7_174880（4.44 ?）显示出较低的Rg值，表明构象更为紧凑且稳定性较弱。在这些复合物中，6GU7_14218027显示出最低的Rg值，表明结构紧凑性最高，暗示了对蛋白质有更强的稳定作用。总体而言，配体结合通常增强了蛋白质的紧凑性，基于Rg分析，6GU7_14218027被认为是最有效的稳定剂，其次是6GU7_174880（图9D）。

3.3.10 模拟后的MMGBSA分析

在100纳秒MD模拟期间，对6GU7_487089、6GU7_174880、6GU7_14218027和参考配体复合物（6GU7_AZD5438）与CDK1（PDB ID: 6GU7）进行了模拟后的MMGBSA分析（图10）。6GU7_487089复合物显示出最有利的结合自由能（?40.29 kcal mol?1），表明与目标蛋白质的结合最强且最稳定。其次是6GU7_174880（?36.06 kcal mol?1），也显示出相对较强的结合亲和力。相比之下，6GU7_14218027（?29.13 kcal mol?1）和对照复合物（?29.53 kcal mol?1）显示出较高的ΔG值，表明结合相互作用较弱。值得注意的是，6GU7_14218027的结合亲和力略低于对照物，表明相对于参考配体有所改进。总体而言，6GU7_487089作为最具潜力的候选者脱颖而出，具有最高的结合亲和力，其次是6GU7_174880。这些发现与在分子动力学模拟中观察到的有利稳定性轮廓一致，进一步支持了它们作为有效抑制剂的潜力。

3.3.11 主成分分析

为了了解每个配体如何影响蛋白质的整体运动，我们进行了主成分分析（PCA）。在研究的复合物中，6GU7_174880在PC1中的贡献最高（44.97%），这与对照复合物6GU7_AZD5438（44.18%）非常接近。这表明在配体CID_174880存在的情况下，蛋白质的大部分运动主要集中在一个主导方向上。这种模式通常表明构象运动更为稳定和明确，类似于参考抑制剂的观察结果。另一方面，6GU7_487089和6GU7_14218027的PC1值较低，分别为31.37%和30.43%。同时，这两个复合物在PC2和PC3中的贡献相对较高，特别是PC2值分别为16.13%和17.88%。这可能表明蛋白质运动分布在多个方向上，表明在模拟期间可能存在更多的构象变化。对于对照复合物，前三个主成分的贡献分别为44.18%（PC1）、10.85%（PC2）和7.21%（PC3），反映了蛋白质-配体复合物的稳定动态行为。有趣的是，6GU7_174880的PCA轮廓与对照物非常相似，这可能支持其作为具有有利结合稳定性的有希望的配体的潜力。相比之下，6GU7_487089和6GU7_14218027较高的PC2和PC3值表明这些复合物可能与目标蛋白质有更灵活的相互作用模式（图S4）。

3.3.12 吉布斯自由能景观（FEL）分析

进行吉布斯自由能景观（FEL）分析，以探讨6GU7蛋白质与选定配体及参考对照复合物的构象稳定性和能量最小值。在FEL图中，深蓝色区域代表最低能量状态，表明在模拟过程中采样的最稳定构象。在研究的复合物中，6GU7_174880（图S5B）和对照复合物6GU7_AZD5438（图S5D）显示出定义明确且能量分布相对紧凑的深能量盆地。这表明这两个复合物在整个模拟过程中保持在更稳定的构象状态，并经历了较少的构象波动。配体CID_174880与对照物之间的相似性表明该化合物可能诱导出与参考抑制剂相当的稳定结合构象。相比之下，6GU7_487089（图S5A）显示出一个较宽的低能量盆地，在构象空间中有较广泛的分布。这可能表明该复合物在模拟过程中探索了多个构象状态，表明其动态行为更为灵活。同样，6GU7_14218027（图S5C）也显示出扩展的能量盆地，低能量区域有明显的扩散。这种模式可能反映了存在多个亚稳态构象和相比于对照复合物更大的结构灵活性。总体而言，FEL结果表明6GU7_174880表现出更稳定的构象景观，与对照物相似，而6GU7_487089和6GU7_14218027显示出相对较宽的能量最小值，表明结合相关的运动更为灵活（图S5）。

3.4 药代动力学研究

3.4.1 物理化学和ADME性质预测

评估了选定植物化学物质的物理化学和ADME性质，并与参考对照化合物AZD5438进行比较，以评估它们的药物相似性和药代动力学适宜性。所有选定化合物的分子量均在口服药物候选物的可接受范围内，范围从340.37到424.49 g mol?1，与对照组（371.46 g mol?1）相当。同样，这些植物化学物质的拓扑极性表面积（TPSA）值（86.99–110.13 ?2）也在有利范围内，表明它们具有良好的膜渗透性和吸收潜力。各种化合物的logP值有所不同，其中CID_487089的亲脂性相对较高（4.35），而CID_174880（1.87）和CID_14218027（3.33）则处于可接受范围内。所有选定的植物化学物质都显示出与对照组相似的高胃肠（GI）吸收率，表明它们适合口服给药。此外，CID_487089和CID_14218027被预测为P-糖蛋白（P-gp）的非底物，而CID_174880被鉴定为P-gp底物，这可能会对其细胞外排行为产生轻微影响。所有化合物均未违反Lipinski五规则、Veber规则或Ghose过滤器规则（所有化合物均无违规），这再次表明它们具有良好的药物相似性。没有化合物被预测具有血脑屏障（BBB）通透性，这可能减少与中枢神经系统相关的副作用风险。所有选定化合物的生物利用度得分为0.55，与对照组相当，表明其具有良好的口服生物利用度。关于结构警报，CID_14218027的表现更为有利，没有PAINS警报，只有一个Brenk警报，而CID_487089和CID_174880的警报数量相对较多。在评估的化合物中，CID_14218027在所有主要CYP亚型（CYP1A2、CYP2C19、CYP2C9、CYP2D6和CYP3A4）上表现出广谱抑制作用，表明其存在较高的代谢干扰风险和药物相互作用风险。相比之下，CID_487089仅对CYP2C9和CYP3A4表现出选择性抑制作用，表明其相互作用风险相对较低。值得注意的是，CID_174880对所评估的任何CYP酶均无抑制活性，显示出较低的代谢相互作用风险。对照组化合物也抑制了CYP2C9和CYP3A4，进一步强调了从药代动力学和安全角度来看，具有最小或无CYP抑制作用的化合物更为理想。所有化合物的详细理化和ADME参数总结在表S10中，其整体ADME概况在图11中可视化。

3.4.2 毒性分析

毒性评估显示，所有三种选定的植物化学物质均表现出良好的安全性特征，未预测出肝毒性、神经毒性、心脏毒性、致突变性或细胞毒性。然而，所有化合物均显示出潜在的呼吸毒性，这与参考药物AZD5438相似，表明可能存在与目标相关的不良反应。值得注意的是，这些植物化学物质表现出免疫毒性，而对照化合物在这方面无活性，这表明所选配体的局限性。重要的是，没有化合物被预测具有致癌性，与对照药物形成对比，突显了其显著的安全优势。根据LD50值和毒性分类，CID_14218027被认为是最安全的候选物，急性毒性较低（第4类），而CID_487089的毒性较高（第2类），使其不太适合进一步开发。因此，进行全面的实验验证对于确认这些化合物的安全性和治疗适用性至关重要。

4. 讨论

肺癌是全球最具破坏性的恶性肿瘤之一。它生长迅速，通常发现时已经晚期，并且对现有治疗方法反应不佳。这些挑战突显了需要新的治疗方法，这些方法需要明确的策略来寻找目标和开发药物。通过转录组分析、网络分析、通路研究和基于计算机的药物发现，本研究系统地确定了肺癌中的一个关键分子靶点，并提出了针对该靶点的药物。转录组学是研究细胞中产生的所有RNA分子的方法，帮助我们理解哪些基因是活跃的以及它们的活性在不同条件下的变化。在癌症研究中，它常用于寻找重要靶点，因为癌细胞通常表现出与正常细胞不同的基因活性模式。然而，单一数据集的结果可能受到技术差异、实验平台或人群差异的影响。为了解决这个问题，分析了四个独立的GEO微阵列数据集（GSE19804、GSE10072、GSE18842和GSE10799），涵盖了不同的患者组（病例组和对照组）。通过这种联合分析，88个基因一致上调，294个基因下调。这些基因的反复变化表明它们与肺癌有真实的关联。在确定了这些差异表达基因后，我们研究了它们的蛋白质如何相互作用。蛋白质-蛋白质相互作用（PPI）网络的强烈富集表明这些连接不是随机的，显示出这些基因在肺癌中的真实生物学关系。然后，基于这个PPI网络，我们使用多种分析方法确定了主要的关键枢纽基因（kHGs），以确保我们的结果稳健且不依赖于单一方法。在这些相互作用中，周期依赖性激酶1（CDK1）显示出最多的相互作用，表明它在调节重要癌症相关过程中可能起核心作用。为了明确kHGs的调控机制，我们分析了它们与转录因子（TFs）和微小RNA（miRNAs）的相互作用。网络分析（图4）确定了FOXC1、GATA2、YY1、E2F1和HINFP为主要TFs，而hsa-miR-192-5p、hsa-miR-92a-3p、hsa-miR-193b-3p、hsa-miR-215-5p和hsa-miR-155-5p被确定为关键的miRNA调节因子。TFs和miRNAs的研究揭示了kHGs在不同水平上的调控机制。转录因子已知能激活基因的表达，而miRNAs则已知在基因转录后调节其活性。这种结合方法使我们能够确定可能导致肺癌中基因表达失调的关键调控因子，并为进一步的功能性和治疗研究提供了清晰的认识。此外，我们使用多个通路数据库进行了功能富集分析（GO和KEGG），以了解这些枢纽基因在肺癌细胞中的作用。这些富集结果表明，肺癌的进展是由细胞周期调控、有丝分裂结构和信号控制的广泛破坏所驱动的，而不是由孤立的基因改变引起的。CDK1在关键有丝分裂事件中的广泛参与、强烈的激酶活性以及与其他细胞周期调节因子的广泛相互作用表明，它在肺癌中作为不受控制的细胞分裂的主要协调者。因此，CDK1成为进一步功能验证和治疗探索的最有希望的候选者。选择LUAD作为验证对象是因为它是最常见的肺癌类型，占肺癌病例的很大比例。在不同独立平台上的多重表达分析表明，CDK1不仅在转录水平上上调，而且在蛋白质水平上也过度表达。此外，CDK1水平较高的患者生存期更短，复发更早，这是一个明显的临床效应。除了表达外，CDK1与免疫细胞浸润的相关性表明，其活性超出了表达范围，影响了肿瘤微环境，尤其是在T细胞亚群和巨噬细胞极化方面。这意味着CDK1可以调节肿瘤生长和免疫反应。最后，在CDK1中观察到的体细胞突变和拷贝数变化增加了这种蛋白质在LUAD患者中的致癌意义。从生物学角度来看，CDK1不仅仅是一个统计上优先的枢纽基因。作为G2/M转换和有丝分裂进入的关键调节器，CDK1在异常激活时促进持续的肿瘤细胞增殖。先前在肺癌中的研究也将CDK1表达升高与较差的生存率、增强的细胞周期和DNA修复信号以及LUAD中的改变的免疫相关通路联系起来。这些观察结果加强了CDK1可能作为肺癌进展的功能性重要驱动因素的观点，使其成为治疗靶向的合理候选者。在这项研究中，我们采用了基于结构的药物发现方法来高效筛选大量植物化学物质库。在这一领域，机器学习（ML）、分子对接和分子动力学（MD）模拟研究的集成改变了新型药物候选物的识别和优化。我们进行了基于化学信息学的筛选程序，筛选了9577种植物化合物的理化和药物相似性特性。在删除重复条目和不需要的子结构后，最终获得了2113种植物化合物，其中1802种化合物没有违反RO5和Veber规则。ML通过减少传统的实验时间和成本，彻底改变了计算药物发现方法。ML是在非常短的时间内高精度计算大型数据集生物活性的最佳选择。为了识别最有希望的候选物，我们应用了基于ML的pIC50预测模型，这帮助我们识别出具有更高生物活性潜力的化合物。这一步将数据集缩小到380种pIC50 > 6.5的植物化合物，以便进行进一步的计算分析，包括分子对接、分子动力学模拟和药代动力学分析。对接结果和对接后的MM-GBSA结合能（ΔG）突出了排名前三的植物化合物（CID_487089、CID_174880和CID_14218027）作为候选药物分子。此外，分子动力学（MD）模拟用于了解所选植物化合物随时间与CDK1的相互作用。虽然apo蛋白和AZD5438结合复合物显示出明显结构波动，但所有植物化合物结合系统在100 ns的模拟过程中都表现出更好的稳定性。值得注意的是，CID_174880始终优于对照化合物，保持较低的RMSD和RMSF值（图9），表明其结合模式更强且更稳定。结构紧凑性分析进一步突出了这一差异。CID_174880复合物显示出最低的SASA和紧密维持的Rg，反映了比AZD5438更深层次的配体埋藏和更紧凑的蛋白质结构。相比之下，对照复合物相对暴露于溶剂中且更灵活，表明CDK1的稳定性较弱。模拟后的MM-GBSA结果进一步支持了这些发现，CID_487089显示出最有利的结合自由能（?40.29 kcal mol?1），其次是CID_174880（?36.06 kcal mol?1），两者都优于对照AZD5438（?29.53 kcal mol?1）。使用PCA进行的动态运动分析显示了对照化合物（AZD5438）和植物化合物之间的明显差异，表明配体与结合位点之间的匹配更好。FEL分析进一步支持了这些发现，展示了配体-蛋白质相互作用的稳定性和构象状态。在研究的化合物中，CID_174880显示出与对照相似的动态和能量特性，表明其结合稳定性良好。对照化合物被预测具有致癌性，而没有任何植物化合物显示出这种风险，表明其总体安全性更好。然而，所有化合物都表现出一定的呼吸毒性，植物化合物还显示出免疫毒性，这可能限制了它们的开发。根据LD50值和毒性分类，CID_14218027是最安全的（第4类），而CID_487089的毒性更高（第2类）。CID_174880的毒性特性接近对照组，但安全性稍有提高。总体而言，CID_174880显示出比对照化合物更好的抑制潜力和更优的安全性。因此，这项研究可以在通过实验（体内和体外）研究验证后，成为识别天然CDK1抑制剂的有用资源。

5. 结论

本文指出CDK1是肺癌的重要驱动因素，它破坏了细胞周期、p53通路和免疫微环境，例如增强的M1巨噬细胞、减少的M2极化以及由于过表达和突变（包括I136N和E57V）导致的不良预后。这些发现表明肺癌存在共同的脆弱性，将有丝分裂失调作为潜在的治疗靶点之一。在评估的植物化学物质中，来自Cichorium intybus的lactupicrin（CID_174880）在基于计算机优先级的候选者中表现最为突出，其对接性能、结构稳定性、动态行为以及预测的药代动力学和毒性特性均优于参考化合物AZD5438。尽管在一些指标上观察到轻微波动，但所有CDK1复合物均保持了整体稳定性和兼容性。对接结果、分子动力学（MD）模拟、分子力场-梯度Brunnelman-Sandau (MM-GBSA) 计算以及 ADMET 分析一致认为 lactupicrin 是一个具有计算优先级的候选药物，值得进一步开发。为了确认其作为 CDK1 抑制剂的功能、抗癌效果及安全性，还需要进行进一步的实验验证，包括体外生化测定、基于细胞的功能研究以及体内实验。总体而言，这项研究为识别有效的 CDK1 靶向化合物提供了明确的工作流程，并为开发广谱、基于机制的肺癌治疗方法指明了方向，将计算预测与潜在的实际应用联系起来。

6. 本研究的局限性

尽管本研究采用了综合性的多组学和计算机模拟方法，但仍需考虑一些局限性。首先，转录组分析仅基于公开可用的微阵列数据集，这些数据集可能会受到批次效应、平台相关偏差以及有限临床信息的影响。虽然使用四个独立 cohort 可以提高可靠性，但缺乏 RNA-seq 数据可能会限制对基因表达细微变异性的捕捉。其次，药物发现流程（包括机器学习预测、分子对接、MM-GBSA 分析、MD 模拟和 AMDET）完全依赖于计算模型。尽管这些方法有助于筛选候选药物，但它们无法完全反映生物系统的复杂性、肿瘤的异质性或真实的药代动力学行为。最后，所发现的候选药物（尤其是 lactupicrin）尚未在体外或体内模型中进行验证。因此，这些候选药物的治疗潜力、最佳剂量及安全性特征仍需通过实验验证来确认。

作者贡献

概念构思 – Md. Ahad Ali 和 Hriddhi Sarker；方法论 – Md. Ahad Ali、Hriddhi Sarker 和 Humaira Sheikh；数据整理 – Hriddhi Sarker、Marguba Kamrun、Bilkis Akter Shifa 和 Sujoy Banik；形式分析 – Md. Ahad Ali、Hriddhi Sarker、Humaira Sheikh 和 Siam Ahmed；可视化 – Hriddhi Sarker、Md. Ahad Ali、Marguba Kamrun 和 Tarikul Islam；验证 – Md. Ahad Ali、Hriddhi Sarker 和 Tarikul Islam；项目管理 – Md. Ahad Ali；软件与资源 – Neeraj Kumar、Sujoy Banik 和 Bilkis Akter Shifa；监督 – Md. Ahad Ali；写作 – 初稿撰写：Hriddhi Sarker、Md. Ahad Ali 和 Humaira Sheikh；修订与编辑：Md. Ahad Ali、Marguba Kamrun 和 Neeraj Kumar。

利益冲突

所有作者声明不存在利益冲突。

数据可获取性

本研究中的原始数据和贡献包含在文章及其补充信息（SI）中（表 S1–S10 和图 S1–S5）。我们选定的化合物的训练集、测试集、chEMBL 数据集（包含共识预测）以及用于预测选定化合物生物活性的 QSAR 模型的相关 Python 代码，都可以在我们的 GitHub 仓库（https://github.com/ahad004/LUAD_ML_QSAR_Modeling）中找到，或者通过以下 DOI 访问：https://doi.org/10.5281/zenodo.19841200。补充信息可获取，详见 DOI：https://doi.org/10.1039/d6dd00045b。

参考文献

?

这些作者作为第一作者对本工作做出了同等贡献。