Gene-interaction perturbation network analysis identifies distinct subtypes and actionable targets for aggressive thyroid cancer: an integrative multi-omics and machine learning approach for precise medicine

Gene-interaction perturbation network analysis identifies distinct subtypes and actionable targets for aggressive thyroid cancer: an integrative multi-omics and machine learning approach for precise medicine | 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 Gene-interaction perturbation network analysis identifies distinct subtypes and actionable targets for aggressive thyroid cancer: an integrative multi-omics and machine learning approach for precise medicine Jing Sun, Jiaxuan Huai, Jianpeng Zhang, Wang Qi, Zhaokai Zhou, and 7 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7834338/v1 This work is licensed under a CC BY 4.0 License Status: Under Revision Version 1 posted 9 You are reading this latest preprint version Abstract High-risk thyroid cancer such as aggressive papillary thyroid cancer (PTC) and anaplastic thyroid cancer (ATC) poses significant clinical challenges due to their tumor heterogeneity and limited therapeutic strategies. We constructed a genetic interaction (GI) perturbation network based on functional interactions from the Reactome database and identified distinct molecular subtypes of thyroid cancer with unique biological characteristics. Furthermore, a GI-based machine learning-assisted prognostic model was established for predicting progression of thyroid cancer patients, and RRM1 was identified as an essential gene that is significantly correlated with progression of thyroid cancer. Drug repositioning and molecular docking analyses demonstrated that gemcitabine targeting RRM1 can serve as a promising therapeutic strategy for ATC. Subsequent in vitro and in vivo experiments validated the therapeutic potential of gemcitabine in targeting RRM1, with significant efficacy in inducing ATC cell apoptosis and inhibiting tumor growth in a xenograft mouse model. Using single-cell RNA-sequencing (scRNA-seq) analysis and multiplex immunohistochemistry (mIHC), we revealed a unique interaction between RRM1-positive epithelial cells and a specific subset of cancer-associated fibroblasts (CAF) in ATC, highlighting the importance of the tumor microenvironment in driving aggressiveness. In summary, our findings highlight the critical role of dynamic gene-interaction network in thyroid cancer progression and offer a promising therapeutic strategy by targeting RRM1 for patients with high-risk thyroid cancer. Aggressive thyroid cancer Genetic interaction perturbation Machine learning Molecular subtyping Drug repositioning Precise medicine Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Introduction Thyroid cancer is one of the most common endocrine malignancies, and its incidence has significantly increased globally in recent decades. Papillary thyroid cancer (PTC) accounts for about 84% of all thyroid cancer cases. Originating from thyroid follicular cells, the papillary, follicular (approximately 4%), and oncocytic (approximately 2%) types are categorized as well-differentiated ones. The more aggressive forms of thyroid cancer include poorly differentiated thyroid cancer (PDTC, about 5%) and anaplastic thyroid cancer (ATC, about 1%). Medullary thyroid cancer (MTC), which develops from parafollicular C cells, accounts for approximately 4% of all cases [ 1 ]. PTC is generally an indolent tumor with a favorable prognosis, exhibiting a 10-year survival rate of approximately 93% [ 2 ]. However, a significant proportion of PTC patients still exhibit a propensity for metastasis and recurrence. Studies have shown that local or regional recurrences occur in 5 to 30% of patients with PTC, which can lead to a poor prognosis [ 3 – 5 ]. In contrast, ATC is a rare but highly aggressive subtype, characterized by its undifferentiated nature and extremely poor prognosis with a 4-month median overall survival [ 6 ]. Therefore, there is an urgent need to identify actionable therapeutic targets and tailor personalized approaches to offer more effective treatment options for high-risk PTC and more aggressive ATC patients. Comprehensive analysis of genetic alterations has become a critical tool for molecular subtyping and facilitates our understanding of the malignant phenotypes of thyroid cancer [ 7 – 9 ]. This approach not only aids in predicting disease progression and prognosis but also guides optimal targeted therapies, which can improve treatment efficacy and control adverse effects. Tailored treatment strategies can be designed for each patient according to the molecular profile of their tumor samples, overcoming the challenges caused by tumor heterogeneity. By integrating molecular typing with precision treatment, clinicians can develop more effective and personalized management strategies, thus improving patient outcomes [ 10 ]. However, many existing molecular subtypes are mainly proposed based on static gene expression patterns, regardless of the dynamic changes in gene regulation [ 9 , 11 ]. Gene expression can be significantly regulated during disease progression, making subtyping classification that is based on transcriptional profiles unstable and hard to reproduce. Conversely, biological networks are more stable and can stably reflect the biological status of tissues [ 12 – 14 ]. Gene network analysis is more reliable and effective than single-gene approaches in handling high-dimensional data [ 15 ]. Nonetheless, most network-based approaches only focus on individual genes and neglect gene interactions (GIs). To address this issue, we proposed gene expression rank-based individual-specific GI perturbation method that incorporates both gene nodes and their interactions. Although gene interactions are highly-conserved in normal tissues, they are significantly disrupted in tumors. The genetic interaction perturbation within the network quantifies the dynamic changes for each gene pair, allowing for a comprehensive reflection of the pathological condition at the individual level. Using this approach, we identified thyroid cancer-specific GI perturbations and GI perturbation-based molecular subtypes with distinct biological characteristics. Furthermore, GI perturbation-derived prognostic models and actionable targets were designed to guide the precise medicine for thyroid cancer patients. Overall, our findings provide a novel perspective on the gene interaction perturbation underlying thyroid cancer and offer promising insights for the development of personalized therapeutic strategies. Materials and methods Transcriptome data collection and processing The transcriptome profiles of 279 normal thyroid specimen in the GTEx database [ 16 ] were used as the benchmark data, and the RNA-seq data of 349 PTC samples from GSE213647 (Illumina NextSeq 500 platform) were used as the training data [ 17 ]. In addition, RNA-seq data of the TCGA-THCA cohort were used for clustering validation and genomic analyses. Another four microarray datasets (GSE60542 [ 18 ], GSE33630 [ 19 ], GSE27155 [ 20 ], and GSE76039 [ 21 ]) containing transcriptomic profiles of diverse thyroid cancer types were obtained from the GEO repository. To investigate c-Myc perturbation effects on RRM1 gene expression, transcriptome data from control siRNA (Ctrl-siRNA) and c-Myc siRNA-treated cells were obtained from GSE126739 [ 22 ] and GSE87693 [ 23 ]. For microarray processing, probe IDs were annotated to gene symbols using platform-specific annotation files. If a gene symbol corresponds to several probes, the probe with the highest measurement value was selected to represent the gene’s expression. A detailed summary of datasets used in this study is shown in Table S1 . Construction of gene interaction-perturbation network A total of 269,769 functional interactions were obtained from the Reactome database and used as the initial background interactome (Table S2 ). These interactions were used to establish the background network for thyroid cancer. Genes that were not present in our cohorts (279 GTEx normal thyroid tissues and 349 PTC samples) were removed, and the remaining genes were incorporated into the background network. For each gene in the background network, its rank was calculated in each sample according to its transcriptional value. The expression matrix was transformed to a rank matrix (the lowest expression means the minimum rank, while the highest expression means the max rank). Subsequently, we generated a delta rank matrix, where rows represent interactions and columns represent samples. The delta rank was calculated by subtracting the ranks of gene pairs connected by an interaction in the background network. Gene interactions (GIs) in normal specimen were considered stable and served as the benchmark network. We calculated the benchmark delta rank vector based on the mean relative ranks of gene pairs in all normal tissues. Finally, we obtained the interaction-perturbation matrix by subtracting the benchmark delta rank vector from the delta rank of each sample. This matrix measures cancer-specific interaction perturbations in the same background network. We performed differential analysis between cancer samples and normal tissues to select the top 3000 interactions that significantly differed (Wilcoxon test) and the top 3000 interactions with the highest SD among all PTC samples. A total of 1,262 interactions were overlapped in these two sets and retained for clustering (Table S3). The consensus clustering algorithm was used to identify distinct subtypes for PTC patients based on PTC-specific GI-perturbation matrix. The nearest template prediction (NTP) algorithm was implemented to identify subtype-specific marker genes for the established subtypes in the training cohort and subsequently used for subtype-prediction in a testing cohort. Multi-omics data processing Multi-omics data including transcriptomic data, copy number variation (CNV), methylation data, and somatic mutation data of the TCGA-THCA cohort were acquired from the UCSC Xena dataset portal ( https://xenabrowser.net/datapages/ ) [ 24 ]. For a standardized comparison of RNA-seq data, fragments per kilobase of exon per million mapped reads (FPKM) were converted to transcripts per million (TPM) values. Somatic variant data was called using MuTect2 and restored in a MAF file. Oncoplots of frequently mutated genes were visualized using the R package “maftools” [ 25 ]. The thresholded gene-level copy number variation (CNV) was estimated using the GISTIC2 method, and different levels represent the status of DNA copy number. DNA methylation profile was measured as beta values in tumor samples, and gene-level methylation was calculated with corresponding probes. Potential targets screening and drug sensitivity analysis The CERES scores and transcriptome profiles of cancer cell lines (CCLs) were collected from the Dependency Map (DepMap) [ 26 ], and their corresponding drug AUC values were collected from the PRISM database. The PRISM database contains AUC data of 1,492 compounds in 480 CCLs. AUC values were used as a measurement of drug sensitivity, and a higher value indicates lower sensitivity to the drug. The K-nearest neighbor (KNN) imputation was applied to impute the missing values. The inbuilt function “calcPhenotype” of the “pRRophetic” R package was used to estimate the sensitivity to different compounds and described as estimated AUC values. In addition, 13,096 drug-target pairs (comprising 4,483 drugs and 2,249 targets) were systematically summarized, and the compounds targeting RRM1 were screen for molecular docking analysis and binding affinity evaluation. Molecular docking analysis Protein structure files (PDB format) were obtained in the Protein Data Bank (RCSB PDB; https://www.rcsb.org/ ). Compound structures in SDF format were acquired from the PubChem database ( https://pubchem.ncbi.nlm.nih.gov/ ), and their lowest-energy conformations were generated using Chem3D Pro 20.0 (PerkinElmer, USA). The optimized structures were converted to MOL format and subsequently prepared as ligand PDBQT files using AutoDockTools 1.5.7 (Scripps Research, USA). The molecular docking simulation was performed with AutoDock Vina 1.2.3, with the docking grid parameters centered on the co-crystallized ligand coordinates. The binding affinity (kcal/mol) and hydrogen bonding interactions between ligands and the target protein were selected as primary evaluation metrics, where lower binding energy and more hydrogen bonds indicate higher binding stability. Post-docking analysis included 3D interaction visualization using PyMOL 2.5.2 (Schrödinger LLC) and 2D interaction diagrams generated by LigPlot + v.2.2.9 (EMBL-EBI). Furthermore, the molecular dynamics simulation is utilized to further evaluate the binding affinity of a specific compound to a specific protein, and the parameters include RMSD (Root Mean Square Deviation), Rg (Radius of Gyration), SASA (Solvent Accessible Surface Area), hydrogen bond number, and RMSF (Root Mean Square Fluctuation). Single-cell RNA-sequencing analysis Five scRNA-seq datasets from thyroid tissues were collected from the GEO database: GSE184362 (7 PTC samples) [ 27 ], GSE148673 (5 ATC samples) [ 28 ], GSE182416 (7 normal thyroid samples, unpublished), GSE193581 (10 ATC, 7 PTC tumors, and 6 normal tissues) [ 29 ], and GSE232237 (7 PTC and 5 ATC samples) [ 17 ]. Data processing and visualization were performed using Seurat (v4.4.1), Scanpy (v1.9.8), and Omicverse (v1.6.5). Quality control included filtering out low-quality cells and removing doublets and contamination using scDblFinder and scCDC. Normalization and batch correction were performed with Seurat and BBKNN, and cell clustering was annotated using classical markers. Four processed datasets (GSE148673, GSE182416, GSE193581, and GSE232237) were integrated. Genes expressed in at least 500 cells were retained, and cell subtypes were reannotated using CellTypist and SCimilarity algorithms. CellHint was used for harmonization and integration, preserving biological information while reducing batch effects. Secondary annotation was performed for thyroid epithelial cells and fibroblasts, and clustering was visualized using UMAP. Differential gene expression analysis was performed with the “FindAllMarkers” function. The Monocle algorithm was used to imitate the developmental trajectory of malignant epithelial cells. To further investigate RRM1-positive thyrocytes, the high dimensional weighted gene co-expression network analysis (hdWGCNA) was performed to identify transcriptional module associated with RRM1 + malignant epithelial cell [ 30 ]. The proportion of RRM1-positive and negative cells were calculated within each epithelial cell type (ATC, PTC, and nTFCs). Enrichment was evaluated using the ratio of observed to expected cell numbers (Ro/e), with values > 1 indicating enrichment [ 31 ]. Cell-cell communication was analyzed using “CellChat” R package [ 32 ]. Immunohistochemistry (IHC) and mIHC Specimens were collected from the Department of Pathology of Shiyan People's Hospital. The study included 5 cases of PTC with different levels of differentiation and 2 cases of ATC. Adjacent thyroid tissues at a distance of 1.0 cm from the tumor were used as control samples. The primary antibody against RRM1 was purchased from Starter Biotech Company, Wuhan, China (Catalog No: S0B0091). The secondary antibody, diaminobenzidine (DAB) chromogen, and other reagents were sourced from Roche, USA. For IHC staining, tissue sections (4-µm thick) were prepared from formalin-fixed, paraffin-embedded blocks and subjected to antigen retrieval. The primary RRM1 antibody was applied and incubated overnight at 4°C. Afterwards, a biotinylated secondary antibody and avidin-biotin-peroxidase complex were incubated for 1 h. DAB was used for color development, and slides were counterstained with hematoxylin. Two independent pathologists, blinded to clinical outcomes, evaluated the staining intensity and percentage of positively stained cells. In addition, mIHC staining was performed on ATC tissues following the protocols described in our previous study [ 33 ]. The images were acquired using the Panoramic MIDI II scanner and subsequently processed with SlideViewer (version 2.6). Dual-luciferase reporter assay A dual-luciferase reporter assay was performed using the Dual-Luciferase® Reporter System (Promega) to investigate the regulatory effect of the transcription factor c-Myc on the RRM1 promoter. The wild-type (pGL4.10-RRM1 promoter-wt) and mutant (pGL4.10-RRM1 promoter-mut) promoter fragments were cloned into the pGL4.10-basic vector, while c-Myc was expressed via pcDNA3.1(+)-c-Myc-3xFLAG, with a control plasmid pcDNA3.1(+)-MCS-3xFLAG. The Renilla luciferase plasmid pRL-CMV served as an internal reference. 293T cells were cultured in DMEM supplemented with 10% FBS at 37°C under 5% CO 2 . Cells were seeded in 96-well plates at 70% confluency and co-transfected with 0.1 µg of firefly luciferase reporter plasmid (wild-type or mutant), 0.1 µg of c-Myc or control plasmid, and 0.005 µg of pRL-CMV using Lipofectamine™ 2000. After a 6-h incubation, the medium was replaced with fresh DMEM with 10% fetal fetal bovine serum (FBS). After 48 h, the cells were lysed with Passive Lysis Buffer, and the luciferase activities were measured using the Dual-Luciferase® System. Firefly luciferase signals were normalized to Renilla luciferase activity, and relative promoter activity was calculated as the ratio of Firefly/Renilla luminescence in c-Myc-transfected cells versus control. Cell viability and colony formation assay Human ATC cell line 8505C was obtained from the Deutsche Sammlung von Mikroorganismen und Zellkulturen (DSMZ). This cell line has undergone STR authentication and is maintained in RPMI-1640 medium supplemented with 10% FBS in a 37°C humidified incubator with 5% CO₂. Gemcitabine (Catalog No: HY-17026) was sourced from MedChemExpress (USA). The cytotoxicity of gemcitabine against 8505C cells was assessed using CCK-8 (40203ES76, Yeasen Biotechnology Shanghai Co., Ltd.). 8505C cells were seeded in 96-well plates (5×10 3 cells per well). Following an overnight incubation, cells were subjected to varying concentrations of gemcitabine (2, 4, 8, 16, 32, and 64 nM) for a duration of 24 hours. Cells treated with PBS served as the control group. Subsequently, the culture medium was removed, and fresh DMEM medium supplemented with CCK-8 reagent was introduced, followed by an additional incubation period of 4 hours. The optical density (OD) values at 450 nm were then measured using a TECAN Infinite M200 Pro microplate reader. The viability of cells in each well was normalized relative to that of the PBS-treated control wells. For the colony formation assay, the 8505C cells (2000 cells per well) were plated in 6-well plates and treated with 10 nM gemcitabine for 2 weeks. Subsequently, the cells were fixed in methanol, stained with crystal violet, and photographed. Flow cytometry analysis The cell cycle and apoptosis were evaluated using flow cytometry. For cell cycle analysis, harvested cells were washed three times with ice-cold PBS. The cells were then fixed with 70% ice-cold ethanol, treated with RNase A, and stained with propidium iodide (PI). For apoptosis analysis, cells were resuspended in an annexin V binding buffer and stained with FITC-conjugated annexin V and PI for 15 minutes at room temperature. Subsequently, the stained samples were analyzed using a FACScan flow cytometer (BD Biosciences, NJ, USA) to determine the cell cycle distribution and the extent of apoptosis, and the data were analyzed using FlowJo software (version 10.8.1, FlowJo LLC) [ 34 ]. Animal study The therapeutic effect of gemcitabine on ATC were evaluated in a xenograft tumor model. Firstly, 5 × 10 6 8505C cells in 100 µL of PBS/Matrigel (1:1) were subcutaneously inoculated into the right flank of BALB/c nude mice (female, 6–8 weeks old). When tumors reached around 80 mm 3 , the mice were randomly divided into two groups: gemcitabine (100 mg/kg, iv, every other day) and PBS control groups (the same volume as above). Moreover, the tumor volumes were calculated by V = (0.5 × (shortest diameter) 2 × (longest diameter)), and body weights were recorded every other day. After 14 days of treatment (6 doses), all mice were sacrificed, and the tumors were harvested and weighed. The isolated tumors were fixed in 4% paraformaldehyde for 24 h at 37°C, then paraffin-embedded and sectioned at 4-µm thick for IHC. All the animal experiments were sanctioned by the ethics committee of the affiliated hospital of Nanjing University of Chinese Medicine (Approval number: 2022NL-KS104). Model selection using machine learning approaches Nine survival-related machine learning approaches including Lasso, Ridge, stepwise Cox, CoxBoost, random survival forest (RSF), elastic net (Enet), partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), and generalized boosted regression modeling (GBM) were incorporated to 95 combinations for model construction based on a 10-fold cross-validation framework. The predictive performance of these models was assessed using the survival C-index, and the model with the highest C-index was selected for the final model. We have referred to a previous study that integrated diverse machine learning approaches for prognostic model construction [ 35 ]. Additional bioinformatic and statistical analyses The normalized Relative Protein Expression (nRPX) data were obtained from the Human Protein Atlas (HPA) database [ 36 ]. Three BIGWIG files of MYC ChIP-seq data were obtained from the GEO repository (GSM822298, GSM3073949, and GSM3360524), and peak signals were depicted using the “EnrichedHeatmap” R package [ 37 ]. The GIs network was plotted using the “igraph” R package. CIBERSORT algorithm was used to quantify the abundance of 22 immune cells with transcriptome profiling data of bulk tissues [ 38 ]. Patient cohort was stratified into three unique immune patterns, designated as TME1, TME2, and TME3, through unsupervised hierarchical clustering according to the ssGSEA z-scores of 29 immune features [ 39 ]. The package “limma” was used to identify differentially expressed genes (DEGs) across different groups of microarray data. Principal coordinates analysis (PCoA) was used to show the dissimilarity of different groups based on the Bray-Curtis distance, and the permutational multivariate analysis of variance (PERMANOVA) test was used for significance testing. The “survminer” R package was used to identify the maximum statistics in Kaplan-Meier analysis and determine the best cut-off value for plotting survival curves. The “DynNom” R package was used to construct a dynamic nomogram model for predicting PFS in PTC patients based on their RSF-derived risk scores (RSFrs) and clinicopathological features including pathological stage, subtypes, gender, etc. The Decision Curve Analysis (DCA) was used to assess the clinical utility of predictive models or variables by comparing their net benefit across threshold probabilities. Student’s t-test or one-way analysis of variance (ANOVA) was used to analyze differences among data subject to normal distribution and variance homogeneity; otherwise, the Mann-Whitney U test or Kruskal-Wallis test was used. Statistical significance was defined as a two-sided p-value or FDR-adjusted q-value below 0.05. The analyses were conducted using GraphPad Prism and R version 4.3.2. Results Construction of the GI-perturbation network and identification of PTC subtypes Figure 1 A shows the workflow and establishment of the cancer GI-perturbation matrix. We initiated the study by retrieving 269,769 functional interactions from the Reactome database to construct the background interactome for thyroid cancer. After integrating these interactions with transcriptome data from 279 normal thyroid tissues (GTEx) and 349 PTC samples (GSE213647), we generated a rank matrix based on gene expression values and subsequently derived a delta rank matrix. By subtracting the benchmark delta rank vector (calculated from normal samples) from the delta rank of each sample, we obtained the cancer-specific GI-perturbation matrix. This matrix was used to quantify cancer-specific interaction perturbations within the same background network. Figure 1 B represents our research theme: from the GIs perturbation to AI-assisted precision medicine. A total of 1,262 gene interactions were identified between tumor and normal samples (Wilcoxon test) and had the highest variability among PTC samples (SD ranking). The relative density of GI-perturbation scores for these interactions showed distinct distributions in normal thyroid tissues and PTC samples (Fig. 1 C). A violin plot revealed a clear separation in the distribution of GI-perturbation scores between normal thyroid and PTC tissues, further highlighting the distinct perturbation patterns in cancer (Fig. 1 D). Three distinct PTC subtypes (C1, C2, and C3) were identified based on the GI-perturbation matrix. Principal component analysis (PCA) demonstrated clear dissimilarity among the three subtypes (Fig. 1 E). In addition, marker genes for each subtype were identified in a 3D volcano plot, and PCA analysis also shows clear dissimilarity among the three subtypes based on these marker genes (Fig. S1 ). We constructed subtype-specific gene interaction networks for C1, C2, and C3 by selecting the top 30 most significantly perturbed interactions in each subtype. These networks revealed unique gene interaction profiles for each subtype, providing insights into their distinct biological characteristics (Fig. 1 F-H). The NTP algorithm was applied to classify the training cohort into the three subtypes based on their specific marker genes, and the classification accuracy was evaluated (Fig. 1 I). NTP-PAM classification, an enhanced version of NTP, was also used to classify the training cohort. The algorithm showed higher accuracy compared to NTP, leading to its selection for subtype prediction in the testing cohort (Fig. 1 J). Comprehensive analysis of immune landscape, genetic mutations, and therapeutic implications across GI-subtypes The TCGA-THCA cohort was stratified into three distinct subtypes with different TME features according to the ssGSEA scores of 29 immune signatures. This classification reflects a significant correlation among the three TME subtypes and the three GI-subtypes (Fig. 2 A). A circos plot depicts the relationship between the three GI-subtypes and the infiltration of 22 immune cell types (Fig. 2 B). PCoA was performed based on the infiltration of 22 immune cell types, revealing obvious dissimilarity among the three GI-subtypes in terms of immune cell infiltration patterns (Fig. 2 C). A ternary plot reflects the infiltrating abundance of 22 immune cells across the three GI-subtypes (Fig. 2 D). The infiltration of B cell naïve in the three subtypes (C1, C2, and C3) was analyzed with a significant difference (KW-test p < 0.001). An oncoplot was generated to display the mutation frequencies of driver genes (BRAF, NRAS, HRAS, and KRAS) in the three subtypes, providing insights into the genetic alterations associated with each subtype (Fig. 2 F). Three pie charts were used to provide an overview of the mutation status of the driver genes in each subtype (Fig. 2 G). Compared to C1 and C3, C2 has the highest proportion of driver gene mutations. We further analyzed the data of expression, methylation, and CNV of various immunomodulators including co-stimulators, co-inhibitors, ligands, receptors, and cell adhesion molecules to assess their heterogeneity across the three subtypes, highlighting the distinct patterns of immunomodulator alterations among the subtypes (Fig. 2 H). Differences in tumor immunogenicity, including CTA score, non-silent mutation burden, intratumor heterogeneity, and SNV neoantigen load, were evaluated among the three subtypes. C3 exhibited the highest CTA score and intratumor heterogeneity, while C2 had the highest non-silent mutation burden and SNV neoantigen load (Fig. 2 I). The ssGSEA score patterns of 7,603 GOBPs were analyzed across the three subtypes using the mfuzz R package, identifying six distinct clusters with varying trends (Fig. 2 J). A heatmap was generated to display the top 20 GOBPs with the greatest variation among the three subtypes as the most representative biological processes for each subtype (Fig. 2 K). Based on the AUC data retrieved from the PRISM database, the sensitivity to various drugs was systematically estimated across the three GI-subtypes (Fig. 2 L). In the CDF curve comparison, GZD824 exhibited the highest sensitivity for C1 (Fig. 2 M), GDC-0879 for C2 (Fig. 2 N), and irinotecan for C3 (Fig. 2 O). A ternary plot was used to show an overview of the sensitivity to PRISM drugs across the three subtypes, with GZD824, GDC-0879, and irinotecan being highlighted. The chemical structures were also displayed (Fig. 2 P). Construction of a prognostic model based on cancer-specific GI-perturbation network A gene interaction network was constructed based on the above-identified 1,262 cancer-specific gene interactions in PTC (Fig. 3 A). Univariate Cox regression analysis was performed on the 1,262 GIs to identify those significantly associated with PFS. A total of 59 survival-related GIs were identified with a threshold of p = 0.05 (26 protective and 33 risk factors, Fig. 3 B; Table S4). Subsequently, the 59 survival-related GIs were subjected to 95 combinations to identify the most effective one for predicting PFS (Fig. 3 C). Each model was evaluated using the survival C-index. The top 10 algorithms are highlighted, with the RSF algorithm ranking first, which achieves the highest C-index of 0.772 (Fig. 3 D; Table S5). The best cutoff of RSFrs was determined using the survminer R package (Fig. 3 E). The RSFrs and clinicopathological features of the TCGA-THCA cohort were summarized in Table S6. Significant difference of PFS for the TCGA-PTC patients was observed between the high-RSFrs and low-RSFrs groups (p < 0.0001, Fig. 3 F). Furthermore, RSFrs retained its prognostic value in each subgroup (Fig. S2 ). Multivariate Cox regression analysis demonstrated RSFrs as an independent risk factor (HR = 1.282, p < 0.001; Fig. 3 G). A dynamic nomogram tool was constructed using the DynNom package, incorporating RSFrs and other clinicopathological features such as age, gender, focus, histology subtype, BRAF status, and pathological stage. This tool allows for interactive prediction of patient outcomes (Fig. 3 H). Figure 3Ia shows the output interface of the nomogram tool. Figure 3Ib displays the predicted survival probabilities for three example patients. Figure 3Ic displays the estimated survival probability curves for the three example patients. DCA analysis was performed to evaluate the clinical application of the nomogram for the prediction of the survival probabilities at different follow-up points. The results showed that the nomogram exhibited the highest net benefit compared to other clinical pathological features, indicating its superior predictive value (Fig. 3 J). Identification of RSF-related druggable targets in thyroid cancer Considering the significant association of RSFrs with recurrence in thyroid cancer patients, we aimed to identify potential therapeutic targets related to RSFrs. Firstly, a list of 2,249 druggable targets with documented drugs were collected (Table S7). We performed a correlation analysis between these 2,249 targets and RSFrs in the TCGA-PTC cohort, and 18 targets were filtered out with a correlation coefficient r > 0.5 and p < 0.001 (left panel of Fig. 4 A). Subsequently, 549 potential targets were identified as significantly upregulated genes in PTC with an FDR-q < 0.001 compared to normal (right panel of Fig. 4 A). Six candidate genes were overlapped in the intersection of the two sets of potential targets: KCNQ3, PI4K2A, IGF2R, ERBB3, LYN, and RRM1 (Fig. 4 B). A landscape distribution of CERES mean scores of all coding genes and dependent levels were shown in Fig. S3. We further evaluated the CERES scores of the six genes in 13 thyroid cancer cells from the DepMap database. RRM1 exhibited a significantly lower CERES score (-2.453 ± 0.077) compared to the other five genes, indicating a high dependency of thyroid cancer cells on RRM1 for survival (Fig. 4 C). We then explored the drug-target landscape network which contains 4,483 drugs and 2,249 targets (Fig. 4 D), and 10 drugs targeting RRM1 protein were identified: caracemide, cladribine, clofarabine, didox, enocitabine, fludarabine, gemcitabine, hydroxyurea, imexon, and triapine (Fig. 4 E). Molecular docking analysis revealed that these drugs show good affinity for the RRM1 protein, with affinity scores all below − 5 kcal/mol (Fig. 4 F). Among them, gemcitabine had the lowest affinity score (-7.727 kcal/mol), suggesting the strongest binding capacity to RRM1 among these compounds (Fig. 4 G). The free energy landscape was plotted, and detailed parameters including RMSD, Rg, SASA, hydrogen bond number, and RMSF values were estimated through the molecular dynamics simulation analysis. These results further revealed the high affinity of gemcitabine binding to the RRM1 protein (Fig. S4). RRM1, regulated by c-MYC, is correlated with malignant progression The scRNA-seq data of seven PTC samples (GSE184362) were processed using the Seurat standard pipeline. Ten major cell populations were identified in the PTC-TME using UMAP reduction (Fig. 5 A; Fig. S5A), and the expression landscape of the marker genes was shown in Fig. S5B and S5C. We then performed higher-resolution re-clustering on the Epithelia subcluster, and five subclusters labeled as C1-C5 were identified based on their marker genes (Fig. 5 B). Figure 5 C provides an overview of RRM1 expression across these five subclusters, while Fig. 5 D shows the proportion of RRM1-positive cells within each subcluster. Figure 5 E illustrates the most enriched pathways in RRM1-positive cells across the subclusters, with C1: proteolysis; C2: protein-RNA complex assembly; C3: organelle assembly; C4: protein import into nucleus; C5: RNA metabolism. Furthermore, we showed the top 10 biological processes enriched in RRM1-positive cells within each of the five subclusters (Fig. 5 F). We further investigated RRM1-positive cells using hdWGCNA. When the soft power threshold reached 9, the scale-free topology fit exceeded 0.8 (Fig. 5 G). Based on this threshold, the cells were divided into eight transcriptional modules, and the gene co-expression network of these modules is depicted (Fig. 5 I). Figure 5 J presents a bubble chart of different harmonized module eigengenes (hMEs) in RRM1-positive and RRM1-negative subclusters, with M3 being particularly notable for its negative correlation with the RRM1-positive feature. A co-expression network of hub genes within the M3 module was further depicted (Fig. 5 K). Subsequently, we performed the monocle algorithm to perform pseudotime analysis on all epithelial cells (Fig. 5 L). The simulated differentiation of RRM1 + epithelial cells showed a similar development trajectory, ultimately divided into three branches (Fig. 5 M). We calculated the ssGSEA scores of “Hallmark_MYC_targets” for all RRM1 + epithelial cells, and found a strong positive correlation with pseudotime progression (r = 0.748, p < 0.001). Three MYC ChIP-seq datasets were analyzed, and we found that MYC peak signals were significantly enriched at the genes TSS in HUVEC, SW1219, and 293T cells (Fig. 5 O). Particularly, a strong enrichment in the promoter region of RRM1 was observed (Fig. 5 P). In addition, RRM1 expression was specifically downregulated in cells treated with siRNA-MYC in the GSE126739 dataset (Fig. 5 Q) and the GSE87693 dataset (Fig. S6). Subsequently, we conducted a dual-luciferase reporter assay in 293T cells (three designed plasmids were shown in Fig. S7), and the relative luciferase intensity in the h-RRM1 (wt) + h-MYC group was significantly increased (Fig. 5 R), demonstrating that MYC has a transcriptional regulatory effect on RRM1. In an integrated scRNA-seq data of 47 samples (including Normal, ANT, PTC, and ATC), MYC activity was highest in ATC cells (Fig. 5 S), and the MYC/RRM1 double-positive cells are enriched in the ATC cells (Fig. 5 T). In the GSE33630 dataset, MYC and RRM1 showed a strong positive correlation in ATC samples only (r = 0.86; Fig. 5 U). mIHC results in two ATC samples demonstrated that RRM1 protein was highly expressed in MYC-positive cells (Fig. 5 V & W). These findings suggest that RRM1 is transcriptionally regulated by MYC. RRM1 acts as a therapeutic target in ATC We further investigated the expression of RRM1 in thyroid cancer and its potential as a therapeutic target. We found that the mRNA expression levels of RRM1 were lowest in normal tissues (or para-tumor tissues) and much higher in cancer samples, especially in PDTC and ATC (Fig. 6 A-E). Subsequently, we collected five PTC tissues with different levels of differentiation and their corresponding para-tumor normal tissues for IHC staining. RRM1 protein was widely expressed in PTCs but almost undetectable in para-tumor normal tissues (Fig. 6 F). Moreover, RRM1 was closely linked to differentiation (Fig. 6 G). We then explored the RRM1 protein levels in seven thyroid cancer cell lines in the HPA database and found the ATC-derived 8505C cell line has the highest RRM1 protein expression (Fig. 6 H). By assessing the viability of 8505C cells under different concentrations of gemcitabine, we determined the IC50 value as 14.66 nM (Fig. 6 I). Using flow cytometry, we examined the cell cycle distribution of 8505C cells in the PBS-treated (ctrl group, Fig. 6 J) and the gemcitabine-treated group (5 nM, Fig. 6 K), and we found that gemcitabine caused a significant S-phase arrest in 8505C cells (p < 0.001; Fig. 6 L). Subsequently, cell apoptosis analysis by flow cytometry was performed on 8505C cells in the ctrl group (Fig. 6 M) and the gemcitabine-treated group (10 nM, Fig. 6 N), respectively. We found that the gemcitabine-treated group (10 nM) exhibited significant apoptosis in 8505C cells (Fig. 6 O), and the LDH levels in the cell culture supernatants of the gemcitabine-treated group were significantly elevated (Fig. 6 P). Colony formation assay also confirmed a significant reduction in the number of colonies in the gemcitabine-treated group (10 nM; Fig. 6 Q). Finally, a xenograft mouse tumor model was used to verify the therapeutic effect of gemcitabine on ATC. Figure 6 R illustrates the changes in tumor volume of each nude mouse in the ctrl group and the gemcitabine-treated group over a two-week period, and the mean values ± SD of the measured tumor volumes for each group were shown in Fig. 6 S. After two weeks, the nude mice were sacrificed, and the tumors were isolated and weighed. The results showed that the tumors in the gemcitabine-treated group were significantly smaller than the control group (p < 0.001; Fig. 6 T). Figure 6 U shows the gross images of the final xenograft tumors from the ctrl group and gemcitabine-treated group. IHC staining showed a significant decrease of Ki67-positive tumor cells in the gemcitabine-treated group compared to the ctrl group (Fig. 6 V). These findings suggest that RRM1 is significantly correlated with differentiation and serves as a therapeutic target in thyroid cancer, especially ATC. The RRM1 + thyrocyte/dividing CAF interaction is unique in ATC To further investigate the characterization of RRM1 + epithelial cells among thyroid tissues and the interplay between RRM1 + epithelial cells with other cell types involved in thyroid cancer TME, we incorporated the scRNA-seq data of 47 samples (including normal, ANT, PTC, and ATC) and deeply identified subpopulations of cell types with higher resolution (Fig. 7 A). The abundance landscape of these identified cell types was also depicted among different tissue types (Fig. 7 B). The epithelial thyrocytes were extracted and identified as three main types: normal thyroid follicular cells (nTFCs), PTC epithelial cells, and ATC epithelial cells (Fig. 7 C). Furthermore, we observed that RRM1 is significantly elevated in ATC thyrocytes, rather than nTFCs and PTC thyrocytes (Fig. 7 D). In addition, RRM1 + thyrocytes were mainly enriched in cancer tissues, especially in ATC thyrocytes (26.5%), while almost few RRM1 + thyrocytes were detected in normal or ANT tissues (Fig. 7 E). The heatmap of Ro/e value illustrates the phenotype preferences of RRM1-positive and RRM1-negative malignant thyrocytes. Among malignant thyrocytes, RRM1-positive cells exhibited the strongest specificity in ATC and PTC cell types, particularly within the ATC cells. In contrast, RRM1-negative thyrocytes displayed relatively weaker specificity in ATC and PTC but showed higher specificity in nTFCs (Fig. 7 F). CellChat algorithm was used to assess the numbers of interactions and strength among these identified cell types based on the integrated scRNA-seq data, and we observed that the interaction between RRM1-positive epithelial cells and CAF4 subpopulation is prominent (Fig. 7 G). Furthermore, we respectively investigated the interactions in PTC and ATC, and we found that the interaction between RRM1-positive epithelial cells and CAF4 is unique in ATC (Fig. 7 H). A bead plot shows the quantified numbers of interactions between RRM1-positive epithelial cells with other cell types in ATC and PTC, respectively (Fig. 7 I). Subsequently, we further extracted all fibroblasts from the integrated scRNA-seq data, and five subtypes of fibroblasts were identified: four CAF types (CAF1-4) and the normal fibroblast (Fig. 7 J). The expression landscape of marker genes for each fibroblast subtype was shown in a dot-plot (Fig. 7 K), and we found that the CAF4 is similar to the dCAF (short for dividing CAF), which is defined by the expression of proliferative genes such as MKI67 and TOP2A, highlighting a more aggressive phenotype of fibroblasts in the TME. CAF4 is significantly elevated in the ATC tissues (Fig. 7 L). To validate the spatial localization and potential interaction between RRM1-positive epithelial cells and CAF4 in thyroid cancer tissues, mIHC was performed on two ATC specimens (previously measured as RRM1-low and RRM1-high). We observed spatial “co-occurrence” or “co-disappearance” of RRM1-positive epithelial cells and CAF4 (Fig. 7 M and 7 N), suggesting that their potential interactions within the TME may provide therapeutic implications for thyroid cancer, particularly for ATC. Discussion Despite the fact most patients with thyroid cancer exhibit favorable outcomes, some subsets such as high-risk PTC and aggressive ATC remain a significant clinical challenge due to poor prognosis and limited therapeutic options [ 40 , 41 ]. In this study, we applied a novel rank-based individual-specific gene-interaction (GI) perturbation approach to molecular subtyping of thyroid cancer, with a focus on identifying actionable therapeutic targets and developing personalized treatment strategies. Our findings not only provide a comprehensive understanding of the GI perturbations in thyroid cancer malignancy and progression, but also highlight the potential of RRM1 (calculated from the GI-perturbation RSF model) as a therapeutic target, particularly in high-risk PTC and ATC. The integration of GI perturbation analysis with network-based approaches has enabled us to reveal the dynamic changes in gene regulation that are often ignored in traditional static transcriptional profiling. By constructing a GI-perturbation network, we identified distinct molecular subtypes of PTC based on their GI perturbation profiles. This approach is more robust than single-gene analysis, as it captures the complex interplay between genes and their interactions, which are often disrupted in cancer. The exploration of GI perturbation network not only aids in predicting disease progression and prognosis but also provides some useful insights for tailoring targeted therapies for high-risk patients. Based on the cancer-specific GI perturbations, an RSF algorithm-derived model was established for the prediction of PFS. Furthermore, RRM1 was identified as a promising therapeutic target in high-risk thyroid cancer, as its expression is significantly associated with model risk score and survival dependency of cancer cells. RRM1, a key enzyme in DNA synthesis and repair, was found to be significantly upregulated in aggressive subtypes of thyroid cancer, especially in PDTC and ATC. High expression of RRM1 predicts worse clinical outcomes and the response to gemcitabine in a variety of malignancies, such as pancreatic cancer [ 42 ], lung cancer [ 43 , 44 ], and adrenal cancer [ 45 ]. In addition, RRM1 acts as a novel therapeutic target in multiple myeloma [ 46 ]. However, few researches reported the biological functions and clinical role of RRM1 in thyroid cancer, let alone ATC. In our study, we first noted the remarkable upregulation of RRM1 in high-risk thyroid cancer and explored the potential of targeting RRM1 as a novel therapeutic strategy for ATC. The transcriptional regulation of RRM1 by MYC further highlights its importance in tumorigenesis, as MYC is a well-known oncogene that drives cancer cell proliferation and survival [ 47 ]. We observed a strong positive correlation between MYC and RRM1 expression in ATC samples, and we also found that RRM1 was significantly downregulated in MYC siRNA-treated cells. These evidences suggest that MYC may directly regulate RRM1 expression in ATC. This hypothesis is further supported by the dual-luciferase reporter assay, which demonstrated that MYC could enhance RRM1 promoter activity and promote its transcription. The therapeutic potential of targeting RRM1 was validated through in vitro and in vivo experiments. Gemcitabine, a nucleoside analog that inhibits RRM1, effectively induced S-phase arrest and apoptosis in RRM1-high ATC cells. Moreover, gemcitabine treatment significantly attenuated tumor growth in a xenograft mouse model, further supporting the notion that RRM1 inhibition could be an effective therapeutic strategy for ATC. Considering the highly aggressive nature and poor prognosis of ATC, combined with the extremely limited and often ineffective treatment options available for ATC, our research findings hold significant promise for advancing the clinical translation of ATC-targeted therapy. Our scRNA-seq analysis revealed a unique interaction between RRM1-positive epithelial cells and a specific subset of CAF in ATC. The CAF4 subpopulation, which exhibits a proliferative phenotype similar to dividing CAFs (dCAFs) [ 48 ], was found to be significantly enriched in ATC tissues. In addition, the spatial co-localization of RRM1-positive epithelial cells and CAF4 suggests that their interplay may greatly contribute to the malignant progression of thyroid cancer. This finding accords with previous studies and highlights the importance of CAFs in shaping the tumor microenvironment (TME), which promotes cancer progression and enhances malignant phenotypes [ 49 – 51 ]. In addition, the interplay between RRM1-positive epithelial cells and CAF4 may provide a novel therapeutic direction for targeting the TME in ATC. The identification of RRM1 as a therapeutic target and the characterization of its interaction with CAF4 subtype have significant implications for precision medicine in thyroid cancer. By integrating molecular profiling data with network-based approaches, we can develop more effective and personalized treatment strategies for high-risk PTC and ATC patients. For instance, patients with high RRM1 expression could benefit from RRM1 inhibitors such as gemcitabine, while those with specific CAF4 interactions may benefit from therapies targeting the TME. Furthermore, the molecular subtyping based on GI perturbation profiles could guide the selection of targeted therapies and facilitate personalized treatment for thyroid cancer. While our study offers some valuable insights into the GI-perturbation network and molecular subtyping in thyroid cancer, several limitations should be acknowledged. First, the sample size of ATC cases in our study was relatively small, which may limit the research robustness. Future studies with a larger sample size using multi-omics analysis are needed to validate the biological role of RRM1 and its interaction with CAFs in ATC. Second, the therapeutic potential of targeting RRM1 is wished to be further validated in clinical trials to assess its efficacy and safety in ATC patients. Additionally, a deeper investigation of the interaction between RRM1-positive epithelial cells and CAF4 is warranted in future studies such as through co-culture experiments of primary cells, as intervening in this interaction could offer a new therapeutic approach for ATC. Our study highlights the significant utility of the genetic interaction perturbation network in revealing the molecular heterogeneity of thyroid cancer, and we also identified actionable targets for aggressive subtypes. Furthermore, RRM1 has emerged as an essential gene for the survival of thyroid cancer cells and exhibits a significant correlation with the progression of aggressive subtypes such as ATC. In particular, targeting RRM1 has shown promise as a therapeutic strategy for ATC. Furthermore, the unique interaction between RRM1-positive epithelial cells and CAF4 (dividing CAF subtype) highlights the importance of the tumor microenvironment in driving aggressiveness. These findings provide some useful insights for the personalized therapies for thyroid cancer, aiming to enhance patient outcomes in the aggressive subtypes. Declarations Ethics approval and consent to participate The animal study was conducted following ethical standards and received approval from the Ethics Committee of the Affiliated Hospital of Nanjing University of Chinese Medicine (approval number: 2022NL-KS104). All the clinical samples were collected in accordance with principles in the Declaration of Helsinki, and approved by the Ethics Committee of Shiyan People's Hospital (approval number: SYRMYY-2024-141). Informed consent was obtained from all patients, and all personal information was kept confidential and used for research purposes only. Declaration of Competing Interest All co-authors declared no competing interests. Funding This work was supported by the National Natural Science Foundation of China (No. 82203850 to J.S.; No. 82274155 to X.W.) and Natural Science Foundation of Jiangsu Province (No. BK20220732 to R.S.). Author Contribution Jing Sun: Conceptualization, Software, Data curation, Investigation, Writing-original draft, Funding acquisition, Formal analysis. Jiaxuan Huai: Resources, Writing-original draft, Visualization. Jianpeng Zhang: Conceptualization, Validation, Visualization. Qi Wang: Resources, Validation. Zhaokai Zhou: Data curation, Validation. Shaohui Xu: Methodology, Investigation, Data curation. Tianyu Zhao: Software, Data curation. Xin Wang: Methodology, Visualization, Investigation. Zhe Cheng: Resources, Methodology, Validation. Xuanbin Wang: Funding acquisition, Project administration, Validation, Supervision. Run Shi: Supervision, Validation, Writing-original draft, Funding acquisition, Project administration. Xiqiao Zhou: Funding acquisition, Supervision. All authors reviewed the manuscript. Acknowledgement We thank World Advanced Science Co., Ltd. for the technical assistance of multiplex immunohistochemistry, Nanjing Advanced Analysis Tech. Co., Ltd. for providing technical support and analysis in flow cytometry. Data availability All presented data and codes in this study can be available from the corresponding author upon reasonable request. References Boucai L, Zafereo M, Cabanillas ME. Thyroid Cancer: A Review. JAMA. 2024;331(5):425–35. Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, Jemal A. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2024;74(3):229–63. Ito Y, Miyauchi A, Kihara M, Fukushima M, Higashiyama T, Miya A. Overall Survival of Papillary Thyroid Carcinoma Patients: A Single-Institution Long-Term Follow-Up of 5897 Patients. World J Surg. 2018;42(3):615–22. Robenshtok E, Fish S, Bach A, Dominguez JM, Shaha A, Tuttle RM. Suspicious cervical lymph nodes detected after thyroidectomy for papillary thyroid cancer usually remain stable over years in properly selected patients. J Clin Endocrinol Metab. 2012;97(8):2706–13. Sitges-Serra A. Local recurrence of papillary thyroid cancer. Expert Rev Endocrinol Metab. 2015;10(4):349–52. Maniakas A, Dadu R, Busaidy NL, Wang JR, Ferrarotto R, Lu C, Williams MD, Gunn GB, Hofmann MC, Cote G, et al. Evaluation of Overall Survival in Patients With Anaplastic Thyroid Carcinoma, 2000–2019. JAMA Oncol. 2020;6(9):1397–404. Cancer Genome Atlas Research N. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676–90. Kim YH, Yoon SJ, Kim M, Kim HH, Song YS, Jung JW, Han D, Cho SW, Kwon SW, Park YJ. Integrative Multi-omics Analysis Reveals Different Metabolic Phenotypes Based on Molecular Characteristics in Thyroid Cancer. Clin Cancer Res. 2024;30(4):883–94. Hong S, Xie Y, Cheng Z, Li J, He W, Guo Z, Zhang Q, Peng S, He M, Yu S, et al. Distinct molecular subtypes of papillary thyroid carcinoma and gene signature with diagnostic capability. Oncogene. 2022;41(47):5121–32. Malone ER, Oliva M, Sabatini PJB, Stockley TL, Siu LL. Molecular profiling for precision cancer therapies. Genome Med. 2020;12(1):8. Li CW, Shi X, Ma B, Wang YL, Lu ZW, Liao T, Wang Y, Ji QH, Wei WJ. A 4 Gene-based Immune Signature Predicts Dedifferentiation and Immune Exhaustion in Thyroid Cancer. J Clin Endocrinol Metab. 2021;106(8):e3208–20. Chen Y, Gu Y, Hu Z, Sun X. Sample-specific perturbation of gene interactions identifies breast cancer subtypes. Brief Bioinform 2021, 22(4). Liu Z, Xu Y, Wang Y, Weng S, Xu H, Ren Y, Guo C, Liu L, Zhang Z, Han X. Immune-related interaction perturbation networks unravel biological peculiars and clinical significance of glioblastoma. Imeta. 2023;2(3):e127. Liu Z, Weng S, Dang Q, Xu H, Ren Y, Guo C, Xing Z, Sun Z, Han X. Gene interaction perturbation network deciphers a high-resolution taxonomy in colorectal cancer. Elife 2022, 11. Yan J, Risacher SL, Shen L, Saykin AJ. Network approaches to systems biology analysis of complex disease: integrative methods for multi-omics data. Brief Bioinform. 2018;19(6):1370–81. Consortium GT. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45(6):580–5. Lee SE, Park S, Yi S, Choi NR, Lim MA, Chang JW, Won HR, Kim JR, Ko HM, Chung EJ, et al. Unraveling the role of the mitochondrial one-carbon pathway in undifferentiated thyroid cancer by multi-omics analyses. Nat Commun. 2024;15(1):1163. Tarabichi M, Saiselet M, Tresallet C, Hoang C, Larsimont D, Andry G, Maenhaut C, Detours V. Revisiting the transcriptional analysis of primary tumours and associated nodal metastases with enhanced biological and statistical controls: application to thyroid cancer. Br J Cancer. 2015;112(10):1665–74. Dom G, Tarabichi M, Unger K, Thomas G, Oczko-Wojciechowska M, Bogdanova T, Jarzab B, Dumont JE, Detours V, Maenhaut C. A gene expression signature distinguishes normal tissues of sporadic and radiation-induced papillary thyroid carcinomas. Br J Cancer. 2012;107(6):994–1000. Giordano TJ, Kuick R, Thomas DG, Misek DE, Vinco M, Sanders D, Zhu Z, Ciampi R, Roh M, Shedden K, et al. Molecular classification of papillary thyroid carcinoma: distinct BRAF, RAS, and RET/PTC mutation-specific gene expression profiles discovered by DNA microarray analysis. Oncogene. 2005;24(44):6646–56. Landa I, Ibrahimpasic T, Boucai L, Sinha R, Knauf JA, Shah RH, Dogan S, Ricarte-Filho JC, Krishnamoorthy GP, Xu B, et al. Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers. J Clin Invest. 2016;126(3):1052–66. Li K, Wang F, Yang ZN, Zhang TT, Yuan YF, Zhao CX, Yeerjiang Z, Cui B, Hua F, Lv XX, et al. TRIB3 promotes MYC-associated lymphoma development through suppression of UBE3B-mediated MYC degradation. Nat Commun. 2020;11(1):6316. Satoh K, Yachida S, Sugimoto M, Oshima M, Nakagawa T, Akamoto S, Tabata S, Saitoh K, Kato K, Sato S, et al. Global metabolic reprogramming of colorectal cancer occurs at adenoma stage and is induced by MYC. Proc Natl Acad Sci U S A. 2017;114(37):E7697–706. Goldman MJ, Craft B, Hastie M, Repecka K, McDade F, Kamath A, Banerjee A, Luo Y, Rogers D, Brooks AN, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38(6):675–8. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56. Tsherniak A, Vazquez F, Montgomery PG, Weir BA, Kryukov G, Cowley GS, Gill S, Harrington WF, Pantel S, Krill-Burger JM, et al. Defining a Cancer Dependency Map. Cell. 2017;170(3):564–e576516. Pu W, Shi X, Yu P, Zhang M, Liu Z, Tan L, Han P, Wang Y, Ji D, Gan H, et al. Single-cell transcriptomic analysis of the tumor ecosystems underlying initiation and progression of papillary thyroid carcinoma. Nat Commun. 2021;12(1):6058. Gao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, Kumar T, Hu M, Sei E, Davis A, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol. 2021;39(5):599–608. Lu L, Wang JR, Henderson YC, Bai S, Yang J, Hu M, Shiau CK, Pan T, Yan Y, Tran TM et al. Anaplastic transformation in thyroid cancer revealed by single-cell transcriptomics. J Clin Invest 2023, 133(11). Morabito S, Reese F, Rahimzadeh N, Miyoshi E, Swarup V. hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data. Cell Rep Methods. 2023;3(6):100498. Gao Y, Li J, Cheng W, Diao T, Liu H, Bo Y, Liu C, Zhou W, Chen M, Zhang Y, et al. Cross-tissue human fibroblast atlas reveals myofibroblast subtypes with distinct roles in immune modulation. Cancer Cell. 2024;42(10):1764–83. e1710. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088. Shi R, Sun J, Zhou Z, Shi M, Wang X, Gao Z, Zhao T, Li M, Shu Y. Integration of multiple machine learning approaches develops a gene mutation-based classifier for accurate immunotherapy outcomes. NPJ Precis Oncol. 2025;9(1):54. Wang X, Wen S, Du X, Zhang Y, Yang X, Zou R, Feng B, Fu X, Jiang F, Zhou G, et al. SAA suppresses alpha-PD-1 induced anti-tumor immunity by driving T(H)2 polarization in lung adenocarcinoma. Cell Death Dis. 2023;14(11):718. Liu Z, Liu L, Weng S, Guo C, Dang Q, Xu H, Wang L, Lu T, Zhang Y, Sun Z, et al. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat Commun. 2022;13(1):816. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson A, Kampf C, Sjostedt E, Asplund A et al. Proteomics. Tissue-based map of the human proteome. Science 2015, 347(6220):1260419. Gu Z, Eils R, Schlesner M, Ishaque N. EnrichedHeatmap: an R/Bioconductor package for comprehensive visualization of genomic signal associations. BMC Genomics. 2018;19(1):234. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7. He Y, Jiang Z, Chen C, Wang X. Classification of triple-negative breast cancers based on Immunogenomic profiling. J Exp Clin Cancer Res. 2018;37(1):327. Bible KC, Kebebew E, Brierley J, Brito JP, Cabanillas ME, Clark TJ Jr., Di Cristofano A, Foote R, Giordano T, Kasperbauer J, et al. 2021 American Thyroid Association Guidelines for Management of Patients with Anaplastic Thyroid Cancer. Thyroid. 2021;31(3):337–86. Zhou H, Wu J, Shi L, Wang Y, Liu B. Analysis of Delayed Surgery and Clinical Outcomes in Intermediate- and High-risk Papillary Thyroid Cancer. J Clin Endocrinol Metab. 2022;107(12):3389–97. Kim R, Tan A, Lai KK, Jiang J, Wang Y, Rybicki LA, Liu X. Prognostic roles of human equilibrative transporter 1 (hENT-1) and ribonucleoside reductase subunit M1 (RRM1) in resected pancreatic cancer. Cancer. 2011;117(14):3126–34. Zhu CM, Lian XY, Bi YH, Hu CC, Liang YW, Li QS. Prognostic value of ribonucleotide reductase subunit M1 (RRM1) in non-small cell lung cancer: A meta-analysis. Clin Chim Acta. 2018;485:67–73. Bepler G, Kusmartseva I, Sharma S, Gautam A, Cantor A, Sharma A, Simon G. RRM1 modulated in vitro and in vivo efficacy of gemcitabine and platinum in non-small-cell lung cancer. J Clin Oncol. 2006;24(29):4731–7. Germano A, Rapa I, Volante M, De Francia S, Migliore C, Berruti A, Papotti M, Terzolo M. RRM1 modulates mitotane activity in adrenal cancer cells interfering with its metabolization. Mol Cell Endocrinol. 2015;401:105–10. Sagawa M, Ohguchi H, Harada T, Samur MK, Tai YT, Munshi NC, Kizaki M, Hideshima T, Anderson KC. Ribonucleotide Reductase Catalytic Subunit M1 (RRM1) as a Novel Therapeutic Target in Multiple Myeloma. Clin Cancer Res. 2017;23(17):5225–37. Dhanasekaran R, Deutzmann A, Mahauad-Fernandez WD, Hansen AS, Gouw AM, Felsher DW. The MYC oncogene - the grand orchestrator of cancer growth and immune evasion. Nat Rev Clin Oncol. 2022;19(1):23–36. Cords L, de Souza N, Bodenmiller B. Classifying cancer-associated fibroblasts-The good, the bad, and the target. Cancer Cell. 2024;42(9):1480–5. Li X, Sun Z, Peng G, Xiao Y, Guo J, Wu B, Li X, Zhou W, Li J, Li Z, et al. Single-cell RNA sequencing reveals a pro-invasive cancer-associated fibroblast subgroup associated with poor clinical outcomes in patients with gastric cancer. Theranostics. 2022;12(2):620–38. Niu N, Shen X, Wang Z, Chen Y, Weng Y, Yu F, Tang Y, Lu P, Liu M, Wang L, et al. Tumor cell-intrinsic epigenetic dysregulation shapes cancer-associated fibroblasts heterogeneity to metabolically support pancreatic cancer. Cancer Cell. 2024;42(5):869–84. e869. Liang T, Tao T, Wu K, Liu L, Xu W, Zhou D, Fang H, Ding Q, Huang G, Wu S. Cancer-Associated Fibroblast-Induced Remodeling of Tumor Microenvironment in Recurrent Bladder Cancer. Adv Sci (Weinh). 2023;10(31):e2303230. Additional Declarations No competing interests reported. Supplementary Files Supplementaryfigureslegends.docx SupplementaryTables17.xlsx Cite Share Download PDF Status: Under Revision Version 1 posted Editorial decision: Revision requested 17 Dec, 2025 Reviews received at journal 14 Dec, 2025 Reviewers agreed at journal 06 Dec, 2025 Reviews received at journal 10 Nov, 2025 Reviewers agreed at journal 03 Nov, 2025 Reviewers invited by journal 15 Oct, 2025 Editor assigned by journal 15 Oct, 2025 Submission checks completed at journal 13 Oct, 2025 First submitted to journal 11 Oct, 2025 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-7834338","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":534980738,"identity":"df434234-f4c0-47c6-bbb3-73a97ba7a108","order_by":0,"name":"Jing Sun","email":"","orcid":"","institution":"Nanjing University of Chinese Medicine","correspondingAuthor":false,"prefix":"","firstName":"Jing","middleName":"","lastName":"Sun","suffix":""},{"id":534980739,"identity":"6f8733e6-dba1-466d-b43a-6b2e4ecaffe0","order_by":1,"name":"Jiaxuan Huai","email":"","orcid":"","institution":"Nanjing University of Chinese Medicine","correspondingAuthor":false,"prefix":"","firstName":"Jiaxuan","middleName":"","lastName":"Huai","suffix":""},{"id":534980740,"identity":"45914c26-053b-4e5e-b384-c346a02d2c34","order_by":2,"name":"Jianpeng Zhang","email":"","orcid":"","institution":"The First Affiliated Hospital of Guangzhou Medical University","correspondingAuthor":false,"prefix":"","firstName":"Jianpeng","middleName":"","lastName":"Zhang","suffix":""},{"id":534980741,"identity":"218f922a-8061-4c8f-8c46-126f1f82869c","order_by":3,"name":"Wang Qi","email":"","orcid":"","institution":"Hubei University of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Wang","middleName":"","lastName":"Qi","suffix":""},{"id":534980742,"identity":"92b16750-f2f2-432a-a4c0-e1bfbeb70dd3","order_by":4,"name":"Zhaokai Zhou","email":"","orcid":"","institution":"Department of Urology, The Second Xiangya Hospital of Central South University","correspondingAuthor":false,"prefix":"","firstName":"Zhaokai","middleName":"","lastName":"Zhou","suffix":""},{"id":534980743,"identity":"b04250af-6572-461a-b450-3b0dfcbe42a8","order_by":5,"name":"Shaohui Xu","email":"","orcid":"","institution":"Soochow University FUNSOM","correspondingAuthor":false,"prefix":"","firstName":"Shaohui","middleName":"","lastName":"Xu","suffix":""},{"id":534980744,"identity":"46998d45-c48b-444b-982d-e6274ec9e94d","order_by":6,"name":"Tianyu Zhao","email":"","orcid":"","institution":"LMU University Hospital Munich","correspondingAuthor":false,"prefix":"","firstName":"Tianyu","middleName":"","lastName":"Zhao","suffix":""},{"id":534980745,"identity":"fa6b0bb7-29a6-4787-bdb8-f2dcabaa07b3","order_by":7,"name":"Xin Wang","email":"","orcid":"","institution":"The Affiliated Cancer Hospital of Nanjing Medical University","correspondingAuthor":false,"prefix":"","firstName":"Xin","middleName":"","lastName":"Wang","suffix":""},{"id":534980746,"identity":"169fd09f-37ba-4913-a3ae-3eb625e74162","order_by":8,"name":"Zhe Cheng","email":"","orcid":"","institution":"Nanjing University of Chinese Medicine","correspondingAuthor":false,"prefix":"","firstName":"Zhe","middleName":"","lastName":"Cheng","suffix":""},{"id":534980747,"identity":"4057d1d3-975e-4a7c-8f88-c4a91806cc69","order_by":9,"name":"Xuanbin Wang","email":"","orcid":"","institution":"Hubei Key Laboratory of Wudang Local Chinese Medicine Research","correspondingAuthor":false,"prefix":"","firstName":"Xuanbin","middleName":"","lastName":"Wang","suffix":""},{"id":534980748,"identity":"887bfe73-c661-490d-89ef-6d91c2409810","order_by":10,"name":"Run Shi","email":"","orcid":"","institution":"The First Affiliated Hospital of Nanjing Medical University","correspondingAuthor":false,"prefix":"","firstName":"Run","middleName":"","lastName":"Shi","suffix":""},{"id":534980749,"identity":"ecb0165a-4f12-4274-9358-fe08fe904c86","order_by":11,"name":"Xiqiao Zhou","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA7klEQVRIiWNgGAWjYDACCRBhwMDDwMB8AMIBcYnUwpZAihYw4IGrxK9FfnbzMWmeAmsZ/vaebw8s2+yiGdibt0kw1NzBqYVxzrE0yRkG6TwSZ85uN5BsS85t4DlWJsFw7BlOLcwSOWYSHwwO8zDcyN0mIdl2ILcBJMLYcBinFjaQggSgFvkbOc8gWuTf4NfCA7PF4EYOG9QWHvxaJCTSki1BfjE8c8xMQuJccm4bT1qxRcIx3FrkZyQfvM3zx9pe7njzM2mJMrvcfvbDG298qMGtBRYK0LAA+Q7ESiCkAaaF8QNhlaNgFIyCUTACAQDaK0rFA77ZcAAAAABJRU5ErkJggg==","orcid":"","institution":"Nanjing University of Chinese Medicine","correspondingAuthor":true,"prefix":"","firstName":"Xiqiao","middleName":"","lastName":"Zhou","suffix":""}],"badges":[],"createdAt":"2025-10-11 10:23:19","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7834338/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7834338/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":94759158,"identity":"ab2fca4d-69a9-41c2-a28f-ee1f23e80925","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":9935900,"visible":true,"origin":"","legend":"","description":"","filename":"Manuscript.docx","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/ebee9dfa8195fddb3d91f97b.docx"},{"id":94759149,"identity":"97cca8eb-959d-4699-a174-03ea2aa8b140","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"json","order_by":1,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":12526,"visible":true,"origin":"","legend":"","description":"","filename":"791078aee51e4c22afb7bbf9400057eb.json","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/9fc2e03149b233532a71a148.json"},{"id":94759161,"identity":"17f8628d-de4a-4783-9dd1-571574f6c097","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10036171,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTables17.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/58028437f2aa57d7f56bf520.xlsx"},{"id":94759155,"identity":"e17c8ea6-20b3-48ac-a352-c233a185aeec","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"docx","order_by":3,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1599324,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementaryfigureslegends.docx","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/2db2953e9bd84546ef4ad2d3.docx"},{"id":94759157,"identity":"f4518300-c8d1-4685-a1ef-935fc74c28af","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"xml","order_by":4,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":180058,"visible":true,"origin":"","legend":"","description":"","filename":"791078aee51e4c22afb7bbf9400057eb1enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/a9830cfadc76c6f235ec6f3d.xml"},{"id":94824442,"identity":"e1adfa06-1cce-40d9-a490-d13e94731c08","added_by":"auto","created_at":"2025-10-31 06:49:00","extension":"jpeg","order_by":5,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1351226,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/97fe535c99674a2a969339df.jpeg"},{"id":94759169,"identity":"39115302-00c1-43a8-813b-e028aaf3ed8c","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"jpeg","order_by":6,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1292832,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/cd7718da74e7b0b3143a3b81.jpeg"},{"id":94824355,"identity":"15a536e2-d646-47de-ae74-94f830ccd375","added_by":"auto","created_at":"2025-10-31 06:48:55","extension":"jpeg","order_by":7,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":885990,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage3.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/bfdd08bd4ae84a2c00069429.jpeg"},{"id":94759159,"identity":"3a6139f1-c414-47c7-a742-afeec4cbe44f","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"jpeg","order_by":8,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1621170,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage4.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/f72b46f720fc3cb1da81a695.jpeg"},{"id":94823697,"identity":"6cdebdee-9e9c-497d-8611-26ee3371a404","added_by":"auto","created_at":"2025-10-31 06:47:51","extension":"jpeg","order_by":9,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1486732,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage5.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/23cbcae6562c6b53014d3269.jpeg"},{"id":94759172,"identity":"19cb753f-a975-4e8b-907e-897b28d43224","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"jpeg","order_by":10,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1769100,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage6.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/069464ae2cc755b9e9077853.jpeg"},{"id":94759162,"identity":"5a3c4355-c1a6-47ee-927d-f48719c1b901","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"jpeg","order_by":11,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1895870,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage7.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/756b579e6b6d4d95023514e7.jpeg"},{"id":94759173,"identity":"8822c5ea-3855-4d62-a1b3-c690fc1273a3","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"png","order_by":12,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":186140,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/7deecfe8dfed132527bb1cb3.png"},{"id":94759177,"identity":"5a282686-9f54-4432-8fbf-e0e6416fd54a","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"png","order_by":13,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":171005,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/6105c2c90e6b85a9dfb90675.png"},{"id":94759167,"identity":"4d9c7ac3-1073-4b3a-8751-65dfa864a909","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"png","order_by":14,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":138160,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/89ce42388ac09d41b8524672.png"},{"id":94759170,"identity":"889fb83c-7658-45b3-bc65-86a6eeabfb70","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"png","order_by":15,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":240199,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/9eb7eb911464ab8f5d8bb54f.png"},{"id":94823755,"identity":"9ef3ad99-5b3d-48ed-9c8e-a37dfa29f5df","added_by":"auto","created_at":"2025-10-31 06:47:57","extension":"png","order_by":16,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":227937,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/04532f4259bf5fb7a84287da.png"},{"id":94759174,"identity":"19158998-1f77-4260-8476-98b2dab7f910","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"png","order_by":17,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":283459,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/696adb2e2a34f1d4197c9fed.png"},{"id":94759181,"identity":"44412ea3-39ec-4c83-a930-215d63a53efa","added_by":"auto","created_at":"2025-10-30 11:53:13","extension":"png","order_by":18,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":303891,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/afe1ea4f12af16a20e16e639.png"},{"id":94759175,"identity":"010a6cb8-a4bd-4bab-aa83-006f6201901f","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"xml","order_by":19,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":176772,"visible":true,"origin":"","legend":"","description":"","filename":"791078aee51e4c22afb7bbf9400057eb1structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/7831963eb4547e592cff8dfe.xml"},{"id":94759179,"identity":"aefc1054-9f4e-4141-8902-ccf8d9a1b07a","added_by":"auto","created_at":"2025-10-30 11:53:13","extension":"html","order_by":20,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":190654,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/81055d4403a8fc72df9e22e3.html"},{"id":94759152,"identity":"548a80ac-8d54-4247-96a5-5711904f360b","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":821639,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eConstruction of the GI-perturbation network and identification of PTC subtypes. \u003c/strong\u003e(A) Workflow for establishment of the cancer GI-perturbation matrix in this study. (B) Illustration of our research theme: from construction of GIs perturbation matrix to AI-assisted precision medicine. (C) Relative density of GI-perturbation scores for interactions between normal thyroid tissues and PTC samples. (D) Violin plot shows the distribution of GI-perturbation scores in normal thyroid and PTC tissues (p \u0026lt; 2.2e-16). (E) PCA demonstrates dissimilarity among the three GI-perturbation subtypes (C1, C2, and C3). (F-H) Subtype-specific gene interaction networks for C1, C2, and C3. (I) NTP classification and its predictive accuracy in the training cohort. (J) NTP-PAM classification shows a higher predictive accuracy in the training cohort.\u003c/p\u003e","description":"","filename":"image1.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/0c83d61a46224996f000d6f5.png"},{"id":94823217,"identity":"b985ab3a-9f8d-4ab0-ba94-a0976e0d7f4e","added_by":"auto","created_at":"2025-10-31 06:46:47","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":802052,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComprehensive analysis of immune landscape, genetic mutations, and therapeutic implications across GI-subtypes. \u003c/strong\u003e(A) A significant correlation was observed among the three TME subtypes and the three GI-subtypes. (B) A circos plot was generated to visualize the relationship between the three GI-subtypes and the infiltration of 22 immune cell types. (C) PCoA reveals strong dissimilarity among the three GI-subtypes in terms of immune cell infiltration patterns. (D) A ternary plot was used to reflect the abundance of 22 immune cell types across the three GI-subtypes. (E) The infiltration of B cell naïve in the three subtypes (C1, C2, and C3) was analyzed with a significant difference (KW-test p \u0026lt; 0.001). (F) An oncoplot was generated to display the mutation frequencies of driver genes (BRAF, NRAS, HRAS, and KRAS) in the three subtypes. (G) Three pie charts were used to provide an overview of the mutation status of the driver genes in each subtype. (H) A landscape of the immunomodulator alterations among the three subtypes. (I) Differences in tumor immunogenicity, including CTA score, non-silent mutation burden, intratumor heterogeneity, and SNV neoantigen load, were evaluated among the three subtypes. (J) The ssGSEA score patterns of 7,603 GOBPs were analyzed across the three subtypes, and six distinct clusters with varying trends were identified. (K) A heatmap was generated to display the top 20 GOBPs with the greatest variation among the three subtypes as the most representative biological processes for each subtype. (L) The sensitivity to various drugs was systematically estimated across the three GI-subtypes. (M-O) In the CDF curve comparison, GZD824 exhibited the highest sensitivity for C1, GDC-0879 for C2, and irinotecan for C3. (P) A ternary plot shows an overview of the sensitivity to PRISM drugs across the three subtypes, and the subtype-specific drugs are highlighted with their chemical structures. **p \u0026lt; 0.01; ***p \u0026lt; 0.001.\u003c/p\u003e","description":"","filename":"image2.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/0bbd851a7006b0602440d0ac.png"},{"id":94824450,"identity":"25e93cc8-ba83-4420-a49e-8985377244eb","added_by":"auto","created_at":"2025-10-31 06:49:00","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":559590,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eConstruction of a prognostic model based on cancer-specific GI-perturbation network. \u003c/strong\u003e(A) Gene interaction network was constructed based on 1,262 cancer-specific gene interactions in PTC. (B) Identification of 59 survival-related GIs (26 protective and 33 risk factors). (C) Evaluation of 95 combinations of nine machine learning algorithms for predicting PFS. (D) Performance of each model evaluated using the survival C-index. (E) Determination of the optimal cutoff of RSFrs. (F) Kaplan-Meier curve shows significant difference of PFS between RSFrs-low and RSFrs-high groups (p \u0026lt; 0.0001). (G) Multivariate Cox regression analysis demonstrated that RSFrs serves as a significant independent risk factor (HR = 1.282, p \u0026lt; 0.001). (H) A dynamic nomogram tool was developed to predict PFS based on RSFrs and other clinicopathological features. (I) Output interface of the nomogram tool with three example patients. (J) DCA analysis showed that the nomogram had the highest net benefit compared to other clinical pathological features, indicating its superior predictive value. *p \u0026lt; 0.05; ***p \u0026lt; 0.001.\u003c/p\u003e","description":"","filename":"image3.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/1ff86ba882dc5f7dc6333131.png"},{"id":94824171,"identity":"c4dc7860-d4c3-41f7-947f-04e0dff17452","added_by":"auto","created_at":"2025-10-31 06:48:36","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":908387,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIdentification of RSF-related druggable targets in thyroid cancer. \u003c/strong\u003e(A) A total of 2,249 druggable targets with documented drugs were collected, and correlation analysis was performed between their expressions and RSFrs in the TCGA-PTC cohort, and 18 targets were filtered out with a correlation coefficient r \u0026gt; 0.5 and p \u0026lt; 0.001 (left panel). Subsequently, a differential expression analysis between PTC and normal tissues, and 549 potential targets that are significantly upregulated in PTC with an FDR-q \u0026lt; 0.001 (right panel). (B) Six candidate genes were overlapped in the intersection of the two sets of potential targets: KCNQ3, PI4K2A, IGF2R, ERBB3, LYN, and RRM1. (C) CERES scores of the six genes in 13 thyroid cancer cells. In particular, RRM1 exhibited the lowest CERES score (-2.453 ± 0.077) compared to the other five genes, indicating a high dependency of thyroid cancer cells on RRM1 for survival. (D) A drug-target landscape network was constructed with 4,483 drugs and 2,249 targets. (E) 10 drugs targeting RRM1 protein were identified: caracemide, cladribine, clofarabine, didox, enocitabine, fludarabine, gemcitabine, hydroxyurea, imexon, and triapine. (F) Molecular docking analysis showed these drugs have good affinity for the RRM1 protein, with affinity scores all below -5 kcal/mol. (G) Gemcitabine exhibits the lowest affinity score (-7.727 kcal/mol), suggesting the strongest binding capacity to RRM1.\u003c/p\u003e","description":"","filename":"image4.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/67a3fd585539cc9938fbb252.png"},{"id":94823369,"identity":"c4396c00-b354-49b4-84e5-a0f7293b06e9","added_by":"auto","created_at":"2025-10-31 06:47:14","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":901728,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eRRM1, regulated by c-MYC, is correlated with malignant progression. \u003c/strong\u003e(A) The scRNA-seq data of seven PTC samples were processed, and 10 major cell populations were identified in the PTC-TME. (B) A higher-resolution re-clustering was performed on the Epithelia subcluster, and five subclusters labeled as C1-C5 were identified based on their marker genes. (C) An overview of RRM1 expression across these five subclusters. (D) The proportion of RRM1-positive cells in each subcluster. (E) The most enriched pathways in RRM1-positive cells across each subclusters. (F) Top 10 biological processes enriched in RRM1-positive cells in the five subclusters. (G) In the hdWGCNA analysis, when the soft power threshold reached 9, the scale-free topology fit exceeded 0.8. (H) RRM1-positive cells were divided into eight transcriptional modules. (I) The gene co-expression network of the 8 modules is depicted. (J) A bubble chart of different harmonized module eigengenes (hMEs) in RRM1-positive and RRM1-negative subclusters. (K) A co-expression network of hub genes within the M3 module. (L) The monocle algorithm was used to perform pseudotime analysis on all epithelial cells. (M) In the development trajectory, the RRM1\u003csup\u003e+\u003c/sup\u003e epithelial cells were finally divided into three branches. (N) A strong positive correlation between pseudotime progression and “MYC targets” activity was observed (r=0.748, p\u0026lt;0.001). (O) MYC peak signals were significantly enriched at the genes TSS in HUVEC, SW1219, and 293T cells. (P) A strong enrichment in the promoter region of RRM1 was observed. (Q) RRM1 expression was specifically downregulated in cells treated with siRNA-MYC in the GSE126739 dataset. (R) A dual-luciferase reporter assay in 293T cells demonstrated that MYC has a transcriptional regulatory effect on RRM1. (S) MYC transcription factor activity was highest in ATC cells. (T) MYC/RRM1 double-positive cells are enriched in the ATC cells. (U) MYC and RRM1 showed a strong positive correlation in ATC samples only (r = 0.86). (V \u0026amp; W) mIHC results in two ATC samples demonstrated that RRM1 protein was highly expressed in MYC-positive cells. Scale bar: 100 µm. **p \u0026lt; 0.01; ***p \u0026lt; 0.001.\u003c/p\u003e","description":"","filename":"image5.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/f5f503ab001fd384f58f9a57.png"},{"id":94823207,"identity":"fad78458-2b27-43d7-b3c3-ff160c88911f","added_by":"auto","created_at":"2025-10-31 06:46:44","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":804876,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eRRM1 acts as a therapeutic target in ATC. \u003c/strong\u003e(A-E) The mRNA expression levels of RRM1 were lowest in normal tissues (or para-tumor tissues) and much higher in cancer samples, especially in PDTC and ATC. (F) IHC staining showed that RRM1 protein was widely expressed in PTC tissues but almost undetectable in para-tumor normal tissues. Scale bar: 100 µm. (G) The expression of RRM1 was significantly correlated with differentiation. (H) The ATC-derived 8505C cell line has the highest RRM1 protein expression. (I) By assessing the viability of 8505C cells under different concentrations of gemcitabine, the IC50 value was determined as 14.66 nM. (J) The cell cycle distribution of 8505C cells in the PBS-treated group and (K) the gemcitabine-treated group (5 nM). (L) Gemcitabine caused a significant S-phase arrest in 8505C cells (p \u0026lt; 0.001). (M) Cell apoptosis analysis by flow cytometry was performed on 8505C cells in the ctrl group and (N) the gemcitabine-treated group (10 nM). (O) Gemcitabine-treated group (10 nM) exhibited significant apoptosis in 8505C cells, and (P) the LDH levels in the cell culture supernatants of the gemcitabine-treated group were significantly elevated. (Q) Colony formation assay confirmed a significant reduction in the number of colonies in the gemcitabine-treated group (10 nM). (R) The changes in tumor volume of each nude mouse in the ctrl group and the gemcitabine-treated group over a two-week period, and (S) the mean values ±SD of the measured tumor volumes for each group were shown. (T) The xenograft tumors in the gemcitabine-treated group were significantly smaller than the control group (p \u0026lt; 0.001). (U) The gross images of the final xenograft tumors from the ctrl group and gemcitabine-treated group. (V) IHC staining showed a significant decrease of Ki67-positive tumor cells in the gemcitabine-treated group compared to the ctrl group. Scale bar: 100 µm. ***p \u0026lt; 0.001.\u003c/p\u003e","description":"","filename":"image6.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/b391c0b6e4ba5aca0d378bf3.png"},{"id":94824371,"identity":"41692999-1382-408c-9e57-beca328bf364","added_by":"auto","created_at":"2025-10-31 06:48:56","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":900678,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe interaction between RRM1\u003c/strong\u003e\u003csup\u003e\u003cstrong\u003e+\u003c/strong\u003e\u003c/sup\u003e\u003cstrong\u003e thyrocyte and dividing CAF is unique in ATC. \u003c/strong\u003e(A) The scRNA-seq data of 47 samples (including normal, ANT, PTC, and ATC) were incorporated and deeply analyzed. (B) The abundance landscape of these identified cell types was also depicted among different tissue types. (C) The epithelial thyrocytes were extracted and identified as three main types: normal thyroid follicular cells (nTFCs), PTC epithelial cells, and ATC epithelial cells. (D) RRM1 is significantly elevated in ATC thyrocytes, rather than nTFCs and PTC thyrocytes. (E) RRM1\u003csup\u003e+\u003c/sup\u003e thyrocytes were mainly enriched in cancer tissues, especially in ATC thyrocytes (26.5%), while almost few RRM1\u003csup\u003e+\u003c/sup\u003e thyrocytes were detected in normal or ANT tissues. (F) The heatmap of Ro/e value illustrates the phenotype preferences of RRM1-positive and RRM1-negative malignant thyrocytes. Among malignant thyrocytes, RRM1-positive cells exhibited the strongest specificity in ATC and PTC cell types, particularly within the ATC cells. (G) CellChat algorithm was used to assess the numbers of interactions and strength among these identified cell types within thyroid cancer TME. (H) The interaction between RRM1-positive epithelial cells and CAF4 is unique in ATC. (I) A bead plot shows the quantified numbers of interactions between RRM1-positive epithelial cells with other cell types in ATC and PTC. (J) Five subtypes of fibroblasts were identified: four CAF types (CAF1-4) and the normal fibroblast. (K) The expression landscape of marker genes for each fibroblast subtype was shown in a dot-plot, and the CAF4 is similar to the dCAF (short for dividing CAF), whchi highlights a more aggressive phenotype of fibroblasts in the TME. (L) The proportion of CAF4 is significantly elevated in the ATC tissues. (M \u0026amp; N) Spatial “co-occurrence” or “co-disappearance” of RRM1-positive epithelial cells and CAF4 in ATC tissues. Scale bar: 100 µm.\u003c/p\u003e","description":"","filename":"image7.png","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/409a38efd1df584b10cfb1a3.png"},{"id":94827193,"identity":"7d3bd102-db1c-499c-a815-3c9c71db9aca","added_by":"auto","created_at":"2025-10-31 06:55:38","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":6017577,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/c170c043-e831-4059-8cbf-3b4f688d3692.pdf"},{"id":94759150,"identity":"41564f2a-dc2e-472d-9057-5fa68fe4ce2e","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":1599324,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementaryfigureslegends.docx","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/ad1cdfe4f4a4120310eee9e2.docx"},{"id":94759166,"identity":"4e5a20b0-f765-4cee-a754-b3ebe37d2823","added_by":"auto","created_at":"2025-10-30 11:53:12","extension":"xlsx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":10036171,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTables17.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7834338/v1/144a99c48225130505fe0cba.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Gene-interaction perturbation network analysis identifies distinct subtypes and actionable targets for aggressive thyroid cancer: an integrative multi-omics and machine learning approach for precise medicine","fulltext":[{"header":"Introduction","content":"\u003cp\u003eThyroid cancer is one of the most common endocrine malignancies, and its incidence has significantly increased globally in recent decades. Papillary thyroid cancer (PTC) accounts for about 84% of all thyroid cancer cases. Originating from thyroid follicular cells, the papillary, follicular (approximately 4%), and oncocytic (approximately 2%) types are categorized as well-differentiated ones. The more aggressive forms of thyroid cancer include poorly differentiated thyroid cancer (PDTC, about 5%) and anaplastic thyroid cancer (ATC, about 1%). Medullary thyroid cancer (MTC), which develops from parafollicular C cells, accounts for approximately 4% of all cases [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e]. PTC is generally an indolent tumor with a favorable prognosis, exhibiting a 10-year survival rate of approximately 93% [\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. However, a significant proportion of PTC patients still exhibit a propensity for metastasis and recurrence. Studies have shown that local or regional recurrences occur in 5 to 30% of patients with PTC, which can lead to a poor prognosis [\u003cspan additionalcitationids=\"CR4\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. In contrast, ATC is a rare but highly aggressive subtype, characterized by its undifferentiated nature and extremely poor prognosis with a 4-month median overall survival [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. Therefore, there is an urgent need to identify actionable therapeutic targets and tailor personalized approaches to offer more effective treatment options for high-risk PTC and more aggressive ATC patients.\u003c/p\u003e\u003cp\u003eComprehensive analysis of genetic alterations has become a critical tool for molecular subtyping and facilitates our understanding of the malignant phenotypes of thyroid cancer [\u003cspan additionalcitationids=\"CR8\" citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]. This approach not only aids in predicting disease progression and prognosis but also guides optimal targeted therapies, which can improve treatment efficacy and control adverse effects. Tailored treatment strategies can be designed for each patient according to the molecular profile of their tumor samples, overcoming the challenges caused by tumor heterogeneity. By integrating molecular typing with precision treatment, clinicians can develop more effective and personalized management strategies, thus improving patient outcomes [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e]. However, many existing molecular subtypes are mainly proposed based on static gene expression patterns, regardless of the dynamic changes in gene regulation [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. Gene expression can be significantly regulated during disease progression, making subtyping classification that is based on transcriptional profiles unstable and hard to reproduce. Conversely, biological networks are more stable and can stably reflect the biological status of tissues [\u003cspan additionalcitationids=\"CR13\" citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. Gene network analysis is more reliable and effective than single-gene approaches in handling high-dimensional data [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. Nonetheless, most network-based approaches only focus on individual genes and neglect gene interactions (GIs).\u003c/p\u003e\u003cp\u003eTo address this issue, we proposed gene expression rank-based individual-specific GI perturbation method that incorporates both gene nodes and their interactions. Although gene interactions are highly-conserved in normal tissues, they are significantly disrupted in tumors. The genetic interaction perturbation within the network quantifies the dynamic changes for each gene pair, allowing for a comprehensive reflection of the pathological condition at the individual level. Using this approach, we identified thyroid cancer-specific GI perturbations and GI perturbation-based molecular subtypes with distinct biological characteristics. Furthermore, GI perturbation-derived prognostic models and actionable targets were designed to guide the precise medicine for thyroid cancer patients. Overall, our findings provide a novel perspective on the gene interaction perturbation underlying thyroid cancer and offer promising insights for the development of personalized therapeutic strategies.\u003c/p\u003e"},{"header":"Materials and methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003eTranscriptome data collection and processing\u003c/h2\u003e\u003cp\u003eThe transcriptome profiles of 279 normal thyroid specimen in the GTEx database [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e] were used as the benchmark data, and the RNA-seq data of 349 PTC samples from GSE213647 (Illumina NextSeq 500 platform) were used as the training data [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. In addition, RNA-seq data of the TCGA-THCA cohort were used for clustering validation and genomic analyses. Another four microarray datasets (GSE60542 [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e], GSE33630 [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], GSE27155 [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e], and GSE76039 [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]) containing transcriptomic profiles of diverse thyroid cancer types were obtained from the GEO repository. To investigate c-Myc perturbation effects on RRM1 gene expression, transcriptome data from control siRNA (Ctrl-siRNA) and c-Myc siRNA-treated cells were obtained from GSE126739 [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e] and GSE87693 [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e]. For microarray processing, probe IDs were annotated to gene symbols using platform-specific annotation files. If a gene symbol corresponds to several probes, the probe with the highest measurement value was selected to represent the gene\u0026rsquo;s expression. A detailed summary of datasets used in this study is shown in Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e.\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003eConstruction of gene interaction-perturbation network\u003c/h3\u003e\n\u003cp\u003eA total of 269,769 functional interactions were obtained from the Reactome database and used as the initial background interactome (Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). These interactions were used to establish the background network for thyroid cancer. Genes that were not present in our cohorts (279 GTEx normal thyroid tissues and 349 PTC samples) were removed, and the remaining genes were incorporated into the background network. For each gene in the background network, its rank was calculated in each sample according to its transcriptional value. The expression matrix was transformed to a rank matrix (the lowest expression means the minimum rank, while the highest expression means the max rank). Subsequently, we generated a delta rank matrix, where rows represent interactions and columns represent samples. The delta rank was calculated by subtracting the ranks of gene pairs connected by an interaction in the background network. Gene interactions (GIs) in normal specimen were considered stable and served as the benchmark network. We calculated the benchmark delta rank vector based on the mean relative ranks of gene pairs in all normal tissues. Finally, we obtained the interaction-perturbation matrix by subtracting the benchmark delta rank vector from the delta rank of each sample. This matrix measures cancer-specific interaction perturbations in the same background network. We performed differential analysis between cancer samples and normal tissues to select the top 3000 interactions that significantly differed (Wilcoxon test) and the top 3000 interactions with the highest SD among all PTC samples. A total of 1,262 interactions were overlapped in these two sets and retained for clustering (Table S3). The consensus clustering algorithm was used to identify distinct subtypes for PTC patients based on PTC-specific GI-perturbation matrix. The nearest template prediction (NTP) algorithm was implemented to identify subtype-specific marker genes for the established subtypes in the training cohort and subsequently used for subtype-prediction in a testing cohort.\u003c/p\u003e\n\u003ch3\u003eMulti-omics data processing\u003c/h3\u003e\n\u003cp\u003eMulti-omics data including transcriptomic data, copy number variation (CNV), methylation data, and somatic mutation data of the TCGA-THCA cohort were acquired from the UCSC Xena dataset portal (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://xenabrowser.net/datapages/\u003c/span\u003e\u003cspan address=\"https://xenabrowser.net/datapages/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. For a standardized comparison of RNA-seq data, fragments per kilobase of exon per million mapped reads (FPKM) were converted to transcripts per million (TPM) values. Somatic variant data was called using MuTect2 and restored in a MAF file. Oncoplots of frequently mutated genes were visualized using the R package \u0026ldquo;maftools\u0026rdquo; [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. The thresholded gene-level copy number variation (CNV) was estimated using the GISTIC2 method, and different levels represent the status of DNA copy number. DNA methylation profile was measured as beta values in tumor samples, and gene-level methylation was calculated with corresponding probes.\u003c/p\u003e\n\u003ch3\u003ePotential targets screening and drug sensitivity analysis\u003c/h3\u003e\n\u003cp\u003eThe CERES scores and transcriptome profiles of cancer cell lines (CCLs) were collected from the Dependency Map (DepMap) [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e], and their corresponding drug AUC values were collected from the PRISM database. The PRISM database contains AUC data of 1,492 compounds in 480 CCLs. AUC values were used as a measurement of drug sensitivity, and a higher value indicates lower sensitivity to the drug. The K-nearest neighbor (KNN) imputation was applied to impute the missing values. The inbuilt function \u0026ldquo;calcPhenotype\u0026rdquo; of the \u0026ldquo;pRRophetic\u0026rdquo; R package was used to estimate the sensitivity to different compounds and described as estimated AUC values. In addition, 13,096 drug-target pairs (comprising 4,483 drugs and 2,249 targets) were systematically summarized, and the compounds targeting RRM1 were screen for molecular docking analysis and binding affinity evaluation.\u003c/p\u003e\n\u003ch3\u003eMolecular docking analysis\u003c/h3\u003e\n\u003cp\u003eProtein structure files (PDB format) were obtained in the Protein Data Bank (RCSB PDB; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.rcsb.org/\u003c/span\u003e\u003cspan address=\"https://www.rcsb.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). Compound structures in SDF format were acquired from the PubChem database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://pubchem.ncbi.nlm.nih.gov/\u003c/span\u003e\u003cspan address=\"https://pubchem.ncbi.nlm.nih.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), and their lowest-energy conformations were generated using Chem3D Pro 20.0 (PerkinElmer, USA). The optimized structures were converted to MOL format and subsequently prepared as ligand PDBQT files using AutoDockTools 1.5.7 (Scripps Research, USA). The molecular docking simulation was performed with AutoDock Vina 1.2.3, with the docking grid parameters centered on the co-crystallized ligand coordinates. The binding affinity (kcal/mol) and hydrogen bonding interactions between ligands and the target protein were selected as primary evaluation metrics, where lower binding energy and more hydrogen bonds indicate higher binding stability. Post-docking analysis included 3D interaction visualization using PyMOL 2.5.2 (Schr\u0026ouml;dinger LLC) and 2D interaction diagrams generated by LigPlot\u0026thinsp;+\u0026thinsp;v.2.2.9 (EMBL-EBI). Furthermore, the molecular dynamics simulation is utilized to further evaluate the binding affinity of a specific compound to a specific protein, and the parameters include RMSD (Root Mean Square Deviation), Rg (Radius of Gyration), SASA (Solvent Accessible Surface Area), hydrogen bond number, and RMSF (Root Mean Square Fluctuation).\u003c/p\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003eSingle-cell RNA-sequencing analysis\u003c/h2\u003e\u003cp\u003eFive scRNA-seq datasets from thyroid tissues were collected from the GEO database: GSE184362 (7 PTC samples) [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e], GSE148673 (5 ATC samples) [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e], GSE182416 (7 normal thyroid samples, unpublished), GSE193581 (10 ATC, 7 PTC tumors, and 6 normal tissues) [\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e], and GSE232237 (7 PTC and 5 ATC samples) [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. Data processing and visualization were performed using Seurat (v4.4.1), Scanpy (v1.9.8), and Omicverse (v1.6.5). Quality control included filtering out low-quality cells and removing doublets and contamination using scDblFinder and scCDC. Normalization and batch correction were performed with Seurat and BBKNN, and cell clustering was annotated using classical markers.\u003c/p\u003e\u003cp\u003eFour processed datasets (GSE148673, GSE182416, GSE193581, and GSE232237) were integrated. Genes expressed in at least 500 cells were retained, and cell subtypes were reannotated using CellTypist and SCimilarity algorithms. CellHint was used for harmonization and integration, preserving biological information while reducing batch effects. Secondary annotation was performed for thyroid epithelial cells and fibroblasts, and clustering was visualized using UMAP. Differential gene expression analysis was performed with the \u0026ldquo;FindAllMarkers\u0026rdquo; function.\u003c/p\u003e\u003cp\u003eThe Monocle algorithm was used to imitate the developmental trajectory of malignant epithelial cells. To further investigate RRM1-positive thyrocytes, the high dimensional weighted gene co-expression network analysis (hdWGCNA) was performed to identify transcriptional module associated with RRM1\u003csup\u003e+\u003c/sup\u003e malignant epithelial cell [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]. The proportion of RRM1-positive and negative cells were calculated within each epithelial cell type (ATC, PTC, and nTFCs). Enrichment was evaluated using the ratio of observed to expected cell numbers (Ro/e), with values\u0026thinsp;\u0026gt;\u0026thinsp;1 indicating enrichment [\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. Cell-cell communication was analyzed using \u0026ldquo;CellChat\u0026rdquo; R package [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e].\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003eImmunohistochemistry (IHC) and mIHC\u003c/h3\u003e\n\u003cp\u003eSpecimens were collected from the Department of Pathology of Shiyan People's Hospital. The study included 5 cases of PTC with different levels of differentiation and 2 cases of ATC. Adjacent thyroid tissues at a distance of 1.0 cm from the tumor were used as control samples. The primary antibody against RRM1 was purchased from Starter Biotech Company, Wuhan, China (Catalog No: S0B0091). The secondary antibody, diaminobenzidine (DAB) chromogen, and other reagents were sourced from Roche, USA. For IHC staining, tissue sections (4-\u0026micro;m thick) were prepared from formalin-fixed, paraffin-embedded blocks and subjected to antigen retrieval. The primary RRM1 antibody was applied and incubated overnight at 4\u0026deg;C. Afterwards, a biotinylated secondary antibody and avidin-biotin-peroxidase complex were incubated for 1 h. DAB was used for color development, and slides were counterstained with hematoxylin. Two independent pathologists, blinded to clinical outcomes, evaluated the staining intensity and percentage of positively stained cells. In addition, mIHC staining was performed on ATC tissues following the protocols described in our previous study [\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e]. The images were acquired using the Panoramic MIDI II scanner and subsequently processed with SlideViewer (version 2.6).\u003c/p\u003e\n\u003ch3\u003eDual-luciferase reporter assay\u003c/h3\u003e\n\u003cp\u003eA dual-luciferase reporter assay was performed using the Dual-Luciferase\u0026reg; Reporter System (Promega) to investigate the regulatory effect of the transcription factor c-Myc on the RRM1 promoter. The wild-type (pGL4.10-RRM1 promoter-wt) and mutant (pGL4.10-RRM1 promoter-mut) promoter fragments were cloned into the pGL4.10-basic vector, while c-Myc was expressed via pcDNA3.1(+)-c-Myc-3xFLAG, with a control plasmid pcDNA3.1(+)-MCS-3xFLAG. The Renilla luciferase plasmid pRL-CMV served as an internal reference. 293T cells were cultured in DMEM supplemented with 10% FBS at 37\u0026deg;C under 5% CO\u003csub\u003e2\u003c/sub\u003e. Cells were seeded in 96-well plates at 70% confluency and co-transfected with 0.1 \u0026micro;g of firefly luciferase reporter plasmid (wild-type or mutant), 0.1 \u0026micro;g of c-Myc or control plasmid, and 0.005 \u0026micro;g of pRL-CMV using Lipofectamine\u0026trade; 2000. After a 6-h incubation, the medium was replaced with fresh DMEM with 10% fetal fetal bovine serum (FBS). After 48 h, the cells were lysed with Passive Lysis Buffer, and the luciferase activities were measured using the Dual-Luciferase\u0026reg; System. Firefly luciferase signals were normalized to Renilla luciferase activity, and relative promoter activity was calculated as the ratio of Firefly/Renilla luminescence in c-Myc-transfected cells versus control.\u003c/p\u003e\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003eCell viability and colony formation assay\u003c/h2\u003e\u003cp\u003eHuman ATC cell line 8505C was obtained from the Deutsche Sammlung von Mikroorganismen und Zellkulturen (DSMZ). This cell line has undergone STR authentication and is maintained in RPMI-1640 medium supplemented with 10% FBS in a 37\u0026deg;C humidified incubator with 5% CO₂. Gemcitabine (Catalog No: HY-17026) was sourced from MedChemExpress (USA). The cytotoxicity of gemcitabine against 8505C cells was assessed using CCK-8 (40203ES76, Yeasen Biotechnology Shanghai Co., Ltd.). 8505C cells were seeded in 96-well plates (5\u0026times;10\u003csup\u003e3\u003c/sup\u003e cells per well). Following an overnight incubation, cells were subjected to varying concentrations of gemcitabine (2, 4, 8, 16, 32, and 64 nM) for a duration of 24 hours. Cells treated with PBS served as the control group. Subsequently, the culture medium was removed, and fresh DMEM medium supplemented with CCK-8 reagent was introduced, followed by an additional incubation period of 4 hours. The optical density (OD) values at 450 nm were then measured using a TECAN Infinite M200 Pro microplate reader. The viability of cells in each well was normalized relative to that of the PBS-treated control wells. For the colony formation assay, the 8505C cells (2000 cells per well) were plated in 6-well plates and treated with 10 nM gemcitabine for 2 weeks. Subsequently, the cells were fixed in methanol, stained with crystal violet, and photographed.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003eFlow cytometry analysis\u003c/h2\u003e\u003cp\u003eThe cell cycle and apoptosis were evaluated using flow cytometry. For cell cycle analysis, harvested cells were washed three times with ice-cold PBS. The cells were then fixed with 70% ice-cold ethanol, treated with RNase A, and stained with propidium iodide (PI). For apoptosis analysis, cells were resuspended in an annexin V binding buffer and stained with FITC-conjugated annexin V and PI for 15 minutes at room temperature. Subsequently, the stained samples were analyzed using a FACScan flow cytometer (BD Biosciences, NJ, USA) to determine the cell cycle distribution and the extent of apoptosis, and the data were analyzed using FlowJo software (version 10.8.1, FlowJo LLC) [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e].\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003eAnimal study\u003c/h2\u003e\u003cp\u003eThe therapeutic effect of gemcitabine on ATC were evaluated in a xenograft tumor model. Firstly, 5 \u0026times; 10\u003csup\u003e6\u003c/sup\u003e 8505C cells in 100 \u0026micro;L of PBS/Matrigel (1:1) were subcutaneously inoculated into the right flank of BALB/c nude mice (female, 6\u0026ndash;8 weeks old). When tumors reached around 80 mm\u003csup\u003e3\u003c/sup\u003e, the mice were randomly divided into two groups: gemcitabine (100 mg/kg, iv, every other day) and PBS control groups (the same volume as above). Moreover, the tumor volumes were calculated by V = (0.5 \u0026times; (shortest diameter)\u003csup\u003e2\u003c/sup\u003e \u0026times; (longest diameter)), and body weights were recorded every other day. After 14 days of treatment (6 doses), all mice were sacrificed, and the tumors were harvested and weighed. The isolated tumors were fixed in 4% paraformaldehyde for 24 h at 37\u0026deg;C, then paraffin-embedded and sectioned at 4-\u0026micro;m thick for IHC. All the animal experiments were sanctioned by the ethics committee of the affiliated hospital of Nanjing University of Chinese Medicine (Approval number: 2022NL-KS104).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\u003ch2\u003eModel selection using machine learning approaches\u003c/h2\u003e\u003cp\u003eNine survival-related machine learning approaches including Lasso, Ridge, stepwise Cox, CoxBoost, random survival forest (RSF), elastic net (Enet), partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), and generalized boosted regression modeling (GBM) were incorporated to 95 combinations for model construction based on a 10-fold cross-validation framework. The predictive performance of these models was assessed using the survival C-index, and the model with the highest C-index was selected for the final model. We have referred to a previous study that integrated diverse machine learning approaches for prognostic model construction [\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e].\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\u003ch2\u003eAdditional bioinformatic and statistical analyses\u003c/h2\u003e\u003cp\u003eThe normalized Relative Protein Expression (nRPX) data were obtained from the Human Protein Atlas (HPA) database [\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e]. Three BIGWIG files of MYC ChIP-seq data were obtained from the GEO repository (GSM822298, GSM3073949, and GSM3360524), and peak signals were depicted using the \u0026ldquo;EnrichedHeatmap\u0026rdquo; R package [\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e]. The GIs network was plotted using the \u0026ldquo;igraph\u0026rdquo; R package. CIBERSORT algorithm was used to quantify the abundance of 22 immune cells with transcriptome profiling data of bulk tissues [\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e]. Patient cohort was stratified into three unique immune patterns, designated as TME1, TME2, and TME3, through unsupervised hierarchical clustering according to the ssGSEA z-scores of 29 immune features [\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e]. The package \u0026ldquo;limma\u0026rdquo; was used to identify differentially expressed genes (DEGs) across different groups of microarray data. Principal coordinates analysis (PCoA) was used to show the dissimilarity of different groups based on the Bray-Curtis distance, and the permutational multivariate analysis of variance (PERMANOVA) test was used for significance testing. The \u0026ldquo;survminer\u0026rdquo; R package was used to identify the maximum statistics in Kaplan-Meier analysis and determine the best cut-off value for plotting survival curves. The \u0026ldquo;DynNom\u0026rdquo; R package was used to construct a dynamic nomogram model for predicting PFS in PTC patients based on their RSF-derived risk scores (RSFrs) and clinicopathological features including pathological stage, subtypes, gender, etc. The Decision Curve Analysis (DCA) was used to assess the clinical utility of predictive models or variables by comparing their net benefit across threshold probabilities. Student\u0026rsquo;s t-test or one-way analysis of variance (ANOVA) was used to analyze differences among data subject to normal distribution and variance homogeneity; otherwise, the Mann-Whitney U test or Kruskal-Wallis test was used. Statistical significance was defined as a two-sided p-value or FDR-adjusted q-value below 0.05. The analyses were conducted using GraphPad Prism and R version 4.3.2.\u003c/p\u003e\u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec17\" class=\"Section2\"\u003e\u003ch2\u003eConstruction of the GI-perturbation network and identification of PTC subtypes\u003c/h2\u003e\u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA shows the workflow and establishment of the cancer GI-perturbation matrix. We initiated the study by retrieving 269,769 functional interactions from the Reactome database to construct the background interactome for thyroid cancer. After integrating these interactions with transcriptome data from 279 normal thyroid tissues (GTEx) and 349 PTC samples (GSE213647), we generated a rank matrix based on gene expression values and subsequently derived a delta rank matrix. By subtracting the benchmark delta rank vector (calculated from normal samples) from the delta rank of each sample, we obtained the cancer-specific GI-perturbation matrix. This matrix was used to quantify cancer-specific interaction perturbations within the same background network. Figure\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB represents our research theme: from the GIs perturbation to AI-assisted precision medicine. A total of 1,262 gene interactions were identified between tumor and normal samples (Wilcoxon test) and had the highest variability among PTC samples (SD ranking). The relative density of GI-perturbation scores for these interactions showed distinct distributions in normal thyroid tissues and PTC samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eC). A violin plot revealed a clear separation in the distribution of GI-perturbation scores between normal thyroid and PTC tissues, further highlighting the distinct perturbation patterns in cancer (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eD). Three distinct PTC subtypes (C1, C2, and C3) were identified based on the GI-perturbation matrix. Principal component analysis (PCA) demonstrated clear dissimilarity among the three subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eE). In addition, marker genes for each subtype were identified in a 3D volcano plot, and PCA analysis also shows clear dissimilarity among the three subtypes based on these marker genes (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). We constructed subtype-specific gene interaction networks for C1, C2, and C3 by selecting the top 30 most significantly perturbed interactions in each subtype. These networks revealed unique gene interaction profiles for each subtype, providing insights into their distinct biological characteristics (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eF-H). The NTP algorithm was applied to classify the training cohort into the three subtypes based on their specific marker genes, and the classification accuracy was evaluated (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eI). NTP-PAM classification, an enhanced version of NTP, was also used to classify the training cohort. The algorithm showed higher accuracy compared to NTP, leading to its selection for subtype prediction in the testing cohort (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eJ).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec18\" class=\"Section2\"\u003e\u003ch2\u003eComprehensive analysis of immune landscape, genetic mutations, and therapeutic implications across GI-subtypes\u003c/h2\u003e\u003cp\u003eThe TCGA-THCA cohort was stratified into three distinct subtypes with different TME features according to the ssGSEA scores of 29 immune signatures. This classification reflects a significant correlation among the three TME subtypes and the three GI-subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA). A circos plot depicts the relationship between the three GI-subtypes and the infiltration of 22 immune cell types (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eB). PCoA was performed based on the infiltration of 22 immune cell types, revealing obvious dissimilarity among the three GI-subtypes in terms of immune cell infiltration patterns (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC). A ternary plot reflects the infiltrating abundance of 22 immune cells across the three GI-subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). The infiltration of B cell na\u0026iuml;ve in the three subtypes (C1, C2, and C3) was analyzed with a significant difference (KW-test p\u0026thinsp;\u0026lt;\u0026thinsp;0.001). An oncoplot was generated to display the mutation frequencies of driver genes (BRAF, NRAS, HRAS, and KRAS) in the three subtypes, providing insights into the genetic alterations associated with each subtype (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eF). Three pie charts were used to provide an overview of the mutation status of the driver genes in each subtype (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eG). Compared to C1 and C3, C2 has the highest proportion of driver gene mutations. We further analyzed the data of expression, methylation, and CNV of various immunomodulators including co-stimulators, co-inhibitors, ligands, receptors, and cell adhesion molecules to assess their heterogeneity across the three subtypes, highlighting the distinct patterns of immunomodulator alterations among the subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eH). Differences in tumor immunogenicity, including CTA score, non-silent mutation burden, intratumor heterogeneity, and SNV neoantigen load, were evaluated among the three subtypes. C3 exhibited the highest CTA score and intratumor heterogeneity, while C2 had the highest non-silent mutation burden and SNV neoantigen load (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eI). The ssGSEA score patterns of 7,603 GOBPs were analyzed across the three subtypes using the mfuzz R package, identifying six distinct clusters with varying trends (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eJ). A heatmap was generated to display the top 20 GOBPs with the greatest variation among the three subtypes as the most representative biological processes for each subtype (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eK). Based on the AUC data retrieved from the PRISM database, the sensitivity to various drugs was systematically estimated across the three GI-subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eL). In the CDF curve comparison, GZD824 exhibited the highest sensitivity for C1 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eM), GDC-0879 for C2 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eN), and irinotecan for C3 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eO). A ternary plot was used to show an overview of the sensitivity to PRISM drugs across the three subtypes, with GZD824, GDC-0879, and irinotecan being highlighted. The chemical structures were also displayed (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eP).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec19\" class=\"Section2\"\u003e\u003ch2\u003eConstruction of a prognostic model based on cancer-specific GI-perturbation network\u003c/h2\u003e\u003cp\u003eA gene interaction network was constructed based on the above-identified 1,262 cancer-specific gene interactions in PTC (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA). Univariate Cox regression analysis was performed on the 1,262 GIs to identify those significantly associated with PFS. A total of 59 survival-related GIs were identified with a threshold of p\u0026thinsp;=\u0026thinsp;0.05 (26 protective and 33 risk factors, Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eB; Table S4). Subsequently, the 59 survival-related GIs were subjected to 95 combinations to identify the most effective one for predicting PFS (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC). Each model was evaluated using the survival C-index. The top 10 algorithms are highlighted, with the RSF algorithm ranking first, which achieves the highest C-index of 0.772 (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD; Table S5). The best cutoff of RSFrs was determined using the survminer R package (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eE). The RSFrs and clinicopathological features of the TCGA-THCA cohort were summarized in Table S6. Significant difference of PFS for the TCGA-PTC patients was observed between the high-RSFrs and low-RSFrs groups (p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001, Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eF). Furthermore, RSFrs retained its prognostic value in each subgroup (Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). Multivariate Cox regression analysis demonstrated RSFrs as an independent risk factor (HR\u0026thinsp;=\u0026thinsp;1.282, p\u0026thinsp;\u0026lt;\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eG). A dynamic nomogram tool was constructed using the DynNom package, incorporating RSFrs and other clinicopathological features such as age, gender, focus, histology subtype, BRAF status, and pathological stage. This tool allows for interactive prediction of patient outcomes (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eH). Figure\u0026nbsp;3Ia shows the output interface of the nomogram tool. Figure\u0026nbsp;3Ib displays the predicted survival probabilities for three example patients. Figure\u0026nbsp;3Ic displays the estimated survival probability curves for the three example patients. DCA analysis was performed to evaluate the clinical application of the nomogram for the prediction of the survival probabilities at different follow-up points. The results showed that the nomogram exhibited the highest net benefit compared to other clinical pathological features, indicating its superior predictive value (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eJ).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec20\" class=\"Section2\"\u003e\u003ch2\u003eIdentification of RSF-related druggable targets in thyroid cancer\u003c/h2\u003e\u003cp\u003eConsidering the significant association of RSFrs with recurrence in thyroid cancer patients, we aimed to identify potential therapeutic targets related to RSFrs. Firstly, a list of 2,249 druggable targets with documented drugs were collected (Table S7). We performed a correlation analysis between these 2,249 targets and RSFrs in the TCGA-PTC cohort, and 18 targets were filtered out with a correlation coefficient r\u0026thinsp;\u0026gt;\u0026thinsp;0.5 and p\u0026thinsp;\u0026lt;\u0026thinsp;0.001 (left panel of Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA). Subsequently, 549 potential targets were identified as significantly upregulated genes in PTC with an FDR-q\u0026thinsp;\u0026lt;\u0026thinsp;0.001 compared to normal (right panel of Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA). Six candidate genes were overlapped in the intersection of the two sets of potential targets: KCNQ3, PI4K2A, IGF2R, ERBB3, LYN, and RRM1 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB). A landscape distribution of CERES mean scores of all coding genes and dependent levels were shown in Fig. S3. We further evaluated the CERES scores of the six genes in 13 thyroid cancer cells from the DepMap database. RRM1 exhibited a significantly lower CERES score (-2.453\u0026thinsp;\u0026plusmn;\u0026thinsp;0.077) compared to the other five genes, indicating a high dependency of thyroid cancer cells on RRM1 for survival (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC). We then explored the drug-target landscape network which contains 4,483 drugs and 2,249 targets (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eD), and 10 drugs targeting RRM1 protein were identified: caracemide, cladribine, clofarabine, didox, enocitabine, fludarabine, gemcitabine, hydroxyurea, imexon, and triapine (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eE). Molecular docking analysis revealed that these drugs show good affinity for the RRM1 protein, with affinity scores all below \u0026minus;\u0026thinsp;5 kcal/mol (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eF). Among them, gemcitabine had the lowest affinity score (-7.727 kcal/mol), suggesting the strongest binding capacity to RRM1 among these compounds (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eG). The free energy landscape was plotted, and detailed parameters including RMSD, Rg, SASA, hydrogen bond number, and RMSF values were estimated through the molecular dynamics simulation analysis. These results further revealed the high affinity of gemcitabine binding to the RRM1 protein (Fig. S4).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec21\" class=\"Section2\"\u003e\u003ch2\u003eRRM1, regulated by c-MYC, is correlated with malignant progression\u003c/h2\u003e\u003cp\u003eThe scRNA-seq data of seven PTC samples (GSE184362) were processed using the Seurat standard pipeline. Ten major cell populations were identified in the PTC-TME using UMAP reduction (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA; Fig. S5A), and the expression landscape of the marker genes was shown in Fig. S5B and S5C. We then performed higher-resolution re-clustering on the Epithelia subcluster, and five subclusters labeled as C1-C5 were identified based on their marker genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB). Figure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eC provides an overview of RRM1 expression across these five subclusters, while Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eD shows the proportion of RRM1-positive cells within each subcluster. Figure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eE illustrates the most enriched pathways in RRM1-positive cells across the subclusters, with C1: proteolysis; C2: protein-RNA complex assembly; C3: organelle assembly; C4: protein import into nucleus; C5: RNA metabolism. Furthermore, we showed the top 10 biological processes enriched in RRM1-positive cells within each of the five subclusters (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eF). We further investigated RRM1-positive cells using hdWGCNA. When the soft power threshold reached 9, the scale-free topology fit exceeded 0.8 (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eG). Based on this threshold, the cells were divided into eight transcriptional modules, and the gene co-expression network of these modules is depicted (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eI). Figure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eJ presents a bubble chart of different harmonized module eigengenes (hMEs) in RRM1-positive and RRM1-negative subclusters, with M3 being particularly notable for its negative correlation with the RRM1-positive feature. A co-expression network of hub genes within the M3 module was further depicted (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eK). Subsequently, we performed the monocle algorithm to perform pseudotime analysis on all epithelial cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eL). The simulated differentiation of RRM1\u003csup\u003e+\u003c/sup\u003e epithelial cells showed a similar development trajectory, ultimately divided into three branches (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eM). We calculated the ssGSEA scores of \u0026ldquo;Hallmark_MYC_targets\u0026rdquo; for all RRM1\u003csup\u003e+\u003c/sup\u003e epithelial cells, and found a strong positive correlation with pseudotime progression (r\u0026thinsp;=\u0026thinsp;0.748, p\u0026thinsp;\u0026lt;\u0026thinsp;0.001). Three MYC ChIP-seq datasets were analyzed, and we found that MYC peak signals were significantly enriched at the genes TSS in HUVEC, SW1219, and 293T cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eO). Particularly, a strong enrichment in the promoter region of RRM1 was observed (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eP). In addition, RRM1 expression was specifically downregulated in cells treated with siRNA-MYC in the GSE126739 dataset (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eQ) and the GSE87693 dataset (Fig. S6). Subsequently, we conducted a dual-luciferase reporter assay in 293T cells (three designed plasmids were shown in Fig. S7), and the relative luciferase intensity in the h-RRM1 (wt)\u0026thinsp;+\u0026thinsp;h-MYC group was significantly increased (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eR), demonstrating that MYC has a transcriptional regulatory effect on RRM1. In an integrated scRNA-seq data of 47 samples (including Normal, ANT, PTC, and ATC), MYC activity was highest in ATC cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eS), and the MYC/RRM1 double-positive cells are enriched in the ATC cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eT). In the GSE33630 dataset, MYC and RRM1 showed a strong positive correlation in ATC samples only (r\u0026thinsp;=\u0026thinsp;0.86; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eU). mIHC results in two ATC samples demonstrated that RRM1 protein was highly expressed in MYC-positive cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eV \u0026amp; W). These findings suggest that RRM1 is transcriptionally regulated by MYC.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec22\" class=\"Section2\"\u003e\u003ch2\u003eRRM1 acts as a therapeutic target in ATC\u003c/h2\u003e\u003cp\u003eWe further investigated the expression of RRM1 in thyroid cancer and its potential as a therapeutic target. We found that the mRNA expression levels of RRM1 were lowest in normal tissues (or para-tumor tissues) and much higher in cancer samples, especially in PDTC and ATC (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA-E). Subsequently, we collected five PTC tissues with different levels of differentiation and their corresponding para-tumor normal tissues for IHC staining. RRM1 protein was widely expressed in PTCs but almost undetectable in para-tumor normal tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eF). Moreover, RRM1 was closely linked to differentiation (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eG). We then explored the RRM1 protein levels in seven thyroid cancer cell lines in the HPA database and found the ATC-derived 8505C cell line has the highest RRM1 protein expression (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eH). By assessing the viability of 8505C cells under different concentrations of gemcitabine, we determined the IC50 value as 14.66 nM (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eI). Using flow cytometry, we examined the cell cycle distribution of 8505C cells in the PBS-treated (ctrl group, Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eJ) and the gemcitabine-treated group (5 nM, Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eK), and we found that gemcitabine caused a significant S-phase arrest in 8505C cells (p\u0026thinsp;\u0026lt;\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eL). Subsequently, cell apoptosis analysis by flow cytometry was performed on 8505C cells in the ctrl group (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eM) and the gemcitabine-treated group (10 nM, Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eN), respectively. We found that the gemcitabine-treated group (10 nM) exhibited significant apoptosis in 8505C cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eO), and the LDH levels in the cell culture supernatants of the gemcitabine-treated group were significantly elevated (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eP). Colony formation assay also confirmed a significant reduction in the number of colonies in the gemcitabine-treated group (10 nM; Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eQ). Finally, a xenograft mouse tumor model was used to verify the therapeutic effect of gemcitabine on ATC. Figure\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eR illustrates the changes in tumor volume of each nude mouse in the ctrl group and the gemcitabine-treated group over a two-week period, and the mean values\u0026thinsp;\u0026plusmn;\u0026thinsp;SD of the measured tumor volumes for each group were shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eS. After two weeks, the nude mice were sacrificed, and the tumors were isolated and weighed. The results showed that the tumors in the gemcitabine-treated group were significantly smaller than the control group (p\u0026thinsp;\u0026lt;\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eT). Figure\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eU shows the gross images of the final xenograft tumors from the ctrl group and gemcitabine-treated group. IHC staining showed a significant decrease of Ki67-positive tumor cells in the gemcitabine-treated group compared to the ctrl group (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eV). These findings suggest that RRM1 is significantly correlated with differentiation and serves as a therapeutic target in thyroid cancer, especially ATC.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cdiv id=\"Sec23\" class=\"Section3\"\u003e\u003ch2\u003eThe RRM1\u003csup\u003e+\u003c/sup\u003e thyrocyte/dividing CAF interaction is unique in ATC\u003c/h2\u003e\u003cp\u003eTo further investigate the characterization of RRM1\u003csup\u003e+\u003c/sup\u003e epithelial cells among thyroid tissues and the interplay between RRM1\u003csup\u003e+\u003c/sup\u003e epithelial cells with other cell types involved in thyroid cancer TME, we incorporated the scRNA-seq data of 47 samples (including normal, ANT, PTC, and ATC) and deeply identified subpopulations of cell types with higher resolution (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA). The abundance landscape of these identified cell types was also depicted among different tissue types (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eB). The epithelial thyrocytes were extracted and identified as three main types: normal thyroid follicular cells (nTFCs), PTC epithelial cells, and ATC epithelial cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eC). Furthermore, we observed that RRM1 is significantly elevated in ATC thyrocytes, rather than nTFCs and PTC thyrocytes (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eD). In addition, RRM1\u003csup\u003e+\u003c/sup\u003e thyrocytes were mainly enriched in cancer tissues, especially in ATC thyrocytes (26.5%), while almost few RRM1\u003csup\u003e+\u003c/sup\u003e thyrocytes were detected in normal or ANT tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE). The heatmap of Ro/e value illustrates the phenotype preferences of RRM1-positive and RRM1-negative malignant thyrocytes. Among malignant thyrocytes, RRM1-positive cells exhibited the strongest specificity in ATC and PTC cell types, particularly within the ATC cells. In contrast, RRM1-negative thyrocytes displayed relatively weaker specificity in ATC and PTC but showed higher specificity in nTFCs (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eF). CellChat algorithm was used to assess the numbers of interactions and strength among these identified cell types based on the integrated scRNA-seq data, and we observed that the interaction between RRM1-positive epithelial cells and CAF4 subpopulation is prominent (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eG). Furthermore, we respectively investigated the interactions in PTC and ATC, and we found that the interaction between RRM1-positive epithelial cells and CAF4 is unique in ATC (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eH). A bead plot shows the quantified numbers of interactions between RRM1-positive epithelial cells with other cell types in ATC and PTC, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eI). Subsequently, we further extracted all fibroblasts from the integrated scRNA-seq data, and five subtypes of fibroblasts were identified: four CAF types (CAF1-4) and the normal fibroblast (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eJ). The expression landscape of marker genes for each fibroblast subtype was shown in a dot-plot (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eK), and we found that the CAF4 is similar to the dCAF (short for dividing CAF), which is defined by the expression of proliferative genes such as MKI67 and TOP2A, highlighting a more aggressive phenotype of fibroblasts in the TME. CAF4 is significantly elevated in the ATC tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eL). To validate the spatial localization and potential interaction between RRM1-positive epithelial cells and CAF4 in thyroid cancer tissues, mIHC was performed on two ATC specimens (previously measured as RRM1-low and RRM1-high). We observed spatial \u0026ldquo;co-occurrence\u0026rdquo; or \u0026ldquo;co-disappearance\u0026rdquo; of RRM1-positive epithelial cells and CAF4 (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eM and \u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eN), suggesting that their potential interactions within the TME may provide therapeutic implications for thyroid cancer, particularly for ATC.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eDespite the fact most patients with thyroid cancer exhibit favorable outcomes, some subsets such as high-risk PTC and aggressive ATC remain a significant clinical challenge due to poor prognosis and limited therapeutic options [\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e, \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e]. In this study, we applied a novel rank-based individual-specific gene-interaction (GI) perturbation approach to molecular subtyping of thyroid cancer, with a focus on identifying actionable therapeutic targets and developing personalized treatment strategies. Our findings not only provide a comprehensive understanding of the GI perturbations in thyroid cancer malignancy and progression, but also highlight the potential of RRM1 (calculated from the GI-perturbation RSF model) as a therapeutic target, particularly in high-risk PTC and ATC.\u003c/p\u003e\u003cp\u003eThe integration of GI perturbation analysis with network-based approaches has enabled us to reveal the dynamic changes in gene regulation that are often ignored in traditional static transcriptional profiling. By constructing a GI-perturbation network, we identified distinct molecular subtypes of PTC based on their GI perturbation profiles. This approach is more robust than single-gene analysis, as it captures the complex interplay between genes and their interactions, which are often disrupted in cancer. The exploration of GI perturbation network not only aids in predicting disease progression and prognosis but also provides some useful insights for tailoring targeted therapies for high-risk patients.\u003c/p\u003e\u003cp\u003eBased on the cancer-specific GI perturbations, an RSF algorithm-derived model was established for the prediction of PFS. Furthermore, RRM1 was identified as a promising therapeutic target in high-risk thyroid cancer, as its expression is significantly associated with model risk score and survival dependency of cancer cells. RRM1, a key enzyme in DNA synthesis and repair, was found to be significantly upregulated in aggressive subtypes of thyroid cancer, especially in PDTC and ATC. High expression of RRM1 predicts worse clinical outcomes and the response to gemcitabine in a variety of malignancies, such as pancreatic cancer [\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e], lung cancer [\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e, \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e], and adrenal cancer [\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e]. In addition, RRM1 acts as a novel therapeutic target in multiple myeloma [\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e]. However, few researches reported the biological functions and clinical role of RRM1 in thyroid cancer, let alone ATC. In our study, we first noted the remarkable upregulation of RRM1 in high-risk thyroid cancer and explored the potential of targeting RRM1 as a novel therapeutic strategy for ATC. The transcriptional regulation of RRM1 by MYC further highlights its importance in tumorigenesis, as MYC is a well-known oncogene that drives cancer cell proliferation and survival [\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e]. We observed a strong positive correlation between MYC and RRM1 expression in ATC samples, and we also found that RRM1 was significantly downregulated in MYC siRNA-treated cells. These evidences suggest that MYC may directly regulate RRM1 expression in ATC. This hypothesis is further supported by the dual-luciferase reporter assay, which demonstrated that MYC could enhance RRM1 promoter activity and promote its transcription.\u003c/p\u003e\u003cp\u003eThe therapeutic potential of targeting RRM1 was validated through in vitro and in vivo experiments. Gemcitabine, a nucleoside analog that inhibits RRM1, effectively induced S-phase arrest and apoptosis in RRM1-high ATC cells. Moreover, gemcitabine treatment significantly attenuated tumor growth in a xenograft mouse model, further supporting the notion that RRM1 inhibition could be an effective therapeutic strategy for ATC. Considering the highly aggressive nature and poor prognosis of ATC, combined with the extremely limited and often ineffective treatment options available for ATC, our research findings hold significant promise for advancing the clinical translation of ATC-targeted therapy.\u003c/p\u003e\u003cp\u003eOur scRNA-seq analysis revealed a unique interaction between RRM1-positive epithelial cells and a specific subset of CAF in ATC. The CAF4 subpopulation, which exhibits a proliferative phenotype similar to dividing CAFs (dCAFs) [\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e], was found to be significantly enriched in ATC tissues. In addition, the spatial co-localization of RRM1-positive epithelial cells and CAF4 suggests that their interplay may greatly contribute to the malignant progression of thyroid cancer. This finding accords with previous studies and highlights the importance of CAFs in shaping the tumor microenvironment (TME), which promotes cancer progression and enhances malignant phenotypes [\u003cspan additionalcitationids=\"CR50\" citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e]. In addition, the interplay between RRM1-positive epithelial cells and CAF4 may provide a novel therapeutic direction for targeting the TME in ATC.\u003c/p\u003e\u003cp\u003eThe identification of RRM1 as a therapeutic target and the characterization of its interaction with CAF4 subtype have significant implications for precision medicine in thyroid cancer. By integrating molecular profiling data with network-based approaches, we can develop more effective and personalized treatment strategies for high-risk PTC and ATC patients. For instance, patients with high RRM1 expression could benefit from RRM1 inhibitors such as gemcitabine, while those with specific CAF4 interactions may benefit from therapies targeting the TME. Furthermore, the molecular subtyping based on GI perturbation profiles could guide the selection of targeted therapies and facilitate personalized treatment for thyroid cancer.\u003c/p\u003e\u003cp\u003eWhile our study offers some valuable insights into the GI-perturbation network and molecular subtyping in thyroid cancer, several limitations should be acknowledged. First, the sample size of ATC cases in our study was relatively small, which may limit the research robustness. Future studies with a larger sample size using multi-omics analysis are needed to validate the biological role of RRM1 and its interaction with CAFs in ATC. Second, the therapeutic potential of targeting RRM1 is wished to be further validated in clinical trials to assess its efficacy and safety in ATC patients. Additionally, a deeper investigation of the interaction between RRM1-positive epithelial cells and CAF4 is warranted in future studies such as through co-culture experiments of primary cells, as intervening in this interaction could offer a new therapeutic approach for ATC.\u003c/p\u003e\u003cp\u003eOur study highlights the significant utility of the genetic interaction perturbation network in revealing the molecular heterogeneity of thyroid cancer, and we also identified actionable targets for aggressive subtypes. Furthermore, RRM1 has emerged as an essential gene for the survival of thyroid cancer cells and exhibits a significant correlation with the progression of aggressive subtypes such as ATC. In particular, targeting RRM1 has shown promise as a therapeutic strategy for ATC. Furthermore, the unique interaction between RRM1-positive epithelial cells and CAF4 (dividing CAF subtype) highlights the importance of the tumor microenvironment in driving aggressiveness. These findings provide some useful insights for the personalized therapies for thyroid cancer, aiming to enhance patient outcomes in the aggressive subtypes.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003eEthics approval and consent to participate\u003c/h2\u003e\u003cp\u003e The animal study was conducted following ethical standards and received approval from the Ethics Committee of the Affiliated Hospital of Nanjing University of Chinese Medicine (approval number: 2022NL-KS104). All the clinical samples were collected in accordance with principles in the Declaration of Helsinki, and approved by the Ethics Committee of Shiyan People's Hospital (approval number: SYRMYY-2024-141). Informed consent was obtained from all patients, and all personal information was kept confidential and used for research purposes only.\u003c/p\u003e\u003ch2\u003eDeclaration of Competing Interest\u003c/h2\u003e\u003cp\u003eAll co-authors declared no competing interests.\u003c/p\u003e\u003ch2\u003eFunding\u003c/h2\u003e\u003cp\u003eThis work was supported by the National Natural Science Foundation of China (No. 82203850 to J.S.; No. 82274155 to X.W.) and Natural Science Foundation of Jiangsu Province (No. BK20220732 to R.S.).\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eJing Sun: Conceptualization, Software, Data curation, Investigation, Writing-original draft, Funding acquisition, Formal analysis. Jiaxuan Huai: Resources, Writing-original draft, Visualization. Jianpeng Zhang: Conceptualization, Validation, Visualization. Qi Wang: Resources, Validation. Zhaokai Zhou: Data curation, Validation. Shaohui Xu: Methodology, Investigation, Data curation. Tianyu Zhao: Software, Data curation. Xin Wang: Methodology, Visualization, Investigation. Zhe Cheng: Resources, Methodology, Validation. Xuanbin Wang: Funding acquisition, Project administration, Validation, Supervision. Run Shi: Supervision, Validation, Writing-original draft, Funding acquisition, Project administration. Xiqiao Zhou: Funding acquisition, Supervision. All authors reviewed the manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe thank World Advanced Science Co., Ltd. for the technical assistance of multiplex immunohistochemistry, Nanjing Advanced Analysis Tech. Co., Ltd. for providing technical support and analysis in flow cytometry.\u003c/p\u003e\u003ch2\u003eData availability\u003c/h2\u003e\u003cp\u003eAll presented data and codes in this study can be available from the corresponding author upon reasonable request.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eBoucai L, Zafereo M, Cabanillas ME. Thyroid Cancer: A Review. JAMA. 2024;331(5):425\u0026ndash;35.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, Jemal A. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2024;74(3):229\u0026ndash;63.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eIto Y, Miyauchi A, Kihara M, Fukushima M, Higashiyama T, Miya A. Overall Survival of Papillary Thyroid Carcinoma Patients: A Single-Institution Long-Term Follow-Up of 5897 Patients. World J Surg. 2018;42(3):615\u0026ndash;22.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eRobenshtok E, Fish S, Bach A, Dominguez JM, Shaha A, Tuttle RM. Suspicious cervical lymph nodes detected after thyroidectomy for papillary thyroid cancer usually remain stable over years in properly selected patients. J Clin Endocrinol Metab. 2012;97(8):2706\u0026ndash;13.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSitges-Serra A. Local recurrence of papillary thyroid cancer. Expert Rev Endocrinol Metab. 2015;10(4):349\u0026ndash;52.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eManiakas A, Dadu R, Busaidy NL, Wang JR, Ferrarotto R, Lu C, Williams MD, Gunn GB, Hofmann MC, Cote G, et al. Evaluation of Overall Survival in Patients With Anaplastic Thyroid Carcinoma, 2000\u0026ndash;2019. JAMA Oncol. 2020;6(9):1397\u0026ndash;404.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eCancer Genome Atlas Research N. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676\u0026ndash;90.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eKim YH, Yoon SJ, Kim M, Kim HH, Song YS, Jung JW, Han D, Cho SW, Kwon SW, Park YJ. Integrative Multi-omics Analysis Reveals Different Metabolic Phenotypes Based on Molecular Characteristics in Thyroid Cancer. Clin Cancer Res. 2024;30(4):883\u0026ndash;94.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHong S, Xie Y, Cheng Z, Li J, He W, Guo Z, Zhang Q, Peng S, He M, Yu S, et al. Distinct molecular subtypes of papillary thyroid carcinoma and gene signature with diagnostic capability. Oncogene. 2022;41(47):5121\u0026ndash;32.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMalone ER, Oliva M, Sabatini PJB, Stockley TL, Siu LL. Molecular profiling for precision cancer therapies. Genome Med. 2020;12(1):8.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi CW, Shi X, Ma B, Wang YL, Lu ZW, Liao T, Wang Y, Ji QH, Wei WJ. A 4 Gene-based Immune Signature Predicts Dedifferentiation and Immune Exhaustion in Thyroid Cancer. J Clin Endocrinol Metab. 2021;106(8):e3208\u0026ndash;20.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eChen Y, Gu Y, Hu Z, Sun X. Sample-specific perturbation of gene interactions identifies breast cancer subtypes. Brief Bioinform 2021, 22(4).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiu Z, Xu Y, Wang Y, Weng S, Xu H, Ren Y, Guo C, Liu L, Zhang Z, Han X. Immune-related interaction perturbation networks unravel biological peculiars and clinical significance of glioblastoma. Imeta. 2023;2(3):e127.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiu Z, Weng S, Dang Q, Xu H, Ren Y, Guo C, Xing Z, Sun Z, Han X. Gene interaction perturbation network deciphers a high-resolution taxonomy in colorectal cancer. Elife 2022, 11.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eYan J, Risacher SL, Shen L, Saykin AJ. Network approaches to systems biology analysis of complex disease: integrative methods for multi-omics data. Brief Bioinform. 2018;19(6):1370\u0026ndash;81.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eConsortium GT. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45(6):580\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLee SE, Park S, Yi S, Choi NR, Lim MA, Chang JW, Won HR, Kim JR, Ko HM, Chung EJ, et al. Unraveling the role of the mitochondrial one-carbon pathway in undifferentiated thyroid cancer by multi-omics analyses. Nat Commun. 2024;15(1):1163.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTarabichi M, Saiselet M, Tresallet C, Hoang C, Larsimont D, Andry G, Maenhaut C, Detours V. Revisiting the transcriptional analysis of primary tumours and associated nodal metastases with enhanced biological and statistical controls: application to thyroid cancer. Br J Cancer. 2015;112(10):1665\u0026ndash;74.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eDom G, Tarabichi M, Unger K, Thomas G, Oczko-Wojciechowska M, Bogdanova T, Jarzab B, Dumont JE, Detours V, Maenhaut C. A gene expression signature distinguishes normal tissues of sporadic and radiation-induced papillary thyroid carcinomas. Br J Cancer. 2012;107(6):994\u0026ndash;1000.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGiordano TJ, Kuick R, Thomas DG, Misek DE, Vinco M, Sanders D, Zhu Z, Ciampi R, Roh M, Shedden K, et al. Molecular classification of papillary thyroid carcinoma: distinct BRAF, RAS, and RET/PTC mutation-specific gene expression profiles discovered by DNA microarray analysis. Oncogene. 2005;24(44):6646\u0026ndash;56.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLanda I, Ibrahimpasic T, Boucai L, Sinha R, Knauf JA, Shah RH, Dogan S, Ricarte-Filho JC, Krishnamoorthy GP, Xu B, et al. Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers. J Clin Invest. 2016;126(3):1052\u0026ndash;66.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi K, Wang F, Yang ZN, Zhang TT, Yuan YF, Zhao CX, Yeerjiang Z, Cui B, Hua F, Lv XX, et al. TRIB3 promotes MYC-associated lymphoma development through suppression of UBE3B-mediated MYC degradation. Nat Commun. 2020;11(1):6316.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSatoh K, Yachida S, Sugimoto M, Oshima M, Nakagawa T, Akamoto S, Tabata S, Saitoh K, Kato K, Sato S, et al. Global metabolic reprogramming of colorectal cancer occurs at adenoma stage and is induced by MYC. Proc Natl Acad Sci U S A. 2017;114(37):E7697\u0026ndash;706.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGoldman MJ, Craft B, Hastie M, Repecka K, McDade F, Kamath A, Banerjee A, Luo Y, Rogers D, Brooks AN, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38(6):675\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747\u0026ndash;56.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTsherniak A, Vazquez F, Montgomery PG, Weir BA, Kryukov G, Cowley GS, Gill S, Harrington WF, Pantel S, Krill-Burger JM, et al. Defining a Cancer Dependency Map. Cell. 2017;170(3):564\u0026ndash;e576516.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003ePu W, Shi X, Yu P, Zhang M, Liu Z, Tan L, Han P, Wang Y, Ji D, Gan H, et al. Single-cell transcriptomic analysis of the tumor ecosystems underlying initiation and progression of papillary thyroid carcinoma. Nat Commun. 2021;12(1):6058.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, Kumar T, Hu M, Sei E, Davis A, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol. 2021;39(5):599\u0026ndash;608.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLu L, Wang JR, Henderson YC, Bai S, Yang J, Hu M, Shiau CK, Pan T, Yan Y, Tran TM et al. Anaplastic transformation in thyroid cancer revealed by single-cell transcriptomics. J Clin Invest 2023, 133(11).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMorabito S, Reese F, Rahimzadeh N, Miyoshi E, Swarup V. hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data. Cell Rep Methods. 2023;3(6):100498.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGao Y, Li J, Cheng W, Diao T, Liu H, Bo Y, Liu C, Zhou W, Chen M, Zhang Y, et al. Cross-tissue human fibroblast atlas reveals myofibroblast subtypes with distinct roles in immune modulation. Cancer Cell. 2024;42(10):1764\u0026ndash;83. e1710.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eJin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShi R, Sun J, Zhou Z, Shi M, Wang X, Gao Z, Zhao T, Li M, Shu Y. Integration of multiple machine learning approaches develops a gene mutation-based classifier for accurate immunotherapy outcomes. NPJ Precis Oncol. 2025;9(1):54.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWang X, Wen S, Du X, Zhang Y, Yang X, Zou R, Feng B, Fu X, Jiang F, Zhou G, et al. SAA suppresses alpha-PD-1 induced anti-tumor immunity by driving T(H)2 polarization in lung adenocarcinoma. Cell Death Dis. 2023;14(11):718.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiu Z, Liu L, Weng S, Guo C, Dang Q, Xu H, Wang L, Lu T, Zhang Y, Sun Z, et al. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat Commun. 2022;13(1):816.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eUhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson A, Kampf C, Sjostedt E, Asplund A et al. Proteomics. Tissue-based map of the human proteome. \u003cem\u003eScience\u003c/em\u003e 2015, 347(6220):1260419.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGu Z, Eils R, Schlesner M, Ishaque N. EnrichedHeatmap: an R/Bioconductor package for comprehensive visualization of genomic signal associations. BMC Genomics. 2018;19(1):234.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eNewman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHe Y, Jiang Z, Chen C, Wang X. Classification of triple-negative breast cancers based on Immunogenomic profiling. J Exp Clin Cancer Res. 2018;37(1):327.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBible KC, Kebebew E, Brierley J, Brito JP, Cabanillas ME, Clark TJ Jr., Di Cristofano A, Foote R, Giordano T, Kasperbauer J, et al. 2021 American Thyroid Association Guidelines for Management of Patients with Anaplastic Thyroid Cancer. Thyroid. 2021;31(3):337\u0026ndash;86.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhou H, Wu J, Shi L, Wang Y, Liu B. Analysis of Delayed Surgery and Clinical Outcomes in Intermediate- and High-risk Papillary Thyroid Cancer. J Clin Endocrinol Metab. 2022;107(12):3389\u0026ndash;97.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eKim R, Tan A, Lai KK, Jiang J, Wang Y, Rybicki LA, Liu X. Prognostic roles of human equilibrative transporter 1 (hENT-1) and ribonucleoside reductase subunit M1 (RRM1) in resected pancreatic cancer. Cancer. 2011;117(14):3126\u0026ndash;34.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhu CM, Lian XY, Bi YH, Hu CC, Liang YW, Li QS. Prognostic value of ribonucleotide reductase subunit M1 (RRM1) in non-small cell lung cancer: A meta-analysis. Clin Chim Acta. 2018;485:67\u0026ndash;73.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBepler G, Kusmartseva I, Sharma S, Gautam A, Cantor A, Sharma A, Simon G. RRM1 modulated in vitro and in vivo efficacy of gemcitabine and platinum in non-small-cell lung cancer. J Clin Oncol. 2006;24(29):4731\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGermano A, Rapa I, Volante M, De Francia S, Migliore C, Berruti A, Papotti M, Terzolo M. RRM1 modulates mitotane activity in adrenal cancer cells interfering with its metabolization. Mol Cell Endocrinol. 2015;401:105\u0026ndash;10.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSagawa M, Ohguchi H, Harada T, Samur MK, Tai YT, Munshi NC, Kizaki M, Hideshima T, Anderson KC. Ribonucleotide Reductase Catalytic Subunit M1 (RRM1) as a Novel Therapeutic Target in Multiple Myeloma. Clin Cancer Res. 2017;23(17):5225\u0026ndash;37.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eDhanasekaran R, Deutzmann A, Mahauad-Fernandez WD, Hansen AS, Gouw AM, Felsher DW. The MYC oncogene - the grand orchestrator of cancer growth and immune evasion. Nat Rev Clin Oncol. 2022;19(1):23\u0026ndash;36.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eCords L, de Souza N, Bodenmiller B. Classifying cancer-associated fibroblasts-The good, the bad, and the target. Cancer Cell. 2024;42(9):1480\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi X, Sun Z, Peng G, Xiao Y, Guo J, Wu B, Li X, Zhou W, Li J, Li Z, et al. Single-cell RNA sequencing reveals a pro-invasive cancer-associated fibroblast subgroup associated with poor clinical outcomes in patients with gastric cancer. Theranostics. 2022;12(2):620\u0026ndash;38.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eNiu N, Shen X, Wang Z, Chen Y, Weng Y, Yu F, Tang Y, Lu P, Liu M, Wang L, et al. Tumor cell-intrinsic epigenetic dysregulation shapes cancer-associated fibroblasts heterogeneity to metabolically support pancreatic cancer. Cancer Cell. 2024;42(5):869\u0026ndash;84. e869.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiang T, Tao T, Wu K, Liu L, Xu W, Zhou D, Fang H, Ding Q, Huang G, Wu S. Cancer-Associated Fibroblast-Induced Remodeling of Tumor Microenvironment in Recurrent Bladder Cancer. Adv Sci (Weinh). 2023;10(31):e2303230.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"cancer-cell-international","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"ccin","sideBox":"Learn more about [Cancer Cell International](http://cancerci.biomedcentral.com/)","snPcode":"12935","submissionUrl":"https://submission.nature.com/new-submission/12935/3","title":"Cancer Cell International","twitterHandle":"@OncoBioMed","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"Aggressive thyroid cancer, Genetic interaction perturbation, Machine learning, Molecular subtyping, Drug repositioning, Precise medicine","lastPublishedDoi":"10.21203/rs.3.rs-7834338/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7834338/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eHigh-risk thyroid cancer such as aggressive papillary thyroid cancer (PTC) and anaplastic thyroid cancer (ATC) poses significant clinical challenges due to their tumor heterogeneity and limited therapeutic strategies. We constructed a genetic interaction (GI) perturbation network based on functional interactions from the Reactome database and identified distinct molecular subtypes of thyroid cancer with unique biological characteristics. Furthermore, a GI-based machine learning-assisted prognostic model was established for predicting progression of thyroid cancer patients, and RRM1 was identified as an essential gene that is significantly correlated with progression of thyroid cancer. Drug repositioning and molecular docking analyses demonstrated that gemcitabine targeting RRM1 can serve as a promising therapeutic strategy for ATC. Subsequent in vitro and in vivo experiments validated the therapeutic potential of gemcitabine in targeting RRM1, with significant efficacy in inducing ATC cell apoptosis and inhibiting tumor growth in a xenograft mouse model. Using single-cell RNA-sequencing (scRNA-seq) analysis and multiplex immunohistochemistry (mIHC), we revealed a unique interaction between RRM1-positive epithelial cells and a specific subset of cancer-associated fibroblasts (CAF) in ATC, highlighting the importance of the tumor microenvironment in driving aggressiveness. In summary, our findings highlight the critical role of dynamic gene-interaction network in thyroid cancer progression and offer a promising therapeutic strategy by targeting RRM1 for patients with high-risk thyroid cancer.\u003c/p\u003e","manuscriptTitle":"Gene-interaction perturbation network analysis identifies distinct subtypes and actionable targets for aggressive thyroid cancer: an integrative multi-omics and machine learning approach for precise medicine","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-10-30 11:53:07","doi":"10.21203/rs.3.rs-7834338/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-12-17T06:08:03+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-12-15T02:38:00+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"157617268388582274717412880591099268805","date":"2025-12-07T03:53:19+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-11-10T07:11:10+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"214763113088402657027598873437357075483","date":"2025-11-03T06:27:19+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-10-16T00:46:57+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-10-16T00:43:34+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-10-13T09:20:44+00:00","index":"","fulltext":""},{"type":"submitted","content":"Cancer Cell International","date":"2025-10-11T10:18:28+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"cancer-cell-international","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"ccin","sideBox":"Learn more about [Cancer Cell International](http://cancerci.biomedcentral.com/)","snPcode":"12935","submissionUrl":"https://submission.nature.com/new-submission/12935/3","title":"Cancer Cell International","twitterHandle":"@OncoBioMed","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"8c7cd736-a0d6-410d-9a80-22633f1b39a5","owner":[],"postedDate":"October 30th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"in-revision","subjectAreas":[],"tags":[],"updatedAt":"2025-12-17T06:23:10+00:00","versionOfRecord":[],"versionCreatedAt":"2025-10-30 11:53:07","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7834338","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7834338","identity":"rs-7834338","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}