Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Oncol., 02 September 2025

Sec. Breast Cancer

Volume 15 - 2025 | https://doi.org/10.3389/fonc.2025.1641195

Multi-disease transcriptomic analysis of sex hormone genes reveals a novel prognostic model for thyroid cancer with breast cancer correlations

Lixue Qiao&#x;Lixue Qiao1†Hao Li&#x;Hao Li2†Keyu YinKeyu Yin3Runsheng Ma,,Runsheng Ma1,4,5Yifei Zhang,Yifei Zhang4,5Yue GuoYue Guo1Detao Yin,,*Detao Yin1,4,5*
  • 1Thyroid Surgery Department, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
  • 2Department of Hepatobiliary and Pancreatic Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
  • 3School of Basic Medical Sciences, Lanzhou University, Lanzhou, Gansu, China
  • 4Engineering Research Center of Multidisciplinary Diagnosis and Treatment of Thyroid Cancer of Henan Province, Zhengzhou,, China
  • 5Key Medicine Laboratory of Thyroid Cancer of Henan Province, Zhengzhou, China

Background: There is a potential bidirectional pathogenicity between thyroid and breast cancers. The association between sex hormones and two types of malignant tumors has emerged as a topic of intense academic debate in recent years. However, the role of sex hormone metabolism-related genes in thyroid cancer still needs to be further explored.

Methods: We obtained thyroid and breast cancer transcriptome data from the TCGA database and sex hormone metabolism-related gene sets from the MSigDB database, thus screening for sex hormone metabolism-related genes linked to the two malignant tumors. Univariate cox regression analysis was used for the screening of disease-free survival (DFS)-associated genes. The TCGA-THCA patients were classified as two categories via a consistent clustering algorithm, and the differential genes between the two categories were subsequently screened. A sex hormone metabolism-related prognostic model (TBSMRPM) of thyroid cancer versus breast cancer consisting of 10 genes was developed by Cox regression analyses and least absolute shrinkage with selection operator (LASSO) cox regression analysis. Finally, we performed clinicopathological subgroup analyses to analyze the correlation between TBSMRPM and clinical characteristics, immune infiltration, tumor mutation burden (TMB), and chemosensitivity, and verified the expression of TBSMRPM signature genes by qRT-PCR.

Results: We identified 2 clusters correlated with sex hormone metabolism, and screened 10 prognostic differential genes related to thyroid cancer, breast cancer and sex hormone metabolism. After establishing the two risk groups for thyroid cancer originated from TBSMRPM, the results showed that the high-risk group exhibited the shorter DFS (P<0.05). In further clinical stratification analysis, immune infiltration analysis, TMB and drug sensitivity analysis, the two TBSMRPM groups showed significant differences. The qRT-PCR results showed that C2CD4A, CERS1, MMP9, SLC5A1, HORMAD2 were highly expressed in the IHH4, KTC-1, and TPC-1 cell lines, while SLITRK2, ARHGEF37, PLP1, RNF223, and F3 were lowly expressed.

Conclusion: The TBSMRPM established in this study has a certain value for the prognosis of thyroid cancer and contributes to refine clinicians’ treatment protocols.

1 Introduction

Thyroid cancer (THCA) and breast cancer (BRCA), as common hormone-dependent malignancies, their bidirectional pathogenic association and the regulatory role of sex hormones therein have become a focus of oncology research. Epidemiological data clearly demonstrate a significant interaction between these two organ malignancies: Nielsen et al. confirmed that individuals with a history of thyroid cancer have a 1.32-fold higher risk of developing breast cancer compared to the general population, while breast cancer patients have a 1.55-fold increased risk of subsequent thyroid cancer (1). This bidirectional risk suggests that, as hormone-dependent organs, the malignant progression of the thyroid and mammary gland is not an isolated event, and their development is influenced by a complex regulatory network (16). Currently, research on their association has involved multiple dimensions, including genetic alterations, hypothalamic-pituitary axis regulation, metabolic abnormalities (e.g., diabetes, obesity), and surveillance bias. Among these, hormone-related mechanisms are particularly critical—various hormone-related factors such as thyroid hormones, sex hormones (with estrogen as the core), endocrine-disrupting chemicals, and adipokines have been confirmed to be involved in the development of both tumors (2). Sex hormones, in particular, have been identified as core regulatory factors due to their direct role in regulating cell proliferation, differentiation, and malignant transformation in thyroid and breast tissues (2). Therefore, systematically elucidating the roles and mechanisms of sex hormone metabolism-related genes associated with both cancers is of great significance for clarifying cross-cancer associations and optimizing clinical diagnosis and treatment.

The metabolic balance of sex hormones is precisely regulated by a network of genes encoding metabolic enzymes and receptors, collectively referred to as sex hormone metabolism-related genes (SMRGs). To date, the roles of SMRGs in breast cancer and thyroid cancer have been partially elucidated. In the field of breast cancer, the clinical value of SMRGs has been verified: for example, aromatase encoded by CYP19A1 is a key rate-limiting enzyme in estrogen synthesis, and its inhibitors have become standard therapeutic agents for ER-positive breast cancer by reducing estrogen production (7); both two-sample Mendelian randomization studies and epidemiological data have shown that elevated levels of total testosterone and bioavailable testosterone increase the risk of ER+ breast cancer, with this association being more pronounced in postmenopausal women (8, 9). In thyroid cancer research, the roles of SMRGs have also been gradually revealed, with the regulatory mechanisms of estrogen and its receptors receiving the most attention: estrogen can enhance the proliferation, migration, and invasion abilities of papillary thyroid carcinoma (PTC) through the ERα/KRT19 signaling axis (10); ERβ is highly expressed in PTC stem cells (PTCSCs), and its knockout can reduce the expression of stemness-related factors, decrease the ALDH+ cell population, and inhibit tumor sphere formation and growth (11); in addition, ERα36, GRP78, and GRP94 are upregulated in primary PTC tissues, and their expression levels can directly affect the malignant phenotype of PTC-derived BCPAP cells (12). Beyond estrogen, androgen receptors (AR), progesterone receptors, and prolactin receptors are also expressed in PTC lesions (13, 14), and AR activation can exert antiproliferative effects (e.g., inducing cellular senescence) in PTC cell models (15, 16). These findings collectively suggest that SMRGs hold promise as novel molecular targets for diagnostic prediction, endocrine therapy, and prognostic evaluation in both cancers.

Despite certain progress in SMRG research, significant limitations remain. Existing studies have shown that Fu et al. identified the COMP gene, which is highly expressed in both breast cancer and thyroid cancer, and confirmed that its overexpression can promote the occurrence and progression of both cancers through the estrogen signaling pathway (17); Jin et al. found that the ratio of ESR1 to ESR can serve as a prognostic marker for predicting survival in female PTC patients and has the potential to be a therapeutic target (18); Zhang et al. constructed a model by screening estrogen-related differential genes associated with THCA and found that this model is closely related to immune infiltration and therapeutic response in THCA (19). However, current research still has three shortcomings: first, the exploration of SMRGs is insufficiently systematic, with most focusing on estrogen metabolism-related genes in PTC, while neglecting the contribution of other sex hormone metabolism genes (such as those involved in androgen and progesterone metabolism) to the development of thyroid cancer; second, most studies are limited to thyroid cancer alone, failing to incorporate its bidirectional risk characteristics with breast cancer, making it difficult to reflect the cross-regulatory association of hormones between the two; third, some studies only focus on the predictive efficacy of SMRGs in PTC, lacking in-depth analysis of their underlying pathogenic mechanisms (e.g., associations with the tumor microenvironment and therapeutic sensitivity). Therefore, there is an urgent need to systematically explore the prognostic efficacy and potential mechanisms of SMRGs associated with both cancers in PTC from a cross-cancer perspective.

