利用转录组分析鉴定乳腺癌中的枢纽基因、单核苷酸多态性及潜在药物靶点

摘要

本研究探讨了与HER2阳性乳腺癌治疗耐药性相关的线粒体氧化应激相关基因。通过转录组数据集、整合生物信息学分析和临床数据（n = 4,929），作者鉴定出MTHFD2和PRDX3为具有潜在预后价值的关键差异表达基因。in silico 分析提示遗传变异可能具有功能影响。

摘要

HER2阳性（HER2⁺）乳腺癌常对拉帕替尼（Lapatinib）等治疗产生耐药性，这一过程可能涉及线粒体代谢和氧化还原重编程。本研究旨在鉴定与线粒体氧化应激相关的差异表达基因（MOS-DEGs）作为耐药性的生物标志物，并表征其功能性非同义单核苷酸多态性（nsSNPs）及结构影响。通过DESeq2分析RNA-seq数据集（GSE231524、GSE231525）以鉴定差异表达基因（DEGs），并将这些基因与线粒体氧化应激相关基因取交集，获得MOS-DEGs。随后进行功能富集分析、蛋白质-蛋白质相互作用（PPI）分析、ROC分析、表达谱分析和生存分析。采用多种预测工具评估nsSNPs，并通过二级结构和三维建模分析其结构影响。整合转录组学与临床分析发现，MTHFD2和PRDX3是在线粒体氧化还原代谢中发挥相反调控作用的核心MOS-DEGs。MTHFD2显著上调，且与不良预后相关（HR = 1.53，p = 1.1 × 10⁻¹⁶），而PRDX3则表现出保护性表达模式（HR = 0.73，p = 7.7×10⁻¹⁰）。ROC分析提示二者具有预测治疗反应的潜力。nsSNP分析揭示了MTHFD2中的5个有害变异和PRDX3中的4个有害变异，其中rs1471336772（MTHFD2）和rs747786383（PRDX3）被鉴定为最具致病性的位点。这些变异在多个计算评分系统中均被预测为有害，包括SIFT ≤ 0.05、PolyPhen-2 ≥ 0.85、有害的CADD评分以及高REVEL评分，且位于催化域或氧化还原活性域内。结构建模表明，这些氨基酸替换可能破坏蛋白构象稳定性，干扰金属离子和辅因子结合位点，并影响NADPH再生或硫氧还蛋白依赖性过氧化物酶活性。分子动力学模拟预测这些变异可能导致结构稳定性下降和灵活性改变，提示其可能损害线粒体氧化还原调控功能。本研究确定MTHFD2和PRDX3为HER2⁺乳腺癌中线粒体氧化应激的关键调控因子。这些基因中的有害nsSNPs可能通过改变氧化还原平衡，进而影响代谢适应（MTHFD2）和抗氧化防御（PRDX3）过程。

引言

乳腺癌是全球女性中最常见的癌症，尽管在早期检测和靶向治疗方面已取得显著进展，但它仍是癌症相关死亡的主要原因之一^1,2,3。乳腺肿瘤内部的分子异质性导致了不同的临床结局和治疗反应，这对精准肿瘤学构成了持续的挑战^4,5。在已知的分子亚型中，人表皮生长因子受体2（HER2）阳性乳腺癌约占所有乳腺癌的15–20%，其特征是HER2基因扩增以及HER2受体酪氨酸激酶的过表达^6,7。尽管诸如曲妥珠单抗（trastuzumab）和拉帕替尼（lapatinib）等HER2靶向治疗已显著改善了患者的预后，但这些药物常出现原发性或获得性耐药，导致肿瘤复发和疾病进展^8,9。因此，阐明驱动HER2治疗耐药的分子机制，仍是临床肿瘤学中尚未满足的重大需求。

越来越多的证据表明，线粒体功能障碍和氧化应激在介导乳腺癌耐药性中发挥关键作用。线粒体除了在ATP生成中的经典作用外，还是细胞凋亡、氧化还原信号传导和代谢可塑性的关键调控者——这些过程均被癌细胞利用以在治疗压力下存活。特别是耐药性HER2⁺癌细胞表现出向氧化磷酸化（OXPHOS）的代谢转变、活性氧（ROS）缓冲能力增强以及线粒体生物发生改变。这些适应性变化使细胞能够在持续抑制HER2信号通路的情况下仍维持生存信号并逃避凋亡。因此，编码线粒体复合物亚基和氧化还原调控酶的基因可能是克服耐药性的潜在生物标志物和治疗靶点。

