Integrated Multi-Omics Profiling Maps Ferroptosis–Cuproptosis Diversity in Cervical Cancer and Identifies a PDGFRB-Driven Monocyte Fibrotic Program Targeted by Sorafenib

Integrated Multi-Omics Profiling Maps Ferroptosis–Cuproptosis Diversity in Cervical Cancer and Identifies a PDGFRB-Driven Monocyte Fibrotic Program Targeted by Sorafenib | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Integrated Multi-Omics Profiling Maps Ferroptosis–Cuproptosis Diversity in Cervical Cancer and Identifies a PDGFRB-Driven Monocyte Fibrotic Program Targeted by Sorafenib Shan Li, Mingyuan Wu, Yuanyuan Meng, Caixia Liang, Xiufang Huang, and 4 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-8811979/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background: Ferroptosis and cuproptosis are two distinct forms of metal-dependent regulated cell death that have emerged as important mechanisms in tumor biology. However, the crosstalk between these pathways and their clinical relevance in cervical cancer remain largely unexplored. Methods: Using the TCGA-CESC cohort (n = 304), we quantified ferroptosis- and cuproptosis-related pathway activities by ssGSEA/GSVA and constructed a Metal Death Activity Score (MDAS). Weighted gene co-expression network analysis (WGCNA) was applied to identify MDAS-associated modules and hub genes. We further integrated single-cell RNA sequencing (scRNA-seq) data comprising 74,454 cells to characterize the cellular heterogeneity and dynamic evolution of metal death phenotypes within the tumor microenvironment. Core genes causally associated with cervical cancer risk were screened using two-sample Mendelian randomization (MR), and the optimal model among 74 machine-learning algorithms was selected to construct a Metal Death Risk Score (MDRS). At the therapeutic level, potential targeted agents were identified via network pharmacology, and PDGFRB–drug interactions were validated using molecular docking, molecular dynamics simulations, and MM/PBSA binding free energy calculations. To bridge molecular-scale evidence with tissue- and single-cell–scale effects, we developed a structure–context–coupled network propagation (SINP) model and performed single-cell pharmacodynamic simulations using scTenifoldKnk, enabling cross-scale mechanistic inference and validation. Results: MDAS was significantly elevated in cervical cancer tissues, while ferroptosis and cuproptosis were largely independent at the global level (r = 0.053, P = 0.355). WGCNA identified 293 hub genes across five MDAS-associated modules. Single-cell analyses revealed higher MDAS activation in myeloid cells and tumor epithelial cells, with a progressive decline along the epithelial malignant transformation trajectory (ρ = −0.231, P = 3.5×10^-44). MR analysis identified eight causal genes, including four risk factors (DACT1, PDGFRB, PRSS23, MYO15B) and four protective factors (MSRB3, CALD1, DAB2, BNC2). The optimal MDRS model (Elastic Net) achieved a C-index of 0.736 after integration with clinical variables and significantly improved risk reclassification (NRI = 0.073, P = 0.044). High-risk patients exhibited enhanced epithelial–mesenchymal transition (EMT), angiogenesis, and suppressed oxidative phosphorylation. Network pharmacology identified sorafenib as a dual-function candidate drug capable of both targeting PDGFRB and inducing ferroptosis. Molecular docking indicated stable binding of sorafenib to the ATP-binding pocket of the PDGFRB kinase domain (Vina score = −10.1 kcal/mol), which was further supported by molecular dynamics simulations, MM/PBSA analysis (ΔG_PB = −34.09 ± 0.26 kcal/mol), and surface plasmon resonance (SPR) validation (KD = 1.95 μM). At the tissue scale, SINP predicted that PDGFRB inhibition markedly suppresses collagen-enriched extracellular matrix (ECM) programs. scTenifoldKnk simulations further demonstrated cell context–dependent effects, revealing that PDGFRB drives pro-fibrotic reprogramming in monocytes and triggers collapse of the collagen network. Cross-scale consistency analysis converged on seven shared ECM core genes (COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, A2M), establishing the PDGFRB–monocyte–collagen axis as a key mechanistic pathway linking the high-MDRS phenotype, stromal stiffening, and EMT. Conclusions: This study establishes an integrative framework encompassing metal death pathway activity scoring, causal gene identification, machine-learning–based risk stratification, and multiscale mechanistic simulation. MDRS enables clinical risk stratification in cervical cancer patients and highlights sorafenib as a potential precision therapeutic candidate for high-MDRS patients, with mechanisms likely involving PDGFRB inhibition, ferroptosis induction, and microenvironmental remodeling through suppression of monocyte-driven pro-fibrotic programs. Oncology cervical cancer ferroptosis cuproptosis machine learning Mendelian randomization single-cell RNA sequencing network pharmacology molecular dynamics systems pharmacology simulation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Introduction Cervical cancer (CESC) remains one of the most common malignancies among women worldwide and continues to impose a substantial public health burden[ 1 , 2 ]. Although the implementation of HPV vaccination programs and improvements in screening systems have significantly reduced incidence and mortality in many countries, pronounced geographic disparities in disease burden persist due to unequal vaccine coverage, limited access to screening, and imbalances in therapeutic resources[ 1 , 2 ]. In particular, for patients with advanced, recurrent, or metastatic disease, the overall efficacy of current treatment strategies—including chemoradiotherapy, anti-angiogenic therapy, and immunotherapy—remains limited, with marked inter-patient heterogeneity in treatment response and the frequent emergence of drug resistance[ 3 – 5 ]. Therefore, elucidating the key molecular mechanisms underlying cervical cancer progression, identifying more reliable risk-driving factors, and developing prognostic models and therapeutic strategies applicable to clinical stratification remain pressing scientific challenges. Dysregulation of regulated cell death (RCD) is a central event in tumor initiation and progression[ 6 – 8 ]. In recent years, metal-dependent forms of cell death—particularly ferroptosis and cuproptosis—have attracted increasing attention because of their close associations with metabolic reprogramming, oxidative stress, and mitochondrial function. Ferroptosis is primarily driven by iron-dependent accumulation of lipid peroxides and is tightly regulated by antioxidant systems such as GPX4 and glutathione metabolism[ 9 ]; In contrast, cuproptosis is linked to copper-mediated aberrant aggregation of lipoylated proteins in the mitochondrial tricarboxylic acid (TCA) cycle, leading to proteotoxic stress[ 10 ]. Although these two pathways differ in their initiating mechanisms and molecular nodes, they share overlapping contexts involving metal homeostasis imbalance, mitochondrial metabolic dependency, and tumor cell adaptation to microenvironmental stress[ 11 , 12 ]. To date, substantial evidence has accumulated regarding the roles of ferroptosis or cuproptosis in various cancers. However, in cervical cancer, whether these two pathways exhibit coordinated regulation or a pattern of “relatively independent coexistence,” how their activities are distributed across different cell types within the tumor microenvironment, and whether such differences can explain heterogeneity in patient prognosis and therapeutic response remain largely unexplored in a systematic, multi-omics manner. Most previous studies have relied primarily on bulk transcriptomic correlations, which, although informative for identifying pathway–phenotype associations, often fail to distinguish true drivers of tumor risk or progression from bystander signals and cannot resolve the differential contributions of distinct cellular lineages within the microenvironment. With advances in single-cell RNA sequencing (scRNA-seq), causal inference methodologies, and machine-learning modeling, multi-scale integrative approaches now offer new avenues to address these challenges. scRNA-seq enables the characterization of tumor microenvironment heterogeneity and evolutionary trajectories at single-cell resolution[ 13 , 14 ]; Mendelian randomization (MR) leverages genetic variants as instrumental variables to provide a relatively robust framework for identifying causally associated candidate genes and potential therapeutic targets[ 15 , 16 ]; and machine-learning models allow the construction of clinically oriented risk scores and stratification tools in multi-gene, high-dimensional settings, provided that rigorous feature selection, cross-validation, and biological interpretability are ensured[ 17 ]. Moreover, static evidence derived solely from network pharmacology or molecular docking is often insufficient to establish a credible causal chain linking drug–target interactions to tissue-level phenotypic alterations. Incorporating molecular dynamics simulations and binding free-energy calculations to provide thermodynamic support, and further connecting these findings to downstream signaling via network propagation or single-cell perturbation simulations, can substantially enhance mechanistic rigor and testability[ 18 , 19 ]. Against this background, the present study aims to establish an integrative framework that spans from “metal-related cell death phenotypes” to “causal target identification–clinical risk stratification–candidate therapeutic strategies.” We integrated ferroptosis- and cuproptosis-related features in the TCGA-CESC cohort to construct a Metal Death Activity Score (MDAS), and applied weighted gene co-expression network analysis (WGCNA) to identify MDAS-associated modules and hub genes. We then leveraged scRNA-seq data to delineate the heterogeneity of MDAS across different cellular lineages and its dynamic changes during disease progression, followed by two-sample MR analysis to screen core genes with causal associations with cervical cancer risk. Based on the MR-identified causal gene set, we further developed and evaluated multiple machine-learning–based survival models to derive a Metal Death Risk Score (MDRS) for clinical risk stratification. At the therapeutic level, we employed network pharmacology to identify candidate drugs and validated key target–drug interactions through molecular docking. We further incorporated molecular dynamics simulations and MM/PBSA binding free-energy calculations to quantitatively characterize the stability of drug–target binding, and combined a structure–context–coupled network propagation (SINP) model with single-cell virtual perturbation analysis (scTenifoldKnk) to achieve multiscale mechanistic inference and validation from the molecular level to tissue and cellular scales. Through this strategy, we aim to provide a more interpretable and verifiable theoretical and analytical foundation for precision stratification and mechanism-driven targeted therapy in cervical cancer. Materials and Methods Data Sources and Preprocessing This study integrated multi‑omics data from several public databases to investigate cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC). Bulk transcriptomic data and corresponding clinical information from The Cancer Genome Atlas (TCGA) were obtained via the Genomic Data Commons (GDC) portal, including 304 primary tumor samples with complete survival information. Gene expression profiles (HTSeq‑FPKM) were log2‑transformed with a pseudocount of 1. In addition, 70 normal cervical tissue samples from the Genotype‑Tissue Expression (GTEx) project were incorporated for comparative analyses. To complement bulk‑level analyses, single‑cell transcriptomic data were retrieved from the Gene Expression Omnibus (GEO; accession GSE208653). This integrated dataset initially comprised 74,454 cells derived from cancer tissues (n = 22,141), high‑grade squamous intraepithelial lesions (HSIL; n = 12,333), and normal tissues (n = 39,980). High‑quality cells were retained after filtering based on the following criteria: detected gene number between 200 and 5,000, mitochondrial read proportion < 20%, and hemoglobin gene expression < 5%. For Mendelian randomization analyses, summary statistics of blood cis‑expression quantitative trait loci (cis‑eQTLs) for the identified hub genes were obtained from the eQTLGen Consortium (n = 31,684; accessed via IEU OpenGWAS) and used as exposure data, while cervical cancer genome‑wide association study (GWAS) summary statistics were obtained from the FinnGen consortium (dataset ID: finn‑b‑C3_CERVIX_UTERI) and used as outcome data. Finally, curated gene sets representing metal‑related regulated cell death pathways were compiled, including 512 ferroptosis‑related genes (drivers, suppressors, and markers) from FerrDb V2 and 43 cuproptosis‑related genes reported in the literature. Construction of the Metal Death Activity Score (MDAS) Based on the TCGA cervical cancer (CESC) cohort, gene set variation analysis (GSVA) was applied to quantify ferroptosis (FAS) and cuproptosis (CAS) pathway activities for each sample. To address scale heterogeneity and construct an integrated metric, FAS and CAS were independently standardized using Z‑score transformation. The standardized values were then summed to define the Metal Death Activity Score (MDAS = Z‑FAS + Z‑CAS). An integrated waterfall plot was used to systematically visualize the distribution of MDAS and the relative contributions of its components. For stratification analyses, in addition to dichotomizing patients into high‑ and low‑MDAS groups based on the median MDAS, patients were further subdivided into four subtypes according to the independent median cutoffs of FAS and CAS (dual‑high, dual‑low, and discordant groups). The prognostic value of these stratification strategies was evaluated using Kaplan–Meier survival analysis and the log‑rank test. Finally, linear models implemented in the limma package were used to identify differentially expressed genes (DEGs) between high‑ and low‑MDAS groups as well as among specific quartile‑based subtypes, thereby elucidating potential underlying molecular mechanisms. Weighted Gene Co‑expression Network Analysis (WGCNA) WGCNA was performed to identify gene modules associated with MDAS and to explore transcriptomic regulatory networks driving metal death activity. First, the top 5,000 genes with the highest variance were selected to reduce noise and enhance analytical robustness, followed by sample clustering to detect and remove outliers. The optimal soft‑thresholding power (β) was determined using the pickSoftThreshold function to construct a signed co‑expression network satisfying scale‑free topology criteria. The topological overlap matrix (TOM) was then calculated, and gene modules were identified using dynamic tree cutting (minModuleSize = 30, mergeCutHeight = 0.25). Module eigengenes (MEs) were computed and subjected to Pearson correlation analysis to assess their relationships with MDAS and other clinical traits; modules with P 0.8 and an absolute gene significance (GS) > 0.2 were defined as MDAS‑associated hub genes. Single‑cell RNA‑seq Analysis Cell Clustering and Hierarchical Annotation To accurately identify cellular subpopulations within the tumor microenvironment, a two‑level hierarchical manual annotation strategy based on canonical marker gene expression was employed. At the first level, cell clusters were classified into three major lineages: immune cells (PTPRC), epithelial cells (EPCAM, KRT8, KRT5), and stromal cells (COL1A1, PECAM1). At the second level, more refined annotations were performed based on lineage‑specific molecular features: Epithelial lineage: subdivided into squamous epithelial cells (KRT5, KRT14, TP63), glandular/columnar epithelial cells (KRT18, KRT19), and proliferative tumor cells (MKI67, TOP2A). Notably, high expression of CDKN2A (p16) was used to identify cell populations with HPV‑associated features. Lymphoid lineage: T cells were identified by CD3D/CD3E and further subdivided into CD4⁺ naïve/memory T cells (IL7R), CD8⁺ cytotoxic T cells (GZMB, GNLY), and CD8⁺ exhausted T cells (PDCD1, HAVCR2, LAG3). NK cells were identified by NKG7 and TYROBP, while B cells and plasma cells were marked by MS4A1 and JCHAIN/MZB1, respectively. Myeloid lineage: including monocytes (CD14, VCAN), macrophages (CD68, CD163, C1QA), and dendritic cells (CD1C, CLEC9A). Stromal lineage: comprising fibroblasts (DCN, LUM) and endothelial cells (VWF, PECAM1). Calculation of MDAS Using UCell To robustly assess ferroptosis and cuproptosis pathway activities at the single‑cell level, the UCell R package was employed. UCell computes gene set enrichment scores based on the Mann–Whitney U statistic, relying on relative gene expression ranks rather than absolute expression values, thereby ensuring robustness to dropout events and batch effects inherent to scRNA‑seq data. Ferroptosis (FAS) and cuproptosis (CAS) gene sets were input into the AddModuleScore_UCell function to calculate initial scores. To ensure comparability, raw UCell scores were Z‑score standardized. The MDAS for each cell was defined as the sum of the standardized scores: MDAS = Z(FAS_UCell) + Z(CAS_UCell). Pseudotime Trajectory Analysis To reconstruct epithelial cell developmental trajectories during cervical cancer progression, the Monocle3 R package was used. Epithelial subsets (including normal squamous epithelial cells, precancerous/HSIL cells, and tumor cells) were extracted from the integrated Seurat object and converted into a cell_data_set object. Data were preprocessed using principal component analysis (PCA), followed by batch‑effect correction based on sample IDs (orig.ident), and then visualized by UMAP for dimensionality reduction. The principal graph was constructed using the learn_graph function with use_partition = FALSE to enforce a single continuous trajectory. The node with the highest proportion of “normal epithelial” cells was defined as the root, and pseudotime values were computed for each cell. Biological relevance was validated by assessing the correspondence between pseudotime and disease stage. Spearman correlation analysis was then used to quantify associations between pseudotime and MDAS, FAS, and CAS. Finally, genes exhibiting significant expression changes along the trajectory were identified using the graph_test function based on Moran’s I statistic (q < 0.05), with particular emphasis on the dynamic expression patterns of hub genes. Mendelian Randomization–Based Causal Inference To further identify reliable therapeutic targets and validate potential causal associations between the identified hub genes and cervical cancer risk, a two‑sample Mendelian randomization (MR) analysis was performed. High‑confidence hub genes identified by WGCNA and scRNA‑seq analyses were used as exposure factors. Gene symbols were converted to Ensembl IDs using the biomaRt package. Cis‑expression quantitative trait loci (cis‑eQTLs) significantly associated with these genes were extracted from the IEU OpenGWAS database and used as IVs. The selection threshold was set at genome‑wide significance (P < 5 × 10⁻⁸). To ensure independence, linkage disequilibrium (LD) pruning was applied with r² < 0.001 within a 10,000‑kb window to remove correlated SNPs. Genetic summary statistics for cervical cancer were obtained from the FinnGen consortium (ID: finn‑b‑C3_CERVIX_UTERI). This dataset is based on a European population, consistent with the exposure data, thereby minimizing population stratification bias. All MR analyses were conducted using the TwoSampleMR package in R. Exposure and outcome datasets were harmonized to align effect alleles. The inverse‑variance weighted (IVW) method was used as the primary approach to estimate causal effects in terms of odds ratios (ORs). Complementary methods, including MR‑Egger regression, the weighted median, and the weighted mode, were applied to assess robustness. Cochran’s Q test was used to evaluate heterogeneity among IVs. Horizontal pleiotropy was assessed using the MR‑Egger intercept test (P > 0.05 indicating no significant pleiotropy). P values were adjusted for multiple testing using the Benjamini–Hochberg false discovery rate (FDR) method. Genes with P < 0.05 and consistent effect directions across multiple MR methods were considered potential causal genes. Machine Learning–Based Construction of the Metal Death Risk Score (MDRS) Based on the eight metal death–related causal genes identified by MR (e.g., DACT1, CALD1), prognostic models were constructed using an integrated framework comprising 74 machine‑learning algorithms derived from 10 base methods, including random survival forests (RSF), LASSO, Ridge, Elastic Net (Enet), Stepwise Cox (StepCox), CoxBoost, partial least squares regression for Cox models (plsRcox), and gradient boosting machines (GBM). Ten‑fold cross‑validation was performed, and the concordance index (C‑index) was calculated to select the optimal algorithm combination with the highest mean C‑index for construction of the Metal Death Risk Score (MDRS). Patients were stratified into high‑ and low‑risk groups based on the median MDRS. Kaplan–Meier survival curves and time‑dependent receiver operating characteristic (ROC) curves were used to evaluate model performance. To assess clinical applicability, a nomogram incorporating MDRS, age, clinical stage, tumor grade, and lymph node status was developed. Calibration curves and decision curve analysis (DCA) were used to evaluate the accuracy and net clinical benefit of the nomogram. In addition, net reclassification improvement (NRI) and integrated discrimination improvement (IDI) were calculated to quantify the incremental prognostic value of MDRS over traditional clinical features. Biological Characteristics and Immune Microenvironment Associated with MDRS Stratification To elucidate the biological mechanisms underlying MDRS‑based stratification, Hallmark gene set activities were first evaluated using GSVA. Linear modeling with limma identified pathways that were significantly different between risk groups (FDR < 0.05). For tumor microenvironment (TME) characterization, the ssGSEA algorithm was applied using ESTIMATE gene sets to calculate stromal and immune scores, and tumor purity was inferred using the Yoshihara formula. Differences in these scores between risk groups were assessed using the Wilcoxon rank‑sum test. In addition, the expression levels of key immune checkpoint genes (e.g., PDCD1, CTLA4, CD274) were compared between groups. Spearman correlation analysis was performed to explore potential regulatory relationships between MDRS model genes and key biological pathways. Network Pharmacology Analysis and Targeted Drug Prediction for MDRS To explore therapeutic strategies targeting MDRS‑associated genes, drug–target interaction data were integrated from DGIdb, DrugBank, and the Comparative Toxicogenomics Database (CTD) to construct a high‑confidence drug screening library. Precision therapeutic strategies were designed for the eight core causal genes (e.g., PDGFRB, DACT1), aiming to “inhibit risk genes” and “activate protective genes.” Identified candidate drugs were cross‑referenced with known ferroptosis and cuproptosis inducers (e.g., sorafenib, erastin) to identify agents with dual mechanisms of “targeting MDRS genes” and “inducing metal‑dependent cell death.” Sankey diagrams were generated using the ggalluvial package to visualize the logical flow of “gene–drug–therapeutic strategy,” and drug–target interaction networks were constructed using the ggraph package with a Kamada–Kawai layout. Molecular Docking To predict the binding affinity and interaction mode between sorafenib and the core target PDGFRB, molecular docking was performed using the CB‑Dock2 online server. PDGFRB comprises extracellular, transmembrane, and intracellular kinase domains, and sorafenib, as a tyrosine kinase inhibitor (TKI), exerts its inhibitory effect by binding to the intracellular kinase domain. Available crystal structures of PDGFRB in the RCSB PDB (e.g., 3MJG, 2L6W) were therefore considered unsuitable because they lack critical kinase regions. Consequently, a high‑confidence full‑length predicted structure of human PDGFRB was obtained from the AlphaFold Protein Structure Database (UniProt ID: P09619; model: AF‑P09619‑F1). The three‑dimensional structure of sorafenib was downloaded from PubChem (CID: 216239). Docking was conducted using the “Auto Blind Docking” mode in CB‑Dock2. The platform automatically preprocessed the protein structure and identified potential binding pockets based on surface curvature using the CurPocket module. Autodock Vina was then employed to perform docking calculations within the identified pockets. The conformation with the lowest Vina score (indicating the strongest predicted affinity) and a binding site consistent with kinase inhibitor characteristics was selected as the optimal binding mode for visualization and downstream analyses. Molecular Dynamics Simulation and MM/PBSA Binding Free Energy Calculation To quantitatively characterize the inhibitory potency of sorafenib toward PDGFRB at atomic resolution and to provide thermodynamic constraints for subsequent structure–network coupled simulations, all‑atom molecular dynamics (MD) simulations were performed, followed by MM/PBSA binding free energy calculations based on equilibrated trajectories. MD simulations of the PDGFRB–sorafenib complex were carried out using GROMACS (v2022.4). The protein structure was preprocessed by removing crystallographic water molecules and completing missing atoms. The protein was parameterized using the Amber99SB force field, while ligand topologies were generated with the GAFF force field. The system was solvated in a TIP3P water box with a minimum distance of 1.0 nm between the solute and the box boundary, and Na⁺/Cl⁻ ions were added to neutralize the system and maintain physiological ionic strength (0.15 M NaCl). After energy minimization using the steepest descent algorithm, the system was equilibrated under the NVT ensemble (300 K, 100 ps) and the NPT ensemble (1 atm, 100 ps), followed by a 10 ns production run with a time step of 2 fs. Bond lengths involving hydrogen atoms were constrained using the LINCS algorithm, and long‑range electrostatic interactions were treated using the particle mesh Ewald (PME) method. Subsequently, MM/PBSA binding free energy calculations were performed using gmx_MMPBSA (v1.6.3) on 301 frames extracted from the equilibrated portion of the trajectory (internal dielectric constant = 2.0, external dielectric constant = 80.0, ionic strength = 0.15 M). Total binding free energies under both PB and GB models, along with individual energy components (ΔE vdW , ΔE ele , ΔG solv ), were obtained and used as thermodynamic driving parameters for downstream network propagation modeling. Surface Plasmon Resonance (SPR) Analysis The binding affinity between sorafenib and PDGFRB was measured using surface plasmon resonance on a Biacore 1K system (Cytiva, Sweden). The recombinant human PDGFRB protein (0.2 mg/mL in PBS) was immobilized on a Series S Sensor Chip CM5 (Cytiva, Cat. No. 29149603) using the Amine Coupling Kit (Cytiva, Cat. No. BR100050). Briefly, the chip surface was activated by injecting a 1:1 mixture of N-hydroxysuccinimide (NHS) and 1-ethyl-3-(3-dimethylaminopropyl) carbodiimide (EDC) at a flow rate of 10 µL/min. PDGFRB protein was diluted to 40 µg/mL in 10 mM sodium acetate buffer (pH 4.0) and injected at 10 µL/min for immobilization, achieving a final immobilization level of approximately 13,000 response units (RU). Remaining active sites were blocked with ethanolamine-HCl. Sorafenib (10 mM stock in DMSO) was serially diluted in PBST running buffer to final concentrations of 0.003, 0.01, 0.04, 0.12, 0.37, 1.1, and 3.3 µM. Binding experiments were performed at 25°C using a flow path of 6 − 5 at a flow rate of 30 µL/min. Each concentration was injected for 90 seconds (association phase), followed by buffer flow for 120 seconds (dissociation phase). Three startup cycles were performed prior to sample injection. Reference-subtracted sensorgrams were analyzed using Biacore Insight Evaluation Software. The equilibrium dissociation constant (KD) was determined by fitting the steady-state response values to a 1:1 Langmuir binding model. Structure‑Informed Network Propagation (SINP) Model To bridge the gap between “atomic‑scale drug–target binding” and “tumor tissue‑scale phenotypic alterations,” we developed a structure‑context–coupled network propagation model termed Structure‑Informed Network Propagation (SINP). This framework integrates thermodynamic parameters derived from MD/MM‑PBSA analyses, transcriptomic co‑regulation relationships (WGCNA), and protein–protein interaction (PPI) networks to enable systems pharmacology simulations with explicit physicochemical constraints. Unlike conventional network models that treat target perturbation as a binary on/off event, SINP converts the sorafenib–PDGFRB binding free energy (ΔG) into a quantitative perturbation strength reflecting inhibitory efficiency. Perturbation efficiency was defined as Ekd = 1/(1 + e^0.2(ΔG + 25)), and was used as the initial perturbation coefficient for network propagation. At the network topology level, to enhance tissue specificity, a generic PPI network was projected onto TCGA‑CESC transcriptomic data, and interaction edges with weak co‑expression were pruned: gene pairs with low correlation (|r| < 0.15) in the cohort were considered inactive in the cervical cancer context and removed. In addition, WGCNA module information was incorporated to dynamically weight propagation: if two genes belonged to the same co‑expression module, their propagation weight was amplified 1.5‑fold to reinforce functionally coupled signal transmission. Network propagation was then initiated from PDGFRB in the cervical cancer–specific weighted network, yielding a perturbation score for each gene along with its predicted direction of regulation, thereby enabling prediction of systemic drug effects at the tissue scale. Single‑cell In Silico Pharmacodynamic Modeling (scTenifoldKnk) To validate the cellular sources of PDGFRB‑mediated effects and delineate downstream regulatory networks at single‑cell resolution, single‑cell in silico pharmacodynamic modeling was performed using scTenifoldKnk. This approach evaluates the impact of PDGFRB inhibition on transcriptional regulatory networks across different cell subpopulations. Based on cell type–specific expression profiles of PDGFRB in the scRNA‑seq dataset, monocytes—exhibiting high PDGFRB expression—were selected as the primary effector cell population for simulation. Proliferative epithelial tumor cells (Epi_Tumor_Prolif) with moderate PDGFRB expression were subjected to the same simulation as a comparative control to assess cell context dependency. Within each target cell population, single‑cell gene regulatory networks (scGRNs) were constructed, and PDGFRB inhibition was simulated by setting its outgoing edge weights to zero. Virtual perturbation scores (Z‑scores) were calculated based on changes in gene distances within the manifold space before and after perturbation. Genes with |Z| > 1.96 and FDR < 0.05 were defined as significantly affected downstream response genes of PDGFRB inhibition and were subsequently used for functional enrichment analyses and cross‑scale consistency validation with SINP predictions. Statistical Analysis All statistical analyses were performed using R (v4.5.0). Continuous variables were compared using the Wilcoxon rank‑sum test or the Kruskal–Wallis test, as appropriate. Correlations were assessed using Spearman’s method. Survival analyses were conducted using Kaplan–Meier curves with the log‑rank test, as well as Cox proportional hazards models. Multiple testing correction was performed using the Benjamini–Hochberg method. A two‑sided P value < 0.05 was considered statistically significant. Results 1. Construction and validation of the metal-dependent death activity score (MDAS) To quantify the activity levels of ferroptosis and cuproptosis in patients with cervical cancer, we constructed a metal-dependent death activity score (MDAS) system based on the TCGA-CESC cohort (n = 304) and normal cervical tissues from GTEx (n = 70). First, the ferroptosis activity score (FAS) and cuproptosis activity score (CAS) were calculated using the ssGSEA algorithm. Tumor-specific analyses showed that both FAS and CAS were significantly elevated in cervical cancer tissues compared with normal tissues (both P < 0.05; Figure 1), indicating aberrant activation of ferroptosis- and cuproptosis-related pathways in cervical cancer. Further correlation analysis revealed only a very weak association between FAS and CAS (Pearson r = 0.053, P = 0.355; Figure 1), suggesting that these two forms of metal-dependent cell death represent relatively independent biological processes in cervical cancer. Based on these findings, the standardized FAS and CAS were integrated into a unified metal-dependent death activity score (MDAS = FAS_z + CAS_z). Quartile-based stratification showed comparable proportions of patients in the high-iron/high-copper, high-iron/low-copper, low-iron/high-copper, and low-iron/low-copper groups (23.7%, 26.3%, 26.3%, and 23.7%, respectively; Figure 1). Notably, the “discordant groups” (high iron/low copper plus low iron/high copper) accounted for 52.6% of the cohort, further supporting the necessity of an integrated analysis. Cox regression analysis was performed to evaluate the prognostic value of FAS and CAS for overall survival (OS). The results showed that CAS was an independent prognostic risk factor (HR = 5.686, 95% CI: 1.009–32.038, P = 0.0488), whereas FAS did not exhibit a significant association with prognosis (HR = 1.368, P = 0.747). These findings suggest that cuproptosis may play a more critical role in the prognosis of cervical cancer and validate the rationale for integrating FAS and CAS to construct the MDAS. Based on the median MDAS value, patients were divided into a high-MDAS group (n = 152) and a low-MDAS group (n = 152) for subsequent multi-omics integrative analyses and classification of metal metabolism–related phenotypes. 2. Identification of MDAS-related gene co-expression modules by WGCNA To identify gene co-expression networks associated with MDAS, WGCNA was performed on transcriptomic data from 304 TCGA-CESC tumor samples. After data preprocessing and removal of low-variance genes, a total of 12,455 genes were included in the subsequent analysis. Soft-threshold power analysis indicated that a power of 9 satisfied the scale-free topology criterion (scale-free R² = 0.915; Figure 1). Using this parameter, a co-expression network was constructed, identifying 13 gene modules (excluding the gray module), with module sizes ranging from 44 genes (salmon module) to 2,352 genes (turquoise module) (Figure 1). Module–trait relationship analysis revealed five gene modules significantly associated with MDAS (P < 0.05; Figure 1). Among them, the magenta module showed the strongest negative correlation with MDAS (r = −0.248, P = 1.19 × 10^-5; 247 genes), whereas the red module (r = 0.213, P = 1.85 × 10^-4; 540 genes) and the pink module (r = 0.211, P = 2.08 × 10^-4; 462 genes) were positively correlated with MDAS. The turquoise module (r = −0.146, P = 0.011; 2,352 genes) and the brown module (r = 0.115, P = 0.046; 1,127 genes) also reached statistical significance. In addition, hub genes within each module were identified based on module membership (MM) and gene significance (GS). Using stringent criteria (|MM| > 0.8 and |GS| > 0.2) or relaxed criteria (|MM| > 0.6 and |GS| > 0.15), a total of 293 hub genes were identified (Figure 1). These genes were selected as candidate genes for subsequent Mendelian randomization analyses and risk score model construction. Notably, the magenta and turquoise modules (negatively correlated with MDAS) may represent molecular features associated with suppression of metal-dependent cell death, whereas the red, pink, and brown modules (positively correlated with MDAS) may reflect biological processes that promote metal-dependent cell death, providing important clues for subsequent functional enrichment analyses. 3. Single-cell analysis 3.1. Single-cell landscape of the metal-dependent death–related prognostic model To elucidate the biological functions of the constructed metal-dependent death activity score (MDAS) in cervical cancer at the single-cell level, we performed an integrated analysis of single-cell transcriptomic data derived from 22,141 cancer tissue cells, 12,333 high-grade squamous intraepithelial lesion (HSIL) cells, and 39,980 normal tissue cells. Through unsupervised clustering combined with marker gene–based annotation, we identified 13 major cell types, including lymphoid lineages (T cells, B cells, and NK cells), myeloid lineages (monocytes and macrophages), stromal lineages (fibroblasts and endothelial cells), and epithelial lineages (normal squamous epithelial cells and tumor epithelial cells) (Figure 2). We first assessed MDAS scores at a global level. The results showed that MDAS scores were significantly higher in HSIL and cancer tissues than in normal tissues, indicating a progressive activation of metal-dependent death pathways during cervical cancer development (Figure 2). At the lineage level, epithelial and myeloid cells exhibited the highest MDAS scores, whereas lymphocytes showed relatively low scores (Figure 2), highlighting marked heterogeneity in MDAS activation across different cellular lineages. 3.2. Activation patterns of MDAS in specific cell subpopulations To more precisely localize the cells responsible for MDAS activation, we ranked and visualized MDAS scores across the 13 identified cell subpopulations. Violin plots showed that CD8+T cells, B cells, and CD4+T cells had the lowest MDAS scores, whereas macrophages, proliferative tumor epithelial cells (Epi_Tumor_Prolif), and tumor epithelial cells (Epithelial_Tumor) exhibited the highest scores (Figure 2). These findings were further supported by UMAP visualizations, in which MDAS-high cells were predominantly enriched in T-cell clusters as well as in subsets of myeloid and epithelial clusters (Figure 2). Percentage-stacked bar plots clearly demonstrated that more than 60% of cells within T-cell subpopulations belonged to the MDAS-high phenotype, while this proportion approached 40% in macrophages (Figure 2). Importantly, we employed split violin plots to investigate changes in MDAS within specific cell types across disease progression. Remarkably, although the overall MDAS scores in T cells were not high, macrophages and tumor epithelial cells displayed a pronounced and progressive increase in MDAS scores from normal tissue to cancer (Figure 4D). These results strongly suggest that macrophages and tumor cells within the tumor microenvironment are the central cellular players mediating MDAS-related functions and may actively drive disease progression. Lineage Characterization of MDAS. (A) UMAP dimensionality reduction clustering plot of the single-cell atlas. Top left: Colored by sample source (Group: GTEx, HSIL, Normal, TCGA). Bottom left: Colored by disease status (Disease). Top right: Colored by single-cell subpopulation annotation (sng.ident), displaying 13 major cell types. (B) Heatmap showing the expression of the top 5 marker genes across cell subpopulations. The horizontal axis represents different cell subpopulations, and the vertical axis represents marker genes; yellow indicates high expression. (C) Dot plot of Level 1 Broad Lineage Markers. Illustrates the expression patterns of canonical markers for immune cells (PTPRC), epithelial cells (EPCAM), and fibroblasts (COL1A1), etc. Dot size represents the percentage of expressing cells, and color intensity represents average expression level. (D) Dot plot of Level 2 Cervical Cancer Specific Markers. Further refines the specific molecular features of each cell subpopulation. (E) Violin plot of MDAS distribution across four broad cell lineages (Epithelial, Lymphoid, Myeloid, Stromal). Results show higher MDAS scores in Epithelial and Myeloid lineages. (F) Violin plot of MDAS distribution across detailed cell subpopulations (ranked by median score). Tumor epithelial cells (Epithelial_Tumor) and Macrophages exhibit the highest metal death activity, while CD8+ T cells and B cells show the lowest scores. (G) Split violin plots showing MDAS changes within cell subpopulations across disease progression. Compares score differences among Normal, High-Grade Squamous Intraepithelial Lesion (HSIL), and Cancer states, revealing a significant increase in MDAS with malignancy in Macrophages and Tumor Epithelial cells. (H) Feature plots projecting Ferroptosis Activity Score (FAS_sc), Cuproptosis Activity Score (CAS_sc), and the integrated score (MDAS_sc) onto the single-cell UMAP. Red areas indicate high scores, visually demonstrating the enrichment of metal death activity in specific cell clusters. 3.3. scRNA-seq validation of prognostic hub genes Next, we validated the expression patterns of key hub genes identified by WGCNA at the single-cell level. We first calculated WGCNA module scores and found that the Brown and Pink modules were expressed across multiple cell types, whereas the Turquoise module showed high expression predominantly in myeloid and epithelial cells (Figure 3). Representative top hub genes from different modules were selected for visualization. Annotated heatmaps clearly illustrated the modular expression patterns: MDK and LCN2 from the Turquoise module were highly and specifically expressed in macrophages and tumor epithelial cells, while COL4A1 and THBS2 from the Brown module were mainly expressed in fibroblasts (Figure 3). DotPlots further confirmed these findings, with prominent high-expression “red dots” precisely mapping hub genes to specific cell types (Figure 3). Finally, UMAP plots intuitively visualized expression hotspots of key genes (e.g., MDK and LCN2), which showed strong spatial overlap with macrophage and tumor cell clusters (Figure 3). Together, these single-cell–level results not only cross-validated the robustness of our WGCNA analysis but also reaffirmed macrophages and tumor cells as the key cellular subpopulations underpinning MDAS-associated functions. 3.4. Pseudotime analysis reveals cell type–specific dynamic evolution of MDAS during malignant transformation Using the Monocle3 algorithm, we analyzed epithelial cell subpopulations and successfully reconstructed a continuous trajectory of malignant transformation from normal epithelium through precancerous lesions (HSIL) to tumor cells. The distribution of pseudotime values strictly followed the order of disease progression (P < 0.0001; Figure 4). Along this evolutionary trajectory, MDAS scores exhibited a significant negative correlation with pseudotime (Spearman ρ = −0.231, P = 3.5 × 10^-44), with a more pronounced decline in ferroptosis activity (FAS) compared with cuproptosis activity (CAS). This suggests that suppression of metal-dependent death pathways—particularly ferroptosis—may be a critical mechanism by which epithelial cells acquire a survival advantage during malignant transformation (Figure 4). In contrast to the dynamic changes observed in epithelial cells, pseudotime analysis of myeloid cells clearly captured the differentiation process from monocytes to macrophages; however, MDAS scores remained relatively stable along this differentiation trajectory (P = 0.185). This indicates that the elevated MDAS observed in the myeloid lineage is primarily driven by disease state rather than by the differentiation process itself (Figure 4). In addition, dynamic expression changes of hub genes such as CLSPN and DDR2 along pseudotime further supported these findings, collectively delineating a lineage-specific and spatiotemporal evolutionary landscape of MDAS across different cell lineages within the cervical cancer microenvironment (Figure 4). 4. Identification of causal genes for cervical cancer by Mendelian randomization To establish causal relationships between MDAS-related hub genes and cervical cancer susceptibility, we performed a two-sample Mendelian randomization (MR) analysis. Blood-derived cis-eQTLs from the eQTLGen Consortium were used as instrumental variables, and cervical cancer GWAS summary statistics from FinnGen were used as the outcome data. Among the 211 high-confidence hub genes identified by WGCNA, 204 were successfully mapped to Ensembl IDs, of which 123 had available cis-eQTL instruments (P < 5 × 10^-8; LD clumping: r^2< 0.001, window = 10,000 kb). In total, 361 independent SNPs were used as instrumental variables. Inverse variance–weighted (IVW) MR analysis identified eight genes with nominally significant causal associations with cervical cancer risk (P 1): DACT1 showed the strongest risk effect (OR = 2.66, 95% CI: 1.50–4.73, P = 0.001), followed by PDGFRB (OR = 1.66), PRSS23 (OR = 1.44), and MYO15B (OR = 1.13). The remaining four genes were protective factors (OR < 1): MSRB3 exhibited the strongest protective effect (OR = 0.48, 95% CI: 0.27–0.83, P = 0.009), followed by CALD1 (OR = 0.84), DAB2 (OR = 0.72), and BNC2 (OR = 0.78). Sensitivity analyses confirmed the robustness of these findings. No significant horizontal pleiotropy was detected for any of the causal genes (MR-Egger intercept P > 0.05), and no significant heterogeneity was observed (Cochran’s Q P > 0.05), indicating that the observed causal effects were not driven by invalid instruments or pleiotropic bias. Notably, DACT1, a negative regulator of the Wnt signaling pathway, exhibited the strongest causal association with cervical cancer risk, suggesting that Wnt dysregulation may play a critical role in tumorigenesis. PDGFRB, which encodes platelet-derived growth factor receptor beta, is involved in angiogenesis and stromal remodeling—hallmarks of tumor progression. Among the protective genes, MSRB3, a methionine sulfoxide reductase, participates in defense against oxidative stress, whereas DAB2 is a well-established tumor suppressor involved in the TGF-β signaling pathway. 5. Construction of a metal-dependent death risk score (MDRS) using 74 machine learning algorithms To develop a robust prognostic risk score based on the eight MR-validated causal genes, we systematically evaluated 74 combinations of machine learning algorithms using a 10-fold cross-validation framework. The algorithm pool comprised 10 core methods, including LASSO-Cox, Ridge-Cox, Elastic Net, Stepwise Cox, random survival forest (RSF), gradient boosting machine (GBM), CoxBoost, and their parameter variants (Figure 6). Among all tested algorithms, Elastic Net (α = 0.5) combined with lambda.min (ElasticNet_a5_lambdamin) achieved the highest cross-validated C-index (0.627) and was therefore selected as the optimal algorithm for MDRS construction. The final MDRS model incorporated all eight causal genes with optimized coefficient weights to calculate individualized risk scores for each patient. 5.1. Prognostic performance and risk stratification of MDRS Using the median MDRS value as the cutoff, patients were stratified into high-risk (n = 146) and low-risk (n = 145) groups. Kaplan–Meier survival analysis demonstrated that patients in the high-risk group had significantly worse overall survival (OS) than those in the low-risk group (log-rank P = 0.0032; Figure 6). Time-dependent ROC analysis showed moderate predictive accuracy, with AUC values of 0.637, 0.647, 0.642, and 0.601 for 1-, 2-, 3-, and 5-year OS, respectively (Figure 6). The distribution of MDRS across patients exhibited a clear gradient, with higher risk scores associated with shorter survival times and increased mortality (Figure 6). These results confirm that MDRS effectively captures prognostic heterogeneity driven by expression patterns of metal-dependent death–related genes. 5.2. Integration of MDRS with clinicopathological features: nomogram construction Given that the standalone MDRS model (C-index = 0.627) still had room for improvement, we investigated whether integrating MDRS with established clinicopathological prognostic factors could enhance predictive performance. A multivariable Cox proportional hazards model was constructed incorporating MDRS together with age, FIGO stage, histological grade, and lymph node status (N stage). The baseline model including only clinical features (age + stage + grade + N stage) yielded a C-index of 0.664. Notably, the addition of MDRS significantly improved model performance, increasing the C-index to 0.736 (ΔC-index = +0.071). This indicates that MDRS provides substantial independent prognostic information beyond traditional clinical parameters, including lymph node metastasis, a strong prognostic factor. A nomogram was constructed to visualize the integrated model and facilitate clinical application (Figure 6). Calibration curves demonstrated good agreement between predicted and observed 3-year survival probabilities, with data points closely aligned along the 45° reference line (Figure 6). 5.3. Incremental value assessment (NRI and IDI analyses) To quantify the prognostic improvement achieved by adding MDRS to the clinical baseline model, we calculated the net reclassification improvement (NRI) and integrated discrimination improvement (IDI) for 3-year survival. Incorporation of MDRS significantly improved risk reclassification, yielding a 3-year NRI of 0.073 (95% CI: 0.001–0.203, P = 0.044). This indicates that a proportion of patients misclassified by the clinical model were correctly reclassified into appropriate risk categories after inclusion of MDRS. In addition, IDI analysis showed a positive trend toward improvement (IDI = 0.294, 95% CI: −0.063–0.566, P = 0.106), further supporting the added clinical value of MDRS. 5.4. Decision curve analysis and subgroup validation Decision curve analysis (DCA) was performed to evaluate the clinical utility of the integrated model across a range of threshold probabilities. The DCA demonstrated that, across most clinically relevant thresholds, the integrated model (clinical features + MDRS; red line) consistently provided greater net benefit than the clinical baseline model alone (blue line) (Figure 6), confirming its superior decision-support capability. Given the clinical challenge of prognostic stratification in early-stage cervical cancer, we further validated the performance of MDRS in patients with FIGO stage I–II disease. Kaplan–Meier analysis showed that MDRS retained significant discriminative ability even within this early-stage subgroup, effectively distinguishing high- and low-risk patients (log-rank P = 0.01; Figure 6). This finding highlights the potential utility of MDRS for early identification of high-risk patients who may benefit from intensified treatment or closer surveillance. 6. Biological characteristics of MDRS-based risk stratification To elucidate the molecular mechanisms underlying the prognostic capability of MDRS, gene set variation analysis (GSVA) was performed to compare pathway activity differences between the MDRS-defined high-risk and low-risk groups. This analysis focused on Hallmark gene sets representing key cancer biological processes, excluding ferroptosis- and cuproptosis-related signatures that had already been analyzed during MDAS construction. 6.1. Differential pathway activity between risk groups GSVA identified 19 pathways with significantly different activity between the MDRS high- and low-risk groups (FDR < 0.05) (Figure 7). Specifically, 16 pathways were significantly upregulated and 3 pathways were downregulated in the high-risk group. The most prominently upregulated pathway in high-risk patients was epithelial–mesenchymal transition (EMT) (logFC = 0.450, FDR = 4.03 × 10^-39), followed by angiogenesis (logFC = 0.337, FDR = 8.79 × 10^-26), Hedgehog signaling (logFC = 0.226, FDR = 7.66 × 10^-17), and TGF-β signaling (logFC = 0.227, FDR = 2.56 × 10^-13) (Figure 7). These findings indicate that MDRS high-risk tumors possess enhanced invasive and metastatic potential. In addition, multiple oncogenic signaling pathways were activated in high-risk patients, including Notch signaling (logFC = 0.159, FDR = 5.46 × 10^-14), Wnt/β-catenin signaling (logFC = 0.143, FDR = 2.50 × 10^-9), and KRAS signaling upregulation (logFC = 0.153, FDR = 4.17 × 10^-9). Inflammatory-related pathways such as hypoxia (logFC = 0.095, FDR = 3.97 × 10^-4), IL2–STAT5 signaling (logFC = 0.096, FDR = 2.89 × 10^-4), and TNF-α/NF-κB signaling (logFC = 0.084, FDR = 0.036) were also significantly elevated. In contrast, pathways associated with normal cellular metabolism were suppressed in high-risk patients, particularly oxidative phosphorylation (logFC = −0.201, FDR = 3.29 × 10^-8), fatty acid metabolism (logFC = −0.088, FDR = 5.08 × 10^-5), and DNA repair (logFC = −0.078, FDR = 4.42 × 10^-3). The downregulation of oxidative phosphorylation concomitant with hypoxia activation suggests a metabolic shift toward the Warburg effect in high-risk tumors. 6.2. Tumor microenvironment characteristics Analysis of tumor microenvironment (TME) scores revealed striking differences in stromal composition between the risk groups (Figure 7). Compared with low-risk patients, high-risk patients exhibited significantly higher stromal scores (P 0.05). Consequently, high-risk tumors displayed significantly lower tumor purity (P < 0.0001), reflecting a more complex and stroma-rich microenvironment. The increased stromal content in high-risk tumors is consistent with the enhanced EMT and angiogenic features observed, as cancer-associated fibroblasts and endothelial cells are major contributors to stromal components. 6.3. Immune checkpoint expression profiles To evaluate the potential responsiveness to immunotherapy across different risk groups, we analyzed the expression of immune checkpoint genes (Figure 7). Among the checkpoints examined, PDCD1LG2 (PD-L2) was significantly upregulated in high-risk patients (P = 0.001), along with CD80 (P = 0.006) and SIGLEC15 (P = 0.006). Interestingly, the key immunosuppressive enzyme IDO1 was significantly downregulated in the high-risk group (P = 0.045). The upregulation of PD-L2 and SIGLEC15 in high-risk patients suggests that these individuals may be responsive to immune checkpoint blockade targeting these molecules, providing a theoretical basis for MDRS-guided immunotherapy stratification. 6.4. Associations between MDRS genes and pathways Correlation analysis between the eight MDRS causal genes and Hallmark pathways revealed distinct functional associations (Figure 7). Risk-promoting genes (PDGFRB, DACT1, and PRSS23) showed strong positive correlations with EMT (r = 0.62–0.68), angiogenesis (r = 0.53–0.68), and TGF-β signaling (r = 0.42–0.64), while exhibiting negative correlations with oxidative phosphorylation and DNA repair pathways. Notably, MYO15B, identified as a risk factor by MR analysis, displayed a unique pattern of negative correlations with most pro-tumorigenic pathways, suggesting a complex regulatory role. In contrast, protective genes (MSRB3, DAB2, BNC2, and CALD1) also showed positive correlations with EMT-related pathways, which may reflect their involvement in stromal remodeling processes that potentially exert tumor-suppressive effects through enhanced immune surveillance. 7. Network pharmacology analysis and molecular docking validation 7.1. Construction of the drug–target interaction network To identify potential therapeutic agents targeting MDRS signature genes, we integrated data from the DGIdb, DrugBank, and CTD databases to perform a network pharmacology analysis. Based on the principle of “inhibiting risk genes while activating or preserving protective genes,” candidate drugs targeting the eight MDRS causal genes were systematically screened. In total, 15 drug–gene interactions were identified (Figure 8). Among them, 12 were classified as therapeutically beneficial, targeting four key genes: PDGFRB (six drugs), DACT1 (two drugs), PRSS23 (two drugs), and MSRB3 (two drugs). The remaining interactions involving CALD1 and DAB2 were categorized as potentially harmful or context-dependent, as inhibition of these protective genes could inadvertently promote tumor progression. 7.2. Target-based therapeutic strategies Based on the distinct functional roles of these targets, specific therapeutic strategies were proposed. PDGFRB (OR = 1.656) was identified as the most druggable risk gene and can be targeted by a range of FDA-approved tyrosine kinase inhibitors (TKIs) with established safety profiles, including imatinib, sunitinib, sorafenib, dasatinib, regorafenib, and lenvatinib. Similarly, other risk factors also provided actionable intervention points: the Wnt signaling regulator DACT1 (OR = 2.664) can be targeted by inhibitors such as ICG-001 and PRI-724, which are currently under clinical investigation; and the serine protease PRSS23 (OR = 1.438) offers opportunities for drug repurposing using FDA-approved protease inhibitors such as nafamostat and camostat. In contrast to these inhibitory strategies, therapeutic intervention for the protective factor MSRB3 (OR = 0.478) involves functional enhancement. As a key gene in oxidative stress defense, its activity may be augmented by antioxidants such as N-acetylcysteine (NAC) and ebselen, thereby restoring cellular resistance to metal-induced cell death. 7.3. Identification of a dual-function drug: sorafenib Notably, our analysis identified sorafenib as a dual-function therapeutic agent of particular relevance for MDRS high-risk patients. Sorafenib not only inhibits the MDRS risk gene PDGFRB, but is also a well-established inducer of ferroptosis by inhibiting System Xc⁻ and subsequently depleting glutathione. This dual mechanism—simultaneously targeting the MDRS molecular signature and inducing metal-dependent cell death—positions sorafenib as an optimal candidate drug for high-risk cervical cancer patients identified by MDRS. 7.4. Molecular docking validation of the sorafenib–PDGFRB interaction To validate the predicted drug–target interaction at the molecular level, molecular docking simulations were performed between sorafenib and the PDGFRB protein structure. Using the AlphaFold-predicted human PDGFRB structure (UniProt ID: P09619), docking results demonstrated that sorafenib binds tightly within the ATP-binding pocket of the intracellular kinase domain of PDGFRB through hydrogen bonds and hydrophobic interactions. Blind docking analysis identified the highest-affinity binding pocket (Pocket C4) with a binding energy (Vina score) of −10.1 kcal/mol, indicating strong binding potential. Binding mode analysis (as shown in the figure) revealed that sorafenib primarily interacts with key residues in the kinase core region, including residues around Val607, Lys634, and Glu651, as well as the DFG motif (Asp822 and Phe823). This binding pattern is characteristic of a type II kinase inhibitor, which suppresses PDGFRB phosphorylation activity by occupying the ATP-binding site and inducing or stabilizing an inactive (DFG-out) conformation. In addition, a high template-matching score of 66.8 based on a homologous structure (PDB ID: 7MGJ) further validated the reliability of this binding conformation. However, docking scores mainly reflect geometric complementarity and potential energy estimates under static conformations and cannot directly quantify the thermodynamic stability of drug–target binding. Therefore, molecular dynamics simulations and MM/PBSA calculations were subsequently performed to validate this interaction at both kinetic and thermodynamic levels. 7.5. Molecular dynamics simulations and MM/PBSA binding free energy calculations confirm strong thermodynamic affinity of PDGFRB–sorafenib To obtain quantitative pharmacodynamic constraints suitable for systems pharmacology modeling and to verify the stability of the docking pose on a dynamic timescale, all-atom molecular dynamics simulations were conducted for the PDGFRB–sorafenib complex, followed by MM/PBSA binding free energy calculations based on equilibrated trajectories (Figure 8). Throughout the production simulations, the overall conformation of the complex remained stable, indicating that sorafenib forms a persistent and robust binding mode within the PDGFRB kinase pocket. MM/PBSA results showed that the binding free energy calculated using the Poisson–Boltzmann (PB) model was ΔG PB = −34.09 ± 0.26 kcal/mol, while the Generalized Born (GB) model yielded ΔG GB = −25.56 kcal/mol. These highly negative binding free energies indicate strong thermodynamic affinity of sorafenib for PDGFRB, corresponding theoretically to nanomolar or even lower dissociation constant ranges. Energy decomposition analysis further revealed that binding was dominated by van der Waals interactions (ΔE vdW = −40.53 kcal/mol), with additional contributions from electrostatic interactions (ΔE ele = −6.59 kcal/mol), whereas desolvation introduced a positive solvation penalty (ΔG solv = +13.03 kcal/mol). Overall, the strong nonbonded interactions were sufficient to overcome the solvation cost, rendering the complex thermodynamically stable. Collectively, these results provide high-confidence molecular-scale evidence supporting sorafenib-mediated inhibition of PDGFRB and offer quantitative parameters for subsequent translation of ΔG values into initial perturbation strengths in network propagation models (SINP). To experimentally validate the computational predictions, surface plasmon resonance (SPR) analysis was performed. The results demonstrated that sorafenib binds to the PDGFRB kinase domain with an equilibrium dissociation constant (KD) of 1.95 μM (Supplementary Figure1), confirming a biophysically relevant interaction consistent with the strong binding affinity predicted by molecular docking (−10.1 kcal/mol) and MM/PBSA calculations (ΔG = −34.09 kcal/mol). 8. Multi-scale mechanistic validation of sorafenib targeting PDGFRB 8.1. SINP predicts that PDGFRB inhibition suppresses the collagen-rich ECM program To dissect the downstream systemic effects of sorafenib-mediated PDGFRB targeting at the tissue scale and to translate molecular-scale thermodynamic evidence into a propagatable pharmacodynamic perturbation, we quantitatively mapped the MM/PBSA-derived binding free energy (ΔG) to the initial inhibition strength of PDGFRB and incorporated it into the SINP model to perform structure–context–coupled network propagation simulations (Figure 8). Driven by thermodynamically informed perturbation coefficients, SINP predicted that inhibition of PDGFRB induces widespread reprogramming of downstream gene expression networks, with the most strongly downregulated gene sets being highly enriched in pathways related to extracellular matrix (ECM) organization, collagen fibril assembly, and matrix remodeling (Figure 8). Among the genes exhibiting the greatest downregulation, collagen family members and canonical stromal molecules were prominently represented, including COL1A1, COL1A2, COL3A1, COL5A2, and COL6A3, suggesting that PDGFRB signaling sustains a “collagen- and matrix-enhanced” pro-fibrotic transcriptional program within the cervical cancer microenvironment. These findings provide tissue-scale mechanistic insights into the stromal enrichment and EMT activation observed in MDRS high-risk tumors. 8.2. Monocyte-specific PDGFRB-driven pro-fibrotic reprogramming revealed by single-cell in silico pharmacodynamic modeling To identify the principal cellular source of the PDGFRB-mediated ECM program, we systematically evaluated the expression distribution of PDGFRB across cell subpopulations in the scRNA-seq dataset. PDGFRB exhibited marked cell type–specific enrichment within the tumor microenvironment, with monocytes serving as the predominant expressing cell population (expression proportion >80%), whereas expression was low or absent in macrophages and other immune cell types (Figure 8). Subsequent single-cell pharmacodynamic modeling using scTenifoldKnk within monocytes revealed that PDGFRB inhibition led to pronounced perturbations in the single-cell regulatory network, with differentially affected genes predominantly enriched in collagen biosynthesis, ECM organization, and matrix remodeling pathways. Specifically, multiple collagen family genes (COL1A1, COL3A1, COL5A1/2, COL6A1/3) showed significant responses, accompanied by marked effects on the matrix remodeling enzyme MMP11 and the collagen cross-linking enzyme LOXL2 (Figure 8). In addition, the immune regulatory gene LAMP5 exhibited the highest Z-score among perturbed genes, suggesting that beyond its pro-fibrotic role, PDGFRB may also contribute to maintenance of an immunosuppressive microenvironment. As a control, the same perturbation analysis performed in the Epi_Tumor_Prolif subpopulation did not reveal significant network disruption (Hits = 0), indicating that the functional impact of PDGFRB is highly context dependent. These results suggest that PDGFRB exerts its key biological effects not within tumor epithelial cells themselves, but primarily by reprogramming monocytes toward a pro-fibrotic phenotype, thereby shaping a collagen-enriched tumor microenvironment. 8.3. Cross-scale convergence identifies seven shared ECM core genes across SINP and single-cell modeling To further assess the concordance between tissue-scale simulations (SINP) and single-cell modeling (scTenifoldKnk), we intersected the top 200 downregulated genes predicted by the SINP model with the significantly perturbed genes identified in monocytes (Figure 8). This cross-scale comparison identified seven fully concordant core genes: COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, and A2M. All of these genes encode key structural components or regulatory molecules of the ECM, collectively constituting a highly convergent collagen–matrix core program. This cross-scale consistency not only validates the tissue-contextual relevance of the SINP model, but also provides single-cell–level confirmation that PDGFRB inhibition can systematically dismantle the collagen-enriched ECM network. Together, these findings define a clear mechanistic axis—PDGFRB–monocyte–collagen—through which sorafenib may exert therapeutic effects in MDRS high-risk cervical cancer. Discussion This study systematically characterized the phenotypic heterogeneity of ferroptosis- and cuproptosis-related pathways in cervical cancer (CESC) by integrating multi-omics data, causal inference, and machine learning modeling, and further established a mechanistic, chain-like linkage between risk stratification results, potentially actionable targets, and candidate therapeutic agents. In contrast to previous studies that focused on a single mode of cell death or a single omics layer, we not only constructed a clinically applicable risk score (MDRS) at the population level, but also established an interpretable, multi-scale mechanistic evidence chain spanning tissue- and single-cell resolutions. Starting from genetically supported causal targets identified by Mendelian randomization (MR), we integrated molecular-scale drug–target thermodynamic evidence (molecular docking and MM/PBSA), tissue-scale signal propagation inference (SINP), and single-cell–level identification of effector cell populations and regulatory network perturbations (scTenifoldKnk). These analyses ultimately converged on a well-defined microenvironmental regulatory axis—the PDGFRB–monocyte–collagen/ECM axis. This axis not only provides a more direct mechanistic explanation for the stromal enrichment and enhanced epithelial–mesenchymal transition (EMT) observed in MDRS high-risk tumors, but also offers a more testable theoretical basis for the potential benefit of sorafenib in specific cervical cancer subtypes. At the pathway level, we observed only a weak overall correlation between ferroptosis and cuproptosis in CESC, with more than half of patients exhibiting discordant phenotypes such as “high ferroptosis/low cuproptosis” or “low ferroptosis/high cuproptosis.” This finding does not imply that the two pathways are independent or unrelated; rather, it likely reflects the strong context dependence and spatial heterogeneity of metal-dependent cell death pathways across patients, cell lineages, and stages of tumor evolution. Single-cell analyses further supported this interpretation: myeloid cells and tumor epithelial cells generally exhibited higher MDAS activation, whereas pseudotime trajectories revealed a progressive decline in MDAS—particularly the ferroptosis-related component—during epithelial malignant transformation. This pattern suggests that, during tumor progression, cancer cells may acquire a survival advantage by reducing their sensitivity to metal-dependent death through metabolic reprogramming. It should be emphasized, however, that these observations provide associative evidence at this level, and the causal direction and underlying molecular mechanisms require further experimental validation. A key strength of this study lies in the use of MR to elevate candidate hub genes from correlation-based signals to high-confidence targets with causal associations, thereby providing genetic support for downstream drug development and mechanistic inference. Among the eight causal genes, factors such as DACT1 implicate Wnt and developmental regulatory networks in the risk phenotype, whereas the importance of PDGFRB stems from its well-established roles in angiogenesis, stromal activation, and fibrotic processes[ 20 – 23 ]. Recent studies have demonstrated that inducing transient target expression can engineer therapeutic vulnerabilities, thereby expanding the scope of targeted cancer therapies[ 24 ]. Importantly, we did not stop at identifying PDGFRB as a “risk gene,” but instead used multi-scale modeling and cross-scale consistency validation to delineate its potential effector pathway. At the tissue scale, SINP simulations predicted that PDGFRB inhibition would markedly suppress collagen-enriched ECM programs. At the single-cell scale, scTenifoldKnk revealed a pronounced cell-context dependency: PDGFRB inhibition triggered substantial collapse of collagen- and matrix remodeling–related networks in monocytes, whereas comparable perturbations were largely absent in tumor epithelial cell controls. These findings suggest that the risk effect of PDGFRB is primarily mediated through myeloid cell–driven pro-fibrotic reprogramming and matrix deposition, rather than solely through intrinsic tumor cell proliferation, thereby mechanistically linking MR-inferred genetic risk with the stromal-enriched/EMT-enhanced phenotype of MDRS high-risk tumors. To minimize methodological bias from any single model, we further performed cross-scale consistency analyses by intersecting downregulated genes predicted by tissue-scale SINP simulations with significantly perturbed genes identified by single-cell scTenifoldKnk modeling. This approach yielded seven fully overlapping ECM core genes (COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, and A2M). These genes are canonical components and regulators of collagen and matrix remodeling, and are closely associated with stromal stiffening, cell migration, and EMT programs in cancer[ 25 – 27 ]. This cross-scale convergence substantially strengthens the robustness of our conclusions, demonstrating that the link between PDGFRB inhibition and ECM program collapse is not an artifact of a single algorithm, but a reproducible biological theme detectable across scales and data structures. Collectively, these findings establish the PDGFRB–monocyte–collagen axis as a key mechanistic pathway connecting risk stratification with microenvironmental remodeling. From a translational perspective, network pharmacology analyses identified sorafenib as a candidate drug with dual potential to both target PDGFRB and induce ferroptosis, and this hypothesis was further reinforced by molecular-scale simulations. Emerging targeted protein degradation technologies may offer additional strategies for targeting previously undruggable components of these pathways[ 28 , 29 ]. Importantly, molecular docking was not treated as standalone evidence; instead, molecular dynamics simulations combined with MM/PBSA binding free energy calculations provided thermodynamically more meaningful quantitative support. The binding free energy of sorafenib to PDGFRB reached ΔG_PB = − 34.09 ± 0.26 kcal/mol under the PB model, indicating strong binding stability. Using this thermodynamic constraint, we quantitatively mapped ΔG to the initial perturbation strength in the SINP model, enabling inference of system-level effects at the tissue scale, which were then validated at the single-cell level as primarily occurring within monocyte-driven pro-fibrotic networks. Thus, the potential role of sorafenib in MDRS high-risk CESC can be more rationally summarized as operating through two complementary pathways: (i) weakening tumor cell survival advantages by inducing ferroptosis at the tumor cell level[ 9 , 30 , 31 ]; and (ii) attenuating myeloid cell–driven collagen deposition by inhibiting PDGFRB signaling within the tumor microenvironment, thereby “softening” the matrix and reducing EMT-supportive conditions. This mechanistic framework not only provides an actionable interpretation of MDRS high-risk phenotypes, but also suggests that the therapeutic benefit of sorafenib in cervical cancer may be closely tied to patient heterogeneity, warranting further investigation in stratified clinical studies. In addition, MDRS, as a machine learning–based risk score derived from a set of causal genes, not only significantly improved predictive performance when integrated with clinical variables, but also enabled a more biologically interpretable stratification. The MDRS high-risk group was characterized by aggressive features such as EMT, angiogenesis, hypoxia, and stromal enrichment, accompanied by upregulation of immune checkpoint molecules (e.g., PD-L2 and SIGLEC15), suggesting a “stroma-enriched/immune-excluded” tumor microenvironment state[ 32 – 34 ]. Novel approaches targeting PD-L1 stability through DHHC3 degradation have shown promise in overcoming immune checkpoint blockade resistance[ 35 , 36 ], suggesting that combination strategies targeting both stromal programs and immune checkpoints may be particularly beneficial in high-risk patients. Accordingly, the value of MDRS extends beyond prognostic prediction, as it may also serve as a mechanism-oriented patient stratification tool to guide future exploration of combination strategies involving stromal/ECM targeting, metal-dependent cell death induction, and immune modulation. Nonetheless, such strategies require rigorous experimental and clinical validation to assess feasibility and safety. Despite establishing a multi-scale framework linking causal genes to mechanistic validation, several limitations should be acknowledged. First, the MDRS model was primarily developed and internally evaluated using the TCGA cohort, and its generalizability requires further validation in independent, multi-center clinical cohorts. Second, the MR analysis relied on blood-derived cis-eQTLs from eQTLGen; although this resource offers large sample size and strong statistical power, the lack of tissue-specific eQTLs may limit the capture of cervix-local regulatory effects. Future availability of large-scale cervical tissue or single-cell eQTL datasets would enhance tissue consistency in causal inference. Third, SINP and scTenifoldKnk are computational simulation frameworks that provide mechanistic inference and consistency evidence; however, key conclusions—such as the role of PDGFRB in monocyte-driven ECM programs—require confirmation through in vitro and in vivo experiments. Finally, the potential benefit of sorafenib in MDRS high-risk subgroups remains to be validated in prospective clinical studies, along with the identification of optimal biomarker combinations for actionable stratified treatment strategies. The SPR analysis revealed that sorafenib binds to the PDGFRB kinase domain with an equilibrium dissociation constant (KD) of 1.95 µM, which falls within the moderate affinity range. This finding is consistent with the pharmacological profile of sorafenib as a multi-target tyrosine kinase inhibitor. As sorafenib was originally developed to target RAF kinases and subsequently found to potently inhibit VEGFR and PDGFR family members, PDGFRB represents a secondary rather than primary target. The observed micromolar affinity is therefore biologically plausible and aligns with previous reports characterizing sorafenib's broad kinase inhibition spectrum. Importantly, this moderate binding affinity, combined with the favorable thermodynamic profile revealed by MM/PBSA calculations (ΔG = − 34.09 kcal/mol) and the strong predicted binding pose from molecular docking (− 10.1 kcal/mol), collectively support PDGFRB as a therapeutically relevant target of sorafenib in cervical cancer. The convergence of computational predictions and experimental validation strengthens the rationale for repurposing sorafenib as a PDGFRB-targeted therapeutic strategy in this disease context. In summary, this study proposes and validates a coherent mechanistic axis linking metal-dependent cell death phenotypes, microenvironmental matrix remodeling, and clinical risk stratification, and establishes an integrative framework spanning pathway phenotypes, causal targets, interpretable prediction, multi-scale mechanistic validation, and candidate therapeutic strategies. This framework offers new insights into heterogeneity-aware management of cervical cancer and the development of mechanism-driven therapeutic approaches. Declarations Data availability statement: The datasets presented in this study are publicly available in online repositories. The bulk transcriptomic data and clinical information (TCGA and GTEx) can be accessed via the GDC Data Portal (https://portal.gdc.cancer.gov/). The single-cell RNA sequencing data generated in this study have been deposited in the Gene Expression Omnibus (GEO) under accession number GSE208653. For Mendelian randomization analysis, summary statistics were obtained from the eQTLGen Consortium (accessed via IEU OpenGWAS) and the FinnGen consortium (Dataset ID: finn-b-C3_CERVIX_UTERI). Funding Statement: This work was supported by the Hubei Provincial Natural Science Foundation Joint Fund (Grant No. 2025AFD319) and the Beijing Science and Technology Innovation Medical Development Foundation (Grant No. KC2023-JX-0288-RQ21). Conflict of interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. References WHO Guidelines Approved by the Guidelines Review Committee , in WHO guideline for screening and treatment of cervical pre-cancer lesions for cervical cancer prevention: Use of dual-stain cytology to triage women after a positive test for human papillomavirus (HPV) . 2024, World Health Organization © World Health Organization 2024.: Geneva. Singh, D., et al., Global estimates of incidence and mortality of cervical cancer in 2020: a baseline analysis of the WHO Global Cervical Cancer Elimination Initiative. Lancet Glob Health, 2023. 11 (2): p. e197-e206. Abu-Rustum, N.R., et al., NCCN Guidelines® Insights: Uterine Neoplasms, Version 3.2025. J Natl Compr Canc Netw, 2025. 23 (8): p. 284-291. Monk, B.J., et al., First-Line Pembrolizumab + Chemotherapy Versus Placebo + Chemotherapy for Persistent, Recurrent, or Metastatic Cervical Cancer: Final Overall Survival Results of KEYNOTE-826. J Clin Oncol, 2023. 41 (36): p. 5505-5511. Oaknin, A., et al., EMPOWER CERVICAL-1: Effects of cemiplimab versus chemotherapy on patient-reported quality of life, functioning and symptoms among women with recurrent cervical cancer. Eur J Cancer, 2022. 174 : p. 299-309. Galluzzi, L., et al., Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ, 2018. 25 (3): p. 486-541. Galluzzi, L., et al., Regulated cell death and adaptive stress responses. Cell Mol Life Sci, 2016. 73 (11-12): p. 2405-10. Peng, F., et al., Regulated cell death (RCD) in cancer: key pathways and targeted therapies. Signal Transduct Target Ther, 2022. 7 (1): p. 286. Jiang, X., B.R. Stockwell, and M. Conrad, Ferroptosis: mechanisms, biology and role in disease. Nat Rev Mol Cell Biol, 2021. 22 (4): p. 266-282. Tsvetkov, P., et al., Copper induces cell death by targeting lipoylated TCA cycle proteins. Science, 2022. 375 (6586): p. 1254-1261. Liu, N. and M. Chen, Crosstalk between ferroptosis and cuproptosis: From mechanism to potential clinical application. Biomed Pharmacother, 2024. 171 : p. 116115. Müller, S., et al., Copper and iron orchestrate cell-state transitions in cancer and immunity. Trends Cell Biol, 2025. 35 (2): p. 105-114. Xia, P., et al., Deciphering the cellular and molecular landscape of cervical cancer progression through single-cell and spatial transcriptomics. NPJ Precis Oncol, 2025. 9 (1): p. 158. Ou, Z., et al., Single-Nucleus RNA Sequencing and Spatial Transcriptomics Reveal the Immunological Microenvironment of Cervical Squamous Cell Carcinoma. Adv Sci (Weinh), 2022. 9 (29): p. e2203040. Hu, X., et al., Benchmarking Mendelian randomization methods for causal inference using genome-wide association study summary statistics. Am J Hum Genet, 2024. 111 (8): p. 1717-1735. Hemani, G., et al., The MR-Base platform supports systematic causal inference across the human phenome. Elife, 2018. 7 . Huang, S., et al., Artificial intelligence in cancer diagnosis and prognosis: Opportunities and challenges. Cancer Lett, 2020. 471 : p. 61-71. Hou, T., et al., Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. J Chem Inf Model, 2011. 51 (1): p. 69-82. Hopkins, A.L., Network pharmacology: the next paradigm in drug discovery. Nat Chem Biol, 2008. 4 (11): p. 682-90. Strell, C., E. Rodríguez-Tomàs, and A. Östman, Functional and clinical roles of stromal PDGF receptors in tumor biology. Cancer Metastasis Rev, 2024. 43 (4): p. 1593-1609. Pandey, P., et al., New insights about the PDGF/PDGFR signaling pathway as a promising target to develop cancer therapeutic strategies. Biomed Pharmacother, 2023. 161 : p. 114491. Shi, J., et al., Pericyte-myofibroblast transition: a novel mechanism in peritoneal fibrosis and the effect of Asiaticoside. J Adv Res, 2025. Carroll, S.H., et al., Genetic requirement of dact1/2 to regulate noncanonical Wnt signaling and calpain 8 during embryonic convergent extension and craniofacial morphogenesis. Elife, 2024. 13 . Yang, J., et al., Targeting an Inducible SALL4-Mediated Cancer Vulnerability with Sequential Therapy. Cancer Res, 2021. 81 (23): p. 6018-6028. Jiang, Y., et al., Targeting extracellular matrix stiffness and mechanotransducers to improve cancer therapy. J Hematol Oncol, 2022. 15 (1): p. 34. Yuan, Z., et al., Collagen remodeling-mediated signaling pathways and their impact on tumor therapy. J Biol Chem, 2025. 301 (3): p. 108330. Suzuki, M., et al., Association of rs6983561 polymorphism at 8q24 with prostate cancer mortality in a Japanese population. Clin Genitourin Cancer, 2011. 9 (1): p. 46-52. Dai, M., et al., Targeted Protein Degradation: An Important Tool for Drug Discovery for "Undruggable" Tumor Transcription Factors. Technol Cancer Res Treat, 2022. 21 : p. 15330338221095950. Hashmi, F., et al., EXPRESS: Phospholipase C gamma mediates endogenous brain-derived neurotrophic factor - regulated calcitonin gene-related peptide expression in colitis - induced visceral pain. Mol Pain, 2016. 12 . Zhang, L., et al., Sorafenib triggers ferroptosis via inhibition of HBXIP/SCD axis in hepatocellular carcinoma. Acta Pharmacol Sin, 2023. 44 (3): p. 622-634. Liu, R., et al., Roles and Mechanisms of Ferroptosis in Sorafenib Resistance for Hepatocellular Carcinoma. J Hepatocell Carcinoma, 2024. 11 : p. 2493-2504. Zheng, S., et al., Tumor battlefield within inflamed, excluded or desert immune phenotypes: the mechanisms and strategies. Exp Hematol Oncol, 2024. 13 (1): p. 80. Chen, D.S. and I. Mellman, Oncology meets immunology: the cancer-immunity cycle. Immunity, 2013. 39 (1): p. 1-10. Li, L., et al., The role of Siglec-15 in tumor immunity: mechanism and therapy. Mol Cancer Ther, 2025. Shi, Y.Y., et al., Treating ICB-resistant cancer by inhibiting PD-L1 via DHHC3 degradation induced by cell penetrating peptide-induced chimera conjugates. Cell Death Dis, 2024. 15 (9): p. 701. Liu, M., et al., A replication study examining three common single-nucleotide polymorphisms and the risk of prostate cancer in a Japanese population. Prostate, 2011. 71 (10): p. 1023-32. Additional Declarations The authors declare no competing interests. Supplementary Files SuppFigure1.tif Supplementary Figure1. Surface Plasmon Resonance (SPR) validation of sorafenib binding to PDGFRB kinase domain. (A) Steady-state affinity analysis. The equilibrium response values were plotted against sorafenib concentrations (0.003–3.3 μM) and fitted to a 1:1 Langmuir binding model. The equilibrium dissociation constant (KD) was determined to be 1.95 μM. (B) SPR sensorgrams showing the real-time binding kinetics of sorafenib to immobilized PDGFRB at seven concentrations ranging from 0.003 to 3.3 μM. The concentration-dependent increase in response units (RU) and the rapid association/dissociation kinetics are consistent with the Type II kinase inhibitor binding mode predicted by molecular docking analysis. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-8811979","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":587279975,"identity":"eb17b8e6-26bb-4877-9a12-d143a4e7e706","order_by":0,"name":"Shan Li","email":"","orcid":"","institution":"Department of Gynecological Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China.","correspondingAuthor":false,"prefix":"","firstName":"Shan","middleName":"","lastName":"Li","suffix":""},{"id":587279976,"identity":"a4c5603d-bfe8-465f-89a6-a0b79c3290ca","order_by":1,"name":"Mingyuan Wu","email":"","orcid":"","institution":"School of Modern Industry for Selenium Science and Engineering, Wuhan Polytechnic University, Wuhan, Hubei, China;","correspondingAuthor":false,"prefix":"","firstName":"Mingyuan","middleName":"","lastName":"Wu","suffix":""},{"id":587279977,"identity":"7e337371-a773-4ff2-85bd-c836dc2f7a13","order_by":2,"name":"Yuanyuan Meng","email":"","orcid":"","institution":"Department of Obstetrics and Gynecology, Huanggang Women and Children Health Hospital, Huanggang, Hubei, China","correspondingAuthor":false,"prefix":"","firstName":"Yuanyuan","middleName":"","lastName":"Meng","suffix":""},{"id":587279978,"identity":"ed1cfd38-12bc-4eb4-828a-deb67978d6c1","order_by":3,"name":"Caixia Liang","email":"","orcid":"","institution":"Department of Gynecological Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China.","correspondingAuthor":false,"prefix":"","firstName":"Caixia","middleName":"","lastName":"Liang","suffix":""},{"id":587279979,"identity":"f8c70918-8412-4f12-898d-5ccefaf0ee7a","order_by":4,"name":"Xiufang Huang","email":"","orcid":"","institution":"Department of Gynecological Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China.","correspondingAuthor":false,"prefix":"","firstName":"Xiufang","middleName":"","lastName":"Huang","suffix":""},{"id":587279980,"identity":"62bc6f59-b55e-46da-a198-da4602f8e700","order_by":5,"name":"Chaoming Li","email":"","orcid":"","institution":"Department of Gynecological Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China.","correspondingAuthor":false,"prefix":"","firstName":"Chaoming","middleName":"","lastName":"Li","suffix":""},{"id":587279981,"identity":"6b967e23-63ec-4ab4-b134-f0ebbf7dd5db","order_by":6,"name":"Chunxu Yang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA20lEQVRIiWNgGAWjYNCCCjDJ+CChooZYLWfAJLPBgzPHiNTB2Aam2CQftjATVm1wvMfwcWXbYTmD473PKhIb2Bj427sT8Gs5c8bY8My5w8YGZ46b3UjcIcMgcebsBvxabuSYSTaUHU7ccCON7UbiGTYGA4lcglrMfzawAbXcf8ZWkNjGTJQWM8aGNpAtbGwMRGmRPHOsWLLhTLqx5Jk0ZomEM8d4CPqF73jzxo8NFdZyfMePMX78UVEjx9/ei1+LwgEw1QwX4MGrHATkG8BUHUGFo2AUjIJRMIIBALp3TrWJeq7tAAAAAElFTkSuQmCC","orcid":"","institution":"Department of Radiation and Medical Oncology, Hubei Key Laboratory of Tumor Biological Behaviors, Hubei Cancer Clinical Study Center, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China","correspondingAuthor":true,"prefix":"","firstName":"Chunxu","middleName":"","lastName":"Yang","suffix":""},{"id":587279982,"identity":"e138616b-a963-4faa-8bc7-9564fd645ccd","order_by":7,"name":"Zhangqian Wang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAxUlEQVRIiWNgGAWjYPACGwYDMM1GvJY00rUcJkGLwfnDmz8Xtp1P3C7dY8DwoewwA//sBgJaDhwrk57Zdjtx55wzBowzzh1mkLhzAL8Ws4M9Zsy8bbdzN9zIMQAygC6USCCg5TCP8WfetnMQLX+J0nKMx0Cat+0ARAsjMVrsz7CVSfOcS67fcOdYwcGec+k8EjcIaJHsB4YYT5mdscHt5o0PfpRZy/HPIKAFCCAxwiDBwHAASPEQVI+iZRSMglEwCkYBVgAAzrFEir+j2M4AAAAASUVORK5CYII=","orcid":"","institution":"School of Modern Industry for Selenium Science and Engineering, Wuhan Polytechnic University, Wuhan, Hubei, China;","correspondingAuthor":true,"prefix":"","firstName":"Zhangqian","middleName":"","lastName":"Wang","suffix":""},{"id":587279983,"identity":"f94ac9fc-446f-4429-ab22-ad83c7bf3fea","order_by":8,"name":"Mengyuan Dai","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA7ElEQVRIiWNgGAWjYFCDA4yNDxgbQKwE4rU0G5CqhYFNgigtuu29h1/zttnl8R0/3FbNu+MwAz97jgHDzx24tZidOZdmzduWXCx5JrHtNu+ZwwySPW8MGHvP4NFyI8fMmLeNOXHDAaCW3LbDDAY3cgyYGdvwaLn/BqSlPnHD+YdtxSAt9gS13OAxfszbdjhxw43ENmawLRKEtJzJMWOcc+544swbD5ul/7al80iceVZwsBefluNnjD+8KatO7Duf/vDjzDZrOf725I0PfuLRAgRsUjxIPDD7AF4NDAzMH38QUDEKRsEoGAUjHAAAY65Z6B0SO+kAAAAASUVORK5CYII=","orcid":"","institution":"Department of Gynecological Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China.","correspondingAuthor":true,"prefix":"","firstName":"Mengyuan","middleName":"","lastName":"Dai","suffix":""}],"badges":[],"createdAt":"2026-02-07 03:40:00","currentVersionCode":1,"declarations":{"humanSubjects":false,"vertebrateSubjects":false,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":false,"humanSubjectConsent":false,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":false},"doi":"10.21203/rs.3.rs-8811979/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-8811979/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":102310422,"identity":"6b06ee3e-985a-42bc-b6fd-67107d11a079","added_by":"auto","created_at":"2026-02-10 11:53:47","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":977780,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eConstruction, Landscape Characterization of Metal Death Activity Score (MDAS), and WGCNA Co-expression Module Identification. \u003c/strong\u003e(A) Scatter plot illustrating the correlation between Ferroptosis Activity Score (FAS) and Cuproptosis Activity Score (CAS). Each point represents a patient, colored according to four subgroups defined by high and low levels of FAS and CAS (Ferro-High/Cupro-High, Ferro-High/Cupro-Low, Ferro-Low/Cupro-High, Ferro-Low/Cupro-Low). The Pearson correlation coefficient indicates no significant correlation between the two scores (r=0.053, P=0.355). (B) Integrated landscape of MDAS. Top: Waterfall plot of patient distribution ranked by MDAS (blue indicates MDAS-Low, red indicates MDAS-High). Middle: Distribution of standardized Z-scores for FAS and CAS corresponding to each patient. Bottom: Annotation bar showing the four-category subgroup classification for each patient. (C) Tumor-specific analysis. Violin plots displaying the differences in CAS (left) and FAS (right) scores between normal tissues (Normal) and primary tumors (Primary Tumor), showing significantly elevated scores for both in tumor tissues. (D) Bar chart showing the number and percentage distribution of patients across the four metal death subgroups. (E) Forest plot of Cox regression analysis evaluating the impact of FAS and CAS on overall survival. Results indicate that CAS is a significant independent risk factor (HR \u0026gt; 1, P=0.049). (F) Dendrogram of Weighted Gene Co-expression Network Analysis (WGCNA). Constructed based on TCGA-CESC transcriptome data, showing gene co-expression modules identified by the dynamic tree cut algorithm, marked by different colors. (G) Module-trait relationship heatmap. Displays the Pearson correlation coefficients and P-values between gene modules (rows) and clinical traits (columns, including MDAS, FAS, CAS, etc.). Red represents positive correlation, and blue represents negative correlation. (H) Scatter plot of Gene Significance (GS) versus Module Membership (MM) within key modules. Shows the distribution of Hub genes in the Magenta (negatively correlated), Pink (positively correlated), and Red (positively correlated) modules significantly associated with MDAS. (I) Bar plot ranking the correlation of each gene module with MDAS. Red indicates positive correlation, and gray indicates negative correlation, highlighting the gene modules most closely related to metal death activity.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/5a9632d7d99d7ffc165f2fad.png"},{"id":102310364,"identity":"704139ea-c47a-4b21-a6c4-337a2fd083c7","added_by":"auto","created_at":"2026-02-10 11:53:30","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1905224,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSingle-Cell Transcriptomic Atlas Mapping of Cervical Cancer and Cell Lineage Characterization of MDAS.\u003c/strong\u003e (A) UMAP dimensionality reduction clustering plot of the single-cell atlas. Top left: Colored by sample source (Group: GTEx, HSIL, Normal, TCGA). Bottom left: Colored by disease status (Disease). Top right: Colored by single-cell subpopulation annotation (sng.ident), displaying 13 major cell types. (B) Heatmap showing the expression of the top 5 marker genes across cell subpopulations. The horizontal axis represents different cell subpopulations, and the vertical axis represents marker genes; yellow indicates high expression. (C) Dot plot of Level 1 Broad Lineage Markers. Illustrates the expression patterns of canonical markers for immune cells (PTPRC), epithelial cells (EPCAM), and fibroblasts (COL1A1), etc. Dot size represents the percentage of expressing cells, and color intensity represents average expression level. (D) Dot plot of Level 2 Cervical Cancer Specific Markers. Further refines the specific molecular features of each cell subpopulation. (E) Violin plot of MDAS distribution across four broad cell lineages (Epithelial, Lymphoid, Myeloid, Stromal). Results show higher MDAS scores in Epithelial and Myeloid lineages. (F) Violin plot of MDAS distribution across detailed cell subpopulations (ranked by median score). Tumor epithelial cells (Epithelial_Tumor) and Macrophages exhibit the highest metal death activity, while CD8+ T cells and B cells show the lowest scores. (G) Split violin plots showing MDAS changes within cell subpopulations across disease progression. Compares score differences among Normal, High-Grade Squamous Intraepithelial Lesion (HSIL), and Cancer states, revealing a significant increase in MDAS with malignancy in Macrophages and Tumor Epithelial cells. (H) Feature plots projecting Ferroptosis Activity Score (FAS_sc), Cuproptosis Activity Score (CAS_sc), and the integrated score (MDAS_sc) onto the single-cell UMAP. Red areas indicate high scores, visually demonstrating the enrichment of metal death activity in specific cell clusters.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/68a61b87dd881acca3542667.png"},{"id":102310413,"identity":"9a6d903a-1940-4563-8b09-db578d942b87","added_by":"auto","created_at":"2026-02-10 11:53:46","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1107641,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSingle-Cell Validation of WGCNA Key Modules and Expression Pattern Analysis of MDAS-Related Hub Genes. \u003c/strong\u003e(A) Dot plot of Hub gene expression levels across cell subpopulations (ordered by WGCNA modules). Validates the specific expression of key genes from the Turquoise (e.g., MDK, LCN2) and Brown (e.g., COL4A1) modules in specific cell types. (B) Feature plots of selected Hub genes on the single-cell atlas. Displays the spatial expression distribution of genes such as ADAMTS5, NID2, and COL4A1, with red indicating high expression. (C) Heatmap of average MDAS scores for each cell subpopulation. Deeper red indicates a higher average MDAS score for that cell type, further confirming the high metal death activity status of Macrophages and Tumor Epithelial cells. (D) Violin plots of activity scores for WGCNA co-expression modules (Brown, Pink, Red, Turquoise) across cell subpopulations. Shows that the Turquoise module is highly expressed mainly in Myeloid and Epithelial cells, while the Brown module is primarily enriched in Fibroblasts. (E) Detailed violin plots of MDAS distinguishing disease states (Cancer vs. HSIL vs. Normal) across different cell types. Provides statistical analysis of MDAS changes with disease progression within each subpopulation. (F) Single-cell landscape of MDAS High/Low groups. Left: Distribution of MDAS-High (red) and MDAS-Low (blue) cells on UMAP. Right: Stacked bar chart showing the proportion of MDAS-High vs. MDAS-Low cells in each subpopulation, indicating a high proportion (\u0026gt;60%) of MDAS-High phenotype cells in T cell subsets and nearly 40% in Macrophages. (G) Grouped heatmap of Top Hub gene expression patterns. The horizontal axis represents cell types, and the vertical axis represents Hub genes; color bars at the top correspond to WGCNA modules and cell lineage classifications, clearly demonstrating the modular and cell-specific expression rules of key genes.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/8793145631a8c46ee4edfc5d.png"},{"id":102310363,"identity":"a80247b5-6482-49c4-94f8-fe4de8ec8989","added_by":"auto","created_at":"2026-02-10 11:53:30","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":1550942,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePseudotime Analysis Reveals Dynamic Evolution Patterns of MDAS During Epithelial Malignant Transformation and Myeloid Differentiation in Cervical Cancer. \u003c/strong\u003e(A) Pseudotime trajectory plot of epithelial malignant transformation (constructed by Monocle3). Left: Trajectory distribution colored by cell type (epithelial subpopulations). Right: Color mapping based on Pseudotime; purple represents early states, and yellow represents late states. (B) Validation boxplot of epithelial pseudotime versus disease stage. Shows a significant increasing trend (P\u0026lt;0.0001) of pseudotime values across Normal Epithelium, Precancerous, and Tumor stages, confirming the trajectory direction aligns with disease progression logic. (C) Dynamic changes of metal death scores along the epithelial-tumor trajectory. Top left: Scatter plot with fitted curve of MDAS score versus pseudotime, showing a significant negative correlation (Spearman ρ=−0.231). Top right: Violin plot distribution of MDAS across different pseudotime quartiles (Q1-Q4), showing significantly higher scores in early stages. Bottom row: Curves showing changes in Ferroptosis Activity Score (FAS, bottom left) and Cuproptosis Activity Score (CAS, bottom right) along pseudotime, with FAS showing a more pronounced declining trend. (D) Scatter plots of expression kinetics for key Hub genes along the epithelial cell pseudotime axis. Displays dynamic expression patterns where some genes (e.g., CLSPN, DDR2) are upregulated with malignant progression, while others (e.g., EFEMP1) are downregulated. (E) Landscape of myeloid cell differentiation trajectory. Colored from left to right by disease source (Normal, HSIL, Cancer), cell type (Monocytes, Macrophages), and pseudotime. (F) Changes in MDAS along the monocyte-macrophage differentiation trajectory. Left: Correlation analysis between overall MDAS and pseudotime (no significant correlation, ρ=−0.018). Right: Distribution of MDAS along pseudotime grouped by disease status. (G) Biological validation of myeloid cell pseudotime. Left: Distribution of pseudotime across different disease stages (Cancer, HSIL, Normal). Right: Comparison of pseudotime differences between Monocytes and Macrophages, confirming the temporal direction of differentiation from monocytes to macrophages (P\u0026lt;0.0001).\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/9c0c143e769ed61e1009f33c.png"},{"id":102310372,"identity":"df428016-c944-46ae-b90e-19e02407173f","added_by":"auto","created_at":"2026-02-10 11:53:41","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":903004,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMendelian Randomization (MR) analysis identifies causal associations between MDAS-related hub genes and cervical cancer risk. \u003c/strong\u003e(A, D) Forest plots of the eight significant causal genes. (A) displays the MR results for BNC2, CALD1, DAB2, and DACT1; (D) displays the results for MSRB3, MYO15B, PDGFRB, and PRSS23. Black dots and error bars represent the effect size (Beta) and standard error for individual SNPs, while red lines represent the pooled causal effect estimates derived from different MR methods (primarily IVW and MR-Egger). Results indicate that DACT1, PDGFRB, PRSS23, and MYO15B are significantly associated with an increased risk of cervical cancer (risk factors, OR \u0026gt; 1), whereas MSRB3, CALD1, DAB2, and BNC2 act as protective factors (OR \u0026lt; 1). (B, E) Funnel plots for heterogeneity testing. (B) corresponds to the four genes in Panel A; (E) corresponds to the four genes in Panel D. The x-axis represents the effect estimate of individual SNPs (β\u003csub\u003eIV\u003c/sub\u003e), and the y-axis represents the inverse of the standard error (1/SE\u003csub\u003eIV\u003c/sub\u003e). The symmetrical distribution of points indicates that the results are not significantly affected by publication bias or heterogeneity, supporting the robustness of the causal inference. (C, F) Scatter plots of SNP effect associations. (C) corresponds to genes in Panel A (BNC2, CALD1, DAB2, DACT1); (F) corresponds to genes in Panel D (MSRB3, MYO15B, PDGFRB, PRSS23). Plots display the correlation between the effect of instrumental variable SNPs on the exposure (gene expression, x-axis) and the outcome (cervical cancer risk, y-axis). Regression lines of different colors represent different MR methods (e.g., light blue for Inverse Variance Weighted, dark blue for MR Egger, green for Weighted Median), where the slope reflects the magnitude and direction of the causal effect.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/77e1b5fe7b9af104fd1b988a.png"},{"id":102310471,"identity":"787bbeff-b6a0-4ed2-bb87-23d336975173","added_by":"auto","created_at":"2026-02-10 11:54:06","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":599251,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eConstruction and evaluation of the Metal-related Death Risk Score (MDRS) and Nomogram based on machine learning. \u003c/strong\u003e(A) Bar chart comparing the C-index of 74 machine learning algorithm combinations. The red bar highlights the optimal algorithm combination (ElasticNet + Cox, alpha=0.5, lambda=min), which achieved the highest cross-validation C-index (0.627). (B) Kaplan-Meier survival curves for MDRS high-risk and low-risk groups in the TCGA-CESC cohort. Patients in the high-risk group exhibited significantly shorter overall survival compared to the low-risk group (Log-rank test, p=0.0032). (C) Time-dependent ROC curves for predicting 1-, 2-, 3-, and 5-year overall survival using MDRS. AUC values at each time point indicate that the model possesses moderate-to-high predictive accuracy. (D) Landscape of MDRS risk score distribution. Upper panel: Patients sorted by risk score from low to high (blue to red); Lower panel: Scatter plot of corresponding survival time and status (red dots indicate death, blue dots indicate survival), showing an increased density of death events as scores rise. (E) Nomogram for predicting 1-, 3-, and 5-year overall survival probabilities in cervical cancer patients. The model integrates MDRS with clinicopathological features including age, FIGO stage, pathological grade, and N stage. (F) Calibration curve for 3-year survival prediction. The gray dashed line represents the ideal prediction (45-degree line), and the black solid line represents the actual performance of the nomogram; the high overlap indicates good model calibration. (G) Decision Curve Analysis (DCA). The red curve represents the integrated model (Clinical + MDRS), the blue curve represents the baseline model (clinical features only), the black dashed line represents the \"Treat-All\" strategy, and the gray solid line represents the \"Treat-None\" strategy. Results show that the integrated model provides a higher net benefit than the baseline model across a wide range of threshold probabilities. (H) Kaplan-Meier survival curves for the early-stage (Stage I-II) patient subgroup. The MDRS significantly distinguished high-risk from low-risk populations even in early-stage patients (Log-rank test, p=0.01), confirming its potential value in early screening.\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/f5aa0e1d0228175b2965892e.png"},{"id":102310460,"identity":"639636e9-4d94-4155-810d-a6ee542e06f8","added_by":"auto","created_at":"2026-02-10 11:54:04","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":205226,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eBiological characteristics, tumor microenvironment landscape, and gene-pathway associations of MDRS risk stratification. \u003c/strong\u003e(A) Bar chart of Hallmark pathway activity differences based on GSVA analysis. Shows signaling pathways significantly differentially expressed between MDRS high-risk and low-risk groups (FDR \u0026lt; 0.05). Red bars indicate pathways significantly activated in the high-risk group (e.g., Epithelial-Mesenchymal Transition, Angiogenesis, TGF-β signaling), while blue bars indicate pathways significantly suppressed in the high-risk group (e.g., Oxidative Phosphorylation, Fatty Acid Metabolism). (B) Box plots comparing the activity of key tumor-related pathways. Detailed display of score differences for specific pathways such as Angiogenesis, DNA Repair, and EMT between high- and low-risk groups (∗p\u0026lt;0.05,∗∗p\u0026lt;0.01,∗∗∗p\u0026lt;0.001,∗∗∗∗p\u0026lt;0.0001, ns: not significant). (C) Assessment of Tumor Microenvironment (TME) characteristics. Immune Score, Stromal Score, and Tumor Purity were calculated and compared between the two groups using the ESTIMATE algorithm. The high-risk group exhibited significantly elevated Stromal Scores and reduced Tumor Purity. (D) Box plots of differential expression of immune checkpoint genes. Shows expression levels of key immune regulatory molecules such as CD274 (PD-L1), PDCD1LG2 (PD-L2), CTLA4, and IDO1 in MDRS high- and low-risk groups. (E) Heatmap of associations between MDRS causal genes and Hallmark pathways. Displays the Spearman correlation between the expression levels of 8 causal genes (y-axis) and the activity scores of various signaling pathways (x-axis). Red indicates positive correlation, and blue indicates negative correlation. Results reveal strong positive correlations between risk-promoting genes (e.g., PDGFRB, DACT1) and EMT/Angiogenesis pathways.\u003c/p\u003e","description":"","filename":"7.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/4bf3f46d2d663d4aee38b8a4.png"},{"id":102310468,"identity":"8f2d724a-8d33-4f1d-bfce-305f18facfe9","added_by":"auto","created_at":"2026-02-10 11:54:05","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":328570,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIdentification and mechanistic validation of PDGFRB-targeted therapeutic strategies based on multi-scale modeling. \u003c/strong\u003e(A) Therapeutic strategy flow Sankey diagram. Illustrates the mapping from MDRS causal genes (left) to pharmacological mechanism categories (middle) and specific candidate drugs (right). The width of the connecting flows represents the strength of association, indicating that tyrosine kinase inhibitors (TKIs) targeting PDGFRB and Wnt pathway inhibitors targeting DACT1 constitute the primary potential therapeutic strategies. (B) Drug–target interaction network. Central circular nodes represent MDRS target genes (red: risk genes; blue: protective genes), while peripheral green square nodes represent potential therapeutic drugs. Edges indicate database-predicted drug–target interactions, highlighting PDGFRB as a core target of multi-target TKIs such as sunitinib and sorafenib. (C) Molecular docking of sorafenib with the PDGFRB kinase domain. The docking binding energy is −10.1 kcal/mol, indicating a very strong binding affinity. The upper panel shows the surface electrostatic potential of the binding pocket, with sorafenib (stick model) stably embedded in the ATP-binding pocket (C4 pocket) of PDGFRB; the lower panel presents a cartoon view of the protein secondary structure, illustrating the position of the drug within the kinase core region. (D) Thermodynamic component analysis of PDGFRB–ligand binding. Energy decomposition histograms calculated using the MM/PBSA (Poisson–Boltzmann Surface Area) method are shown. The total binding free energy is −34.09 kcal/mol, demonstrating a very strong thermodynamic affinity between sorafenib and PDGFRB. The decomposition reveals that van der Waals interactions (−40.53 kcal/mol) are the dominant driving force, effectively compensating for the unfavorable desolvation penalty (+13.03 kcal/mol). (E) Per-residue energy decomposition. Displays the top 20 amino acid residues contributing most favorably to ligand binding (Top 20 favorable residues). Red bars indicate favorable energetic contributions, with LYS758, ARG785, and LYS387 showing particularly strong contributions, revealing key anchoring residues within the binding pocket. (F) Visualization of the drug-induced perturbation response network. Downstream gene regulatory networks reconstructed using the SINP model. Node color represents response magnitude (blue indicates downregulation), and node size reflects node type. The network illustrates how inhibition of PDGFRB propagates through the network topology, leading to pronounced downregulation of collagen family genes (e.g., COL1A1, COL1A2) and matrix remodeling–related genes (MMP2, VCAN). (G) Cell type–specific expression and co-expression profile of PDGFRB. The left violin plot shows PDGFRB expression across tumor microenvironment cell populations, confirming its highly specific expression in monocytes (\u0026gt;80%). The right dot plot demonstrates strong co-expression between PDGFRB and core extracellular matrix genes (COL1A1, COL3A1, LUM) in monocytes, suggesting cell type–specific functional coupling. (H) Topological relationship between network distance and perturbation strength. Box plots depict the relationship between a gene’s topological distance from PDGFRB and the magnitude of its downregulation. Genes classified as “Very Close” exhibit stronger suppression, validating the SINP model’s structure- and context-coupled propagation of pharmacological perturbations. (I) Ranked list of suppressed genes predicted by bulk SINP simulation. Shows the top genes most strongly downregulated following sorafenib treatment as predicted by the SINP algorithm. These genes are highly enriched for collagen family members (COLs) and ECM organization–related molecules, indicating that PDGFRB inhibition primarily disrupts the pro-fibrotic stromal network. (J) Monocyte-specific single-cell virtual knockout (Virtual KO) simulation. PDGFRB perturbation was simulated in the monocyte subset using the scTenifoldKnk algorithm. The scatter plot displays gene expression changes (perturbation Z-scores), highlighting significant downregulation of LAMP5 and collagen genes such as COL3A1. This result confirms, at the single-cell scale, that PDGFRB is a key driver of pro-fibrotic reprogramming in monocytes. (K) Consensus stromal signature across scales. Cross-validation of tissue-scale (SINP; y-axis: bulk differential score) and single-cell–scale (scTenifoldKnk; color/size: single-cell Z-score) simulations identifies seven consistently and strongly suppressed core genes (COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, A2M). These genes constitute a highly convergent collagen–stromal core program, establishing a cross-scale therapeutic axis of PDGFRB–monocyte–collagen matrix.\u003c/p\u003e","description":"","filename":"8.png","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/1a4b7cecf3e6443aa8020ea2.png"},{"id":102311605,"identity":"85340c2d-994e-48c3-a00e-2c92389a686e","added_by":"auto","created_at":"2026-02-10 11:58:26","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":8431087,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/20cb1112-c5cb-4f53-887f-7e3214a272e5.pdf"},{"id":102310367,"identity":"f6a727a0-7fc4-4855-a3a5-969e32a3261b","added_by":"auto","created_at":"2026-02-10 11:53:37","extension":"tif","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":4378612,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSupplementary Figure1. Surface Plasmon Resonance (SPR) validation of sorafenib binding to PDGFRB kinase domain.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e(A) Steady-state affinity analysis. The equilibrium response values were plotted against sorafenib concentrations (0.003–3.3 μM) and fitted to a 1:1 Langmuir binding model. The equilibrium dissociation constant (KD) was determined to be 1.95 μM.\u003c/p\u003e\n\u003cp\u003e(B) SPR sensorgrams showing the real-time binding kinetics of sorafenib to immobilized PDGFRB at seven concentrations ranging from 0.003 to 3.3 μM. The concentration-dependent increase in response units (RU) and the rapid association/dissociation kinetics are consistent with the Type II kinase inhibitor binding mode predicted by molecular docking analysis.\u003c/p\u003e","description":"","filename":"SuppFigure1.tif","url":"https://assets-eu.researchsquare.com/files/rs-8811979/v1/a6e6ea6249c59708fc451371.tif"}],"financialInterests":"The authors declare no competing interests.","formattedTitle":"\u003cp\u003e\u003cstrong\u003eIntegrated Multi-Omics Profiling Maps Ferroptosis–Cuproptosis Diversity in Cervical Cancer and Identifies a PDGFRB-Driven Monocyte Fibrotic Program Targeted by Sorafenib\u003c/strong\u003e\u003c/p\u003e","fulltext":[{"header":"Introduction","content":"\u003cp\u003eCervical cancer (CESC) remains one of the most common malignancies among women worldwide and continues to impose a substantial public health burden[\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Although the implementation of HPV vaccination programs and improvements in screening systems have significantly reduced incidence and mortality in many countries, pronounced geographic disparities in disease burden persist due to unequal vaccine coverage, limited access to screening, and imbalances in therapeutic resources[\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. In particular, for patients with advanced, recurrent, or metastatic disease, the overall efficacy of current treatment strategies\u0026mdash;including chemoradiotherapy, anti-angiogenic therapy, and immunotherapy\u0026mdash;remains limited, with marked inter-patient heterogeneity in treatment response and the frequent emergence of drug resistance[\u003cspan additionalcitationids=\"CR4\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. Therefore, elucidating the key molecular mechanisms underlying cervical cancer progression, identifying more reliable risk-driving factors, and developing prognostic models and therapeutic strategies applicable to clinical stratification remain pressing scientific challenges.\u003c/p\u003e \u003cp\u003eDysregulation of regulated cell death (RCD) is a central event in tumor initiation and progression[\u003cspan additionalcitationids=\"CR7\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]. In recent years, metal-dependent forms of cell death\u0026mdash;particularly ferroptosis and cuproptosis\u0026mdash;have attracted increasing attention because of their close associations with metabolic reprogramming, oxidative stress, and mitochondrial function. Ferroptosis is primarily driven by iron-dependent accumulation of lipid peroxides and is tightly regulated by antioxidant systems such as GPX4 and glutathione metabolism[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]; In contrast, cuproptosis is linked to copper-mediated aberrant aggregation of lipoylated proteins in the mitochondrial tricarboxylic acid (TCA) cycle, leading to proteotoxic stress[\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e]. Although these two pathways differ in their initiating mechanisms and molecular nodes, they share overlapping contexts involving metal homeostasis imbalance, mitochondrial metabolic dependency, and tumor cell adaptation to microenvironmental stress[\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e, \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eTo date, substantial evidence has accumulated regarding the roles of ferroptosis or cuproptosis in various cancers. However, in cervical cancer, whether these two pathways exhibit coordinated regulation or a pattern of \u0026ldquo;relatively independent coexistence,\u0026rdquo; how their activities are distributed across different cell types within the tumor microenvironment, and whether such differences can explain heterogeneity in patient prognosis and therapeutic response remain largely unexplored in a systematic, multi-omics manner. Most previous studies have relied primarily on bulk transcriptomic correlations, which, although informative for identifying pathway\u0026ndash;phenotype associations, often fail to distinguish true drivers of tumor risk or progression from bystander signals and cannot resolve the differential contributions of distinct cellular lineages within the microenvironment. With advances in single-cell RNA sequencing (scRNA-seq), causal inference methodologies, and machine-learning modeling, multi-scale integrative approaches now offer new avenues to address these challenges. scRNA-seq enables the characterization of tumor microenvironment heterogeneity and evolutionary trajectories at single-cell resolution[\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]; Mendelian randomization (MR) leverages genetic variants as instrumental variables to provide a relatively robust framework for identifying causally associated candidate genes and potential therapeutic targets[\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]; and machine-learning models allow the construction of clinically oriented risk scores and stratification tools in multi-gene, high-dimensional settings, provided that rigorous feature selection, cross-validation, and biological interpretability are ensured[\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. Moreover, static evidence derived solely from network pharmacology or molecular docking is often insufficient to establish a credible causal chain linking drug\u0026ndash;target interactions to tissue-level phenotypic alterations. Incorporating molecular dynamics simulations and binding free-energy calculations to provide thermodynamic support, and further connecting these findings to downstream signaling via network propagation or single-cell perturbation simulations, can substantially enhance mechanistic rigor and testability[\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eAgainst this background, the present study aims to establish an integrative framework that spans from \u0026ldquo;metal-related cell death phenotypes\u0026rdquo; to \u0026ldquo;causal target identification\u0026ndash;clinical risk stratification\u0026ndash;candidate therapeutic strategies.\u0026rdquo; We integrated ferroptosis- and cuproptosis-related features in the TCGA-CESC cohort to construct a Metal Death Activity Score (MDAS), and applied weighted gene co-expression network analysis (WGCNA) to identify MDAS-associated modules and hub genes. We then leveraged scRNA-seq data to delineate the heterogeneity of MDAS across different cellular lineages and its dynamic changes during disease progression, followed by two-sample MR analysis to screen core genes with causal associations with cervical cancer risk. Based on the MR-identified causal gene set, we further developed and evaluated multiple machine-learning\u0026ndash;based survival models to derive a Metal Death Risk Score (MDRS) for clinical risk stratification. At the therapeutic level, we employed network pharmacology to identify candidate drugs and validated key target\u0026ndash;drug interactions through molecular docking. We further incorporated molecular dynamics simulations and MM/PBSA binding free-energy calculations to quantitatively characterize the stability of drug\u0026ndash;target binding, and combined a structure\u0026ndash;context\u0026ndash;coupled network propagation (SINP) model with single-cell virtual perturbation analysis (scTenifoldKnk) to achieve multiscale mechanistic inference and validation from the molecular level to tissue and cellular scales. Through this strategy, we aim to provide a more interpretable and verifiable theoretical and analytical foundation for precision stratification and mechanism-driven targeted therapy in cervical cancer.\u003c/p\u003e"},{"header":"Materials and Methods","content":"\u003cp\u003eData Sources and Preprocessing\u003c/p\u003e \u003cp\u003eThis study integrated multi‑omics data from several public databases to investigate cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC). Bulk transcriptomic data and corresponding clinical information from The Cancer Genome Atlas (TCGA) were obtained via the Genomic Data Commons (GDC) portal, including 304 primary tumor samples with complete survival information. Gene expression profiles (HTSeq‑FPKM) were log2‑transformed with a pseudocount of 1. In addition, 70 normal cervical tissue samples from the Genotype‑Tissue Expression (GTEx) project were incorporated for comparative analyses. To complement bulk‑level analyses, single‑cell transcriptomic data were retrieved from the Gene Expression Omnibus (GEO; accession GSE208653). This integrated dataset initially comprised 74,454 cells derived from cancer tissues (n\u0026thinsp;=\u0026thinsp;22,141), high‑grade squamous intraepithelial lesions (HSIL; n\u0026thinsp;=\u0026thinsp;12,333), and normal tissues (n\u0026thinsp;=\u0026thinsp;39,980). High‑quality cells were retained after filtering based on the following criteria: detected gene number between 200 and 5,000, mitochondrial read proportion\u0026thinsp;\u0026lt;\u0026thinsp;20%, and hemoglobin gene expression\u0026thinsp;\u0026lt;\u0026thinsp;5%. For Mendelian randomization analyses, summary statistics of blood cis‑expression quantitative trait loci (cis‑eQTLs) for the identified hub genes were obtained from the eQTLGen Consortium (n\u0026thinsp;=\u0026thinsp;31,684; accessed via IEU OpenGWAS) and used as exposure data, while cervical cancer genome‑wide association study (GWAS) summary statistics were obtained from the FinnGen consortium (dataset ID: finn‑b‑C3_CERVIX_UTERI) and used as outcome data. Finally, curated gene sets representing metal‑related regulated cell death pathways were compiled, including 512 ferroptosis‑related genes (drivers, suppressors, and markers) from FerrDb V2 and 43 cuproptosis‑related genes reported in the literature.\u003c/p\u003e \u003cp\u003eConstruction of the Metal Death Activity Score (MDAS)\u003c/p\u003e \u003cp\u003eBased on the TCGA cervical cancer (CESC) cohort, gene set variation analysis (GSVA) was applied to quantify ferroptosis (FAS) and cuproptosis (CAS) pathway activities for each sample. To address scale heterogeneity and construct an integrated metric, FAS and CAS were independently standardized using Z‑score transformation. The standardized values were then summed to define the Metal Death Activity Score (MDAS\u0026thinsp;=\u0026thinsp;Z‑FAS\u0026thinsp;+\u0026thinsp;Z‑CAS). An integrated waterfall plot was used to systematically visualize the distribution of MDAS and the relative contributions of its components.\u003c/p\u003e \u003cp\u003eFor stratification analyses, in addition to dichotomizing patients into high‑ and low‑MDAS groups based on the median MDAS, patients were further subdivided into four subtypes according to the independent median cutoffs of FAS and CAS (dual‑high, dual‑low, and discordant groups). The prognostic value of these stratification strategies was evaluated using Kaplan\u0026ndash;Meier survival analysis and the log‑rank test. Finally, linear models implemented in the limma package were used to identify differentially expressed genes (DEGs) between high‑ and low‑MDAS groups as well as among specific quartile‑based subtypes, thereby elucidating potential underlying molecular mechanisms.\u003c/p\u003e \u003cp\u003eWeighted Gene Co‑expression Network Analysis (WGCNA)\u003c/p\u003e \u003cp\u003eWGCNA was performed to identify gene modules associated with MDAS and to explore transcriptomic regulatory networks driving metal death activity. First, the top 5,000 genes with the highest variance were selected to reduce noise and enhance analytical robustness, followed by sample clustering to detect and remove outliers. The optimal soft‑thresholding power (β) was determined using the pickSoftThreshold function to construct a signed co‑expression network satisfying scale‑free topology criteria. The topological overlap matrix (TOM) was then calculated, and gene modules were identified using dynamic tree cutting (minModuleSize\u0026thinsp;=\u0026thinsp;30, mergeCutHeight\u0026thinsp;=\u0026thinsp;0.25).\u003c/p\u003e \u003cp\u003eModule eigengenes (MEs) were computed and subjected to Pearson correlation analysis to assess their relationships with MDAS and other clinical traits; modules with P\u0026thinsp;\u0026lt;\u0026thinsp;0.05 were considered significant. Finally, genes with an absolute module membership (MM)\u0026thinsp;\u0026gt;\u0026thinsp;0.8 and an absolute gene significance (GS)\u0026thinsp;\u0026gt;\u0026thinsp;0.2 were defined as MDAS‑associated hub genes.\u003c/p\u003e \u003cp\u003eSingle‑cell RNA‑seq Analysis\u003c/p\u003e \u003cp\u003eCell Clustering and Hierarchical Annotation\u003c/p\u003e \u003cp\u003eTo accurately identify cellular subpopulations within the tumor microenvironment, a two‑level hierarchical manual annotation strategy based on canonical marker gene expression was employed. At the first level, cell clusters were classified into three major lineages: immune cells (PTPRC), epithelial cells (EPCAM, KRT8, KRT5), and stromal cells (COL1A1, PECAM1). At the second level, more refined annotations were performed based on lineage‑specific molecular features:\u003c/p\u003e \u003cp\u003eEpithelial lineage: subdivided into squamous epithelial cells (KRT5, KRT14, TP63), glandular/columnar epithelial cells (KRT18, KRT19), and proliferative tumor cells (MKI67, TOP2A). Notably, high expression of CDKN2A (p16) was used to identify cell populations with HPV‑associated features.\u003c/p\u003e \u003cp\u003eLymphoid lineage: T cells were identified by CD3D/CD3E and further subdivided into CD4⁺ na\u0026iuml;ve/memory T cells (IL7R), CD8⁺ cytotoxic T cells (GZMB, GNLY), and CD8⁺ exhausted T cells (PDCD1, HAVCR2, LAG3). NK cells were identified by NKG7 and TYROBP, while B cells and plasma cells were marked by MS4A1 and JCHAIN/MZB1, respectively.\u003c/p\u003e \u003cp\u003eMyeloid lineage: including monocytes (CD14, VCAN), macrophages (CD68, CD163, C1QA), and dendritic cells (CD1C, CLEC9A).\u003c/p\u003e \u003cp\u003eStromal lineage: comprising fibroblasts (DCN, LUM) and endothelial cells (VWF, PECAM1).\u003c/p\u003e \u003cp\u003eCalculation of MDAS Using UCell\u003c/p\u003e \u003cp\u003eTo robustly assess ferroptosis and cuproptosis pathway activities at the single‑cell level, the UCell R package was employed. UCell computes gene set enrichment scores based on the Mann\u0026ndash;Whitney U statistic, relying on relative gene expression ranks rather than absolute expression values, thereby ensuring robustness to dropout events and batch effects inherent to scRNA‑seq data. Ferroptosis (FAS) and cuproptosis (CAS) gene sets were input into the AddModuleScore_UCell function to calculate initial scores. To ensure comparability, raw UCell scores were Z‑score standardized. The MDAS for each cell was defined as the sum of the standardized scores: MDAS\u0026thinsp;=\u0026thinsp;Z(FAS_UCell)\u0026thinsp;+\u0026thinsp;Z(CAS_UCell).\u003c/p\u003e \u003cp\u003ePseudotime Trajectory Analysis\u003c/p\u003e \u003cp\u003eTo reconstruct epithelial cell developmental trajectories during cervical cancer progression, the Monocle3 R package was used. Epithelial subsets (including normal squamous epithelial cells, precancerous/HSIL cells, and tumor cells) were extracted from the integrated Seurat object and converted into a cell_data_set object. Data were preprocessed using principal component analysis (PCA), followed by batch‑effect correction based on sample IDs (orig.ident), and then visualized by UMAP for dimensionality reduction. The principal graph was constructed using the learn_graph function with use_partition\u0026thinsp;=\u0026thinsp;FALSE to enforce a single continuous trajectory.\u003c/p\u003e \u003cp\u003eThe node with the highest proportion of \u0026ldquo;normal epithelial\u0026rdquo; cells was defined as the root, and pseudotime values were computed for each cell. Biological relevance was validated by assessing the correspondence between pseudotime and disease stage. Spearman correlation analysis was then used to quantify associations between pseudotime and MDAS, FAS, and CAS. Finally, genes exhibiting significant expression changes along the trajectory were identified using the graph_test function based on Moran\u0026rsquo;s I statistic (q\u0026thinsp;\u0026lt;\u0026thinsp;0.05), with particular emphasis on the dynamic expression patterns of hub genes.\u003c/p\u003e \u003cp\u003eMendelian Randomization\u0026ndash;Based Causal Inference\u003c/p\u003e \u003cp\u003eTo further identify reliable therapeutic targets and validate potential causal associations between the identified hub genes and cervical cancer risk, a two‑sample Mendelian randomization (MR) analysis was performed. High‑confidence hub genes identified by WGCNA and scRNA‑seq analyses were used as exposure factors. Gene symbols were converted to Ensembl IDs using the biomaRt package. Cis‑expression quantitative trait loci (cis‑eQTLs) significantly associated with these genes were extracted from the IEU OpenGWAS database and used as IVs. The selection threshold was set at genome‑wide significance (P\u0026thinsp;\u0026lt;\u0026thinsp;5 \u0026times; 10⁻⁸). To ensure independence, linkage disequilibrium (LD) pruning was applied with r\u0026sup2; \u0026lt; 0.001 within a 10,000‑kb window to remove correlated SNPs. Genetic summary statistics for cervical cancer were obtained from the FinnGen consortium (ID: finn‑b‑C3_CERVIX_UTERI). This dataset is based on a European population, consistent with the exposure data, thereby minimizing population stratification bias. All MR analyses were conducted using the TwoSampleMR package in R. Exposure and outcome datasets were harmonized to align effect alleles. The inverse‑variance weighted (IVW) method was used as the primary approach to estimate causal effects in terms of odds ratios (ORs). Complementary methods, including MR‑Egger regression, the weighted median, and the weighted mode, were applied to assess robustness. Cochran\u0026rsquo;s Q test was used to evaluate heterogeneity among IVs. Horizontal pleiotropy was assessed using the MR‑Egger intercept test (P\u0026thinsp;\u0026gt;\u0026thinsp;0.05 indicating no significant pleiotropy). P values were adjusted for multiple testing using the Benjamini\u0026ndash;Hochberg false discovery rate (FDR) method. Genes with P\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and consistent effect directions across multiple MR methods were considered potential causal genes.\u003c/p\u003e \u003cp\u003eMachine Learning\u0026ndash;Based Construction of the Metal Death Risk Score (MDRS)\u003c/p\u003e \u003cp\u003eBased on the eight metal death\u0026ndash;related causal genes identified by MR (e.g., DACT1, CALD1), prognostic models were constructed using an integrated framework comprising 74 machine‑learning algorithms derived from 10 base methods, including random survival forests (RSF), LASSO, Ridge, Elastic Net (Enet), Stepwise Cox (StepCox), CoxBoost, partial least squares regression for Cox models (plsRcox), and gradient boosting machines (GBM). Ten‑fold cross‑validation was performed, and the concordance index (C‑index) was calculated to select the optimal algorithm combination with the highest mean C‑index for construction of the Metal Death Risk Score (MDRS). Patients were stratified into high‑ and low‑risk groups based on the median MDRS. Kaplan\u0026ndash;Meier survival curves and time‑dependent receiver operating characteristic (ROC) curves were used to evaluate model performance. To assess clinical applicability, a nomogram incorporating MDRS, age, clinical stage, tumor grade, and lymph node status was developed. Calibration curves and decision curve analysis (DCA) were used to evaluate the accuracy and net clinical benefit of the nomogram. In addition, net reclassification improvement (NRI) and integrated discrimination improvement (IDI) were calculated to quantify the incremental prognostic value of MDRS over traditional clinical features.\u003c/p\u003e \u003cp\u003eBiological Characteristics and Immune Microenvironment Associated with MDRS Stratification\u003c/p\u003e \u003cp\u003eTo elucidate the biological mechanisms underlying MDRS‑based stratification, Hallmark gene set activities were first evaluated using GSVA. Linear modeling with limma identified pathways that were significantly different between risk groups (FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05). For tumor microenvironment (TME) characterization, the ssGSEA algorithm was applied using ESTIMATE gene sets to calculate stromal and immune scores, and tumor purity was inferred using the Yoshihara formula. Differences in these scores between risk groups were assessed using the Wilcoxon rank‑sum test. In addition, the expression levels of key immune checkpoint genes (e.g., PDCD1, CTLA4, CD274) were compared between groups. Spearman correlation analysis was performed to explore potential regulatory relationships between MDRS model genes and key biological pathways.\u003c/p\u003e \u003cp\u003eNetwork Pharmacology Analysis and Targeted Drug Prediction for MDRS\u003c/p\u003e \u003cp\u003eTo explore therapeutic strategies targeting MDRS‑associated genes, drug\u0026ndash;target interaction data were integrated from DGIdb, DrugBank, and the Comparative Toxicogenomics Database (CTD) to construct a high‑confidence drug screening library. Precision therapeutic strategies were designed for the eight core causal genes (e.g., PDGFRB, DACT1), aiming to \u0026ldquo;inhibit risk genes\u0026rdquo; and \u0026ldquo;activate protective genes.\u0026rdquo; Identified candidate drugs were cross‑referenced with known ferroptosis and cuproptosis inducers (e.g., sorafenib, erastin) to identify agents with dual mechanisms of \u0026ldquo;targeting MDRS genes\u0026rdquo; and \u0026ldquo;inducing metal‑dependent cell death.\u0026rdquo; Sankey diagrams were generated using the ggalluvial package to visualize the logical flow of \u0026ldquo;gene\u0026ndash;drug\u0026ndash;therapeutic strategy,\u0026rdquo; and drug\u0026ndash;target interaction networks were constructed using the ggraph package with a Kamada\u0026ndash;Kawai layout.\u003c/p\u003e \u003cp\u003eMolecular Docking\u003c/p\u003e \u003cp\u003eTo predict the binding affinity and interaction mode between sorafenib and the core target PDGFRB, molecular docking was performed using the CB‑Dock2 online server. PDGFRB comprises extracellular, transmembrane, and intracellular kinase domains, and sorafenib, as a tyrosine kinase inhibitor (TKI), exerts its inhibitory effect by binding to the intracellular kinase domain. Available crystal structures of PDGFRB in the RCSB PDB (e.g., 3MJG, 2L6W) were therefore considered unsuitable because they lack critical kinase regions. Consequently, a high‑confidence full‑length predicted structure of human PDGFRB was obtained from the AlphaFold Protein Structure Database (UniProt ID: P09619; model: AF‑P09619‑F1). The three‑dimensional structure of sorafenib was downloaded from PubChem (CID: 216239).\u003c/p\u003e \u003cp\u003eDocking was conducted using the \u0026ldquo;Auto Blind Docking\u0026rdquo; mode in CB‑Dock2. The platform automatically preprocessed the protein structure and identified potential binding pockets based on surface curvature using the CurPocket module. Autodock Vina was then employed to perform docking calculations within the identified pockets. The conformation with the lowest Vina score (indicating the strongest predicted affinity) and a binding site consistent with kinase inhibitor characteristics was selected as the optimal binding mode for visualization and downstream analyses.\u003c/p\u003e \u003cp\u003eMolecular Dynamics Simulation and MM/PBSA Binding Free Energy Calculation\u003c/p\u003e \u003cp\u003eTo quantitatively characterize the inhibitory potency of sorafenib toward PDGFRB at atomic resolution and to provide thermodynamic constraints for subsequent structure\u0026ndash;network coupled simulations, all‑atom molecular dynamics (MD) simulations were performed, followed by MM/PBSA binding free energy calculations based on equilibrated trajectories. MD simulations of the PDGFRB\u0026ndash;sorafenib complex were carried out using GROMACS (v2022.4). The protein structure was preprocessed by removing crystallographic water molecules and completing missing atoms. The protein was parameterized using the Amber99SB force field, while ligand topologies were generated with the GAFF force field. The system was solvated in a TIP3P water box with a minimum distance of 1.0 nm between the solute and the box boundary, and Na⁺/Cl⁻ ions were added to neutralize the system and maintain physiological ionic strength (0.15 M NaCl).\u003c/p\u003e \u003cp\u003eAfter energy minimization using the steepest descent algorithm, the system was equilibrated under the NVT ensemble (300 K, 100 ps) and the NPT ensemble (1 atm, 100 ps), followed by a 10 ns production run with a time step of 2 fs. Bond lengths involving hydrogen atoms were constrained using the LINCS algorithm, and long‑range electrostatic interactions were treated using the particle mesh Ewald (PME) method. Subsequently, MM/PBSA binding free energy calculations were performed using gmx_MMPBSA (v1.6.3) on 301 frames extracted from the equilibrated portion of the trajectory (internal dielectric constant\u0026thinsp;=\u0026thinsp;2.0, external dielectric constant\u0026thinsp;=\u0026thinsp;80.0, ionic strength\u0026thinsp;=\u0026thinsp;0.15 M). Total binding free energies under both PB and GB models, along with individual energy components (ΔE\u003csub\u003evdW\u003c/sub\u003e, ΔE\u003csub\u003eele\u003c/sub\u003e, ΔG\u003csub\u003esolv\u003c/sub\u003e), were obtained and used as thermodynamic driving parameters for downstream network propagation modeling.\u003c/p\u003e \u003cp\u003eSurface Plasmon Resonance (SPR) Analysis\u003c/p\u003e \u003cp\u003eThe binding affinity between sorafenib and PDGFRB was measured using surface plasmon resonance on a Biacore 1K system (Cytiva, Sweden). The recombinant human PDGFRB protein (0.2 mg/mL in PBS) was immobilized on a Series S Sensor Chip CM5 (Cytiva, Cat. No. 29149603) using the Amine Coupling Kit (Cytiva, Cat. No. BR100050). Briefly, the chip surface was activated by injecting a 1:1 mixture of N-hydroxysuccinimide (NHS) and 1-ethyl-3-(3-dimethylaminopropyl) carbodiimide (EDC) at a flow rate of 10 \u0026micro;L/min. PDGFRB protein was diluted to 40 \u0026micro;g/mL in 10 mM sodium acetate buffer (pH 4.0) and injected at 10 \u0026micro;L/min for immobilization, achieving a final immobilization level of approximately 13,000 response units (RU). Remaining active sites were blocked with\u003c/p\u003e \u003cp\u003eethanolamine-HCl.\u003c/p\u003e \u003cp\u003eSorafenib (10 mM stock in DMSO) was serially diluted in PBST running buffer to final concentrations of 0.003, 0.01, 0.04, 0.12, 0.37, 1.1, and 3.3 \u0026micro;M. Binding experiments were performed at 25\u0026deg;C using a flow path of 6\u0026thinsp;\u0026minus;\u0026thinsp;5 at a flow rate of 30 \u0026micro;L/min. Each concentration was injected for 90 seconds (association phase), followed by buffer flow for 120 seconds (dissociation phase). Three startup cycles were performed prior to sample injection. Reference-subtracted sensorgrams were analyzed using Biacore Insight Evaluation Software. The equilibrium dissociation constant (KD) was determined by fitting the steady-state response values to a 1:1 Langmuir binding model.\u003c/p\u003e \u003cp\u003eStructure‑Informed Network Propagation (SINP) Model\u003c/p\u003e \u003cp\u003eTo bridge the gap between \u0026ldquo;atomic‑scale drug\u0026ndash;target binding\u0026rdquo; and \u0026ldquo;tumor tissue‑scale phenotypic alterations,\u0026rdquo; we developed a structure‑context\u0026ndash;coupled network propagation model termed Structure‑Informed Network Propagation (SINP). This framework integrates thermodynamic parameters derived from MD/MM‑PBSA analyses, transcriptomic co‑regulation relationships (WGCNA), and protein\u0026ndash;protein interaction (PPI) networks to enable systems pharmacology simulations with explicit physicochemical constraints.\u003c/p\u003e \u003cp\u003eUnlike conventional network models that treat target perturbation as a binary on/off event, SINP converts the sorafenib\u0026ndash;PDGFRB binding free energy (ΔG) into a quantitative perturbation strength reflecting inhibitory efficiency. Perturbation efficiency was defined as Ekd\u0026thinsp;=\u0026thinsp;1/(1\u0026thinsp;+\u0026thinsp;e^0.2(ΔG\u0026thinsp;+\u0026thinsp;25)), and was used as the initial perturbation coefficient for network propagation.\u003c/p\u003e \u003cp\u003eAt the network topology level, to enhance tissue specificity, a generic PPI network was projected onto TCGA‑CESC transcriptomic data, and interaction edges with weak co‑expression were pruned: gene pairs with low correlation (|r| \u0026lt; 0.15) in the cohort were considered inactive in the cervical cancer context and removed. In addition, WGCNA module information was incorporated to dynamically weight propagation: if two genes belonged to the same co‑expression module, their propagation weight was amplified 1.5‑fold to reinforce functionally coupled signal transmission. Network propagation was then initiated from PDGFRB in the cervical cancer\u0026ndash;specific weighted network, yielding a perturbation score for each gene along with its predicted direction of regulation, thereby enabling prediction of systemic drug effects at the tissue scale.\u003c/p\u003e \u003cp\u003eSingle‑cell In Silico Pharmacodynamic Modeling (scTenifoldKnk)\u003c/p\u003e \u003cp\u003eTo validate the cellular sources of PDGFRB‑mediated effects and delineate downstream regulatory networks at single‑cell resolution, single‑cell in silico pharmacodynamic modeling was performed using scTenifoldKnk. This approach evaluates the impact of PDGFRB inhibition on transcriptional regulatory networks across different cell subpopulations.\u003c/p\u003e \u003cp\u003eBased on cell type\u0026ndash;specific expression profiles of PDGFRB in the scRNA‑seq dataset, monocytes\u0026mdash;exhibiting high PDGFRB expression\u0026mdash;were selected as the primary effector cell population for simulation. Proliferative epithelial tumor cells (Epi_Tumor_Prolif) with moderate PDGFRB expression were subjected to the same simulation as a comparative control to assess cell context dependency. Within each target cell population, single‑cell gene regulatory networks (scGRNs) were constructed, and PDGFRB inhibition was simulated by setting its outgoing edge weights to zero. Virtual perturbation scores (Z‑scores) were calculated based on changes in gene distances within the manifold space before and after perturbation. Genes with |Z| \u0026gt; 1.96 and FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05 were defined as significantly affected downstream response genes of PDGFRB inhibition and were subsequently used for functional enrichment analyses and cross‑scale consistency validation with SINP predictions.\u003c/p\u003e \u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eStatistical Analysis\u003c/h2\u003e \u003cp\u003eAll statistical analyses were performed using R (v4.5.0). Continuous variables were compared using the Wilcoxon rank‑sum test or the Kruskal\u0026ndash;Wallis test, as appropriate. Correlations were assessed using Spearman\u0026rsquo;s method. Survival analyses were conducted using Kaplan\u0026ndash;Meier curves with the log‑rank test, as well as Cox proportional hazards models. Multiple testing correction was performed using the Benjamini\u0026ndash;Hochberg method. A two‑sided P value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 was considered statistically significant.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cp\u003e1. Construction and validation of the metal-dependent death activity score (MDAS)\u003c/p\u003e\n\u003cp\u003eTo quantify the activity levels of ferroptosis and cuproptosis in patients with cervical cancer, we constructed a metal-dependent death activity score (MDAS) system based on the TCGA-CESC cohort (n = 304) and normal cervical tissues from GTEx (n = 70). First, the ferroptosis activity score (FAS) and cuproptosis activity score (CAS) were calculated using the ssGSEA algorithm. Tumor-specific analyses showed that both FAS and CAS were significantly elevated in cervical cancer tissues compared with normal tissues (both P \u0026lt; 0.05; Figure 1), indicating aberrant activation of ferroptosis- and cuproptosis-related pathways in cervical cancer. Further correlation analysis revealed only a very weak association between FAS and CAS (Pearson r = 0.053, P = 0.355; Figure 1), suggesting that these two forms of metal-dependent cell death represent relatively independent biological processes in cervical cancer.\u003c/p\u003e\n\u003cp\u003eBased on these findings, the standardized FAS and CAS were integrated into a unified metal-dependent death activity score (MDAS = FAS_z + CAS_z). Quartile-based stratification showed comparable proportions of patients in the high-iron/high-copper, high-iron/low-copper, low-iron/high-copper, and low-iron/low-copper groups (23.7%, 26.3%, 26.3%, and 23.7%, respectively; Figure 1). Notably, the \u0026ldquo;discordant groups\u0026rdquo; (high iron/low copper plus low iron/high copper) accounted for 52.6% of the cohort, further supporting the necessity of an integrated analysis. Cox regression analysis was performed to evaluate the prognostic value of FAS and CAS for overall survival (OS). The results showed that CAS was an independent prognostic risk factor (HR = 5.686, 95% CI: 1.009\u0026ndash;32.038, P = 0.0488), whereas FAS did not exhibit a significant association with prognosis (HR = 1.368, P = 0.747). These findings suggest that cuproptosis may play a more critical role in the prognosis of cervical cancer and validate the rationale for integrating FAS and CAS to construct the MDAS. Based on the median MDAS value, patients were divided into a high-MDAS group (n = 152) and a low-MDAS group (n = 152) for subsequent multi-omics integrative analyses and classification of metal metabolism\u0026ndash;related phenotypes.\u003c/p\u003e\n\u003cp\u003e2. Identification of MDAS-related gene co-expression modules by WGCNA\u003c/p\u003e\n\u003cp\u003eTo identify gene co-expression networks associated with MDAS, WGCNA was performed on transcriptomic data from 304 TCGA-CESC tumor samples. After data preprocessing and removal of low-variance genes, a total of 12,455 genes were included in the subsequent analysis. Soft-threshold power analysis indicated that a power of 9 satisfied the scale-free topology criterion (scale-free R\u0026sup2; = 0.915; Figure 1). Using this parameter, a co-expression network was constructed, identifying 13 gene modules (excluding the gray module), with module sizes ranging from 44 genes (salmon module) to 2,352 genes (turquoise module) (Figure 1).\u003c/p\u003e\n\u003cp\u003eModule\u0026ndash;trait relationship analysis revealed five gene modules significantly associated with MDAS (P \u0026lt; 0.05; Figure 1). Among them, the magenta module showed the strongest negative correlation with MDAS (r = \u0026minus;0.248, P = 1.19 \u0026times; 10^-5; 247 genes), whereas the red module (r = 0.213, P = 1.85 \u0026times; 10^-4; 540 genes) and the pink module (r = 0.211, P = 2.08 \u0026times; 10^-4; 462 genes) were positively correlated with MDAS. The turquoise module (r = \u0026minus;0.146, P = 0.011; 2,352 genes) and the brown module (r = 0.115, P = 0.046; 1,127 genes) also reached statistical significance. In addition, hub genes within each module were identified based on module membership (MM) and gene significance (GS). Using stringent criteria (|MM| \u0026gt; 0.8 and |GS| \u0026gt; 0.2) or relaxed criteria (|MM| \u0026gt; 0.6 and |GS| \u0026gt; 0.15), a total of 293 hub genes were identified (Figure 1). These genes were selected as candidate genes for subsequent Mendelian randomization analyses and risk score model construction. Notably, the magenta and turquoise modules (negatively correlated with MDAS) may represent molecular features associated with suppression of metal-dependent cell death, whereas the red, pink, and brown modules (positively correlated with MDAS) may reflect biological processes that promote metal-dependent cell death, providing important clues for subsequent functional enrichment analyses.\u003c/p\u003e\n\u003cp\u003e3. Single-cell analysis\u003c/p\u003e\n\u003cp\u003e3.1. Single-cell landscape of the metal-dependent death\u0026ndash;related prognostic model\u003c/p\u003e\n\u003cp\u003eTo elucidate the biological functions of the constructed metal-dependent death activity score (MDAS) in cervical cancer at the single-cell level, we performed an integrated analysis of single-cell transcriptomic data derived from 22,141 cancer tissue cells, 12,333 high-grade squamous intraepithelial lesion (HSIL) cells, and 39,980 normal tissue cells. Through unsupervised clustering combined with marker gene\u0026ndash;based annotation, we identified 13 major cell types, including lymphoid lineages (T cells, B cells, and NK cells), myeloid lineages (monocytes and macrophages), stromal lineages (fibroblasts and endothelial cells), and epithelial lineages (normal squamous epithelial cells and tumor epithelial cells) (Figure 2). We first assessed MDAS scores at a global level. The results showed that MDAS scores were significantly higher in HSIL and cancer tissues than in normal tissues, indicating a progressive activation of metal-dependent death pathways during cervical cancer development (Figure 2). At the lineage level, epithelial and myeloid cells exhibited the highest MDAS scores, whereas lymphocytes showed relatively low scores (Figure 2), highlighting marked heterogeneity in MDAS activation across different cellular lineages.\u003c/p\u003e\n\u003cp\u003e3.2. Activation patterns of MDAS in specific cell subpopulations\u003c/p\u003e\n\u003cp\u003eTo more precisely localize the cells responsible for MDAS activation, we ranked and visualized MDAS scores across the 13 identified cell subpopulations. Violin plots showed that CD8+T cells, B cells, and CD4+T cells had the lowest MDAS scores, whereas macrophages, proliferative tumor epithelial cells (Epi_Tumor_Prolif), and tumor epithelial cells (Epithelial_Tumor) exhibited the highest scores (Figure 2). These findings were further supported by UMAP visualizations, in which MDAS-high cells were predominantly enriched in T-cell clusters as well as in subsets of myeloid and epithelial clusters (Figure 2). Percentage-stacked bar plots clearly demonstrated that more than 60% of cells within T-cell subpopulations belonged to the MDAS-high phenotype, while this proportion approached 40% in macrophages (Figure 2). Importantly, we employed split violin plots to investigate changes in MDAS within specific cell types across disease progression. Remarkably, although the overall MDAS scores in T cells were not high, macrophages and tumor epithelial cells displayed a pronounced and progressive increase in MDAS scores from normal tissue to cancer (Figure 4D). These results strongly suggest that macrophages and tumor cells within the tumor microenvironment are the central cellular players mediating MDAS-related functions and may actively drive disease progression.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eLineage Characterization of MDAS.\u003c/strong\u003e (A) UMAP dimensionality reduction clustering plot of the single-cell atlas. Top left: Colored by sample source (Group: GTEx, HSIL, Normal, TCGA). Bottom left: Colored by disease status (Disease). Top right: Colored by single-cell subpopulation annotation (sng.ident), displaying 13 major cell types. (B) Heatmap showing the expression of the top 5 marker genes across cell subpopulations. The horizontal axis represents different cell subpopulations, and the vertical axis represents marker genes; yellow indicates high expression. (C) Dot plot of Level 1 Broad Lineage Markers. Illustrates the expression patterns of canonical markers for immune cells (PTPRC), epithelial cells (EPCAM), and fibroblasts (COL1A1), etc. Dot size represents the percentage of expressing cells, and color intensity represents average expression level. (D) Dot plot of Level 2 Cervical Cancer Specific Markers. Further refines the specific molecular features of each cell subpopulation. (E) Violin plot of MDAS distribution across four broad cell lineages (Epithelial, Lymphoid, Myeloid, Stromal). Results show higher MDAS scores in Epithelial and Myeloid lineages. (F) Violin plot of MDAS distribution across detailed cell subpopulations (ranked by median score). Tumor epithelial cells (Epithelial_Tumor) and Macrophages exhibit the highest metal death activity, while CD8+ T cells and B cells show the lowest scores. (G) Split violin plots showing MDAS changes within cell subpopulations across disease progression. Compares score differences among Normal, High-Grade Squamous Intraepithelial Lesion (HSIL), and Cancer states, revealing a significant increase in MDAS with malignancy in Macrophages and Tumor Epithelial cells. (H) Feature plots projecting Ferroptosis Activity Score (FAS_sc), Cuproptosis Activity Score (CAS_sc), and the integrated score (MDAS_sc) onto the single-cell UMAP. Red areas indicate high scores, visually demonstrating the enrichment of metal death activity in specific cell clusters.\u003c/p\u003e\n\u003cp\u003e3.3. scRNA-seq validation of prognostic hub genes\u003c/p\u003e\n\u003cp\u003eNext, we validated the expression patterns of key hub genes identified by WGCNA at the single-cell level. We first calculated WGCNA module scores and found that the Brown and Pink modules were expressed across multiple cell types, whereas the Turquoise module showed high expression predominantly in myeloid and epithelial cells (Figure 3). Representative top hub genes from different modules were selected for visualization. Annotated heatmaps clearly illustrated the modular expression patterns: MDK and LCN2 from the Turquoise module were highly and specifically expressed in macrophages and tumor epithelial cells, while COL4A1 and THBS2 from the Brown module were mainly expressed in fibroblasts (Figure 3). DotPlots further confirmed these findings, with prominent high-expression \u0026ldquo;red dots\u0026rdquo; precisely mapping hub genes to specific cell types (Figure 3). Finally, UMAP plots intuitively visualized expression hotspots of key genes (e.g., MDK and LCN2), which showed strong spatial overlap with macrophage and tumor cell clusters (Figure 3). Together, these single-cell\u0026ndash;level results not only cross-validated the robustness of our WGCNA analysis but also reaffirmed macrophages and tumor cells as the key cellular subpopulations underpinning MDAS-associated functions.\u003c/p\u003e\n\u003cp\u003e3.4. Pseudotime analysis reveals cell type\u0026ndash;specific dynamic evolution of MDAS during malignant transformation\u003c/p\u003e\n\u003cp\u003eUsing the Monocle3 algorithm, we analyzed epithelial cell subpopulations and successfully reconstructed a continuous trajectory of malignant transformation from normal epithelium through precancerous lesions (HSIL) to tumor cells. The distribution of pseudotime values strictly followed the order of disease progression (P \u0026lt; 0.0001; Figure 4). Along this evolutionary trajectory, MDAS scores exhibited a significant negative correlation with pseudotime (Spearman \u0026rho; = \u0026minus;0.231, P = 3.5 \u0026times; 10^-44), with a more pronounced decline in ferroptosis activity (FAS) compared with cuproptosis activity (CAS). This suggests that suppression of metal-dependent death pathways\u0026mdash;particularly ferroptosis\u0026mdash;may be a critical mechanism by which epithelial cells acquire a survival advantage during malignant transformation (Figure 4). In contrast to the dynamic changes observed in epithelial cells, pseudotime analysis of myeloid cells clearly captured the differentiation process from monocytes to macrophages; however, MDAS scores remained relatively stable along this differentiation trajectory (P = 0.185). This indicates that the elevated MDAS observed in the myeloid lineage is primarily driven by disease state rather than by the differentiation process itself (Figure 4). In addition, dynamic expression changes of hub genes such as CLSPN and DDR2 along pseudotime further supported these findings, collectively delineating a lineage-specific and spatiotemporal evolutionary landscape of MDAS across different cell lineages within the cervical cancer microenvironment (Figure 4).\u003c/p\u003e\n\u003cp\u003e4. Identification of causal genes for cervical cancer by Mendelian randomization\u003c/p\u003e\n\u003cp\u003eTo establish causal relationships between MDAS-related hub genes and cervical cancer susceptibility, we performed a two-sample Mendelian randomization (MR) analysis. Blood-derived cis-eQTLs from the eQTLGen Consortium were used as instrumental variables, and cervical cancer GWAS summary statistics from FinnGen were used as the outcome data. Among the 211 high-confidence hub genes identified by WGCNA, 204 were successfully mapped to Ensembl IDs, of which 123 had available cis-eQTL instruments (P \u0026lt; 5 \u0026times; 10^-8; LD clumping: r^2\u0026lt; 0.001, window = 10,000 kb). In total, 361 independent SNPs were used as instrumental variables. Inverse variance\u0026ndash;weighted (IVW) MR analysis identified eight genes with nominally significant causal associations with cervical cancer risk (P \u0026lt; 0.05) (Figure 5). Among these, four genes were risk factors (OR \u0026gt; 1): DACT1 showed the strongest risk effect (OR = 2.66, 95% CI: 1.50\u0026ndash;4.73, P = 0.001), followed by PDGFRB (OR = 1.66), PRSS23 (OR = 1.44), and MYO15B (OR = 1.13). The remaining four genes were protective factors (OR \u0026lt; 1): MSRB3 exhibited the strongest protective effect (OR = 0.48, 95% CI: 0.27\u0026ndash;0.83, P = 0.009), followed by CALD1 (OR = 0.84), DAB2 (OR = 0.72), and BNC2 (OR = 0.78). Sensitivity analyses confirmed the robustness of these findings. No significant horizontal pleiotropy was detected for any of the causal genes (MR-Egger intercept P \u0026gt; 0.05), and no significant heterogeneity was observed (Cochran\u0026rsquo;s Q P \u0026gt; 0.05), indicating that the observed causal effects were not driven by invalid instruments or pleiotropic bias. Notably, DACT1, a negative regulator of the Wnt signaling pathway, exhibited the strongest causal association with cervical cancer risk, suggesting that Wnt dysregulation may play a critical role in tumorigenesis. PDGFRB, which encodes platelet-derived growth factor receptor beta, is involved in angiogenesis and stromal remodeling\u0026mdash;hallmarks of tumor progression. Among the protective genes, MSRB3, a methionine sulfoxide reductase, participates in defense against oxidative stress, whereas DAB2 is a well-established tumor suppressor involved in the TGF-\u0026beta; signaling pathway.\u003c/p\u003e\n\u003cp\u003e5. Construction of a metal-dependent death risk score (MDRS) using 74 machine learning algorithms\u003c/p\u003e\n\u003cp\u003eTo develop a robust prognostic risk score based on the eight MR-validated causal genes, we systematically evaluated 74 combinations of machine learning algorithms using a 10-fold cross-validation framework. The algorithm pool comprised 10 core methods, including LASSO-Cox, Ridge-Cox, Elastic Net, Stepwise Cox, random survival forest (RSF), gradient boosting machine (GBM), CoxBoost, and their parameter variants (Figure 6). Among all tested algorithms, Elastic Net (\u0026alpha; = 0.5) combined with lambda.min (ElasticNet_a5_lambdamin) achieved the highest cross-validated C-index (0.627) and was therefore selected as the optimal algorithm for MDRS construction. The final MDRS model incorporated all eight causal genes with optimized coefficient weights to calculate individualized risk scores for each patient.\u003c/p\u003e\n\u003cp\u003e5.1. Prognostic performance and risk stratification of MDRS\u003c/p\u003e\n\u003cp\u003eUsing the median MDRS value as the cutoff, patients were stratified into high-risk (n = 146) and low-risk (n = 145) groups. Kaplan\u0026ndash;Meier survival analysis demonstrated that patients in the high-risk group had significantly worse overall survival (OS) than those in the low-risk group (log-rank P = 0.0032; Figure 6). Time-dependent ROC analysis showed moderate predictive accuracy, with AUC values of 0.637, 0.647, 0.642, and 0.601 for 1-, 2-, 3-, and 5-year OS, respectively (Figure 6). The distribution of MDRS across patients exhibited a clear gradient, with higher risk scores associated with shorter survival times and increased mortality (Figure 6). These results confirm that MDRS effectively captures prognostic heterogeneity driven by expression patterns of metal-dependent death\u0026ndash;related genes.\u003c/p\u003e\n\u003cp\u003e5.2. Integration of MDRS with clinicopathological features: nomogram construction\u003c/p\u003e\n\u003cp\u003eGiven that the standalone MDRS model (C-index = 0.627) still had room for improvement, we investigated whether integrating MDRS with established clinicopathological prognostic factors could enhance predictive performance. A multivariable Cox proportional hazards model was constructed incorporating MDRS together with age, FIGO stage, histological grade, and lymph node status (N stage). The baseline model including only clinical features (age + stage + grade + N stage) yielded a C-index of 0.664. Notably, the addition of MDRS significantly improved model performance, increasing the C-index to 0.736 (\u0026Delta;C-index = +0.071). This indicates that MDRS provides substantial independent prognostic information beyond traditional clinical parameters, including lymph node metastasis, a strong prognostic factor. A nomogram was constructed to visualize the integrated model and facilitate clinical application (Figure 6). Calibration curves demonstrated good agreement between predicted and observed 3-year survival probabilities, with data points closely aligned along the 45\u0026deg; reference line (Figure 6).\u003c/p\u003e\n\u003cp\u003e5.3. Incremental value assessment (NRI and IDI analyses)\u003c/p\u003e\n\u003cp\u003eTo quantify the prognostic improvement achieved by adding MDRS to the clinical baseline model, we calculated the net reclassification improvement (NRI) and integrated discrimination improvement (IDI) for 3-year survival. Incorporation of MDRS significantly improved risk reclassification, yielding a 3-year NRI of 0.073 (95% CI: 0.001\u0026ndash;0.203, P = 0.044). This indicates that a proportion of patients misclassified by the clinical model were correctly reclassified into appropriate risk categories after inclusion of MDRS. In addition, IDI analysis showed a positive trend toward improvement (IDI = 0.294, 95% CI: \u0026minus;0.063\u0026ndash;0.566, P = 0.106), further supporting the added clinical value of MDRS.\u003c/p\u003e\n\u003cp\u003e5.4. Decision curve analysis and subgroup validation\u003c/p\u003e\n\u003cp\u003eDecision curve analysis (DCA) was performed to evaluate the clinical utility of the integrated model across a range of threshold probabilities. The DCA demonstrated that, across most clinically relevant thresholds, the integrated model (clinical features + MDRS; red line) consistently provided greater net benefit than the clinical baseline model alone (blue line) (Figure 6), confirming its superior decision-support capability. Given the clinical challenge of prognostic stratification in early-stage cervical cancer, we further validated the performance of MDRS in patients with FIGO stage I\u0026ndash;II disease. Kaplan\u0026ndash;Meier analysis showed that MDRS retained significant discriminative ability even within this early-stage subgroup, effectively distinguishing high- and low-risk patients (log-rank P = 0.01; Figure 6). This finding highlights the potential utility of MDRS for early identification of high-risk patients who may benefit from intensified treatment or closer surveillance.\u003c/p\u003e\n\u003cp\u003e6. Biological characteristics of MDRS-based risk stratification\u003c/p\u003e\n\u003cp\u003eTo elucidate the molecular mechanisms underlying the prognostic capability of MDRS, gene set variation analysis (GSVA) was performed to compare pathway activity differences between the MDRS-defined high-risk and low-risk groups. This analysis focused on Hallmark gene sets representing key cancer biological processes, excluding ferroptosis- and cuproptosis-related signatures that had already been analyzed during MDAS construction.\u003c/p\u003e\n\u003cp\u003e6.1. Differential pathway activity between risk groups\u003c/p\u003e\n\u003cp\u003eGSVA identified 19 pathways with significantly different activity between the MDRS high- and low-risk groups (FDR \u0026lt; 0.05) (Figure 7). Specifically, 16 pathways were significantly upregulated and 3 pathways were downregulated in the high-risk group. The most prominently upregulated pathway in high-risk patients was epithelial\u0026ndash;mesenchymal transition (EMT) (logFC = 0.450, FDR = 4.03 \u0026times; 10^-39), followed by angiogenesis (logFC = 0.337, FDR = 8.79 \u0026times; 10^-26), Hedgehog signaling (logFC = 0.226, FDR = 7.66 \u0026times; 10^-17), and TGF-\u0026beta; signaling (logFC = 0.227, FDR = 2.56 \u0026times; 10^-13) (Figure 7). These findings indicate that MDRS high-risk tumors possess enhanced invasive and metastatic potential. In addition, multiple oncogenic signaling pathways were activated in high-risk patients, including Notch signaling (logFC = 0.159, FDR = 5.46 \u0026times; 10^-14), Wnt/\u0026beta;-catenin signaling (logFC = 0.143, FDR = 2.50 \u0026times; 10^-9), and KRAS signaling upregulation (logFC = 0.153, FDR = 4.17 \u0026times; 10^-9). Inflammatory-related pathways such as hypoxia (logFC = 0.095, FDR = 3.97 \u0026times; 10^-4), IL2\u0026ndash;STAT5 signaling (logFC = 0.096, FDR = 2.89 \u0026times; 10^-4), and TNF-\u0026alpha;/NF-\u0026kappa;B signaling (logFC = 0.084, FDR = 0.036) were also significantly elevated. In contrast, pathways associated with normal cellular metabolism were suppressed in high-risk patients, particularly oxidative phosphorylation (logFC = \u0026minus;0.201, FDR = 3.29 \u0026times; 10^-8), fatty acid metabolism (logFC = \u0026minus;0.088, FDR = 5.08 \u0026times; 10^-5), and DNA repair (logFC = \u0026minus;0.078, FDR = 4.42 \u0026times; 10^-3). The downregulation of oxidative phosphorylation concomitant with hypoxia activation suggests a metabolic shift toward the Warburg effect in high-risk tumors.\u003c/p\u003e\n\u003cp\u003e6.2. Tumor microenvironment characteristics\u003c/p\u003e\n\u003cp\u003eAnalysis of tumor microenvironment (TME) scores revealed striking differences in stromal composition between the risk groups (Figure 7). Compared with low-risk patients, high-risk patients exhibited significantly higher stromal scores (P \u0026lt; 0.0001), indicating increased stromal cell infiltration. However, immune scores did not differ significantly between the two groups (P \u0026gt; 0.05). Consequently, high-risk tumors displayed significantly lower tumor purity (P \u0026lt; 0.0001), reflecting a more complex and stroma-rich microenvironment. The increased stromal content in high-risk tumors is consistent with the enhanced EMT and angiogenic features observed, as cancer-associated fibroblasts and endothelial cells are major contributors to stromal components.\u003c/p\u003e\n\u003cp\u003e6.3. Immune checkpoint expression profiles\u003c/p\u003e\n\u003cp\u003eTo evaluate the potential responsiveness to immunotherapy across different risk groups, we analyzed the expression of immune checkpoint genes (Figure 7). Among the checkpoints examined, PDCD1LG2 (PD-L2) was significantly upregulated in high-risk patients (P = 0.001), along with CD80 (P = 0.006) and SIGLEC15 (P = 0.006). Interestingly, the key immunosuppressive enzyme IDO1 was significantly downregulated in the high-risk group (P = 0.045). The upregulation of PD-L2 and SIGLEC15 in high-risk patients suggests that these individuals may be responsive to immune checkpoint blockade targeting these molecules, providing a theoretical basis for MDRS-guided immunotherapy stratification.\u003c/p\u003e\n\u003cp\u003e6.4. Associations between MDRS genes and pathways\u003c/p\u003e\n\u003cp\u003eCorrelation analysis between the eight MDRS causal genes and Hallmark pathways revealed distinct functional associations (Figure 7). Risk-promoting genes (PDGFRB, DACT1, and PRSS23) showed strong positive correlations with EMT (r = 0.62\u0026ndash;0.68), angiogenesis (r = 0.53\u0026ndash;0.68), and TGF-\u0026beta; signaling (r = 0.42\u0026ndash;0.64), while exhibiting negative correlations with oxidative phosphorylation and DNA repair pathways. Notably, MYO15B, identified as a risk factor by MR analysis, displayed a unique pattern of negative correlations with most pro-tumorigenic pathways, suggesting a complex regulatory role. In contrast, protective genes (MSRB3, DAB2, BNC2, and CALD1) also showed positive correlations with EMT-related pathways, which may reflect their involvement in stromal remodeling processes that potentially exert tumor-suppressive effects through enhanced immune surveillance.\u003c/p\u003e\n\u003cp\u003e7. Network pharmacology analysis and molecular docking validation\u003c/p\u003e\n\u003cp\u003e7.1. Construction of the drug\u0026ndash;target interaction network\u003c/p\u003e\n\u003cp\u003eTo identify potential therapeutic agents targeting MDRS signature genes, we integrated data from the DGIdb, DrugBank, and CTD databases to perform a network pharmacology analysis. Based on the principle of \u0026ldquo;inhibiting risk genes while activating or preserving protective genes,\u0026rdquo; candidate drugs targeting the eight MDRS causal genes were systematically screened. In total, 15 drug\u0026ndash;gene interactions were identified (Figure 8). Among them, 12 were classified as therapeutically beneficial, targeting four key genes: PDGFRB (six drugs), DACT1 (two drugs), PRSS23 (two drugs), and MSRB3 (two drugs). The remaining interactions involving CALD1 and DAB2 were categorized as potentially harmful or context-dependent, as inhibition of these protective genes could inadvertently promote tumor progression.\u003c/p\u003e\n\u003cp\u003e7.2. Target-based therapeutic strategies\u003c/p\u003e\n\u003cp\u003eBased on the distinct functional roles of these targets, specific therapeutic strategies were proposed. PDGFRB (OR = 1.656) was identified as the most druggable risk gene and can be targeted by a range of FDA-approved tyrosine kinase inhibitors (TKIs) with established safety profiles, including imatinib, sunitinib, sorafenib, dasatinib, regorafenib, and lenvatinib. Similarly, other risk factors also provided actionable intervention points: the Wnt signaling regulator DACT1 (OR = 2.664) can be targeted by inhibitors such as ICG-001 and PRI-724, which are currently under clinical investigation; and the serine protease PRSS23 (OR = 1.438) offers opportunities for drug repurposing using FDA-approved protease inhibitors such as nafamostat and camostat. In contrast to these inhibitory strategies, therapeutic intervention for the protective factor MSRB3 (OR = 0.478) involves functional enhancement. As a key gene in oxidative stress defense, its activity may be augmented by antioxidants such as N-acetylcysteine (NAC) and ebselen, thereby restoring cellular resistance to metal-induced cell death.\u003c/p\u003e\n\u003cp\u003e7.3. Identification of a dual-function drug: sorafenib\u003c/p\u003e\n\u003cp\u003eNotably, our analysis identified sorafenib as a dual-function therapeutic agent of particular relevance for MDRS high-risk patients. Sorafenib not only inhibits the MDRS risk gene PDGFRB, but is also a well-established inducer of ferroptosis by inhibiting System Xc⁻ and subsequently depleting glutathione. This dual mechanism\u0026mdash;simultaneously targeting the MDRS molecular signature and inducing metal-dependent cell death\u0026mdash;positions sorafenib as an optimal candidate drug for high-risk cervical cancer patients identified by MDRS.\u003c/p\u003e\n\u003cp\u003e7.4. Molecular docking validation of the sorafenib\u0026ndash;PDGFRB interaction\u003c/p\u003e\n\u003cp\u003eTo validate the predicted drug\u0026ndash;target interaction at the molecular level, molecular docking simulations were performed between sorafenib and the PDGFRB protein structure. Using the AlphaFold-predicted human PDGFRB structure (UniProt ID: P09619), docking results demonstrated that sorafenib binds tightly within the ATP-binding pocket of the intracellular kinase domain of PDGFRB through hydrogen bonds and hydrophobic interactions. Blind docking analysis identified the highest-affinity binding pocket (Pocket C4) with a binding energy (Vina score) of \u0026minus;10.1 kcal/mol, indicating strong binding potential. Binding mode analysis (as shown in the figure) revealed that sorafenib primarily interacts with key residues in the kinase core region, including residues around Val607, Lys634, and Glu651, as well as the DFG motif (Asp822 and Phe823). This binding pattern is characteristic of a type II kinase inhibitor, which suppresses PDGFRB phosphorylation activity by occupying the ATP-binding site and inducing or stabilizing an inactive (DFG-out) conformation. In addition, a high template-matching score of 66.8 based on a homologous structure (PDB ID: 7MGJ) further validated the reliability of this binding conformation.\u003c/p\u003e\n\u003cp\u003eHowever, docking scores mainly reflect geometric complementarity and potential energy estimates under static conformations and cannot directly quantify the thermodynamic stability of drug\u0026ndash;target binding. Therefore, molecular dynamics simulations and MM/PBSA calculations were subsequently performed to validate this interaction at both kinetic and thermodynamic levels.\u003c/p\u003e\n\u003cp\u003e7.5. Molecular dynamics simulations and MM/PBSA binding free energy calculations confirm strong thermodynamic affinity of PDGFRB\u0026ndash;sorafenib\u003c/p\u003e\n\u003cp\u003eTo obtain quantitative pharmacodynamic constraints suitable for systems pharmacology modeling and to verify the stability of the docking pose on a dynamic timescale, all-atom molecular dynamics simulations were conducted for the PDGFRB\u0026ndash;sorafenib complex, followed by MM/PBSA binding free energy calculations based on equilibrated trajectories (Figure 8). Throughout the production simulations, the overall conformation of the complex remained stable, indicating that sorafenib forms a persistent and robust binding mode within the PDGFRB kinase pocket. MM/PBSA results showed that the binding free energy calculated using the Poisson\u0026ndash;Boltzmann (PB) model was \u0026Delta;G\u003csub\u003ePB\u003c/sub\u003e = \u0026minus;34.09 \u0026plusmn; 0.26 kcal/mol, while the Generalized Born (GB) model yielded \u0026Delta;G\u003csub\u003eGB\u003c/sub\u003e = \u0026minus;25.56 kcal/mol. These highly negative binding free energies indicate strong thermodynamic affinity of sorafenib for PDGFRB, corresponding theoretically to nanomolar or even lower dissociation constant ranges.\u003c/p\u003e\n\u003cp\u003eEnergy decomposition analysis further revealed that binding was dominated by van der Waals interactions (\u0026Delta;E\u003csub\u003evdW\u003c/sub\u003e = \u0026minus;40.53 kcal/mol), with additional contributions from electrostatic interactions (\u0026Delta;E\u003csub\u003eele\u003c/sub\u003e = \u0026minus;6.59 kcal/mol), whereas desolvation introduced a positive solvation penalty (\u0026Delta;G\u003csub\u003esolv\u003c/sub\u003e = +13.03 kcal/mol). Overall, the strong nonbonded interactions were sufficient to overcome the solvation cost, rendering the complex thermodynamically stable. Collectively, these results provide high-confidence molecular-scale evidence supporting sorafenib-mediated inhibition of PDGFRB and offer quantitative parameters for subsequent translation of \u0026Delta;G values into initial perturbation strengths in network propagation models (SINP).\u003c/p\u003e\n\u003cp\u003eTo experimentally validate the computational predictions, surface plasmon resonance (SPR) analysis was performed. The results demonstrated that sorafenib binds to the PDGFRB kinase domain with an equilibrium dissociation constant (KD) of 1.95 \u0026mu;M (Supplementary Figure1), confirming a biophysically relevant interaction consistent with the strong binding affinity predicted by molecular docking (\u0026minus;10.1 kcal/mol) and MM/PBSA calculations (\u0026Delta;G = \u0026minus;34.09 kcal/mol).\u003c/p\u003e\n\u003cp\u003e8. Multi-scale mechanistic validation of sorafenib targeting PDGFRB\u003c/p\u003e\n\u003cp\u003e8.1. SINP predicts that PDGFRB inhibition suppresses the collagen-rich ECM program\u003c/p\u003e\n\u003cp\u003eTo dissect the downstream systemic effects of sorafenib-mediated PDGFRB targeting at the tissue scale and to translate molecular-scale thermodynamic evidence into a propagatable pharmacodynamic perturbation, we quantitatively mapped the MM/PBSA-derived binding free energy (\u0026Delta;G) to the initial inhibition strength of PDGFRB and incorporated it into the SINP model to perform structure\u0026ndash;context\u0026ndash;coupled network propagation simulations (Figure 8). Driven by thermodynamically informed perturbation coefficients, SINP predicted that inhibition of PDGFRB induces widespread reprogramming of downstream gene expression networks, with the most strongly downregulated gene sets being highly enriched in pathways related to extracellular matrix (ECM) organization, collagen fibril assembly, and matrix remodeling (Figure 8).\u003c/p\u003e\n\u003cp\u003eAmong the genes exhibiting the greatest downregulation, collagen family members and canonical stromal molecules were prominently represented, including COL1A1, COL1A2, COL3A1, COL5A2, and COL6A3, suggesting that PDGFRB signaling sustains a \u0026ldquo;collagen- and matrix-enhanced\u0026rdquo; pro-fibrotic transcriptional program within the cervical cancer microenvironment. These findings provide tissue-scale mechanistic insights into the stromal enrichment and EMT activation observed in MDRS high-risk tumors.\u003c/p\u003e\n\u003cp\u003e8.2. Monocyte-specific PDGFRB-driven pro-fibrotic reprogramming revealed by single-cell in silico pharmacodynamic modeling\u003c/p\u003e\n\u003cp\u003eTo identify the principal cellular source of the PDGFRB-mediated ECM program, we systematically evaluated the expression distribution of PDGFRB across cell subpopulations in the scRNA-seq dataset. PDGFRB exhibited marked cell type\u0026ndash;specific enrichment within the tumor microenvironment, with monocytes serving as the predominant expressing cell population (expression proportion \u0026gt;80%), whereas expression was low or absent in macrophages and other immune cell types (Figure 8).\u003c/p\u003e\n\u003cp\u003eSubsequent single-cell pharmacodynamic modeling using scTenifoldKnk within monocytes revealed that PDGFRB inhibition led to pronounced perturbations in the single-cell regulatory network, with differentially affected genes predominantly enriched in collagen biosynthesis, ECM organization, and matrix remodeling pathways. Specifically, multiple collagen family genes (COL1A1, COL3A1, COL5A1/2, COL6A1/3) showed significant responses, accompanied by marked effects on the matrix remodeling enzyme MMP11 and the collagen cross-linking enzyme LOXL2 (Figure 8). In addition, the immune regulatory gene LAMP5 exhibited the highest Z-score among perturbed genes, suggesting that beyond its pro-fibrotic role, PDGFRB may also contribute to maintenance of an immunosuppressive microenvironment.\u003c/p\u003e\n\u003cp\u003eAs a control, the same perturbation analysis performed in the Epi_Tumor_Prolif subpopulation did not reveal significant network disruption (Hits = 0), indicating that the functional impact of PDGFRB is highly context dependent. These results suggest that PDGFRB exerts its key biological effects not within tumor epithelial cells themselves, but primarily by reprogramming monocytes toward a pro-fibrotic phenotype, thereby shaping a collagen-enriched tumor microenvironment.\u003c/p\u003e\n\u003cp\u003e8.3. Cross-scale convergence identifies seven shared ECM core genes across SINP and single-cell modeling\u003c/p\u003e\n\u003cp\u003eTo further assess the concordance between tissue-scale simulations (SINP) and single-cell modeling (scTenifoldKnk), we intersected the top 200 downregulated genes predicted by the SINP model with the significantly perturbed genes identified in monocytes (Figure 8). This cross-scale comparison identified seven fully concordant core genes: COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, and A2M. All of these genes encode key structural components or regulatory molecules of the ECM, collectively constituting a highly convergent collagen\u0026ndash;matrix core program.\u003c/p\u003e\n\u003cp\u003eThis cross-scale consistency not only validates the tissue-contextual relevance of the SINP model, but also provides single-cell\u0026ndash;level confirmation that PDGFRB inhibition can systematically dismantle the collagen-enriched ECM network. Together, these findings define a clear mechanistic axis\u0026mdash;PDGFRB\u0026ndash;monocyte\u0026ndash;collagen\u0026mdash;through which sorafenib may exert therapeutic effects in MDRS high-risk cervical cancer.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eThis study systematically characterized the phenotypic heterogeneity of ferroptosis- and cuproptosis-related pathways in cervical cancer (CESC) by integrating multi-omics data, causal inference, and machine learning modeling, and further established a mechanistic, chain-like linkage between risk stratification results, potentially actionable targets, and candidate therapeutic agents. In contrast to previous studies that focused on a single mode of cell death or a single omics layer, we not only constructed a clinically applicable risk score (MDRS) at the population level, but also established an interpretable, multi-scale mechanistic evidence chain spanning tissue- and single-cell resolutions. Starting from genetically supported causal targets identified by Mendelian randomization (MR), we integrated molecular-scale drug\u0026ndash;target thermodynamic evidence (molecular docking and MM/PBSA), tissue-scale signal propagation inference (SINP), and single-cell\u0026ndash;level identification of effector cell populations and regulatory network perturbations (scTenifoldKnk). These analyses ultimately converged on a well-defined microenvironmental regulatory axis\u0026mdash;the PDGFRB\u0026ndash;monocyte\u0026ndash;collagen/ECM axis. This axis not only provides a more direct mechanistic explanation for the stromal enrichment and enhanced epithelial\u0026ndash;mesenchymal transition (EMT) observed in MDRS high-risk tumors, but also offers a more testable theoretical basis for the potential benefit of sorafenib in specific cervical cancer subtypes.\u003c/p\u003e \u003cp\u003eAt the pathway level, we observed only a weak overall correlation between ferroptosis and cuproptosis in CESC, with more than half of patients exhibiting discordant phenotypes such as \u0026ldquo;high ferroptosis/low cuproptosis\u0026rdquo; or \u0026ldquo;low ferroptosis/high cuproptosis.\u0026rdquo; This finding does not imply that the two pathways are independent or unrelated; rather, it likely reflects the strong context dependence and spatial heterogeneity of metal-dependent cell death pathways across patients, cell lineages, and stages of tumor evolution. Single-cell analyses further supported this interpretation: myeloid cells and tumor epithelial cells generally exhibited higher MDAS activation, whereas pseudotime trajectories revealed a progressive decline in MDAS\u0026mdash;particularly the ferroptosis-related component\u0026mdash;during epithelial malignant transformation. This pattern suggests that, during tumor progression, cancer cells may acquire a survival advantage by reducing their sensitivity to metal-dependent death through metabolic reprogramming. It should be emphasized, however, that these observations provide associative evidence at this level, and the causal direction and underlying molecular mechanisms require further experimental validation.\u003c/p\u003e \u003cp\u003eA key strength of this study lies in the use of MR to elevate candidate hub genes from correlation-based signals to high-confidence targets with causal associations, thereby providing genetic support for downstream drug development and mechanistic inference. Among the eight causal genes, factors such as DACT1 implicate Wnt and developmental regulatory networks in the risk phenotype, whereas the importance of PDGFRB stems from its well-established roles in angiogenesis, stromal activation, and fibrotic processes[\u003cspan additionalcitationids=\"CR21 CR22\" citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e]. Recent studies have demonstrated that inducing transient target expression can engineer therapeutic vulnerabilities, thereby expanding the scope of targeted cancer therapies[\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. Importantly, we did not stop at identifying PDGFRB as a \u0026ldquo;risk gene,\u0026rdquo; but instead used multi-scale modeling and cross-scale consistency validation to delineate its potential effector pathway. At the tissue scale, SINP simulations predicted that PDGFRB inhibition would markedly suppress collagen-enriched ECM programs. At the single-cell scale, scTenifoldKnk revealed a pronounced cell-context dependency: PDGFRB inhibition triggered substantial collapse of collagen- and matrix remodeling\u0026ndash;related networks in monocytes, whereas comparable perturbations were largely absent in tumor epithelial cell controls. These findings suggest that the risk effect of PDGFRB is primarily mediated through myeloid cell\u0026ndash;driven pro-fibrotic reprogramming and matrix deposition, rather than solely through intrinsic tumor cell proliferation, thereby mechanistically linking MR-inferred genetic risk with the stromal-enriched/EMT-enhanced phenotype of MDRS high-risk tumors.\u003c/p\u003e \u003cp\u003eTo minimize methodological bias from any single model, we further performed cross-scale consistency analyses by intersecting downregulated genes predicted by tissue-scale SINP simulations with significantly perturbed genes identified by single-cell scTenifoldKnk modeling. This approach yielded seven fully overlapping ECM core genes (COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, and A2M). These genes are canonical components and regulators of collagen and matrix remodeling, and are closely associated with stromal stiffening, cell migration, and EMT programs in cancer[\u003cspan additionalcitationids=\"CR26\" citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. This cross-scale convergence substantially strengthens the robustness of our conclusions, demonstrating that the link between PDGFRB inhibition and ECM program collapse is not an artifact of a single algorithm, but a reproducible biological theme detectable across scales and data structures. Collectively, these findings establish the PDGFRB\u0026ndash;monocyte\u0026ndash;collagen axis as a key mechanistic pathway connecting risk stratification with microenvironmental remodeling.\u003c/p\u003e \u003cp\u003eFrom a translational perspective, network pharmacology analyses identified sorafenib as a candidate drug with dual potential to both target PDGFRB and induce ferroptosis, and this hypothesis was further reinforced by molecular-scale simulations. Emerging targeted protein degradation technologies may offer additional strategies for targeting previously undruggable components of these pathways[\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e, \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]. Importantly, molecular docking was not treated as standalone evidence; instead, molecular dynamics simulations combined with MM/PBSA binding free energy calculations provided thermodynamically more meaningful quantitative support. The binding free energy of sorafenib to PDGFRB reached ΔG_PB\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;34.09\u0026thinsp;\u0026plusmn;\u0026thinsp;0.26 kcal/mol under the PB model, indicating strong binding stability. Using this thermodynamic constraint, we quantitatively mapped ΔG to the initial perturbation strength in the SINP model, enabling inference of system-level effects at the tissue scale, which were then validated at the single-cell level as primarily occurring within monocyte-driven pro-fibrotic networks. Thus, the potential role of sorafenib in MDRS high-risk CESC can be more rationally summarized as operating through two complementary pathways: (i) weakening tumor cell survival advantages by inducing ferroptosis at the tumor cell level[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e, \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]; and (ii) attenuating myeloid cell\u0026ndash;driven collagen deposition by inhibiting PDGFRB signaling within the tumor microenvironment, thereby \u0026ldquo;softening\u0026rdquo; the matrix and reducing EMT-supportive conditions. This mechanistic framework not only provides an actionable interpretation of MDRS high-risk phenotypes, but also suggests that the therapeutic benefit of sorafenib in cervical cancer may be closely tied to patient heterogeneity, warranting further investigation in stratified clinical studies.\u003c/p\u003e \u003cp\u003eIn addition, MDRS, as a machine learning\u0026ndash;based risk score derived from a set of causal genes, not only significantly improved predictive performance when integrated with clinical variables, but also enabled a more biologically interpretable stratification. The MDRS high-risk group was characterized by aggressive features such as EMT, angiogenesis, hypoxia, and stromal enrichment, accompanied by upregulation of immune checkpoint molecules (e.g., PD-L2 and SIGLEC15), suggesting a \u0026ldquo;stroma-enriched/immune-excluded\u0026rdquo; tumor microenvironment state[\u003cspan additionalcitationids=\"CR33\" citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e]. Novel approaches targeting PD-L1 stability through DHHC3 degradation have shown promise in overcoming immune checkpoint blockade resistance[\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e, \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e], suggesting that combination strategies targeting both stromal programs and immune checkpoints may be particularly beneficial in high-risk patients. Accordingly, the value of MDRS extends beyond prognostic prediction, as it may also serve as a mechanism-oriented patient stratification tool to guide future exploration of combination strategies involving stromal/ECM targeting, metal-dependent cell death induction, and immune modulation. Nonetheless, such strategies require rigorous experimental and clinical validation to assess feasibility and safety.\u003c/p\u003e \u003cp\u003eDespite establishing a multi-scale framework linking causal genes to mechanistic validation, several limitations should be acknowledged. First, the MDRS model was primarily developed and internally evaluated using the TCGA cohort, and its generalizability requires further validation in independent, multi-center clinical cohorts. Second, the MR analysis relied on blood-derived cis-eQTLs from eQTLGen; although this resource offers large sample size and strong statistical power, the lack of tissue-specific eQTLs may limit the capture of cervix-local regulatory effects. Future availability of large-scale cervical tissue or single-cell eQTL datasets would enhance tissue consistency in causal inference. Third, SINP and scTenifoldKnk are computational simulation frameworks that provide mechanistic inference and consistency evidence; however, key conclusions\u0026mdash;such as the role of PDGFRB in monocyte-driven ECM programs\u0026mdash;require confirmation through in vitro and in vivo experiments. Finally, the potential benefit of sorafenib in MDRS high-risk subgroups remains to be validated in prospective clinical studies, along with the identification of optimal biomarker combinations for actionable stratified treatment strategies.\u003c/p\u003e \u003cp\u003eThe SPR analysis revealed that sorafenib binds to the PDGFRB kinase domain with an equilibrium dissociation constant (KD) of 1.95 \u0026micro;M, which falls within the moderate affinity range. This finding is consistent with the pharmacological profile of sorafenib as a multi-target tyrosine kinase inhibitor. As sorafenib was originally developed to target RAF kinases and subsequently found to potently inhibit VEGFR and PDGFR family members, PDGFRB represents a secondary rather than primary target. The observed\u003c/p\u003e \u003cp\u003emicromolar affinity is therefore biologically plausible and aligns with previous reports characterizing sorafenib's broad kinase inhibition spectrum. Importantly, this moderate binding affinity, combined with the favorable thermodynamic profile revealed by MM/PBSA calculations (ΔG\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;34.09 kcal/mol) and the strong predicted binding pose from molecular docking (\u0026minus;\u0026thinsp;10.1 kcal/mol), collectively support PDGFRB as a therapeutically relevant target of sorafenib in cervical cancer. The convergence of computational predictions and experimental validation strengthens the rationale for repurposing sorafenib as a PDGFRB-targeted therapeutic strategy in this disease context.\u003c/p\u003e \u003cp\u003eIn summary, this study proposes and validates a coherent mechanistic axis linking metal-dependent cell death phenotypes, microenvironmental matrix remodeling, and clinical risk stratification, and establishes an integrative framework spanning pathway phenotypes, causal targets, interpretable prediction, multi-scale mechanistic validation, and candidate therapeutic strategies. This framework offers new insights into heterogeneity-aware management of cervical cancer and the development of mechanism-driven therapeutic approaches.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eData availability statement:\u003c/strong\u003e The datasets presented in this study are publicly available in online repositories. The bulk transcriptomic data and clinical information (TCGA and GTEx) can be accessed via the GDC Data Portal (https://portal.gdc.cancer.gov/). The single-cell RNA sequencing data generated in this study have been deposited in the Gene Expression Omnibus (GEO) under accession number GSE208653. For Mendelian randomization analysis, summary statistics were obtained from the eQTLGen Consortium (accessed via IEU OpenGWAS) and the FinnGen consortium (Dataset ID: finn-b-C3_CERVIX_UTERI).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding Statement:\u0026nbsp;\u003c/strong\u003eThis work was supported by the Hubei Provincial Natural Science Foundation Joint Fund (Grant No. 2025AFD319) and the Beijing Science and Technology Innovation Medical Development Foundation (Grant No. KC2023-JX-0288-RQ21).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConflict of interest:\u0026nbsp;\u003c/strong\u003eThe authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n \u003cli\u003e\u003cem\u003eWHO Guidelines Approved by the Guidelines Review Committee\u003c/em\u003e, in \u003cem\u003eWHO guideline for screening and treatment of cervical pre-cancer lesions for cervical cancer prevention: Use of dual-stain cytology to triage women after a positive test for human papillomavirus (HPV)\u003c/em\u003e. 2024, World Health Organization \u0026copy; World Health Organization 2024.: Geneva.\u003c/li\u003e\n \u003cli\u003eSingh, D., et al., \u003cem\u003eGlobal estimates of incidence and mortality of cervical cancer in 2020: a baseline analysis of the WHO Global Cervical Cancer Elimination Initiative.\u003c/em\u003e Lancet Glob Health, 2023. \u003cstrong\u003e11\u003c/strong\u003e(2): p. e197-e206.\u003c/li\u003e\n \u003cli\u003eAbu-Rustum, N.R., et al., \u003cem\u003eNCCN Guidelines\u0026reg; Insights: Uterine Neoplasms, Version 3.2025.\u003c/em\u003e J Natl Compr Canc Netw, 2025. \u003cstrong\u003e23\u003c/strong\u003e(8): p. 284-291.\u003c/li\u003e\n \u003cli\u003eMonk, B.J., et al., \u003cem\u003eFirst-Line Pembrolizumab + Chemotherapy Versus Placebo + Chemotherapy for Persistent, Recurrent, or Metastatic Cervical Cancer: Final Overall Survival Results of KEYNOTE-826.\u003c/em\u003e J Clin Oncol, 2023. \u003cstrong\u003e41\u003c/strong\u003e(36): p. 5505-5511.\u003c/li\u003e\n \u003cli\u003eOaknin, A., et al., \u003cem\u003eEMPOWER CERVICAL-1: Effects of cemiplimab versus chemotherapy on patient-reported quality of life, functioning and symptoms among women with recurrent cervical cancer.\u003c/em\u003e Eur J Cancer, 2022. \u003cstrong\u003e174\u003c/strong\u003e: p. 299-309.\u003c/li\u003e\n \u003cli\u003eGalluzzi, L., et al., \u003cem\u003eMolecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018.\u003c/em\u003e Cell Death Differ, 2018. \u003cstrong\u003e25\u003c/strong\u003e(3): p. 486-541.\u003c/li\u003e\n \u003cli\u003eGalluzzi, L., et al., \u003cem\u003eRegulated cell death and adaptive stress responses.\u003c/em\u003e Cell Mol Life Sci, 2016. \u003cstrong\u003e73\u003c/strong\u003e(11-12): p. 2405-10.\u003c/li\u003e\n \u003cli\u003ePeng, F., et al., \u003cem\u003eRegulated cell death (RCD) in cancer: key pathways and targeted therapies.\u003c/em\u003e Signal Transduct Target Ther, 2022. \u003cstrong\u003e7\u003c/strong\u003e(1): p. 286.\u003c/li\u003e\n \u003cli\u003eJiang, X., B.R. Stockwell, and M. Conrad, \u003cem\u003eFerroptosis: mechanisms, biology and role in disease.\u003c/em\u003e Nat Rev Mol Cell Biol, 2021. \u003cstrong\u003e22\u003c/strong\u003e(4): p. 266-282.\u003c/li\u003e\n \u003cli\u003eTsvetkov, P., et al., \u003cem\u003eCopper induces cell death by targeting lipoylated TCA cycle proteins.\u003c/em\u003e Science, 2022. \u003cstrong\u003e375\u003c/strong\u003e(6586): p. 1254-1261.\u003c/li\u003e\n \u003cli\u003eLiu, N. and M. Chen, \u003cem\u003eCrosstalk between ferroptosis and cuproptosis: From mechanism to potential clinical application.\u003c/em\u003e Biomed Pharmacother, 2024. \u003cstrong\u003e171\u003c/strong\u003e: p. 116115.\u003c/li\u003e\n \u003cli\u003eM\u0026uuml;ller, S., et al., \u003cem\u003eCopper and iron orchestrate cell-state transitions in cancer and immunity.\u003c/em\u003e Trends Cell Biol, 2025. \u003cstrong\u003e35\u003c/strong\u003e(2): p. 105-114.\u003c/li\u003e\n \u003cli\u003eXia, P., et al., \u003cem\u003eDeciphering the cellular and molecular landscape of cervical cancer progression through single-cell and spatial transcriptomics.\u003c/em\u003e NPJ Precis Oncol, 2025. \u003cstrong\u003e9\u003c/strong\u003e(1): p. 158.\u003c/li\u003e\n \u003cli\u003eOu, Z., et al., \u003cem\u003eSingle-Nucleus RNA Sequencing and Spatial Transcriptomics Reveal the Immunological Microenvironment of Cervical Squamous Cell Carcinoma.\u003c/em\u003e Adv Sci (Weinh), 2022. \u003cstrong\u003e9\u003c/strong\u003e(29): p. e2203040.\u003c/li\u003e\n \u003cli\u003eHu, X., et al., \u003cem\u003eBenchmarking Mendelian randomization methods for causal inference using genome-wide association study summary statistics.\u003c/em\u003e Am J Hum Genet, 2024. \u003cstrong\u003e111\u003c/strong\u003e(8): p. 1717-1735.\u003c/li\u003e\n \u003cli\u003eHemani, G., et al., \u003cem\u003eThe MR-Base platform supports systematic causal inference across the human phenome.\u003c/em\u003e Elife, 2018. \u003cstrong\u003e7\u003c/strong\u003e.\u003c/li\u003e\n \u003cli\u003eHuang, S., et al., \u003cem\u003eArtificial intelligence in cancer diagnosis and prognosis: Opportunities and challenges.\u003c/em\u003e Cancer Lett, 2020. \u003cstrong\u003e471\u003c/strong\u003e: p. 61-71.\u003c/li\u003e\n \u003cli\u003eHou, T., et al., \u003cem\u003eAssessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations.\u003c/em\u003e J Chem Inf Model, 2011. \u003cstrong\u003e51\u003c/strong\u003e(1): p. 69-82.\u003c/li\u003e\n \u003cli\u003eHopkins, A.L., \u003cem\u003eNetwork pharmacology: the next paradigm in drug discovery.\u003c/em\u003e Nat Chem Biol, 2008. \u003cstrong\u003e4\u003c/strong\u003e(11): p. 682-90.\u003c/li\u003e\n \u003cli\u003eStrell, C., E. Rodr\u0026iacute;guez-Tom\u0026agrave;s, and A. \u0026Ouml;stman, \u003cem\u003eFunctional and clinical roles of stromal PDGF receptors in tumor biology.\u003c/em\u003e Cancer Metastasis Rev, 2024. \u003cstrong\u003e43\u003c/strong\u003e(4): p. 1593-1609.\u003c/li\u003e\n \u003cli\u003ePandey, P., et al., \u003cem\u003eNew insights about the PDGF/PDGFR signaling pathway as a promising target to develop cancer therapeutic strategies.\u003c/em\u003e Biomed Pharmacother, 2023. \u003cstrong\u003e161\u003c/strong\u003e: p. 114491.\u003c/li\u003e\n \u003cli\u003eShi, J., et al., \u003cem\u003ePericyte-myofibroblast transition: a novel mechanism in peritoneal fibrosis and the effect of Asiaticoside.\u003c/em\u003e J Adv Res, 2025.\u003c/li\u003e\n \u003cli\u003eCarroll, S.H., et al., \u003cem\u003eGenetic requirement of dact1/2 to regulate noncanonical Wnt signaling and calpain 8 during embryonic convergent extension and craniofacial morphogenesis.\u003c/em\u003e Elife, 2024. \u003cstrong\u003e13\u003c/strong\u003e.\u003c/li\u003e\n \u003cli\u003eYang, J., et al., \u003cem\u003eTargeting an Inducible SALL4-Mediated Cancer Vulnerability with Sequential Therapy.\u003c/em\u003e Cancer Res, 2021. \u003cstrong\u003e81\u003c/strong\u003e(23): p. 6018-6028.\u003c/li\u003e\n \u003cli\u003eJiang, Y., et al., \u003cem\u003eTargeting extracellular matrix stiffness and mechanotransducers to improve cancer therapy.\u003c/em\u003e J Hematol Oncol, 2022. \u003cstrong\u003e15\u003c/strong\u003e(1): p. 34.\u003c/li\u003e\n \u003cli\u003eYuan, Z., et al., \u003cem\u003eCollagen remodeling-mediated signaling pathways and their impact on tumor therapy.\u003c/em\u003e J Biol Chem, 2025. \u003cstrong\u003e301\u003c/strong\u003e(3): p. 108330.\u003c/li\u003e\n \u003cli\u003eSuzuki, M., et al., \u003cem\u003eAssociation of rs6983561 polymorphism at 8q24 with prostate cancer mortality in a Japanese population.\u003c/em\u003e Clin Genitourin Cancer, 2011. \u003cstrong\u003e9\u003c/strong\u003e(1): p. 46-52.\u003c/li\u003e\n \u003cli\u003eDai, M., et al., \u003cem\u003eTargeted Protein Degradation: An Important Tool for Drug Discovery for \u0026quot;Undruggable\u0026quot; Tumor Transcription Factors.\u003c/em\u003e Technol Cancer Res Treat, 2022. \u003cstrong\u003e21\u003c/strong\u003e: p. 15330338221095950.\u003c/li\u003e\n \u003cli\u003eHashmi, F., et al., \u003cem\u003eEXPRESS: Phospholipase C gamma mediates endogenous brain-derived neurotrophic factor - regulated calcitonin gene-related peptide expression in colitis - induced visceral pain.\u003c/em\u003e Mol Pain, 2016. \u003cstrong\u003e12\u003c/strong\u003e.\u003c/li\u003e\n \u003cli\u003eZhang, L., et al., \u003cem\u003eSorafenib triggers ferroptosis via inhibition of HBXIP/SCD axis in hepatocellular carcinoma.\u003c/em\u003e Acta Pharmacol Sin, 2023. \u003cstrong\u003e44\u003c/strong\u003e(3): p. 622-634.\u003c/li\u003e\n \u003cli\u003eLiu, R., et al., \u003cem\u003eRoles and Mechanisms of Ferroptosis in Sorafenib Resistance for Hepatocellular Carcinoma.\u003c/em\u003e J Hepatocell Carcinoma, 2024. \u003cstrong\u003e11\u003c/strong\u003e: p. 2493-2504.\u003c/li\u003e\n \u003cli\u003eZheng, S., et al., \u003cem\u003eTumor battlefield within inflamed, excluded or desert immune phenotypes: the mechanisms and strategies.\u003c/em\u003e Exp Hematol Oncol, 2024. \u003cstrong\u003e13\u003c/strong\u003e(1): p. 80.\u003c/li\u003e\n \u003cli\u003eChen, D.S. and I. Mellman, \u003cem\u003eOncology meets immunology: the cancer-immunity cycle.\u003c/em\u003e Immunity, 2013. \u003cstrong\u003e39\u003c/strong\u003e(1): p. 1-10.\u003c/li\u003e\n \u003cli\u003eLi, L., et al., \u003cem\u003eThe role of Siglec-15 in tumor immunity: mechanism and therapy.\u003c/em\u003e Mol Cancer Ther, 2025.\u003c/li\u003e\n \u003cli\u003eShi, Y.Y., et al., \u003cem\u003eTreating ICB-resistant cancer by inhibiting PD-L1 via DHHC3 degradation induced by cell penetrating peptide-induced chimera conjugates.\u003c/em\u003e Cell Death Dis, 2024. \u003cstrong\u003e15\u003c/strong\u003e(9): p. 701.\u003c/li\u003e\n \u003cli\u003eLiu, M., et al., \u003cem\u003eA replication study examining three common single-nucleotide polymorphisms and the risk of prostate cancer in a Japanese population.\u003c/em\u003e Prostate, 2011. \u003cstrong\u003e71\u003c/strong\u003e(10): p. 1023-32.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[{"identity":"e567747b-3e0c-4858-a05e-da1eb6cd60d3","identifier":"10.13039/501100003819","name":"Natural Science Foundation of Hubei Province","awardNumber":"2025AFD319","order_by":0}],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"Department of Gynecological Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei, China","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"cervical cancer, ferroptosis, cuproptosis, machine learning, Mendelian randomization, single-cell RNA sequencing, network pharmacology, molecular dynamics, systems pharmacology simulation","lastPublishedDoi":"10.21203/rs.3.rs-8811979/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-8811979/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eBackground: Ferroptosis and cuproptosis are two distinct forms of metal-dependent regulated cell death that have emerged as important mechanisms in tumor biology. However, the crosstalk between these pathways and their clinical relevance in cervical cancer remain largely unexplored.\u003c/p\u003e\n\u003cp\u003eMethods: Using the TCGA-CESC cohort (n = 304), we quantified ferroptosis- and cuproptosis-related pathway activities by ssGSEA/GSVA and constructed a Metal Death Activity Score (MDAS). Weighted gene co-expression network analysis (WGCNA) was applied to identify MDAS-associated modules and hub genes. We further integrated single-cell RNA sequencing (scRNA-seq) data comprising 74,454 cells to characterize the cellular heterogeneity and dynamic evolution of metal death phenotypes within the tumor microenvironment. Core genes causally associated with cervical cancer risk were screened using two-sample Mendelian randomization (MR), and the optimal model among 74 machine-learning algorithms was selected to construct a Metal Death Risk Score (MDRS). At the therapeutic level, potential targeted agents were identified via network pharmacology, and PDGFRB–drug interactions were validated using molecular docking, molecular dynamics simulations, and MM/PBSA binding free energy calculations. To bridge molecular-scale evidence with tissue- and single-cell–scale effects, we developed a structure–context–coupled network propagation (SINP) model and performed single-cell pharmacodynamic simulations using scTenifoldKnk, enabling cross-scale mechanistic inference and validation.\u003c/p\u003e\n\u003cp\u003eResults: MDAS was significantly elevated in cervical cancer tissues, while ferroptosis and cuproptosis were largely independent at the global level (r = 0.053, P = 0.355). WGCNA identified 293 hub genes across five MDAS-associated modules. Single-cell analyses revealed higher MDAS activation in myeloid cells and tumor epithelial cells, with a progressive decline along the epithelial malignant transformation trajectory (ρ = −0.231, P = 3.5×10^-44). MR analysis identified eight causal genes, including four risk factors (DACT1, PDGFRB, PRSS23, MYO15B) and four protective factors (MSRB3, CALD1, DAB2, BNC2). The optimal MDRS model (Elastic Net) achieved a C-index of 0.736 after integration with clinical variables and significantly improved risk reclassification (NRI = 0.073, P = 0.044). High-risk patients exhibited enhanced epithelial–mesenchymal transition (EMT), angiogenesis, and suppressed oxidative phosphorylation. Network pharmacology identified sorafenib as a dual-function candidate drug capable of both targeting PDGFRB and inducing ferroptosis. Molecular docking indicated stable binding of sorafenib to the ATP-binding pocket of the PDGFRB kinase domain (Vina score = −10.1 kcal/mol), which was further supported by molecular dynamics simulations, MM/PBSA analysis (ΔG_PB = −34.09 ± 0.26 kcal/mol), and surface plasmon resonance (SPR) validation (KD = 1.95 μM). At the tissue scale, SINP predicted that PDGFRB inhibition markedly suppresses collagen-enriched extracellular matrix (ECM) programs. scTenifoldKnk simulations further demonstrated cell context–dependent effects, revealing that PDGFRB drives pro-fibrotic reprogramming in monocytes and triggers collapse of the collagen network. Cross-scale consistency analysis converged on seven shared ECM core genes (COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, LUM, A2M), establishing the PDGFRB–monocyte–collagen axis as a key mechanistic pathway linking the high-MDRS phenotype, stromal stiffening, and EMT.\u003c/p\u003e\n\u003cp\u003eConclusions: This study establishes an integrative framework encompassing metal death pathway activity scoring, causal gene identification, machine-learning–based risk stratification, and multiscale mechanistic simulation. MDRS enables clinical risk stratification in cervical cancer patients and highlights sorafenib as a potential precision therapeutic candidate for high-MDRS patients, with mechanisms likely involving PDGFRB inhibition, ferroptosis induction, and microenvironmental remodeling through suppression of monocyte-driven pro-fibrotic programs.\u003c/p\u003e","manuscriptTitle":"Integrated Multi-Omics Profiling Maps Ferroptosis–Cuproptosis Diversity in Cervical Cancer and Identifies a PDGFRB-Driven Monocyte Fibrotic Program Targeted by Sorafenib","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-02-10 11:36:01","doi":"10.21203/rs.3.rs-8811979/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"651d6d76-e912-493c-beaa-705c54cc003f","owner":[],"postedDate":"February 10th, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":62488879,"name":"Oncology"}],"tags":[],"updatedAt":"2026-02-10T11:36:02+00:00","versionOfRecord":[],"versionCreatedAt":"2026-02-10 11:36:01","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-8811979","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-8811979","identity":"rs-8811979","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}