To address the above research gaps, this study integrates transcriptomic data of THCA and BRCA from The Cancer Genome Atlas (TCGA) database, and screens SMRGs directly or potentially associated with both cancers through cross-analysis, providing a molecular basis for deciphering cross-cancer metabolic associations. On this basis, consensus clustering and LASSO Cox regression methods are used to construct a thyroid cancer prognostic model (TBSMRPM). The innovation of this study lies in: analyzing the role of SMRGs for the first time from a cross-disease perspective, overcoming the limitations of single-cancer research; and systematically associating the model with immune infiltration, tumor mutational burden (TMB), and drug sensitivity to provide multidimensional references for clinical treatment. Ultimately, it aims to provide new insights for the prognostic evaluation and treatment strategy optimization of thyroid cancer, while laying a foundation for elucidating the role of SMRGs in cross-cancer regulation.

2 Materials and methods

2.1 Acquisition and processing of information

We obtained transcriptomic data for thyroid and breast cancers from the TCGA database and also downloaded clinical and survival information for thyroid cancer. The THCA data lacking DFS survival information were excluded, and the final THCA cohort used for the study included 498 tumor samples and 59 normal samples. The BRCA cohort included 1,118 tumor samples and 113 normal samples. The baseline characteristics table of the patient cohort is shown in Supplementary Table 1. We also logged on to MSigDB (Molecular Signatures Database) to obtain 401 sex hormone metabolism-related genes (20). This research was authorized by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University [approval number: 2020-KY-0075-002]. The overall flowchart of this study is shown in Figure 1.

Figure 1
Flowchart depicting a research process. It starts with the TCGA-BRCA, TCGA-THCA cohorts, and sex hormone metabolism-related genes. These lead to identifying sex hormone-related genes linked to thyroid or breast cancer. The process involves screening prognosis-related genes, clinical subtyping using consensus clustering, and screening SHMRGs from subtypes. Subsequent steps include the TBSMRPM identification, nomogram analysis, immune infiltration analysis, GSEA, TMB analysis, and drug sensitivity prediction. It concludes with verifying key gene expression using qPCR.

Figure 1. Sex hormone-related gene analysis workflow.

2.2 Consistent clustering based on gene features

We firstly extracted the different genes of thyroid cancer and breast cancer from the acquired transcriptome data, and then intersected with the sex hormone metabolism-related genes to extract the sex hormone metabolism-related genes associated with thyroid cancer or breast cancer. We further obtained thyroid cancer prognosis-related genes via univariate Cox analysis. By utilizing the “ConsensusClusterPlus” R package, we accomplished consensus clustering on the basis of genetic characteristics for the THCA patients (patients were classified as two subclasses, C1 and C2) (21). The specific steps are as follows: First, z-score standardization was performed on the gene expression profiles to eliminate differences in dimensions. The key parameters were set as follows: the number of repeated sampling (reps) was 500, with 80% of samples randomly selected in each sampling, while all features (pFeature=1) were retained for clustering analysis to ensure the robustness of results against data perturbations. The “km” (k-means clustering) algorithm was used for clustering, with Euclidean distance as the distance measurement method, and the range of potential clustering numbers evaluated was k=2 to k=10. The optimal number of clusters was determined through cumulative distribution function (CDF) curve analysis. And we plotted KM curves to determine the DFS difference between the clusters. After that, we also assessed the differences in immune infiltration between clusters and evaluated the stromal and immune scores utilizing the ESTIMATE algorithm.

2.3 Identification of differential genes and establishment of TBSMRPM

We identified differential genes in different clusters using the “Deseq2” package (22). A total of 498 thyroid cancer (THCA) patients were randomly divided into a training set and a test set at a ratio of 5:5, with 249 patients in each of the training set and the test set. We used the training set to screen the optimal differential genes by univariate Cox, LASSO cox and stepwise cox regression analyses, respectively, in which the “glmnet” R package was called (23). The median risk scores established by multivariate Cox regression was utilized for the differentiating of the risk groups. We also plotted risk factor linkage charts and KM curves to show the differences in population proportions, survival (DFS) status, and gene expression across risk groups. The sensitivity and specificity of the above differential genomes are reflected by receiver operating characteristic curves (ROCs).

2.4 Exploration of model-clinical correlations

To explore the correlation between TBSMRPM and clinical features, we plotted box line plots of the relationship between risk scores and age, gender, stage, extra-glandular invasion, T stage, N stage, M stage, and tumor burden. Stratified analysis further demonstrated the variability of the risk groups within each clinical subgroup.

2.5 Modeling and evaluation of the clinical factor-related nomogram

Our Cox regression analysis included TBSMRPM and clinical factors. The univariate Cox regression analysis was utilized for the screening of model components (24). Furthermore, we plotted ROC curves, calibration curves, and decision curves to fully assess the predictive efficacy, calibration, and clinical utility of the model, respectively.

2.6 Immunocorrelation analysis and gene set enrichment analysis

We assessed the correlation of immune cells with risk scores by means of seven algorithms (2530). The immune checkpoints, immune pathways and immune cells were evaluated between the two risk groups (31). The ESTIMATE algorithm was used for the calculating of immune and stromal scores (32). In addition, changes in pathway activity between different risk groups were analyzed by GSEA (33).

2.7 Tumor mutational burden

TMB is the counting of somatic mutation sites in the tumor genome, usually described as mutations per million bases (mut/Mb). This indicator may provide some indication of a tumor’s ability to generate neoantigens and predict the efficacy of tumor immunotherapy (34). After integrating the data, we analyzed the somatic mutations associated with risk groupings using the maftools R package (35).

2.8 Drug sensitivity analysis

We screened potential therapeutic agents by correlating drug sensitivity with risk subgroups and invoked the “oncoPredict” R package to complete the calculating of drug half-maximal inhibitory concentration (IC50) (36).

2.9 Procedures of quantitative real-time PCR

Our research purchased the human papillary thyroid cancer cell lines (IHH4, KTC-1, TPC-1) and normal thyroid cell line (Nthy ori-3-1) from the Shanghai Cell Biochemical Institute (Shanghai, China). These cell lines were cultured in RPMI-1640 medium (Gemini) containing 10% fetal bovine serum (Gemini). Then the medium was placed in an incubator at 5% CO2 and 37°C for further incubation. The extraction of total RNA was performed utilizing TRIzol reagent (Invitrogen). After that, the Prime Script RT reagent kit with gDNA Eraser (Takara) was utilized for the reverse transcription process. We referred to the instructions to complete the PCR reaction by using the SYBR Green Detection kit (Takara) and utilized the 2-ΔΔCt method to complete the calculation of relative gene expression levels. Detailed primer sequences are shown in Supplementary Table 4.

2.10 Statistical analysis

All data analyses were finished by R software (version 4.3.1) and GraphPad Prism (version 8.0.2). Differences between the groups were compared using the Wilcoxon test. Kaplan-Meier curves and log-rank tests were utilized for assessment of the DFS discrepancy between the groups. P<0.05 was considered statistically significant.

3 Results

3.1 Construction of the consistent clustering