同时，线粒体或氧化应激相关基因中的非同义单核苷酸多态性（nsSNP）可改变蛋白质的结构与功能，从而影响酶活性、药物反应性及疾病易感性¹⁰。利用SIFT、PolyPhen-2和CADD等in silico工具对nsSNP的有害效应进行计算预测，可快速识别可能参与肿瘤异质性和治疗耐药的功能性变异。因此，整合转录组与突变数据可从转录和基因组两个层面，提供基因失调与结构改变的多维视角。

生物信息学和系统生物学的最新进展使得通过大规模表达谱分析和基于网络的分析方法高通量地鉴定潜在驱动基因成为可能。诸如基因表达综合数据库（Gene Expression Omnibus, GEO）等公共存储库提供了宝贵的RNA测序数据集，能够捕捉不同实验模型和患者队列中的分子变化。结合DESeq2、clusterProfiler以及STRING–Cytoscape网络分析等分析工具，这些资源可实现差异表达基因（DEGs）的全面表征、通路富集分析以及枢纽基因的发现。尤为重要的是，将线粒体氧化应激基因特征与差异表达基因整合分析，有助于揭示HER2靶向治疗耐药性的新分子机制，并鉴定可用于干预治疗的候选生物标志物。

HER3驱动的耐药轴最近受到关注，被认为是HER2抑制后的一条重要代偿通路。HER3的活化可恢复下游PI3K/AKT信号传导，促进细胞存活和药物耐受性（图1）。DUSP6（双特异性磷酸酶6）是ERK信号通路的负调控因子，已被发现参与调控这种适应性反应。抑制DUSP6可重新激活ERK信号通路，可能逆转HER3介导的耐药性，但其与线粒体氧化应激及代谢适应之间的相互作用仍缺乏充分理解。因此，在DUSP6抑制和长期拉帕替尼暴露条件下，对HER2⁺细胞系（BT474和MDA-MB-453）进行比较转录组学分析，为解析药物耐药性的分子决定因素提供了一个理想模型。

图1。线粒体功能障碍、氧化应激与乳腺癌进展之间分子机制的示意图。该图展示了乳腺癌中由线粒体氧化应激介导的肿瘤发生所涉及的整合性分子机制。线粒体电子传递链复合物的失调，特别是涉及 DUFS3、UQCRC1、COX4I1、SDHA 和 ATP5PO 的复合物，导致活性氧（ROS）过度生成，进而引起氧化性DNA损伤、脂质过氧化以及线粒体膜电位受损。这些氧化性紊乱可激活PI3K/AKT、MAPK和NF-κB等促癌信号通路，同时抑制p53和BAX等凋亡调控因子，从而促进肿瘤细胞的存活、增殖及免疫逃逸。该示意图还强调了线粒体氧化应激对癌症代谢和肿瘤微环境的双重影响，突出了其在HER2+乳腺癌条件下对治疗耐药性和免疫调节的作用。这一整合模型为后续基于转录组和单核苷酸多态性（SNP）的分析提供了机制基础，旨在识别线粒体驱动基因及潜在的治疗靶点。 请点击此处查看该图的放大版本。

本研究旨在系统性地鉴定与HER2⁺乳腺癌中线粒体氧化应激相关的枢纽基因、功能性非同义单核苷酸多态性（nsSNP）以及潜在的药物靶点。通过整合来自GSE231524和GSE231525的转录组数据与经过人工整理的线粒体及氧化应激相关基因集，本研究旨在确定与线粒体氧化应激相关的差异表达基因（MOS-DEGs），这些基因可能参与治疗耐药性的形成。后续分析包括基因本体（GO）和KEGG通路富集分析、蛋白质-蛋白质相互作用（PPI）网络构建、枢纽基因优先排序，以及利用Kaplan–Meier免疫治疗数据库进行生存相关性分析。此外，还进一步探究了关键基因的nsSNP变异及其潜在的结构影响，以阐明由突变驱动的功能性后果。

通过整合转录组学分析、网络生物学和in silico突变分析，本研究为发现HER2⁺乳腺癌中潜在的线粒体生物标志物和治疗靶点提供了一个综合框架。氧化磷酸化相关枢纽基因及其功能性单核苷酸多态性的鉴定，可能为开发精准治疗策略铺平道路，以克服HER2靶向药物耐药性并改善患者预后。

方案

数据获取与预处理

