This study integrates transcriptomic and single-cell sequencing data from BC and TC, and employs a multi-faceted, multi-method analytical approach to identify shared molecular targets for the treatment of BC and TC (Figure 1).

Transcriptomic data (RNA-seq) for BC and TC were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Specifically, GSE124646 (which includes 90 breast cancer samples and 10 normal samples) and GSE126698 (comprising 22 thyroid cancer samples and 6 normal samples) were used for the primary data analysis. For validation of hub genes, GSE109169 (containing 25 breast cancer samples and 25 normal samples) and GSE140109 (which includes 6 thyroid cancer samples and 4 normal samples) were utilized.

Single-cell transcriptomic data (scRNA-seq) for BC and TC were also sourced from the GEO database. The datasets included GSE161529, which contains 6 breast cancer samples and 13 normal samples, and GSE191288, which comprises 6 thyroid cancer samples and 1 normal sample.

Differentially expressed genes (DEGs) between cancer and normal samples were identified using the limma package in R, with the thresholds set to adjusted P < 0.05 and |log₂FC| > 0.585. Subsequently, the ClusterProfiler package was employed to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) enrichment analyses on the identified DEGs, followed by visual representation of the results. The enrichment analysis parameters were set to pvalueCutoff = 0.05 and qvalueCutoff = 0.5.

Quality control of the single-cell transcriptomic datasets was performed using the Seurat R package (v4.1.2). Low-quality cells and low-expression genes were excluded based on the following criteria: the number of detected features per cell ranged from 200 to 5,000, the number of transcripts per cell ranged from 1,000 to 20,000, and the proportion of mitochondrial gene expression per cell was less than 20%. Data normalization was carried out using the NormalizeData function, and highly variable genes were identified using the FindVariableFeatures function with nfeatures = 2000.

To correct batch effects, the Harmony R package (version 0.1.1) was applied. Subsequently, linear transformation was performed using the ScaleData function. Classification was carried out through Principal Component Analysis (PCA) and Partial Least Squares Discriminant Analysis (PLS-DA) models (Chen et al., 2021). The optimal number of principal components (PCs) was assessed using the ElbowPlot function, and 50 PCs were selected for further analysis. Cell clustering was performed using the FindNeighbors and FindClusters functions (dims = 1:30, resolution = 2).

Cell clusters were annotated based on canonical marker genes reported in the literature. DEGs among cell types were identified using the FindAllMarkers function in Seurat, with the parameters min. pct = 0.1 and logfc. threshold = 0.25, retaining only genes with p < 0.05.

Additionally, to estimate the cellular composition from bulk RNA-seq data, the CIBERSORT function in the IOBR R package was applied for deconvolution analysis, with the permutation count (perm) set to 100.

Weighted Gene Co-expression Network Analysis (WGCNA) was performed using the WGCNA R package (v1.72-1) to identify gene modules correlated with BC and TC phenotypes. Differentially expressed genes from the GSE124646 (BC) and GSE126698 (TC) datasets were used to construct gene expression matrices for network analysis. To determine the appropriate soft-thresholding power (β) required for scale-free topology, the pickSoftThreshold function was applied to both datasets, with optimal β values ranging from 10 to 12. Subsequently, an adjacency matrix was computed using the formula (aij = |Sij|β). The adjacency matrix was then transformed into a topological overlap matrix (TOM), and a dissimilarity matrix (1−TOM) was calculated. Hierarchical clustering based on this dissimilarity matrix was performed to identify distinct gene modules. Modules with strong correlations to clinical phenotypes were selected as candidate modules for downstream analyses, including biomarker discovery and functional annotation.

The “randomForest” R software package was used to use the random forest (RF) machine learning algorithm to screen central genes that are highly related to thyroid cancer and breast cancer. The classification accuracy of different numbers of RF feature genes was determined, and the feature genes with the highest classification accuracy were retained to determine the final Hub genes. In order to test the diagnostic efficacy of Hub genes, the receiver operating characteristic (ROC) curve and the corresponding area under the ROC curve (AUC) of each Hub gene were calculated based on the normalized expression level of each Hub gene.

The CIBERSORT algorithm was applied to the GSE126698 (TC) and GSE124646 (BC) datasets to estimate the relative proportions of 22 immune cell types within each sample. Differences in immune cell composition between cancer and normal tissues were assessed using the Wilcoxon rank-sum test, with P values calculated for statistical significance. To further evaluate tumor microenvironment (TME) characteristics at the sample level, the ESTIMATE R package was employed to calculate the immune infiltration score (ImmuneScore), stromal cell content (StromalScore), composite microenvironment score (ESTIMATEScore), and tumor purity (TumorPurity). Spearman’s rank correlation analysis was conducted to assess the association between hub gene expression and immune cell infiltration levels. A P value <0.05 was considered statistically significant.