By screening tumor and normal tissue differential genes from the TCGA-THCA cohort and TCGA-BRCA cohort, and then taking the intersection with sex hormone metabolism-related genes, we obtained a total of 159 sex hormone-related genes associated with thyroid or breast cancer (Supplementary Table 2, Supplementary Figure 1C). We then extracted 26 differentially expressed genes associated with DFS by univariate Cox regression. According to the gene characteristics, the consistent clustering results showed that the curves showed a clear flattening trend and the best clustering effect when k = 2 (Figures 2A, B,Supplementary Figures 1D, E). KM curves of the clusters also demonstrated marked differences (Figure 2C, Supplementary Table 3). Enrichment scores for immune cells were assessed in the two clusters using ssGSEA analysis (Figure 2D), which showed that the cells with significant immune infiltration in the two clusters included T helper cells, aDCs, and so on. And the result demonstrated between-cluster variability in infiltration of all immune cells (P<0.05). The immune pathway result illustrated that the differential genes of the two clusters were clearly correlated with pathways of MHC class-I, HLA, and parainflammation (Figure 2E). The enrichment of the 15 immunization pathways also differed significantly across clusters (P<0.05). C1 corresponded to higher stromal and immune scores based on ESTIMATE results (Figure 2F). In addition, the immune checkpoints that were significantly expressed in C2 over C1 included TGFB1, KDR, and CD27, and those that were significantly expressed in C1 over C2 included CD274, VTCN1, and LGALS9 (Figure 2G).

Figure 2
Graphical summary of clustering analysis with seven panels labeled A to G. Panel A shows a consensus cumulative distribution function plot with colored lines for clusters one to ten. Panel B is a consensus matrix heatmap for k equals two. Panel C displays a Kaplan-Meier survival plot comparing clusters C1 and C2, with statistical significance. Panels D and E are box plots showing immune infiltration across various immune cells and pathways for clusters C1 and C2. Panel F includes a box plot for immune scores, and panel G shows expression levels of various genes, all compared between clusters C1 and C2.

Figure 2. Construction of the consensus clustering. (A) The result of consensus cumulative distribution function. (B) Determination of the ideal number of clusters. (C) K-M curves with between-cluster differences. (D) The ssGSEA scores of 16 kinds of immune cells. (E) Boxplot presenting significant differences in immune pathway enrichment between the two clusters. (F) Correlations between the tumor microenvironment score and the two clusters. (G) The presence of differentially expressed immune checkpoints. ns, no significance, *P<0.05; **P<0.01, ***P<0.001.

3.2 Identification and evaluation of the TBSMRPM

After consistent clustering, we performed differential gene analysis for both clusters and screened a total of 838 differential genes (Supplementary Table 4). Univariate cox regression analysis, LASSO cox regression analysis, and stepwise cox regression analysis was used to screen for the prognostically strong related genes and to fit the optimal model. Ultimately, we screened out a total of 10 differential genes correlated with DFS to establish the TBSMRPM (Figure 3C). The formula is as follows: Risk score = (-5.548expression of SLITRK2) + (-0.450* expression of ARHGEF37) + (-0.296* expression of F3) + (0.240* expression of MMP9) + (0.313* expression of C2CD4A) + (0.449* expression of CERS1) + (0.537* expression of PLP1) + (0.635* expression of RNF223) + (0.641* expression of SLC5A1) + (1.143* expression of HORMAD2). A total of 498 THCA patients were randomly divided into a training set and a test set at a ratio of 5:5, with 249 patients in each of the training set and the test set. The expression distributions, survival status and risk scores of the 10 genes are displayed in Supplementary Figure 2A. The KM curves of the DFS status of the two groups in the training, test, and total sets showed significant differences (Supplementary Figure 2B, Supplementary Table 3), and the ROC curves demonstrated the good predictive efficacy of the model (Supplementary Figure 2C) (33), the area under the ROC curves (AUCs) of the training-focused model to predict 1-, 3-, and 5-year DFS were 0.914, 0.859, and 0.772, respectively, and in the validation-focused model, they were 0.750, 0.795, and 0.770.

Figure 3
Panel A shows a graph of coefficients versus Log Lambda, illustrating how coefficients vary across lambda values. Panel B displays a plot of partial likelihood deviance against Log Lambda, with red dots indicating data points. Panel C is a table listing variables with hazard ratios, confidence intervals, P-values, and coefficients, alongside a forest plot for visual comparison.

Figure 3. Establishment of sex hormone metabolism-related genes prognostic model. (A, B) Parameterization for LASSO analysis. (C) Forest plots for stepwise cox analysis.

3.3 Correlation of TBSMRPM with clinical features

To determine the value of the model in depth, our research evaluated the correlation of TBSMRPM and important clinical features and found that the risk score was significantly correlated with stage, extra-glandular invasion, T stage, N stage, and tumor burden (Figure 4). Then, our study performed stratified analyses of clinical factors, age<60 years, age≥60 years, female, male, absence of extra-glandular invasion, presence of extra-glandular invasion, clinical stage I/II, clinical stage III/IV, T1/T2, T3/T4, N0, N1, M0, M1, absence of tumor burden, and presence of tumor burden, and other subgroup analyses showed significant variability in the DFS status-KM curves among different risk groups (Figure 5, Supplementary Table 3). These results suggest that TBSMRPM has excellent anticipation ability.

Figure 4
Box plots show risk scores across different clinical features: A) Age below and above 60, B) Female vs. male, C) Stage I+II vs. III+IV, D) Tumor stages T1+T2 vs. T3+T4, E) Lymph node involvement N0 vs. N1, F) Metastasis M0 vs. M1, G) No extrathyroidal extension vs. extension, H) Tumor-free vs. with tumor. Significant differences are marked with p-values.

Figure 4. Risk scores of the patients are classified by common clinical features.

Figure 5
Survival probability graphs show various medical analysis results. Each graph represents subgroups categorized by age, gender, cancer stage, tumor extension, and node involvement. The x-axis shows time in years, and the y-axis shows survival probability. High and low-risk groups are displayed with distinct lines. Significant differences are indicated by log-rank p-values.

Figure 5. Subgroup KM analysis on the basis of various clinical features.

3.4 Validation and application of the prognostic model

3.4.1 Integration of clinical factors and evaluation of nomogram performance

Our analysis results indicated that stage, T stage, N stage, M stage, and risk score were prognosis-related factors (Figure 6A). The risk score was established as an independent predictor of prognosis for DFS by multivariate Cox regression analysis (Figure 6B). Subsequently, we incorporated stage, T stage, N stage, M stage, and risk score into the model to predict DFS status of thyroid cancer patients (Figure 6C). ROC results indicated that the AUCs of the nomogram for predicting DFS status at 1, 3, and 5 years were 0.807 (Figure 6D), 0.776 (Figure 6E), and 0.752 (Figure 6F), respectively. Both decision curves (Figures 6G, I) and calibration curves (Figures 6J–L) indicated the ideality of the predictive modeling.

Figure 6
Composite image showing statistical analyses and model evaluations across different panels. Panel A presents a forest plot with hazard ratios for various variables. Panel B shows a similar plot focusing on different variables. Panel C illustrates a nomogram for predicted probabilities. Panels D, E, and F feature ROC curves, indicating model performance. Panels G, H, and I display decision curve analyses depicting net benefits across risk thresholds. Panels J, K, and L compare predicted probabilities of survival with actual survival proportions, using calibration plots with confidence intervals.

Figure 6. A nomogram model with clinical factors. (A, B) Univariate and multivariate Cox analysis results. (C) The nomogram for predicting DFS of thyroid cancer in total set. ROCs (D-F), decision curves (G-I) and calibration curves (J-L) for the nomogram and other features used to predict DFS.