本研究使用公开的转录组学和遗传学数据进行，未直接涉及人类或动物受试者。与HER2⁺乳腺癌治疗耐药性相关的转录组数据集从NCBI基因表达综合数据库（GEO）（https://www.ncbi.nlm.nih.gov/geo/）¹¹中获取。选择了两个RNA-seq数据集GSE231524和GSE231525，因其特别关注HER3驱动的耐药性以及在HER2+乳腺癌细胞系（BT474和MDA-MB-453）中DUSP6的抑制作用。这些数据集包括经拉帕替尼（1 µM）处理和DUSP6敲低后获得的亲本型、药物耐受型和药物耐药型表型。使用RStudio（v4.3.2）中的GEOquery（v2.70.0）和Biobase（v2.62.0）软件包访问原始计数矩阵及相应的元数据文件¹²。对元数据进行整理，为每个数据集定义两个主要对比：GSE231524将对照组（BT474亲本型，第0天）与药物耐受和药物耐药样本（第9天至第9个月）进行比较，而GSE231525将对照组（随机siRNA）与DUSP6敲低组（DUSP6-KD）进行比较。使用DESeq2（v1.42.0）框架进行质量控制和数据标准化，该框架采用方差稳定变换（VST）以减少异方差性，并确保样本间的可比性。使用ggplot2（v3.5.0）和pheatmap（v1.0.12）对数据分布和聚类模式进行可视化评估，以确认数据的一致性，并在差异表达分析之前识别潜在的离群值^13,14。

本研究中，严格区分了来自细胞系转录组数据集的发现与来自患者来源临床数据集的发现。细胞系数据主要用于探索性分析，包括在受控实验模型中鉴定差异表达基因以及获得初步的机制性见解。相比之下，患者来源的数据集则用于验证基因表达模式的外部一致性，并评估其临床相关性，包括预后评估。因此，细胞系模型的结果与临床队列的结果被分别解读，以避免过度泛化，并确保所有发现均具有恰当的转化研究背景。

差异基因表达分析

通过差异表达分析，鉴定在对照组与处理组之间显著变化的基因。使用整合在DESeq2中的校正回归模型（ARM）对标准化计数进行处理，以准确估计log₂倍数变化及统计显著性。实验设计公式定义为 ~condition，代表对照组与处理组之间的比较。筛选标准为：校正p值（FDR）< 0.05，且绝对log₂倍数变化 ≥ 1的基因被视为显著差异表达基因。采用apeglm方法对log₂倍数变化进行压缩调整，以提高效应量估计的稳健性。分析结果通过EnhancedVolcano（v1.22.0）¹⁵和ggplot2¹⁶进行可视化，生成火山图和MA图，展示表达幅度与统计置信度之间的关系。同时在DESeq2中评估了离散估计值，以确保在生物学重复样本间实现准确的方差建模和一致的标准化¹⁷。

线粒体氧化应激相关差异表达基因（MOS-DEGs）的获取与鉴定

为了研究能量代谢、氧化应激与药物耐药性之间的关联，从多个数据库中整合编制了一份线粒体及氧化应激相关基因的综合列表，包括 Human MitoCarta3.0¹⁸（https://personal.broadinstitute.org/scalvo/MitoCarta3.0/human.mitocarta3.0.html）、基因本体数据库（GO:0006979，氧化应激响应）（http://geneontology.org/）、京都基因与基因组百科全书（KEGG）氧化磷酸化通路（https://www.genome.jp/kegg/）以及人类氧化应激基因数据库（HOSGDB）（http://hosgdb.com/）。所有获取的基因均使用 org.Hs.eg.db（v3.18.0）和 AnnotationDbi（v1.64.0）统一转换为 HGNC 批准的基因符号，并剔除重复条目、假基因和非编码 RNA，以确保注释的准确性。最终构建的经人工审编的线粒体氧化应激基因集（MOS 基因）随后被用作参考基因集，与两个转录组数据集中鉴定出的差异表达基因进行整合分析。

使用 R 语言中的 dplyr（v1.1.3）¹⁹ 和基础 R 的 intersect() 函数，对精选的 MOS 基因列表与从 GSE231524 和 GSE231525 获得的差异表达基因（DEGs）进行交集分析。该整合性方法用于识别与线粒体代谢、氧化还原调控及氧化应激适应功能相关的 MOS-DEGs。通过 R 中的 VennDiagram（v1.7.3）包对数据集之间的重叠情况进行可视化，以展示不同实验模型间共有和特有的基因²⁰。获得的精炼 MOS-DEGs 基因列表用于后续分析，揭示 HER2 靶向治疗耐药背后的转录与代谢重编程机制。

MOS-DEGs 的表达谱分析与可视化