Tumor-related microRNAs (miRNAs) were retrieved from the Human miRNA Disease Database (HMDD) (http://www.cuilab.cn/hmdd). Hub gene-associated mRNA–miRNA interaction pairs were obtained from the miRWalk database (http://mirwalk.umm.uni-heidelberg.de/), and intersected with the tumor-related miRNAs to identify relevant mRNA–miRNA regulatory relationships. Only pairs with a target score >80 were retained for further analysis. To identify potential long non-coding RNAs (lncRNAs) interacting with the tumor-associated miRNAs, predictions were performed using the ENCORI database, and corresponding lncRNA–miRNA interaction pairs were collected.

For transcriptional regulatory analysis, the RcisTarget R package was used to predict transcription factors (TFs) targeting the hub genes (Aibar and González-Blas, 2017). Motif enrichment analysis was conducted to identify significant TF-binding motifs, which were then used to construct the TF regulatory network.

Fresh tumor tissues and corresponding adjacent normal tissues were collected from six patients diagnosed with both TC and BC, ensuring the integrity and freshness of all specimens. Each sample was rinsed with physiological saline and homogenized into a tissue suspension. Total RNA was extracted using TRIzol reagent according to the manufacturer’s protocol, followed by DNase treatment to eliminate genomic DNA contamination. Reverse transcription was performed using a commercial reverse transcription kit to synthesize cDNA. For qRT-PCR, 2 μL of cDNA was placed into an EP tube and amplified using a SYBR Green pre-mix. GAPDH was used as the internal reference gene. The threshold cycle (Ct) values were determined, and relative mRNA expression levels were calculated using the 2^−ΔΔCT method.

Tissues from six patients with concurrent TC and BC, including tumours and adjacent normal tissues, were fixed with 4% paraformaldehyde, then embedded in paraffin and sectioned. The sections were stained according to the manufacturer’s instructions to observe the expression levels of relevant proteins in the sections.

All calculations and statistical analyses in this study were performed using R (https://www.r-project.org/, version 4.1.2). To determine the statistical significance of differences between two groups of normally distributed data, an independent Student’s t-test was used, while the Mann-Whitney U test (i.e., Wilcoxon rank-sum test) was employed to assess differences between non-normally distributed variables. All p-values were calculated based on two samples, and p-values less than 0.05 were considered statistically significant. Additionally, Spearman correlation analysis was used in this study to obtain the correlation coefficients between variables. P-values were calculated on both sides of the equation, and values less than 0.05 were considered statistically significant.

Immunohistochemistry and RT-qPCR data were analysed using SPSS 26.00 statistical software. All data are presented as mean ± standard deviation. A non-parametric t-test was used to compare the two groups of data. A p-value <0.05 was considered statistically significant.

Differential expression analysis of the BC dataset identified a total of 1,242 dysregulated genes, including 626 upregulated and 616 downregulated genes (Figure 2A). In the TC dataset, 101 genes were upregulated and 342 were downregulated (Figure 2B). Cross-comparison of the two datasets revealed 76 overlapping genes (Co-DEGs) shared between BC and TC (Figures 2C,D). Among them, 50 genes were consistently downregulated, 22 consistently upregulated, and 4 displayed divergent expression trends between the two cancers. GO and KEGG pathway enrichment analyses indicated that the Co-DEGs were predominantly involved in ECM-receptor interaction, ABC transporter regulation, and cortisol synthesis and secretion (Figures 2E–H).

To explore shared cellular components between BC and TC, we integrated scRNA-seq data and deconvolution analysis. scRNA-seq datasets GSE161529 (BC) and GSE191288 (TC) were analyzed separately. Following quality control (Figure 3A) and batch effect correction (Figure 3C), 75,069 cells from the BC dataset were classified into 11 distinct cell types (Figure 3D), and deconvolution analysis confirmed significant variation in these cell types across samples (Figure 3E). In the TC dataset, 32,735 cells were classified into 9 cell types (Figures 3B,F,G), and similarly, transcriptome deconvolution revealed significant heterogeneity across TC samples (Figure 3H). Notably, both BC and TC patients exhibited increased monocyte infiltration compared to healthy controls (Figure 3I). This observation led to a focus on monocyte infiltration in subsequent analyses.

To identify gene expression modules shared by BC and TC, WGCNA was performed. For the BC dataset, a soft-threshold power of 10 was selected to ensure scale-free topology (Figures 4A,B), resulting in the identification of 14 co-expression modules, each represented by a unique color (Figure 4C). Module–trait relationships were assessed to identify modules significantly associated with BC phenotypes (Figure 4D). Similarly, WGCNA was applied to the TC dataset using a soft-threshold power of 12 (Figures 4E,F), which yielded 8 distinct modules (Figure 4G). The correlation between each module and TC phenotypes was also analyzed (Figure 4H). To assess cross-cancer module conservation, we evaluated the overlap between BC-specific and TC-specific modules. The analysis revealed substantial overlap, indicating shared transcriptional networks between BC and TC (Figure 4I).

Moreover, consistent with the single-cell analysis findings, WGCNA revealed a significant positive correlation between monocyte infiltration and both the turquoise module in BC and the blue module in TC. Therefore, these two modules were designated as key consensus gene modules shared between breast and thyroid cancer.

We next evaluated whether the gene signatures from key consensus modules could distinguish patients with BC and TC from healthy individuals. Analyses were performed using Co-DEGs from the blue module of the BC consensus network and the blue module of the TC consensus network. Heatmap visualization revealed two major findings: (1) Co-DEGs were consistently upregulated in both BC and TC patients compared to controls. (2) Hierarchical clustering of key consensus module genes effectively separated cancer patients from healthy controls (Figures 5A,E). Principal PCA based on DEGs from the consensus gene modules of BC and TC demonstrated clear separation between the disease and control groups (Figures 5B,F), with distinct clustering patterns observed for patients with the two cancer types. Additionally, further validation using PLS-DA also revealed different classification patterns between patients with the two cancer types (Figures 5C,G).

To explore potential functional relevance, KEGG pathway enrichment analysis was performed on DEGs from the consensus module of BC and TC. Notably, the glycosaminoglycan biosynthesis pathway was significantly enriched in both cancer types (Figures 5D,H), suggesting a potentially important role for this pathway in the pathogenesis of both BC and TC.

To identify hub genes within the key consensus modules, DEGs from the BC turquoise module (n = 16) and TC blue module (n = 13) were analyzed using a random forest classifier. The optimal mtry parameter was selected based on the lowest classification error (Figures 6A,E). For BC, the error rate stabilized when the number of decision trees reached 1,500 (Figure 6B), while for TC, stabilization occurred at 500 trees (Figure 6F). These values were used for subsequent analyses, with all other parameters set to default. Based on the Gini coefficient, the top 10 most informative genes for each cancer type were identified (Figures 6C,D,G,H). Intersection of the top-ranked genes from BC (n = 15) and TC (n = 12) yielded seven shared hub genes: MKI67, TAP1, UBE2C, CENPF, PILRA, SPP1, and TMEM51. Receiver operating characteristic (ROC) curve analysis showed that each hub gene exhibited moderate discriminatory ability between cancer patients and healthy controls (Figures 6I,J).

We further assessed immune infiltration patterns using single-cell transcriptomic data and tumor microenvironment (TME) scoring at the sample level. CIBERSORT analysis revealed significant differences in CD4 memory resting T cells and regulatory T cells (Tregs) between BC and TC samples (Figures 7A,D). ESTIMATE analysis also indicated significant variation in ImmuneScore between BC and TC (Figures 7B,E). These findings suggest that T cell populations may play a central role in shaping the immune landscape in both cancers.

Spearman correlation analysis between the seven hub genes and immune cell proportions showed that, in both BC and TC: Hub genes were negatively correlated with CD4 memory resting T cells and monocytes. Hub genes were positively correlated with Tregs and M1 macrophages (Figures 7C,F). These results suggest that hub genes may influence BC and TC progression through modulation of immune cell infiltration.

To investigate potential regulatory mechanisms underlying the expression of the identified hub genes, we analyzed both TF and miRNA regulatory networks. Tumor-associated miRNAs were obtained from the Human miRNA Disease Database (HMDD). mRNA–miRNA interaction pairs involving hub genes were extracted from the miRWalk database and intersected with 657 tumor-related miRNAs, resulting in 12 validated mRNA–miRNA interaction pairs (Figure 8A). Using the ENCORI database, we then predicted lncRNAs interacting with these tumor-associated miRNAs, and corresponding interaction pairs were identified (Figure 8B). In parallel, motif enrichment analysis was performed to identify key TFs potentially regulating the hub genes. The results suggested that hub gene expression may be regulated by transcription factors such as NFYB, NR4A2, and BPTF (Figure 8C).

To validate hub gene expression, we performed RT-qPCR on tumor and adjacent normal tissue samples. In TC tissues, the expression levels of SPP1, TAP1, PILRA, TMEM51, UBE2C, and MKI67 were significantly upregulated compared to adjacent normal tissues, whereas CENPF did not show a statistically significant difference (Figures 9A–G). In BC tissues, TAP1, PILRA, UBE2C, TMEM51, and MKI67 were significantly upregulated, while the expression levels of SPP1 and CENPF remained unchanged between cancer and adjacent tissues (Figures 9H–N). Subsequently, we selected PILRA, MKI67, and UBE2C, which were significantly upregulated in both cancer types, for further validation using immunohistochemistry. The IHC results confirmed that the protein levels of these three genes were significantly elevated in both BC and TC tissues compared to adjacent normal tissues (Figures 8O–R).