《Computational Biology and Chemistry》：In-Silico Transcriptomic Characterization Analysis Identifies SMAD5, USP2, and CUL4A as Key Biomarkers and Therapeutic Targets in Colorectal Cancer
整合转录组学与网络分析鉴定出SMAD5、USP2和CUL4A为候选的结直肠癌（colorectal cancer, CRC）生物标志物。针对约25,000种化合物进行的虚拟筛选（virtual screening）鉴定出靶向上述三个优先蛋白的有前景配体。在100纳秒（ns）分子动力学（molecular dynamics, MD）模拟过程中，CUL4A–橙皮苷（Hesperidin）与USP2–R428复合物表现出最持久的相互作用。吸收、分布、代谢、排泄及毒性（ADMET）与结构分析突出显示MK3207、R428和橙皮苷为需要进一步评估的独特先导候选化合物。该整合的转录组学–分子对接–分子动力学工作流程为优先筛选潜在的结直肠癌治疗靶点提供了一个计算框架。
一、研究背景与意义
结直肠癌（colorectal cancer, CRC）是全球最常见且致死率最高的癌症之一，对人群健康构成重大威胁。尽管近年来在筛查、诊断和治疗方面取得了进展，CRC仍具有高度的临床和生物学复杂性，表现为显著的分子异质性、异质性的疾病进展过程以及对治疗反应的差异。结直肠肿瘤的发病机制是一个多阶段过程，涉及遗传、表观遗传和转录组学改变的持续积累，推动正常结肠黏膜向腺瘤乃至浸润性癌的进展。这种腺瘤-癌序列代表了一个动态的、生物学上复杂的疾病谱系，其中初始分子缺陷可能影响后续的恶性进展。因此，发现有效的生物标志物和具有生物学重要性的治疗靶点对于增强CRC的早期诊断、预后分层和靶向干预仍然至关重要。
临床上，CRC是全球第三大常见确诊癌症和第二大癌症相关死亡原因，且50岁以下成年人的发病率正在上升。治疗策略依分期而定：局限性疾病通常采用内镜或手术切除治疗，而FOLFOX、FOLFIRI或CAPOX方案仍是晚期疾病全身治疗的基石。在转移性CRC中，生物标志物指导的治疗包括：针对RAS野生型左侧肿瘤的抗EGFR抗体；抗血管生成药物；BRAF V600E导向的encorafenib联合cetuximab；HER2靶向的基于trastuzumab的方案；KRAS G12C抑制剂联合抗EGFR治疗；以及针对错配修复缺陷/微卫星高度不稳定肿瘤的免疫检查点抑制剂。尽管取得了这些进展，大多数转移性CRC仍然难以治疗，特别是微卫星稳定、RAS突变型疾病，而获得性耐药很常见，五年生存率仍低于20%。这些挑战凸显了识别在独立CRC队列中持续失调、与已确立的疾病驱动机制相关联、且适用于小分子靶向的分子靶点的必要性。
转录组学分析已被证明是阐明CRC分子景观和识别哪些基因及通路促进肿瘤发展和进展的有用方法。来自人类结直肠组织的基因表达数据集提供了研究正常黏膜、腺瘤和癌等各种病理阶段转录变化的潜力。大型转录组数据在公共数据库（如美国国家生物技术信息中心基因表达综合数据库，NCBI GEO）中已变得更加可获取，使得二次分析和跨队列比较以发现生物标志物成为可能。然而，单一数据集的研究结果可能受到队列特异性变异、样本异质性和平台特异性偏倚的影响。在这方面，两个或更多独立数据集的整合分析越来越被认为是识别CRC中可重复且具有生物学意义的候选基因的更强大方法。
尽管差异表达分析在识别失调基因方面有效，但基因列表并不总能提供癌症候选生物标志物的生物学和机制性见解。整合性生物信息学方法可以通过将转录变化置于更广泛的功能背景中来克服这一缺点。蛋白质-蛋白质相互作用（protein-protein interaction, PPI）网络分析、模块检测和功能富集分析提供了互补的信息层，有助于识别功能相关的基因组、中心调控节点以及与CRC进展相关的富集生物学过程和信号通路。这些方法通过利用转录组可重复性、网络拓扑和通路水平解释，能够更有信心地优先筛选具有机制意义的候选基因。
已有众多CRC生物标志物被报道，包括癌胚抗原（CEA）和CA19-9、KRAS/NRAS/BRAF/PIK3CA改变、错配修复/微卫星不稳定状态以及ERBB2扩增，以及转录组学候选基因如CDH3、CEMIP（KIAA1199）、ETV4、FOXQ1、ASCL2、CD44和AXIN2。然而，许多转录水平候选基因缺乏独立的跨队列验证，而其治疗可靶向性也很少在同一发现框架内进行评估。为解决这些局限性，本研究应用了一个可重复性门控的多队列工作流程，在靶点选择前评估方向一致性，并整合了网络中心性、模块成员、通路相关性、基于结构的虚拟筛选、分子对接和计算机模拟ADMET评估。该方法优先筛选出SMAD5、USP2和CUL4A作为可重复支持的CRC相关靶点，并鉴定出MK3207、R428和橙皮苷作为候选配体。据研究人员所知，在统一的转录组-结构工作流程中对这些靶点进行优先筛选代表了本研究新颖性的一个独特方面。
二、研究内容与主要发现
研究人员通过整合多队列转录组学、网络生物学和基于结构的药物发现，旨在识别可重复的、潜在结构可靶向的CRC生物标志物和治疗靶点。研究分析了三个GEO数据集（GSE8671、GSE14333和GSE89076），进行差异表达、数据集一致性、重叠和聚类、PPI网络、模块检测和功能富集分析。随后使用虚拟筛选、分子对接、ADMET分析和分子动力学模拟对优先靶点进行了测试。这一整合工作流程提示SMAD5、USP2和CUL4A是可能与CRC相关的靶点，并鉴定出MK3207、R428和橙皮苷为先导分子。该研究的重要性在于它将可重复的基因失调与网络相关性和动态结构验证相关联，使得能够以更强的计算基础选择靶点和化合物用于后续实验。
三、关键技术方法
研究人员从NCBI GEO数据库获取了三个公开可用的CRC转录组数据集：GSE8671（包含32对配对正常黏膜和腺瘤样本，基于GPL570平台）、GSE14333（包含290例原发性CRC样本，基于GPL570平台）和GSE89076（包含配对正常和肿瘤结直肠组织，包括腺瘤样本，基于GPL16699平台）。在R环境中进行数据预处理，包括log2转换、探针到基因的注释、多探针基因保留最高表达探针以及去除未注释探针。使用limma包进行基因水平线性模型拟合和经验贝叶斯方差调节，以Benjamini-Hochberg错误发现率（false discovery rate, FDR）校正多重检验，以校正p值<0.05且|log2差异倍数（fold change, FC）|≥1为差异表达基因（differentially expressed genes, DEGs）的筛选标准。跨数据集一致性通过基因水平log2FC值的Pearson相关性分析评估。使用STRING数据库（综合评分≥0.4）构建PPI网络，并在Cytoscape 3.10.4中进行拓扑分析，使用MCODE插件进行模块检测（degree cutoff=2, node score cutoff=0.2, k-core=2, maximum depth=100）。功能富集分析使用clusterProfiler进行基因本体论（Gene Ontology, GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）通路富集。结构生物学方面，USP2（PDB ID: 5XU8, chain A, 1.81 ?）和CUL4A（PDB ID: 2HYE, chain C, 3.10 ?）的三维结构来自RCSB蛋白质数据银行，而全长人SMAD5结构来自AlphaFold蛋白质结构数据库（AF-Q99717-F1）。使用MTiOpenScreen网络服务器对约25,000种化合物进行基于结构的虚拟筛选，采用AutoDock Vina工作流程。分子对接使用AutoDock Vina 1.2.3进行，分子动力学模拟使用Schr?dinger套件中的Desmond引擎进行约100纳秒的全原子模拟，并使用Prime MM-GBSA计算结合能。ADMET分析使用SwissADME和pkCSM平台进行。
四、研究结果
差异表达分析：在GSE8671中，评估了20,150个探针水平条目，其中11,116个条目满足校正p<0.05的阈值，6,041个具有正logFC值，5,075个具有负logFC值，表明广泛的转录组重编程且略微偏向于上调。最显著上调的基因包括ADAMDEC1（logFC=3.78）、LPAR1（logFC=2.27）、PAG1（logFC=2.32）、ABCG2（logFC=4.43）和PRDX6（logFC=1.56）；最显著下调的基因包括AXIN2（logFC=?2.59）、CDH3（logFC=?6.37）、CD44（logFC=?2.14）、CEMIP（logFC=?6.01）和SORD（logFC=?1.98）。在GSE14333中，评估了20,845个基因，其中926个上调、945个下调、18,974个不显著，表明只有相对较小比例的转录组在该标准下显示差异表达。值得注意的是，CUL4A（logFC=?0.60，校正p=0.725）和SMAD5（logFC=0.28，校正p=0.925）在该数据集中未达到DEG筛选标准。在GSE89076中，包含22,698个基因，其中2,317个上调、1,456个下调、18,925个不显著。最上调的基因包括SULT1A2、SCGN、UGP2、HIGD1A和GLP2R；最下调的基因包括KRT80、KIAA1199、ETV4、FOXQ1和ASCL2。
跨数据集一致性分析：GSE8671与GSE14333的比较显示弱负相关（r=?0.11, n=16,655），表明两个数据集间基因表达差异的方向和幅度一致性差。GSE8671与GSE89076的相关系数高度正相关（r=0.70, n=16,655），表明两个数据集在转录组反应上具有强一致性。GSE14333与GSE89076的相关性不显著（r=?0.01, n=16,655）。总体而言，GSE8671和GSE89076表现出最高水平的跨数据集转录一致性，而GSE14333与两个比较数据集的一致性均弱至可忽略。
重叠分析：上调DEG的重叠模式显示，GSE8671与GSE89076之间的成对重叠最大，共有1,040个共同上调基因；GSE8671与GSE14333之间仅共享31个上调基因，GSE14333与GSE89076之间仅共享44个。仅37个基因在所有三个数据集中一致上调。下调DEG方面，GSE8671与GSE89076之间共享500个下调基因，GSE8671与GSE14333之间仅14个，GSE14333与GSE89076之间52个，仅16个基因在所有三个数据集中共同下调。每个数据集还包含大量独特基因，表明尽管GSE8671和GSE89076之间存在坚实的共享特征，但整体基因图景受数据特异性变异的影响很大。
无监督层次聚类：GSE8671的聚类显示UR和DR基因集之间的清晰互斥行为，支持腺瘤表型相关的有序转录重塑。GSE14333的聚类显示较不明显的分组，SMAD5和CKAP4等基因在特定样本亚群中表达更明显，而STARD5、APOA5、CUL4A、NUDT22和CDC27等基因表现出中间、样本限制的表达行为，被描述为分期变异型基因。GSE89076显示三个数据集中最强的双向表达模式，KRT80、ETV4、ENC1、ASCL2、FOXQ1、CDH3和KIAA1199等基因在单一主要样本块中聚类，而SULT1A2、GLP2R、SPIB、SCGN和HIGD1A等UR基因在相同样本块中显示一致降低的标准化表达，形成清晰的互斥条带。
网络与功能富集分析：弦图分析显示CRC生物标志物集形成结构良好的非随机分子网络，具有不同但相互连接的基因模块。红色扇区包含ADAMDEC1、LPAR1、PAG1、ABCG2和PRDX6等基因；蓝色扇区包含CDH3、CD44、CEMIP、SORD、GIMAP1、GAPT、SCIN、ERICH5、BEST4和LEAP2等基因。功能富集分析显示GO生物学过程（biological process, BP）以对刺激的反应、对应激的反应、多细胞生物过程、生物过程的正调控、细胞对刺激的反应、细胞群增殖、发育过程、解剖结构发育、运动和对外源化学物的反应为主。GO细胞组分（cellular component, CC）富集集中在细胞外空间、细胞外周、细胞外囊泡、细胞外膜结合细胞器、细胞外细胞器、细胞外外泌体、细胞外区域、细胞顶端部分、质膜和顶端质膜。GO分子功能（molecular function, MF）分析鉴定出蛋白质结合、小分子结合、离子结合、转运蛋白活性、金属离子跨膜转运蛋白活性、阴离子结合、跨膜转运蛋白活性和钠离子跨膜转运蛋白活性的富集。通路富集进一步突出了矿物质摄取、金属硫蛋白-金属相互作用、金属反应过程、基质金属蛋白酶激活和锌稳态。
靶点优先筛选：结合跨队列转录组可重复性、网络拓扑、功能富集和结构可靶向性，鉴定出三个最有前景的CRC靶点：SMAD5、USP2和CUL4A。这三个基因均位于可重复的转录核心中，主要来自GSE8671和GSE89076的比较，且均存在于整合表达分析中。网络分析显示这三个基因在STRING PPI网络中相互连接，并与Wnt/黏附相关信号的MCODE衍生模块相关联。结构评估显示三个蛋白质均具有适合基于结构配体评估的区域：USP2具有明确的催化/底物结合裂隙，CUL4A具有cullin重复相互作用区域，SMAD5具有功能重要的MH2相互作用结构域。
虚拟筛选：针对CUL4A、SMAD5和USP2的化合物库虚拟筛选鉴定出几个具有理想结合评分的高亲和力化合物。对于CUL4A，橙皮苷（Hesperidin_ZINC000103939005）具有最高的初始结合亲和力（?10.8 kcal/mol）和精炼对接评分（?8.7 kcal/mol）。对于SMAD5，Metergotamine（ZINC000072266819）在初始筛选中能量最佳（?10.0 kcal/mol），而精炼对接能量最佳的是MK3207（ZINC000103760984，?10.6 kcal/mol）。对于USP2，初始筛选最佳命中为Ravidasvir（ZINC000150607150，?12.5 kcal/mol），但精炼对接评分最佳的是R428（ZINC000051951668，?12.1 kcal/mol）。值得注意的是，USP2在所有靶点中表现出最佳的总体结合谱，具有最低的筛选和对接评分。
分子对接分析：MK3207与SMAD5的对接亲和力预测为?10.6 kcal/mol，通过GLN169和HIS367的常规氢键以及GLY365、PRO368、LEU160、GLU164、ASN159、ASN170和SER176的范德华接触稳定复合物。R428与USP2的对接能量为?12.1 kcal/mol，预测结合位点包括与SER5–N1（3.16 ?）、SER190–N22（3.24 ?）和GLU360–N14（3.45 ?）的三个氢键，以及涉及PHE364、TYR186和PHE505的芳香族稳定化相互作用。重要的是，R428的预测姿态位于经典USP2催化三联体（CYS276–HIS557–ASP574）之外，与CYS276的距离约12.9 ?，与HIS557的距离约11.6 ?，表明与非催化表面沟槽而非催化裂隙相关联。橙皮苷与CUL4A的对接亲和力为?8.7 kcal/mol，最短氢键接触为THR684（2.71 ?）、ARG687（2.80 ?）、GLN673（2.83 ?）、THR537（2.89 ?）、THR535（2.90 ?）和THR677（2.94 ?）。
分子动力学模拟：对SMAD5–MK3207、USP2–R428和CUL4A–橙皮苷复合物进行了约100 ns的全原子MD模拟。系统分别包含64,403、121,451和99,062个原子以及18,980、37,219和28,789个水分子。通过比较蛋白质/配体RMSD、Cα-RMSF、配体回旋半径（Rg）、溶剂可及表面积（SASA）、分子表面积（MolSA）、极性表面积（PSA）、二级结构变化、分子内氢键和蛋白质-配体接触占有率来检查结构稳定性和相互作用的持续性。CUL4A–橙皮苷和USP2–R428复合物在100 ns模拟过程中表现出最持久的相互作用。使用主成分分析（PCA）表征集体运动，并使用动态交叉相关矩阵（DCCM）分析评估相关和反相关残基运动。Prime MM-GBSA计算在0和100 ns结构上进行，作为结合能端点估计。
ADMET分析：使用SwissADME和pkCSM对候选化合物进行计算机模拟ADMET分析，评估了分子量、亲脂性、氢键供体和受体数量、拓扑极性表面积以及Lipinski五规则符合性等物理化学参数。预测的药代动力学参数包括胃肠吸收、血脑屏障通透性、P-糖蛋白底物状态和细胞色素P450相互作用潜力。毒理学终点如Ames致突变性、肝毒性和口服毒性相关参数也被分析。分析突出显示MK3207、R428和橙皮苷为需要进一步评估的独特先导候选化合物。
五、讨论与结论
该研究通过整合多队列转录组学、网络生物学和基于结构的药物发现，成功鉴定出SMAD5、USP2和CUL4A作为CRC的可重复、结构可靶向的生物标志物和治疗靶点，并确定了MK3207、R428和橙皮苷作为候选先导分子。该研究的意义在于将可重复的基因失调与网络相关性和动态结构验证相关联，使得能够以更强的计算基础选择靶点和化合物用于后续实验。这一整合的转录组学-对接-分子动力学工作流程为优先筛选潜在的CRC治疗靶点提供了一个计算框架。尽管这些发现需要在体外和体内模型中进一步验证，但该研究为CRC的精准医疗提供了新的候选靶点和先导化合物，并为整合计算生物学在癌症药物发现中的应用提供了范例。