使用 RStudio 中的 ComplexHeatmap（v2.18.0）²¹ 和 pheatmap（v1.0.12）软件包对鉴定出的 MOS-DEGs 进行表达谱分析，以可视化亲本、药物耐受及耐药条件下全局表达模式。对标准化计数数据采用 z-score 标准化方法，以使不同样本间的基因表达矩阵标准化。采用欧氏距离和完全连接法进行聚类分析，以识别共表达模式并区分不同条件特异的转录谱。使用 ggplot2 生成热图和聚类图，确保不同条件之间具有清晰的视觉区分。该可视化方法有助于识别与线粒体活性、氧化应激调控以及耐药状态下代谢重编程相关的基因群组。

功能富集与通路注释

为探究鉴定出的MOS-DEGs的生物学意义及其调控机制，使用R Studio（版本4.3.1）进行了基因本体（Gene Ontology, GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）富集分析。分析在tidyverse环境中进行，利用多个Bioconductor软件包实现可重复的计算与可视化。基因注释和标识符映射基于org.Hs.eg.db数据库（https://bioconductor.org/packages/org.Hs.eg.db/）完成 Homo sapiens 参考基因组（GRCh38）。使用 clusterProfiler 软件包（版本 4.8.1；https://bioconductor.org/packages/clusterProfiler/）进行基因本体（GO）富集分析，该软件将基因分为三大本体类别——生物过程（Biological Process, BP）、细胞组分（Cellular Component, CC）和分子功能（Molecular Function, MF）。enrichGO²² 功能被用于参数设置为 p值 < 0.05 和 校正 p 值 （FDR） < 0.05，采用Benjamini–Hochberg校正方法。可视化图表（包括柱状图、点图和弦图）使用enrichplot（https://bioconductor.org/packages/enrichplot/）和ggplot2生成。¹³ （https://cran.r-project.org/web/packages/ggplot2/）和 GOplot（https://cran.r-project.org/web/packages/GOplot/）。这些工具提供了对富集的 GO 条目及其基因关联的结构化视图。

使用 clusterProfiler 软件包中的 enrichKEGG() 函数进行 KEGG 通路富集分析，并参考人类 KEGG 数据库（https://www.genome.jp/kegg/）。采用 KEGGREST 软件包（https://bioconductor.org/packages/KEGGREST/）进行通路数据的获取与注释。调整后的 p 值（q 值）< 0.05 的通路被视为具有显著性。通路的可视化与映射分析通过 pathview（https://bioconductor.org/packages/pathview/）、ggplot2 和 enrichplot 完成，而 igraph 与 ggraph 用于构建网络图谱²³。所有富集分析与可视化均在 R Studio（v4.3.1）中通过可重复的代码和标准化的 Bioconductor 工作流程实现，确保可靠地识别与 MOS-DEGs 相关的富集功能类别及生物学通路。

基于ROC的乳腺癌预测性生物标志物验证

为了验证MOS-DEGs的临床预测能力，使用ROCplotter在线工具（https://www.rocplot.org/）²⁴进行了受试者工作特征（ROC）曲线分析。ROCplotter是一个集成的基于网络的平台，整合了基因表达数据与来自3,104例乳腺癌患者的临床注释治疗反应数据集，其中包括接受化疗、激素治疗或抗HER2药物治疗的患者。

分析采用“病理学完全缓解”作为结果变量，“任何化疗”作为治疗类别。Affymetrix 芯片数据集得出的基因表达值根据平台内的临床注释自动分为应答组和非应答组。

采用受试者工作特征（ROC）曲线下面积（AUC）、Mann-Whitney U检验、倍数变化和卡方检验来评估每个基因区分应答者与非应答者的能力。曲线下面积（AUC）被用作评估区分性能的主要指标。AUC值大于0.55且ROC p值<0.05被视为具有显著性，代表转录组学生物标志物典型的中等预测性能，同时应用错误发现率（FDR）校正以保持分析的严谨性。

使用相应的Affymetrix探针ID对所有选定的MOS-DEG进行查询。通过临床乳腺癌队列评估每个基因的判别潜力，其中表达数据被分为应答者和非应答者两组。ROC曲线、箱形图及相关统计结果由ROCplotter平台直接生成，并导出用于后续可视化和比较。该分析量化了参与线粒体氧化应激的氧化还原与代谢调控因子的预测价值。在临床样本中表现出一致预测意义的基因被保留，用于最终预测面板的构建。

肿瘤、正常及转移组织中的差异表达分析（TNMplot 分析）