The results of the seven immune algorithms indicated that the risk scores were negatively correlated with cells such as endothelial cells, CD8 T cells, and hemapoietic stem cells (coefficient < -0.3), and positively correlated with NKT cells, sebocytes, and macrophages M0 relationship (coefficient > 0.2) (Supplementary Figure 3). The ESTIMATE result indicated that the high-risk group corresponded to higher ESTIMATE scores (Figure 7A), which suggested a higher level of immune infiltration. The immune pathway analysis showed that the high-risk group was significantly enriched in Parainflammation, MHC class-I and HLA pathways (Figure 7B). Meanwhile, 12 immune pathways such as co-inhibition, APC co-stimulation, and CCR were significantly different between the TBSMRPM groups. The CD4 T cells, CD56 dim tural killer cells, etc., also showed significant differences between different risk groups, and a total of 26 immune cell infiltrations were significantly different between the two risk groups (Figure 7C). Differential expression analysis of immune checkpoints in different risk groups showed that the immune checkpoints more significantly expressed in the low-risk group included CD160, KDR, CD27, and so on. The more significantly expressed immune checkpoints in the high-risk group included LAG3, TGFBR1, LGALS9 and so on (Figure 7D).

Figure 7
Box plot panels labeled A, B, C, and D compare immune infiltration and gene expression levels. Red boxes indicate high risk; blue boxes indicate low risk. Significant differences are marked by asterisks and “ns” for non-significant. The data categories are stromal score, immune score, and various immune cell and gene expression metrics.

Figure 7. Immunocorrelation results of the two risk groups. (A) Correlations between the tumor microenvironment score and different risk groups. (B) Boxplot for differential analysis of immune pathway enrichment. (C) Boxplot presenting the significant differences in 28 types of immune cell enrichment between different risk groups. (D) Boxplot for differential analysis of immune checkpoint expression. ns, no significance, *P<0.05; **P<0.01, ***P<0.001.

3.4.2 Correlation with tumor immune microenvironment and functional pathway enrichment

To reveal gene expression patterns associated with specific biological processes, pathways or functions, our study performed GSEA analysis and it showed that apical junction, estrogen response late, antigen processing and presentation signaling pathways and other pathways had important roles in the high-risk group (Figures 8A–H), which are important in the development and progression of malignant tumors. While signaling pathways such as fatty acid metabolism, glycerolipid metabolism, phosphatidylinositol signaling system, and steroid hormone biosynthesis were significant in the low-risk group (Figures 8I–L).

Figure 8
Twelve line graphs display running enrichment scores for various biological pathways. Each panel (A-L) represents different pathways such as HALLMARK and KEGG categories. The x-axes show the rank in ordered datasets, and the y-axes show enrichment scores. Colored bars and lines depict different genetic expressions and pathways, highlighting patterns and peak scores across the datasets.

Figure 8. An in-depth exploration of the role of molecules. Significantly enriched signaling pathways in high-risk group (A-G) and low-risk group (I-L) in accordance to GSEA.

3.4.3 Association with tumor mutational burden and prognostic implications

The waterfall plot indicated that the rank order of frequency of somatic mutations in the low-risk group was BRAF (42%), NRAS (13%), TG (5%), HRAS (4%), and other genes were mutated in less than 4% of the genes (Figure 9A). The rank order of frequency of somatic mutations occurring in the high-risk group was BRAF (77%), TTN (5%), NRAS (4%), and other genes had a mutation frequency of less than 4% (Figure 9B). The somatic mutations in both groups were predominantly missense mutations (Supplementary Figures 4A, B). The box line plot visualized the significant variability of somatic mutations between the risk groups, and the high-risk group had higher TMB values (Figure 9C). Distinguishing the high and low TMB groups by the median somatic mutation scores of all patients showed a significant difference in KM analysis of DFS status between the two groups (Figure 9D, Supplementary Table 3). The KM curve analysis of the combined TMB group and risk group still showed a significant difference between the four groups (Figure 9E, Supplementary Table 3). In addition, Supplementary Figure 4C demonstrates that risk score and TMB were significantly correlated (R=0.25, P<0.0001). Our findings revealed that poor prognosis in the high-risk group may be associated with higher TMB.

Figure 9
Panel A and B: Oncoprints showing gene mutations across two sample groups, highlighting the percentage of samples altered and specific gene mutations. Panel C: Box plot comparing tumor mutation burden between high and low groups, showing significant difference (p=5.4e-06). Panel D: Survival plot stratified by high and low tumor mutation burden, with a log-rank p-value of 0.0067. Panel E: Kaplan-Meier survival curves for combinations of tumor mutation burden and risk levels, showing significant survival differences (p<0.0001).

Figure 9. Differential analysis of TMB. The mutations landscape of low-risk group (A) and high-risk group (B). (C) Relevance of TMB to risk groups. (D) Differences in DFS of TMB groups. (E) Comparison of KM curves for the combined TMB and risk groups.

3.4.4 Prediction of drug sensitivity based on risk stratification

The drug sensitivity analysis revealed that Axitinib (Figure 10A), Cyclophosphamide (Figure 10B), Lapatinib (Figure 10D), fulvestrant (Figure 10G), nilotinib (Figure 10H) were more sensitive in the low-risk group, and Erlotinib (Figure 10C), ERK-6640 (Figure 10E), AZD7762 (Figure 10F), and Sapitinib (Figure 10I) were more sensitive in the high-risk group.

Figure 10
Violin plots labeled A to I compare estimated IC50 values for various drugs between high and low groups. Each plot shows red (high) and blue (low) distributions with associated p-values, indicating statistical significance in the differences between the groups for each drug.

Figure 10. Drug sensitivity in different risk groups.

3.5 PCR for determining gene expression

We examined the expression of 10 genes in IHH4, KTC-1, TPC-1 human papillary thyroid carcinoma cell lines and human thyroid normal cell line Nthy ori-3-1 by PCR analysis. The results showed that C2CD4A, CERS1, MMP9, SLC5A1, and HORMAD2 were highly expressed in IHH4, KTC-1, and TPC-1 cell lines, and SLITRK2, ARHGEF37, PLP1, RNF223, and F3 were lowly expressed in IHH4, KTC-1, and TPC-1 cell lines (Figure 11). The associations between the 10 genes and sex hormone metabolism are detailed in Supplementary Table 8.

Figure 11
Bar graphs illustrate the relative mRNA expression levels of various genes (CERS1, C2CD4A, HORMAD2, SLC5A1, MMP9, PLP1, RNF223, SLITRK2, ARHGEF37, F3) across different cell types (NTHY, IHH4, KTC-1, TPC-1). Each graph shows expression levels with statistical significance indicated by asterisks and symbols above the bars.

Figure 11. Expression of 10 genes in IHH4, KTC-1, TPC-1 human papillary thyroid cancer cell lines and human thyroid normal cell line Nthy ori-3-1. * indicates P<0.05, ** indicates P < 0.01, *** indicates P < 0.001; # indicates P<0.05, ## indicates P<0.01, ### indicates P<0.001, & indicates P<0.05, && indicates P<0.01, &&& indicates P<0.001.

4 Discussion

Numerous studies have shown that there is a potential, bidirectional pathogenic relationship between thyroid cancer and breast cancer, and that the development of one malignant tumor may increase the risk of developing the other (3739). Of interest, although both malignancies are highly prevalent in the female population, the simultaneous occurrence of both diseases is more common in men (40). This suggests that sex hormones may serve as a vital role in the tumorigenesis and progression of thyroid cancer. However, studies on the prognostic relationship between sex hormones and thyroid cancer are still very limited (2, 41). In our research, we constructed a prognostic model for the first time on the basis of sex hormone metabolism-related genes in thyroid cancer and breast cancer, which was able to better predict the prognostic (DFS) profile of thyroid cancer patients.