使用 TNMplot 网络工具 v2（https://tnmplot.com/analysis/）²⁵，分析了前导 MOS-DEG 在正常、肿瘤及转移性乳腺组织中的表达模式。通过 RNA-Seq（TCGA + GTEx + MET500）和基因芯片数据集进行分析，以确保跨平台验证。采用“多基因分析”模块，并选择浸润性乳腺癌作为目标组织类型²⁶。对肿瘤 vs. 正常（TvsN）、转移 vs. 肿瘤（MvsT）以及转移 vs. 正常（MvsN）三组之间的表达值进行 log₂ 转换后比较。TNMplot 使用 Mann–Whitney U 检验自动计算倍数变化（FC）和 p 值，以评估统计学显著性。表达分布以箱线图和密度图形式可视化，图形直接由 TNMplot 界面生成，其中绿色、红色和灰色分别代表正常组织、肿瘤组织和转移组织。所有图像均以高分辨率导出，用于整合至结果部分。该双平台分析有助于稳健地识别和验证与乳腺癌进展相关的关键线粒体氧化还原-代谢调控因子²⁷。

使用Kaplan–Meier绘图仪进行生存与预后分析

为评估MOS-DEGs在乳腺癌中的预后相关性，使用Kaplan–Meier Plotter在线工具（https://kmplot.com/analysis/）进行生存分析²⁸该数据库整合了来自超过4,900名乳腺癌患者的基因表达与生存数据，这些数据源自多个GEO、EGA和TCGA数据集。分析采用与优先基因对应的个体Affymetrix探针ID进行无复发生存期（RFS）分析：225609_at（GSR）、201761_at（MTHFD2）、201619_at（PRDX3/AOP1）和201128_s_at（ACLY）。患者被分为 高- 根据中位表达值将样本分为高表达组和低表达组，采用Kaplan-Meier法估计生存概率。使用对数秩检验（log-rank test）评估生存曲线之间的统计学显著性，风险比（HR）及其95%置信区间（CI）由工具自动计算。所有分析均以无复发生存期（RFS）为终点，未对激素受体或HER2状态（ER、PR、HER2 = 全部）进行限制。去除重复样本，并验证比例风险假设以确保统计稳健性。质量控制过滤排除了存在偏倚的微阵列数据。根据KM Plotter的默认设置，未进行人工探针筛选或对多重检验进行p值校正。统计学显著性定义为 p. < 0.05。通过高分辨率生存曲线图对每位候选MOS基因的高表达组与低表达组之间的预后差异进行可视化分析，并下载图像以供进一步解读^29,30.

本研究分析最初采用完全由HER2⁺乳腺癌样本组成的数据库，用于鉴定差异表达基因（DEGs）和核心基因（hub genes）。随后，生存分析未限制于HER2状态（ER、PR、HER2 = 全部），以评估所鉴定基因更广泛的预后相关性与普适性。该策略作为二次验证步骤而采用，而非重新定义研究重点。因此，对所鉴定核心基因的预后意义进行谨慎解读，主要结论仍特异性地针对HER2⁺乳腺癌。

MTHFD2-201（ENST00000394053.7）和 PRDX3-201（ENST00000298510.4）的典型转录本序列来自 Ensembl 基因组浏览器（https://www.ensembl.org）^31,32。使用 Ensembl 变异效应预测工具（VEP）（https://www.ensembl.org/vep）进行变异注释与分类，该工具为每个已识别的变异提供了详细的基因组背景、密码子改变及氨基酸替换信息。下游分析仅选择错义变异（非同义单核苷酸多态性）。

致病性预测与变异位点优先排序

每个非同义单核苷酸多态性（nsSNP）的功能后果通过多种计算预测工具联合评估。使用SIFT（https://sift.bii.a-star.edu.sg）评估氨基酸保守性，将得分≤0.05的变异归类为有害³³。PolyPhen-2（http://genetics.bwh.harvard.edu/pph2）用于估计氨基酸替换对蛋白质结构和进化的影响，得分≥0.85提示可能具有破坏性³⁴。CADD（https://cadd.gs.washington.edu）提供整合多种注释信息的综合有害性评分，评分≥20表示具有较高的致病潜力³⁵。此外，通过VEP分析界面整合了MetaLR³⁶、Mutation Assessor和REVEL等互补性预测指标，以提高预测的可靠性³⁷。满足以下阈值的变异——MetaLR≥0.70、Mutation Assessor≥3.5、REVEL≥0.75——被优先认定为可能致病。

结构与机制影响预测

为了评估氨基酸替换对结构完整性和生化功能的影响，使用 MutPred2（http://mutpred.mutdb.org）³⁸ 和 DynaMut（http://biosig.unimelb.edu.au/dynamut）³⁹ 对每个排名靠前的非同义单核苷酸多态性（nsSNP）进行进一步分析。MutPred2 估算了功能破坏的概率，包括催化活性改变、金属结合残基的获得或丧失、溶剂可及性的变化以及别构调节，评分 ≥ 0.80 被归类为高度致病性。DynaMut 计算了野生型与突变蛋白之间的吉布斯自由能变化（ΔΔG），评估了稳定性改变的方向和程度，并生成了原子位移和氢键重排的可视化图谱。

二级与三级结构建模及溶剂可及性分析

MTHFD2 和 PRDX3 的实验解析晶体结构取自蛋白质数据库（Protein Data Bank, PDB），并使用 PyMOL V:3.1（https://pymol.org）⁴⁰ 进行处理，以可视化有害残基的空间分布。通过引入相应的氨基酸替换构建突变体模型，随后进行结构优化和能量最小化。三维结构的比较分析揭示了二级结构元件的位移、原子间相互作用的改变，以及非同义单核苷酸多态性（nsSNP）与催化位点及辅因子结合结构域之间的空间邻近性，提示其可能破坏氧化还原和代谢功能。

使用 PSIPRED V: 3.2（http://bioinf.cs.ucl.ac.uk/psipred）^41,42 和 NetSurfP 3.0（https://services.healthtech.dtu.dk/service.php?NetSurfP-2.0）⁴³ 进行了二级结构和溶剂可及性分析。这些工具可预测α-螺旋、β-折叠、卷曲以及无序区域，并提供相对溶剂可及性（RSA）评分。通过映射具有中等到高 RSA 值且结构有序的残基，以识别溶剂暴露且功能关键的位点。受影响的位点在二维拓扑图中进行可视化，以判断有害突变是否发生在刚性的催化核心或柔性的环状区域，从而预测其对蛋白质折叠动力学和酶催化效率的可能影响。

数据库交叉验证、功能整合与稳定性验证

每个优先筛选的非同义单核苷酸多态性（nsSNP）均与群体水平的基因组数据库（包括dbSNP、1000基因组计划、ExAC和gnomAD）进行交叉比对，以确认其变异频率、全球等位基因分布以及先前报道的临床关联。通过整合进化保守性分析、结构建模和基于机器学习的功能预测，鉴定了MTHFD2和PRDX3中高置信度的有害变异。这些高影响突变随后被定位至功能结构域，以阐明其在乳腺癌线粒体氧化应激失衡、代谢信号通路改变及治疗耐药性中的潜在作用。为进一步验证每个有害替换的热力学效应，使用iMutant 3.0（https://folding.biofold.org/i-mutant/i-mutant3.0.html）⁴⁴基于序列和结构数据预测突变对蛋白质稳定性的影响。该分析计算了ΔΔG值（kcal/mol），代表野生型与突变型蛋白质之间自由能的变化。ΔΔG值为负的变异被归类为破坏性突变，表明蛋白质稳定性降低且解折叠概率增加。将iMutant的预测结果与DynaMut和MutPred2结果整合，实现了对可能影响氧化还原功能、催化完整性及整体蛋白质构象稳定性的关键结构残基的交叉验证。

整合的功能解释与治疗相关性

所有鉴定出的有害非同义单核苷酸多态性（nsSNP）均通过与dbSNP、gnomAD和ExAC群体数据库的交叉比对进行验证，以确认其次要等位基因频率及先前报道的与癌症表型的关联。综合进化保守性分析、结构建模和稳定性数据的解读表明，高影响突变rs1471336772（MTHFD2）和rs747786383（PRDX3）对蛋白质构象和催化效率具有最强的有害效应。计算分析结果共同提示，MTHFD2基因的突变会破坏NADPH依赖性氧化还原代谢，而PRDX3基因的突变则损害过氧化物酶介导的氧化应激防御功能，从而导致线粒体功能障碍和肿瘤侵袭性增强。本基于nsSNP的结构与功能分析为未来的治疗靶点筛选和突变验证提供了计算基础，并强调MTHFD2和PRDX3可作为针对氧化还原通路的乳腺癌精准治疗的生物标志物。为提高清晰度并全面概述分析策略，图2展示了总结本研究主要步骤的示意图工作流程。该流程整合了差异基因表达分析、线粒体基因筛选、蛋白质-蛋白质相互作用网络构建、基于ROC分析的临床验证以及基于nsSNP的结构表征。这一逐步推进的分析框架凸显了从转录组数据处理到生物标志物识别及功能解读的逻辑递进过程。