We first screened for genes correlated with sex hormone metabolism. Notably, considering the potential bidirectional pathogenic relationship between thyroid and breast cancers, we could not ignore the possible indirect effects of genes with altered expression only in a single malignancy on the other tumor. Thus, we included genes related to sex hormone metabolism associated with thyroid or breast cancer together in our study. This idea is different from the study of Duan et al. (42). Subsequently, we classified thyroid cancer patients into two subtypes using consensus clustering. Since the core basis of consensus clustering lies in the similarity of expression patterns of “sex hormone metabolism-related genes”, the algorithm automatically clusters samples with more similar expression patterns into the same group. Therefore, the division of C1 and C2 inherently implies “distinguishable differences in the overall expression profiles of sex hormone-related genes”—this serves as the premise for the stable output of 2 subtypes by clustering. In addition, we utilized existing analyses (pathway enrichment, immune microenvironment, and prognosis association) to cross-validate the expression differences of sex hormone-related genes at the functional level. Furthermore, the significant differences in prognosis and immune infiltration between different subtypes can provide references for clinical staging. To further screen the more desirable prognostic genes and fit an excellent prognostic model, we performed univariate Cox, LASSO Cox and stepwise cox regression analyses after differential gene analysis of subtypes. The ROC results of both training and validation sets suggested a better predictive performance of the model. And the significant survival difference results in the clinical correlation study indicated the rationality of the TBSMRPM risk grouping based on TBSMRPM. We further established a nomogram model in accordance with clinical factors for the prognostic prediction.

In addition, understanding the tumor immune microenvironment contributes to the diagnosis and treatment of thyroid cancer (43). Various cells in the tumor microenvironment such as: endothelial cells, fibroblasts, and infiltrating immune cells serve as vital roles in the processes of tumorigenesis, invasion, and metastasis. In our results, activated dendritic cells, activated CD8+ T cells, Th1 cells, and macrophages were significantly infiltrated in the high-risk group, which may be linked to the adverse prognosis. Relevant research has shown that dendritic cells may result in the tumor progression by affecting the tumor angiogenesis and causing immune dysregulation (44), whereas significant infiltration of activated CD8+ T cells, Th1 cells suggests a higher risk of recurrence (4547), and significant infiltration of macrophages in thyroid cancer is strongly linked to tumor invasion and short survival (48). The immune pathway analysis illustrated that the risk grouping was strongly associated with the MHC class I pathway. Angell et al. demonstrated that the absence of MHC class I expression was a vital mechanism of immune escape in thyroid malignancies, which highly overlapped with our findings (49). The immune checkpoint results demonstrated the strong correlation of risk groupings with some immune checkpoints. The study by Park et al. also confirmed that these differential genes associated with immune escape signaling were overexpressed in patients with advanced thyroid cancer, which could help to predict the recurrence of advanced thyroid cancer (47). In conclusion, our immune correlation analysis demonstrated the variability of immune function between different subtypes and risk groups in the consistent clustering. The different risk subgroups should be treated with different therapeutic strategies, and the rational formulation of immunotherapy can help to improve the efficacy of the treatment and prevent the development of drug resistance (50).

In the subsequent GSEA analysis, we found that signaling pathways such as fatty acid metabolism, glycerolipid metabolism, phosphatidylinositol signaling system, and steroid hormone biosynthesis were significantly enriched in the low-risk group. It is well known that sex hormones are important members of the steroid family, and their synthesis and metabolism are inextricably linked to these signaling pathways. Sun et al. found that increased fatty acid metabolism was linked to the decreased infiltration of T- and B-cells in male breast cancer and other cancer types, and that targeting fatty acid metabolism pathways may alleviate the immunosuppressive microenvironment in a variety of cancers (51). Lu et al. demonstrated the importance of lipid metabolism in thyroid cancer (52). Liu et al. experimentally verified that pyruvate carboxylase is closely associated with tumor aggressiveness in thyroid cancer by stimulating fatty acid synthesis (53). The phosphatidylinositol (PI) signaling system is a key cellular signaling pathway concerning the cellular processes such as cell growth, differentiation, survival, and intracellular communication, and dysregulation of the related signaling pathway is associated with thyroid cancer and breast cancer (54, 55). Furthermore, a small number of studies examining the association of other sex hormones with thyroid cancer, in addition to estrogen-related studies, have also yielded positive findings (56, 57).

After that, we performed somatic mutation and drug sensitivity analysis. And we found the high-risk group had a significantly higher frequency of BRAF gene mutations. The high-risk group had a higher TMB value and corresponded to the poor prognosis. Studies have confirmed the importance of the detection of BRAF gene mutations in the diagnosis, treatment and determination of prognosis of thyroid cancer (58). This is consistent with the results of our study. There are fewer studies related to other genes with higher mutation frequencies such as NRAS in our findings (59), which found that NRAS mutations were linked to the high risk of distant metastasis of thyroid cancer (60, 61). These findings suggest, to some extent, the potential role of this gene in thyroid cancer. Furthermore, we gained more insight into the correlation of drug sensitivity between different risk groups.

Ultimately, we screened 10 ideal prognosis-related genes. Among them, MMP9 has been shown to promote thyroid cancer incidence and progression in a large number of literatures (62, 63), which is consistent with our findings. HORMAD2 has been shown to inhibit the incidence and progression of thyroid cancer and its hypermethylation promotes the progression of thyroid cancer in a study by Lin et al. (64). Han et al. mentioned that the mRNA level of C2CD4A was dysregulated in the lung cancer (65). Rong et al. confirmed that C2CD4A was highly expressed in colorectal cancer tissues and contributed to tumor growth by inhibiting the P53 signaling pathway (66). Chen et al. verified that CERS1 was highly expressed in colorectal cancer by PCR (67). Xu et al. found that CERS1 was significantly downregulated in non-small-cell lung cancer cell lines and brain metastatic tissue and its upregulation corresponded to better prognosis (68). Lei et al. proposed that RNF223 could promote pancreatic cancer growth and migration, and identified potential protein targets and metabolism-related pathways (69). Zhang et al. demonstrated that ARHGEF37 promotes hepatocellular carcinoma outgrowth and metastasis through activation of Cdc42 (70). SLC5A1 has been confirmed to be highly expressed in many malignant tumors by several studies, including our research, and a study on pancreatic cancer demonstrated that SLC5A1 regulates cancer cell growth through AMPK/mTOR signaling (71). However, studies on the association of these genes with thyroid cancer have rarely been reported and need to be further explored.

There are some limitations to this study. First, due to limited sample acquisition, we performed PCR validation only in papillary thyroid carcinoma. Second, the model established in this study lacks validation in the external cohorts. Third, the actual role of the screened genes in thyroid cancer lacks sufficient experimental validation. In conclusion, by screening differential genes and establishing a prognostic model, we mined potential novel biomarkers related to thyroid and breast cancers, which are highly suggestive for prognostic prediction of thyroid cancer. However, the actual clinical value of the model needs to be further validated in large cohort studies, and the actual roles and specific mechanisms of these potential biomarkers in malignant tumors need to be further determined by a large number of in vivo and ex vivo experiments.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

LQ: Writing – original draft, Writing – review & editing. HL: Writing – review & editing, Writing – original draft. KY: Writing – original draft, Writing – review & editing. RM: Writing – review & editing, Writing – original draft. YZ: Writing – original draft, Writing – review & editing. YG: Writing – original draft, Writing – review & editing. DY: Writing – review & editing, Writing – original draft.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Acknowledgments