图2. HER2⁺乳腺癌中线粒体氧化应激相关生物标志物识别与验证的整合多步工作流程。该示意图总结了本研究采用的分析流程。首先，对RNA-seq数据集（GSE231524和GSE231525）进行差异基因表达（DEG）分析，以鉴定显著改变的基因。将这些差异表达基因与人工整理的线粒体氧化应激相关基因取交集，获得MOS-DEGs。接着，利用STRING和Cytoscape进行蛋白质-蛋白质相互作用（PPI）网络分析，以识别枢纽基因和功能模块。随后，使用ROCplotter平台进行受试者工作特征（ROC）曲线分析，评估所选基因在临床队列中的预测性能。最后，对优先基因（MTHFD2和PRDX3）中的非同义SNP（nsSNP）进行分析并开展结构建模，以评估关键变异可能产生的功能和结构影响。该整合性工作流程结合了转录组学、网络分析、临床数据和结构分析，用于识别潜在的生物标志物和治疗靶点。请点击此处查看该图的放大版本。

结果

差异基因表达分析

为研究与HER2⁺乳腺癌进展及治疗耐药性相关的转录组改变，采用两个独立的RNA-seq数据集GSE231524和GSE231525进行了差异表达谱分析。在第一个数据集（GSE231524）中，最初共定量了19,727个基因。经过标准化、过滤以及应用校正回归模型（ARM）进行差异表达分析后，保留了6,604个显著表达的基因（补充表S1）。在完成基因符号标准化、去除重复项及注释优化后，最终确定6,300个唯一基因用于后续分析。火山图、MA图和离散度估计结果（图3A-C）清晰展示了显著上调和下调基因的分布情况（校正后p值< 0.05，|log₂FC| ≥ 1）。离散度估计显示，在平均标准化计数范围内，方差模式拟合良好，表明DESeq2分析流程中的标准化效果可靠且模型拟合恰当⁴⁵。

同样，对于第二个数据集（GSE23...

讨论

线粒体氧化应激是肿瘤适应性的显著特征，调控细胞存活、代谢可塑性及治疗耐药性。本研究通过整合转录组数据、功能富集分析、临床验证以及详细的in silico突变和结构分析，系统探讨了HER2⁺乳腺癌中与线粒体氧化应激相关的差异表达基因（MOS-DEGs）。差异表达与突变谱的整合分析，从多维度揭示了氧化还原-代谢调控与基因组不稳定性之间的关联。在所鉴定的MOS-DEGs中，MTHFD2和PRDX3作为关键的线粒体决定因子脱颖而出，它们在线粒体氧化应激适应、代谢调控及治疗反应中表现出独特但互补的功能⁴⁸。在这些基因中发现的多个有害非同义单核苷酸多态性（nsSNPs）进一步揭示了单个氨基酸替换如何影响酶活性、蛋白质折叠以及乳腺癌细胞的线粒体韧性。需要特别强调的是，本研究中所呈现的结果来源于整合性的计算和in silico分析，因此仍需进一步的实验验证以确认其生物学和临床相关性。

MTHFD2 作为一种代谢性癌基因及预后驱动因子

材料

本文使用的材料清单

姓名	公司	目录编号	评论
AnnotationDbi	Bioconductor	v1.64.0	用于 HGNC 基因符号标准化与注释的 R 注释包。
apeglm	Bioconductor	DESeq2 中的方法	用于差异表达分析中 log2 倍数变化收缩的方法。
Biobase	Bioconductor	v2.62.0	用于访问和管理 GEO 衍生的表达数据与元数据对象的 R 包。
CADD	华盛顿大学 / Kircher 实验室	网络工具	用于预测非同义单核苷酸多态性（nsSNP）致病性的综合注释依赖性耗竭评分。
clusterProfiler	Bioconductor	v4.8.1	用于 GO 和 KEGG 富集分析的 R 包。
ComplexHeatmap	Bioconductor	v2.18.0	用于生成热图及表达聚类的 R 包。
Cytoscape	Cytoscape 联盟	软件	用于可视化由 STRING 衍生的蛋白质-蛋白质相互作用网络并结合 CytoHubba 进行关键基因优先排序的网络可视化平台。
dbSNP	NCBI	数据库	用于交叉验证优先级 nsSNP 的群体变异数据库。
DESeq2	Bioconductor	v1.42.0	用于标准化、方差建模和差异表达分析的 R 包。
dplyr	CRAN / tidyverse	v1.1.3	用于数据操作及基因集交集分析的 R 包。
DynaMut	墨尔本大学 / BioSig	网络服务器	用于估算蛋白质中突变相关稳定性与柔韧性变化的工具。
EnhancedVolcano	Bioconductor	v1.22.0	用于将差异表达结果以火山图形式可视化的 R 包。
Ensembl 基因组浏览器	EMBL-EBI / Ensembl	数据库	MTHFD2 和 PRDX3 的经典转录本序列来源。
Ensembl 变异效应预测器（VEP）	EMBL-EBI / Ensembl	网络工具	用于注释错义变异并整合预测评分的工具。
ExAC	布罗德研究所	数据库	用于 nsSNP 交叉验证的群体水平外显子组数据库。
基因表达综合数据库（Gene Expression Omnibus, GEO）	NCBI	数据库	用于获取 HER2 阳性乳腺癌分析的转录组数据集的存储库。
基因本体（Gene Ontology）	基因本体联盟	GO:0006979	用于获取氧化应激相关基因并进行富集分析的本体资源。
GEOquery	Bioconductor	v2.70.0	用于访问 GEO 计数矩阵和元数据的 R 包。
ggplot2	CRAN / tidyverse	v3.5.0	用于数据可视化、聚类图及图形输出的 R 包。
ggraph	CRAN	R 包	用于富集分析和通路表示过程中的图结构与网络可视化的 R 包。
gnomAD	布罗德研究所	数据库	用于检查优先级 nsSNP 频率与分布的群体变异数据库。
GOplot	CRAN	R 包	用于 GO 富集可视化（包括弦图和汇总图）的 R 包。
GSE231524	NCBI GEO	登录号	用于分析亲本、药物耐受及耐药 BT474 表型的 RNA-seq 数据集。
GSE231525	NCBI GEO	登录号	用于分析 DUSP6 沉默在 HER2 阳性乳腺癌细胞中作用的 RNA-seq 数据集。
人类线粒体基因集 MitoCarta3.0	布罗德研究所	数据库	用于定义线粒体基因集的经人工整理的线粒体基因资源。
人类氧化应激基因数据库（HOSGDB）	HOSGDB	数据库	用于获取氧化应激相关基因的数据库。
igraph	CRAN	R 包	用于网络表示和通路可视化的 R 包。
iMutant 3.0	博洛尼亚大学 / Biofold	网络服务器	用于预测突变对蛋白质稳定性影响的工具。
Kaplan–Meier Plotter	KMplot	网络工具	用于乳腺癌队列中无复发生存分析的在线平台。
KEGG	京都大学	数据库	用于氧化磷酸化及通路富集分析的通路数据库。
KEGGREST	Bioconductor	R 包	用于获取和注释 KEGG 通路信息的包。
Lapatinib	文中未指明	1 μM 处理条件	在原始实验数据集中用于诱导药物耐受/耐药表型的 HER2 靶向抑制剂。
MetaLR	通过 Ensembl VEP 集成	评分	用于变异优先排序的计算型致病性指标。
MutPred2	MutPred	网络服务器	用于预测氨基酸替换功能后果的工具。
NetSurfP 3.0	丹麦技术大学	网络服务器	用于二级结构和溶剂可及性预测的工具。
org.Hs.eg.db	Bioconductor	v3.18.0	用于基因标识符映射的人类基因组注释包。
pathview	Bioconductor	R 包	用于将富集基因映射到 KEGG 通路的包。
pheatmap	CRAN	v1.0.12	用于热图绘制和聚类可视化的 R 包。
PolyPhen-2	哈佛大学 / 布莱根妇女医院	网络工具	用于估算氨基酸替换对结构/功能影响的计算预测工具。
蛋白质数据库（Protein Data Bank, PDB）	RCSB PDB	数据库	MTHFD2 和 PRDX3 的实验解析蛋白质结构来源。
PSIPRED	伦敦大学学院	v3.2	用于二维拓扑分析的蛋白质二级结构预测服务器。
PyMOL	Schrödinger / PyMOL	v3.1	用于可视化蛋白质结构和突变位点的分子图形软件。
REVEL	通过 Ensembl VEP 集成	评分	用于错义变异优先排序的集成致病性评分。
ROCplotter	ROCplot.org	网络工具	用于基于 ROC 曲线验证乳腺癌应答队列中基因表达的在线平台。
RStudio	Posit	v4.3.1/v4.3.2 环境	用于转录组分析、富集分析和可视化工作流的统计计算环境。
SIFT	A*STAR	网络工具	用于将氨基酸替换分类为可耐受或有害的预测工具。
STRING	STRING 联盟	数据库	用于网络分析和枢纽基因识别的蛋白质–蛋白质相互作用数据库。
TMNplot	TMNplot.com	网络工具	用于比较正常、肿瘤和转移性乳腺组织中基因表达的平台。
VennDiagram	CRAN	v1.7.3	用于可视化差异表达基因与人工整理线粒体基因集之间重叠的 R 包。