We are very grateful to the staff in TCGA Program for their kind work in data collection and delivery. This study was funded by the General Project of Natural Science Foundation of Henan Province (No.222300420568), the Key Medical Science and Technology Project of Henan Province (No. SBGJ202101014), the Major Scientific Research Projects of Traditional Chinese Medicine in Henan Province (No.20-21ZYZD14) and the Cultivation of Young and Middle-aged Health Science and Technology Innovation Leading Talents in Henan Province (No. YXKC2020015). Thanks to the above funds for supporting this research.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1641195/full#supplementary-material

References

1. Nielsen SM, White MG, Hong S, Aschebrook-Kilfoy B, Kaplan EL, Angelos P, et al. The breast-thyroid cancer link: A systematic review and meta-analysis. Cancer Epidemiol Biomarkers Prev. (2016) 25:231–8. doi: 10.1158/1055-9965.EPI-15-0833

PubMed Abstract | Crossref Full Text | Google Scholar

2. Halada S, Casado-Medrano V, Baran JA, Lee J, Chinmay P, Bauer AJ, et al. Hormonal crosstalk between thyroid and breast cancer. Endocrinology. (2022) 163. doi: 10.1210/endocr/bqac075.Cited in: Pubmed

PubMed Abstract | Crossref Full Text | Google Scholar

3. Bolf EL, Sprague BL, and Carr FE. A Linkage Between Thyroid and Breast Cancer: A Common Etiology? Cancer Epidemiol Biomarkers Prev. (2019) 28:643–9. doi: 10.1158/1055-9965.EPI-18-0877

PubMed Abstract | Crossref Full Text | Google Scholar

4. Bakos B, Kiss A, Árvai K, Szili B, Deák-Kocsis B, Tobiás B, et al. Co-occurrence of thyroid and breast cancer is associated with an increased oncogenic SNP burden. BMC Cancer. (2021) 21:706. doi: 10.1186/s12885-021-08377-4

PubMed Abstract | Crossref Full Text | Google Scholar

5. Avgerinos KI, Spyrou N, Mantzoros CS, and Dalamaga M. Obesity and cancer risk: Emerging biological mechanisms and perspectives. Metabolism. (2019) 92:121–35. doi: 10.1016/j.metabol.2018.11.001

PubMed Abstract | Crossref Full Text | Google Scholar

6. Yin D-T, He H, Yu K, Xie J, Lei M, Ma R, et al. The association between thyroid cancer and insulin resistance, metabolic syndrome and its components: A systematic review and meta-analysis. Int J Surg. (2018) 57:66–75. doi: 10.1016/j.ijsu.2018.07.013

PubMed Abstract | Crossref Full Text | Google Scholar

7. Hamadeh IS, Patel JN, Rusin S, and Tan AR. Personalizing aromatase inhibitor therapy in patients with breast cancer. Cancer Treat Rev. (2018) 70:47–55. doi: 10.1016/j.ctrv.2018.07.014

PubMed Abstract | Crossref Full Text | Google Scholar

8. Nounu A, Kar SP, Relton CL, and Richmond RC. Sex steroid hormones and risk of breast cancer: a two-sample Mendelian randomization study. Breast Cancer Res. (2022) 24:66. doi: 10.1186/s13058-022-01553-9

PubMed Abstract | Crossref Full Text | Google Scholar

9. Coradini D and Oriana S. Impact of sex hormones dysregulation and adiposity on the outcome of postmenopausal breast cancer patients. Clin Obes. (2021) 11:e12423. doi: 10.1111/cob.12423

PubMed Abstract | Crossref Full Text | Google Scholar

10. Song ZM, Wang YD, Chai F, Zhang J, Lv S, Wang JX, et al. Estrogen enhances the proliferation, migration, and invasion of papillary thyroid carcinoma via the ERα/KRT19 signaling axis. J Endocrinol Invest. (2025) 48:653–70. doi: 10.1007/s40618-024-02473-5

PubMed Abstract | Crossref Full Text | Google Scholar

11. Li M, Chai H-F, Peng F, Meng Y-T, Zhang L-Z, Zhang L, et al. Estrogen receptor β upregulated by lncRNA-H19 to promote cancer stem-like properties in papillary thyroid carcinoma. Cell Death Dis. (2018) 9:1120. doi: 10.1038/s41419-018-1077-9

PubMed Abstract | Crossref Full Text | Google Scholar

12. Dai Y-J, Qiu Y-B, Jiang R, Xu M, Liao L-Y, Chen GG, et al. Concomitant high expression of ERα36, GRP78 and GRP94 is associated with aggressive papillary thyroid cancer behavior. Cell Oncol (Dordr). (2018) 41:269–82. doi: 10.1007/s13402-017-0368-y

PubMed Abstract | Crossref Full Text | Google Scholar

13. Ahn HY, Song R-Y, Ahn HS, and Kim HS. Expression of estrogen and progesterone receptors in papillary thyroid carcinoma in korea. Cancer Res Treat. (2021) 53:1204–12. doi: 10.4143/crt.2020.1201

PubMed Abstract | Crossref Full Text | Google Scholar

14. Costa P, Catarino AL, Silva F, Sobrinho LG, and Bugalho MJ. Expression of prolactin receptor and prolactin in normal and Malignant thyroid: a tissue microarray study. Endocr Pathol. (2006) 17:377–86. doi: 10.1007/s12022-006-0009-x

PubMed Abstract | Crossref Full Text | Google Scholar

15. Banu KS, Govindarajulu P, and Aruldhas MM. Testosterone and estradiol have specific differential modulatory effect on the proliferation of human thyroid papillary and follicular carcinoma cell lines independent of TSH action. Endocr Pathol. (2001) 12:315–27. doi: 10.1385/EP:12:3:315

PubMed Abstract | Crossref Full Text | Google Scholar

16. Gupta A, Carnazza M, Jones M, Darzynkiewicz Z, Halicka D, O’Connell T, et al. Androgen receptor activation induces senescence in thyroid cancer cells. Cancers (Basel). (2023) 15. doi: 10.3390/cancers15082198

PubMed Abstract | Crossref Full Text | Google Scholar

17. Fu J, He M, Wu Q, Zhang X, Qi X, Shen K, et al. The clinical and genetic features in patients coexisting primary breast and thyroid cancers. Front Endocrinol (Lausanne). (2023) 14:1136120. doi: 10.3389/fendo.2023.1136120

PubMed Abstract | Crossref Full Text | Google Scholar

18. Yi JW, Kim S-J, Kim JK, Seong CY, Yu HW, Chai YJ, et al. Upregulation of the ESR1 gene and ESR ratio (ESR1/ESR2) is associated with a worse prognosis in papillary thyroid carcinoma: the impact of the estrogen receptor α/β Expression on clinical outcomes in papillary thyroid carcinoma patients. Ann Surg Oncol. (2017) 24:3754–62. doi: 10.1245/s10434-017-5780-z

PubMed Abstract | Crossref Full Text | Google Scholar

19. Zhang L, Zhou M, Gao X, Xie Y, Xiao J, Liu T, et al. Estrogen-related genes for thyroid cancer prognosis, immune infiltration, staging, and drug sensitivity. BMC Cancer. (2023) 23:1048. doi: 10.1186/s12885-023-11556-0

PubMed Abstract | Crossref Full Text | Google Scholar

20. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, and Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | Crossref Full Text | Google Scholar

21. Wilkerson MD and Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | Crossref Full Text | Google Scholar

22. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

23. Friedman J, Hastie T, and Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software. (2010) 33. doi: 10.18637/jss.v033.i01

PubMed Abstract | Crossref Full Text | Google Scholar

24. Iasonos A, Schrag D, Raj GV, and Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol. (2008) 26:1364–70. doi: 10.1200/JCO.2007.12.9791

PubMed Abstract | Crossref Full Text | Google Scholar

25. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | Crossref Full Text | Google Scholar

26. Chen B, Khodadoust MS, Liu CL, Newman AM, and Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12

PubMed Abstract | Crossref Full Text | Google Scholar

27. Plattner C, Finotello F, and Rieder D. Deconvoluting tumor-infiltrating immune cells from RNA-seq data using quanTIseq. Methods Enzymol. (2020) 636:261–85. doi: 10.1016/bs.mie.2019.05.056

PubMed Abstract | Crossref Full Text | Google Scholar

28. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. (2016) 17:218. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | Crossref Full Text | Google Scholar

29. Hong K, Cen K, Chen Q, Dai Y, Mai Y, and Guo Y. Identification and validation of a novel senescence-related biomarker for thyroid cancer to predict the prognosis and immunotherapy. Front Immunol. (2023) 14:1128390. doi: 10.3389/fimmu.2023.1128390

PubMed Abstract | Crossref Full Text | Google Scholar

30. Aran D, Hu Z, and Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | Crossref Full Text | Google Scholar

31. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | Crossref Full Text | Google Scholar

32. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | Crossref Full Text | Google Scholar

33. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. (2005) 102:15545–50.

PubMed Abstract | Google Scholar

34. Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. (2017) 9:34. doi: 10.1186/s13073-017-0424-2

PubMed Abstract | Crossref Full Text | Google Scholar

35. Mayakonda A, Lin D-C, Assenov Y, Plass C, and Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | Crossref Full Text | Google Scholar

36. Maeser D, Gruener RF, and Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. (2021) 22. doi: 10.1093/bib/bbab260

PubMed Abstract | Crossref Full Text | Google Scholar

37. Van Fossen VL, Wilhelm SM, Eaton JL, and McHenry CR. Association of thyroid, breast and renal cell cancer: a population-based study of the prevalence of second Malignancies. Ann Surg Oncol. (2013) 20:1341–7. doi: 10.1245/s10434-012-2718-3

PubMed Abstract | Crossref Full Text | Google Scholar

38. Li S, Yang J, Shen Y, Zhao X, Zhang L, Wang B, et al. Clinicopathological features, survival and risk in breast cancer survivors with thyroid cancer: an analysis of the SEER database. BMC Public Health. (2019) 19:1592. doi: 10.1186/s12889-019-7947-y

PubMed Abstract | Crossref Full Text | Google Scholar

39. Joseph KR, Edirimanne S, and Eslick GD. The association between breast cancer and thyroid cancer: a meta-analysis. Breast Cancer Res Treat. (2015) 152:173–81. doi: 10.1007/s10549-015-3456-6

PubMed Abstract | Crossref Full Text | Google Scholar

40. Zafon C, Obiols G, and Mesa J. Second primary cancer in patients with papillary thyroid carcinoma. Anticancer Res. (2013) 33:337–40.

Google Scholar

41. Tan H, Wang S, Huang F, and Tong Z. Association between breast cancer and thyroid cancer risk: a two-sample Mendelian randomization study. Front Endocrinol (Lausanne). (2023) 14:1138149. doi: 10.3389/fendo.2023.1138149

PubMed Abstract | Crossref Full Text | Google Scholar

42. Duan J, Liu C, Yi J, and Wang Y. Shared sex hormone metabolism-related gene prognostic index between breast and endometrial cancers. Front Endocrinol (Lausanne). (2023) 14:1126862. doi: 10.3389/fendo.2023.1126862

PubMed Abstract | Crossref Full Text | Google Scholar

43. Menicali E, Guzzetti M, Morelli S, Moretti S, and Puxeddu E. Immune landscape of thyroid cancers: new insights. Front Endocrinol (Lausanne). (2020) 11:637826. doi: 10.3389/fendo.2020.637826

PubMed Abstract | Crossref Full Text | Google Scholar

44. Qiu H, Hu X, He C, Yu B, Li Y, and Li J. Identification and validation of an individualized prognostic signature of bladder cancer based on seven immune related genes. Front Genet. (2020) 11:12. doi: 10.3389/fgene.2020.00012

PubMed Abstract | Crossref Full Text | Google Scholar

45. Cunha LL, Marcello MA, Nonogaki S, Morari EC, Soares FA, Vassallo J, et al. CD8+ tumour-infiltrating lymphocytes and COX2 expression may predict relapse in differentiated thyroid cancer. Clin Endocrinol (Oxf). (2015) 83:246–53. doi: 10.1111/cen.12586

PubMed Abstract | Crossref Full Text | Google Scholar

46. Chen Z, Guo M-L, Li Y-Y, Yan K, Li L, Shen F, et al. Immune profiling identifies CD8+ T-cell subset signatures as prognostic markers for recurrence in papillary thyroid cancer. Front Immunol. (2022) 13:894919. doi: 10.3389/fimmu.2022.894919

PubMed Abstract | Crossref Full Text | Google Scholar

47. Park S-J, Kang YE, Kim J-H, Park J-L, Kim S-K, Baek S-W, et al. Transcriptomic analysis of papillary thyroid cancer: A focus on immune-subtyping, oncogenic fusion, and recurrence. Clin Exp Otorhinolaryngol. (2022) 15:183–93. doi: 10.21053/ceo.2021.02215

PubMed Abstract | Crossref Full Text | Google Scholar

48. Yin H, Tang Y, Guo Y, and Wen S. Immune microenvironment of thyroid cancer. J Cancer. (2020) 11:4884–96. doi: 10.7150/jca.44506

PubMed Abstract | Crossref Full Text | Google Scholar

49. Angell TE, Lechner MG, Jang JK, LoPresti JS, and Epstein AL. MHC class I loss is a frequent mechanism of immune escape in papillary thyroid cancer that is reversed by interferon and selumetinib treatment in vitro. Clin Cancer Res. (2014) 20:6034–44. doi: 10.1158/1078-0432.CCR-14-0879

PubMed Abstract | Crossref Full Text | Google Scholar

50. Mould RC, van Vloten JP, AuYeung AWK, Karimi K, and Bridle BW. Immune responses in the thyroid cancer microenvironment: making immunotherapy a possible mission. Endocr Relat Cancer. (2017) 24:T311–29. doi: 10.1530/ERC-17-0316

PubMed Abstract | Crossref Full Text | Google Scholar

51. Sun H, Zhang L, Wang Z, Gu D, Zhu M, Cai Y, et al. Single-cell transcriptome analysis indicates fatty acid metabolism-mediated metastasis and immunosuppression in male breast cancer. Nat Commun. (2023) 14:5590. doi: 10.1038/s41467-023-41318-2

PubMed Abstract | Crossref Full Text | Google Scholar

52. Lu J, Zhang Y, Sun M, Ding C, Zhang L, Kong Y, et al. Multi-omics analysis of fatty acid metabolism in thyroid carcinoma. Front Oncol. (2021) 11:737127. doi: 10.3389/fonc.2021.737127

PubMed Abstract | Crossref Full Text | Google Scholar

53. Liu C, Zhou X, Pan Y, Liu Y, and Zhang Y. Pyruvate carboxylase promotes thyroid cancer aggressiveness through fatty acid synthesis. BMC Cancer. (2021) 21:722. doi: 10.1186/s12885-021-08499-9

PubMed Abstract | Crossref Full Text | Google Scholar

54. Santarpia L, El-Naggar AK, Cote GJ, Myers JN, and Sherman SI. Phosphatidylinositol 3-kinase/akt and ras/raf-mitogen-activated protein kinase pathway mutations in anaplastic thyroid cancer. J Clin Endocrinol Metab. (2008) 93:278–84. doi: 10.1210/jc.2007-1076

PubMed Abstract | Crossref Full Text | Google Scholar

55. Amirani E, Hallajzadeh J, Asemi Z, Mansournia MA, and Yousefi B. Effects of chitosan and oligochitosans on the phosphatidylinositol 3-kinase-AKT pathway in cancer therapy. Int J Biol Macromol. (2020) 164:456–67. doi: 10.1016/j.ijbiomac.2020.07.137

PubMed Abstract | Crossref Full Text | Google Scholar

56. Zhang LJ, Xiong Y, Nilubol N, He M, Bommareddi S, Zhu X, et al. Testosterone regulates thyroid cancer progression by modifying tumor suppressor genes and tumor immunity. Carcinogenesis. (2015) 36:420–8. doi: 10.1093/carcin/bgv001

PubMed Abstract | Crossref Full Text | Google Scholar

57. O’Connell TJ, Dadafarin S, Jones M, Rodríguez T, Gupta A, Shin E, et al. Androgen activity is associated with PD-L1 downregulation in thyroid cancer. Front Cell Dev Biol. (2021) 9:663130. doi: 10.3389/fcell.2021.663130

PubMed Abstract | Crossref Full Text | Google Scholar

58. Caronia LM, Phay JE, and Shah MH. Role of BRAF in thyroid oncogenesis. Clin Cancer Res. (2011) 17:7511–7. doi: 10.1158/1078-0432.CCR-11-1155

PubMed Abstract | Crossref Full Text | Google Scholar

59. Yin D-T, Yu K, Lu R-Q, Li X, Xu J, Lei M, et al. Clinicopathological significance of TERT promoter mutation in papillary thyroid carcinomas: a systematic review and meta-analysis. Clin Endocrinol (Oxf). (2016) 85:299–305. doi: 10.1111/cen.13017

PubMed Abstract | Crossref Full Text | Google Scholar

60. Melo M, Gaspar da Rocha A, Batista R, Vinagre J, Martins MJ, Costa G, et al. BRAF, and NRAS in primary thyroid cancer and metastatic disease. J Clin Endocrinol Metab. (2017) 102:1898–907. doi: 10.1210/jc.2016-2785

PubMed Abstract | Crossref Full Text | Google Scholar

61. Jang EK, Song DE, Sim SY, Kwon H, Choi YM, Jeon MJ, et al. NRAS codon 61 mutation is associated with distant metastasis in patients with follicular thyroid carcinoma. Thyroid. (2014) 24:1275–81. doi: 10.1089/thy.2014.0053

PubMed Abstract | Crossref Full Text | Google Scholar

62. Guan H, Guo Y, Liu L, Ye R, Liang W, Li H, et al. INAVA promotes aggressiveness of papillary thyroid cancer by upregulating MMP9 expression. Cell Biosci. (2018) 8:26. doi: 10.1186/s13578-018-0224-4

PubMed Abstract | Crossref Full Text | Google Scholar

63. Ouyang X, Feng L, Yao L, Xiao Y, Hu X, Zhang G, et al. Testicular orphan receptor 4 (TR4) promotes papillary thyroid cancer invasion via activating circ-FNLA/miR-149-5p/MMP9 signaling. Mol Ther Nucleic Acids. (2021) 24:755–67. doi: 10.1016/j.omtn.2021.03.021

PubMed Abstract | Crossref Full Text | Google Scholar

64. Lin Q, Hou S, Guan F, and Lin C. HORMAD2 methylation-mediated epigenetic regulation of gene expression in thyroid cancer. J Cell Mol Med. (2018) 22:4640–52. doi: 10.1111/jcmm.13680

PubMed Abstract | Crossref Full Text | Google Scholar

65. Han F, Qian L, Zhang Y, Liu P, Li R, and Chen M. C2CD4A/B variants in the predisposition of lung cancer in the Chinese Han population. Funct Integr Genomics. (2022) 22:331–40. doi: 10.1007/s10142-022-00827-x

PubMed Abstract | Crossref Full Text | Google Scholar

66. Rong Z, Luo Z, Fu Z, Zhang P, Li T, Zhang J, et al. The novel circSLC6A6/miR-1265/C2CD4A axis promotes colorectal cancer growth by suppressing p53 signaling pathway. J Exp Clin Cancer Res. (2021) 40:324. doi: 10.1186/s13046-021-02126-y

PubMed Abstract | Crossref Full Text | Google Scholar

67. Chen L, Chen H, Li Y, Li L, Qiu Y, and Ren J. Endocannabinoid and ceramide levels are altered in patients with colorectal cancer. Oncol Rep. (2015) 34:447–54. doi: 10.3892/or.2015.3973

PubMed Abstract | Crossref Full Text | Google Scholar

68. Xu Y, Pan J, Lin Y, Wu Y, Chen Y, and Li H. Ceramide synthase 1 inhibits brain metastasis of non-small cell lung cancer by interacting with USP14 and downregulating the PI3K/AKT/mTOR signaling pathway. Cancers (Basel). (2023) 15. doi: 10.3390/cancers15071994

PubMed Abstract | Crossref Full Text | Google Scholar

69. Feng L, Wang J, Zhang J, Diao J, He L, Fu C, et al. Comprehensive analysis of E3 ubiquitin ligases reveals ring finger protein 223 as a novel oncogene activated by KLF4 in pancreatic cancer. Front Cell Dev Biol. (2021) 9:738709. doi: 10.3389/fcell.2021.738709

PubMed Abstract | Crossref Full Text | Google Scholar

70. Zhang X, Ren L, Wu J, Feng R, Chen Y, Li R, et al. ARHGEF37 overexpression promotes extravasation and metastasis of hepatocellular carcinoma via directly activating Cdc42. J Exp Clin Cancer Res. (2022) 41:230. doi: 10.1186/s13046-022-02441-y

PubMed Abstract | Crossref Full Text | Google Scholar

71. Gao H-F, Chen L-Y, Cheng C-S, Chen H, Meng Z-Q, and Chen Z. SLC5A1 promotes growth and proliferation of pancreatic carcinoma via glucose-dependent AMPK/mTOR signaling. Cancer Manag Res. (2019) 11:3171–85. doi: 10.2147/CMAR.S195424

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: thyroid cancer, breast cancer, sex hormone metabolism-related gene, prognostic model, immune infiltrate

Citation: Qiao L, Li H, Yin K, Ma R, Zhang Y, Guo Y and Yin D (2025) Multi-disease transcriptomic analysis of sex hormone genes reveals a novel prognostic model for thyroid cancer with breast cancer correlations. Front. Oncol. 15:1641195. doi: 10.3389/fonc.2025.1641195

Received: 10 June 2025; Accepted: 18 August 2025;
Published: 02 September 2025.

Edited by:

Mamunur Rahaman, University of New South Wales, Australia

Reviewed by:

Tomohiro Chiba, Cancer Institute Hospital of Japanese Foundation for Cancer Research, Japan
Hao Wang, Shenzhen University General Hospital, China

Copyright © 2025 Qiao, Li, Yin, Ma, Zhang, Guo and Yin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Detao Yin, ZGV0YW95aW5Aenp1LmVkdS5jbg==

These authors have contributed equally to this work and share first authorship

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.