The Comprehensive Analysis of Weighted Gene Co-Expression Network Analysis and Machine Learning Revealed Diagnostic Biomarkers for Breast Implant Illness Complicated with Breast Cancer

Introduction

Breast cancer (BC) is the most frequently diagnosed type of cancer, with approximately 300,590 new cases anticipated.¹ In numerous instances, mastectomy serves as the primary treatment for malignant tumours and is thus regarded as a therapeutic option. The removal of a breast can significantly affect these women’s overall quality of life, resulting in feelings of psychological inadequacy and reduced social integration.^2,3 Breast reconstruction is a crucial component of holistic treatment for breast cancer, as it improves the aesthetic appearance and restores the natural shape of the breast.⁴ Prosthetic implants, also known as implant-based breast reconstruction (IBR), are the most commonly used method for breast reconstruction.⁵ A longitudinal trend analysis of the NIS database covering 1998–2008 revealed that the annual rate of implant reconstruction increased by 11%,⁶ indicating that medical staff should choose personalized approaches for each individual, considering issues related to implants and their psychological impact.

Breast implant illness (BII) refers to a range of systemic symptoms that patients attribute to theirthe silicone breast implants. These symptoms typically include fatigue, fever, muscle pain, and joint pain.^7,8 Many females have likewise been afflicted with autoimmune inflammatory conditions, such as arthritic inflammation, Sjögren’s condition, systemic sclerosis, and syndromes resembling lupus, and numerous patients indicate that their symptoms improve after their implants are removed. The evidence associating breast implants with BII has been limited to potential links, without a proven mechanism or established cause. This has led to several hypotheses concerning the aetiology of breast implant illness, including those related to autoimmune responses, biofilm formation, and psychosomatic influences.^9,10 There is much misinformation online about breast implant illness and suggested treatments. This misinformation should be addressed to enhance our understanding of this condition.¹¹ The exploration of relevant diagnostic markers of BII has become an urgent clinical problem.

Over the past decade, genomic expression profiling has proven to be valuable in delivering crucial information for predicting survival in patients with cancer.¹² Furthermore, genetic biomarkers are regarded as significant tools for forecasting prognostic outcomes and treatment responses in cancer patients.

To our knowledge, there is a lack of research aimed at identifying immune-related biomarkers for BII diagnosis. The relationship between BC and BII remains uncertain. Therefore, we initially utilized an array of bioinformatics techniques to explore the link between BC and BII through high-throughput sequencing and to identify the potential cellular and molecular pathways. In our research, we pinpointed possible crosstalk genes (CGs) that may link BC and BII and evaluated the interactions between these CGs and infiltrating immune cells by employing sophisticated statistical methods. This approach was taken to enhance our comprehension of the pathophysiological mechanisms that could connect BC and BII. Furthermore, the possible utility of CGs in disease diagnosis was assessed and validated across multiple patient groups. By uncovering these common gene signatures and molecular pathways, our objective is to enhance our comprehension of the development of BC-BII and to identify prospective therapeutic targets to improve the prognosis of individuals afflicted by this ailment.

Materials and Methods

Data Retrieval

RNA sequencing information pertaining to individuals with BC was sourced from the Cancer Genome Atlas (TCGA) repository. The unprocessed datasets, specifically GSE178425, which included gene expression information for both patients with BII and healthy subjects, were accessed from the Gene Expression Omnibus (GEO) repository. To assess the diagnostic efficacy, we also used the GSE31448 dataset as a validation dataset. Supplementary Figure 1 presents a flowchart detailing the research methodology.

Identification of Differentially Expressed Genes

The “limma” R package was employed to detect genes that exhibited significant differences in expression (DEGs) within the TCGA-BC dataset. RNA-seq data processing and normalization for GSE178425 were performed using the DESeq2 package. The criteria for identifying DEGs were established at a log fold change (log FC) of at least ±1.2 and a p value of less than 0.05, and the results were visualized in volcano plots. Additionally, a comprehensive analysis of DEGs from both the TCGA-BC cohort and the GSE178425 dataset was performed, with overlaps highlighted using Venn diagrams.

Weighted Gene Co-Expression Network Analysis and Module Gene Selection

The Weighted Gene Co-Expression Network Analysis (WGCNA) package in R was used to construct the coexpression network for the GSE178425 dataset. Initially, the median absolute deviation (MAD) for each gene was computed, and the genes with the lowest 50% MAD values were excluded. The DEG expression matrix was subsequently filtered using the goodSamplesGenes function to remove any substandard genes or samples, resulting in the creation of a scale-free coexpression network. The appropriate soft threshold for forming a scale-free network was identified using the “pickSoftThreshold” function. A hierarchical clustering dendrogram was then generated to categorize similar genes into distinct modules, ensuring that each module contained no fewer than 150 genes. Finally, Pearson correlation analysis was performed to evaluate the associations between each module. The modules with the strongest positive correlations with BII were chosen for further analysis.

Enrichment Analysis

Next, we pinpointed the intersection of DEGs that are common to both BC and the principal module genes associated with BII, as well as the intersection of DEGs within BII and the key module genes of BII. We performed GO and KEGG analyses on these shared gene sets. The “clusterProfiler” R package was used to conduct enrichment analysis for both the GO and KEGG pathways of the overlapping genes, setting the threshold at P < 0.05. We generated histograms, bubble plots, and circle plots of the results using the “enrichplot” and “ggplot2” R packages.

Protein–Protein Interaction Network Construction

To investigate the interplay among proteins, pathways, and coexpression, we employed the STRING database (https://cn.string-db.org/) to construct a protein‒protein interaction (PPI) network for the DEGs associated with BC and the pivotal module genes of BII, setting the minimum interaction score threshold at 0.400. We identified hub genes employing CytoHubba, an add-on for the Cytoscape application. CytoHubba is a plugin that evaluates nodes based on their network properties, allowing for the detection of pivotal nodes within biological networks.¹³ Furthermore, MCODE, a plugin for Cytoscape, was applied to filter the key modules of the core genes from the PPI network.

Machine Learning

Three machine learning algorithms were utilized to refine the selection of candidate genes for BC diagnosis: the least absolute shrinkage and selection operator (LASSO) regression method,¹⁴ support vector machine (SVM),¹⁵ and random forest (RF).¹⁶ We utilized Lasso Cox proportional hazards regression to pinpoint changes in the regression coefficients of pivotal genes. The most effective parameter, λ, was determined through 10-fold cross-validation with the R package glmnet, and genes were chosen based on lambda.min. We employed R package visualizations to depict the shrinkage of coefficients in LASSO Cox regression. SVM is a supervised machine learning method applied to regression or classification tasks, necessitating a labelled training dataset. SVM-RFE, a machine learning technique, trains a subset of features from various groups to reduce the feature set and identify the most influential features. Finally, the RF algorithm employs the R package “Random Forest” to evaluate the importance of genes. The intersection of the outcomes from the three methods can be regarded as potential hub genes for diagnosis.

Development of a Nomograph and Assessment of the Receiver Operating Characteristic Curve

To assess the significance of the candidate genes in diagnosing BC in the context of BII, we constructed a nomograph utilizing the “rms” R package, in which “Points” denotes the score assigned to the candidate genes. The receiver operating characteristic (ROC) curve was plotted to evaluate the diagnostic impact of the candidate genes, and a nomograph was constructed in the context of BC diagnosis. The area under the curve (AUC) and 95% confidence interval (CI) were computed to measure its effectiveness. The ROC analysis yielded the AUC and 95% confidence interval (CI) as measures of the effectiveness of the nomograph, where an AUC value exceeding 0.7 was considered to indicate robust diagnostic capability.

Immune Infiltration Analysis

To enhance the understanding of immune cell characteristics in the tissues of the BC versus normal group and the BII versus normal group, we analysed the variations in immune cell subsets among the samples. The infiltration of 28 immune cell types within tumour samples was assessed through single-sample gene set enrichment analysis (ssGSEA). The average values between the BC and normal groups, as well as the BII and normal groups, were compared using the t_test function in the R package rstatix.

Single-Cell RNA Sequencing (scRNA-Seq) Analysis

Six single-cell RNA sequencing samples of breast cancer (GSM4909289, GSM4909296, and GSM4909281) and normal breast tissue (GSM4909253, GSM4909254, and GSM4909257) were obtained from GSE161529. The Seurat package was then employed to enhance the scRNA-seq data by selecting higher-quality cells. Normalization was performed using the “normalizedata” function in the Seurat R package. Principal component analysis (PCA) was executed with the “RunPCA” function in Seurat, emphasizing the top 2000 genes to minimize data dimensionality. The cells were classified and characterized using the “FindNeighbors” and “FindClusters” functions with a resolution of 1. The resulting clusters were displayed in two dimensions using the “RunTSNE” and “RunUMAP” functions.

Cell Culture

The human breast cancer cell lines MCF-7, T47D, SUM-159PT, and UACC-812, along with the normal human epithelial breast cell line MCF-10A, were obtained from Procell Life Science & Technology Co., Ltd. (Wuhan, China). MCF-7, T47D, and MCF-10A cells were maintained in DMEM (Gibco) supplemented with 10% FBS (Gibco) and 10 U/mL penicillin‒streptomycin (Gibco). Moreover, SUM-159PT and UACC-812 cells were cultured in RPMI 1640 medium (Sigma‒Aldrich, St. Louis, MO, USA). All the cells were incubated at 37 °C in an atmosphere of 5% CO2.

Quantitative RT‒PCR Analysis

Total RNA was harvested from both breast cancer cell lines and a normal human epithelial breast cell line using TRIzol reagent (Invitrogen). The cDNA was subsequently converted into cDNA through reverse transcription using a Transcriptor First Strand cDNA Synthesis Kit (Roche, Penzberg, Germany). FastStart Universal SYBR Green Master Mix (Roche) was used to amplify each sample within a 20 µL reaction mixture. The fold changes were determined using the 2-DDCt method. The expression levels were quantified and normalized to those of the internal control GAPDH. The primer sequences are listed in Table 1.

Table 1 The Sequences of Primers

Immunohistochemistry

We obtained formalin-fixed, paraffin-embedded tissue sections from 25 patients with breast cancer who were treated surgically at Harbin Medical University Cancer Hospital between March 2013 and June 2022. The patients had undergone mastectomy, followed by either immediate or delayed implant-based breast reconstruction. Individuals were identified with BII based on a clinical assessment procedure that comprised an in-depth medical history consultation and an examination of a thorough symptom checklist as outlined by Imran¹⁷ and a physical examination following the protocols established in previous BII studies.^18–20 The paraffin-embedded sections fixed with formalin were subjected to xylene and ethanol to eliminate the paraffin, followed by a rinse with distilled water. The sections were pretreated with EDTA Target Retrieval Solution at 120 °C and pH 8.0 for 3 minutes in a pressure cooker. To block endogenous peroxidase activity, a 3% H2O2 solution in PBS was applied for 10 minutes. Nonspecific binding was reduced by incubating the sections with goat serum for 1 hour. The sections were subsequently treated with primary antibodies overnight at 4 °C and then incubated with secondary antibodies for 30 minutes at 37 °C. The primary antibody utilized was KRT14 (1:200), which was obtained from Proteintech (Wuhan, CN), and the secondary antibody used was goat anti-rabbit IgG. Colour development was performed using diaminobenzidine (DAB) staining. Two pathologists evaluated all the specimens independently and without bias, assessing the percentage of positively stained membranes. Written informed consent was obtained from each patient.

Statistical Analysis

Statistical analyses were performed using R (version 3.6.1). To examine multiple data groups, we used analysis of variance (ANOVA) followed by either Student’s t test or the Wilcoxon rank sum test for pairwise comparisons. A p value of less than 0.05 was considered statistically significant.

Results

Detection of Genes with Altered Expression

In total, 10,643 DEGs were detected in the TCGA-BRCA dataset, meeting the criteria of a p value < 0.05 and a |log2FC| > 1.2, including 5249 genes whose expression levels were decreased and 5394 genes whose expression levels were increased in breast cancer tissue compared with normal breast tissue. The volcano plot and heatmap depicted in Figure 1A and B, respectively, display the differential expression patterns of these DEGs. Similarly, for BII, a complete set of 2466 DEGs was identified (consisting of 772 upregulated and 1694 downregulated genes), applying identical criteria. The patterns of differential expression for these DEGs in the context of BII are displayed in Figure 2A and B.

Figure 1 Differentially expressed genes (DEGs) between the BC and normal groups. (A) In this figure, red and green colours indicate DEGs with notably higher and lower expression levels in the BC groups, respectively. (B) Heatmap presenting the significantly expressed genes in both the BC and control groups. Abbreviation: BC, breast cancer.

Figure 2 Identification of DEGs in BII using Limma and WGCNA module genes. (A) Volcano plot displaying DEGs, with red and green triangles highlighting significant genes. (B) Heatmap displaying the upregulated and downregulated DEGs identified from the BII dataset. (C and D) The soft threshold of b = 9 was selected based on scale Independence and average connectivity. (E) Clustering dendrogram illustrating the separation between BII and non-BII samples. (F) Gene coexpression modules are represented by different colours under the gene tree. (G) Heatmap illustrating eigengene adjacency. (H) Heatmap depicting the correlation between module genes and BII indicating that the magenta module is most strongly associated with BII. For each pair, the top left triangle represents the correlation coefficient, whereas the bottom right triangle indicates the p value. (I) Correlation plot between module membership and the significance of genes within the magenta module. Abbreviations: BII, breast implant illness; DEGs, Differentially expressed genes; WGCNA, weighted gene co-expression network analysis.

Analysis of Weighted Gene Coexpression Networks and Identification of Key Modules

We established a scale-free coexpression network using weighted gene coexpression network analysis (WGCNA) to pinpoint the module most significantly associated with BII. We chose b = 9 (scale-free R² = 0.86) as the “soft” threshold, considering both scale independence and average connectivity (Figure 2C and D). Figure 2E displays the clustering dendrogram for BII and the control group. Using this threshold, 19 gene coexpression modules of distinct colours were detected, with a module merging threshold of 0.25 and a minimum module size of 150, as shown in Figure 2F–H. Figure 2I shows the correlation between BII and the coexpression modules, with the magenta module (consisting of 1137 genes) showing the strongest correlation with BII and identified as the principal module for subsequent analysis.

Functional Enrichment Analysis of BII

GSE178425 represents a novel dataset for BII that has not been previously investigated. To determine whether this dataset can reliably mirror the pathogenesis of BII, we conducted functional enrichment analysis using the intersecting genes from Deseq2 and the WGCNA module genes. A total of 386 common genes (CGs) were identified by intersecting the 2466 DEGs with the 1137 genes in the magenta module (Figure 3A).

Figure 3 Functional enrichment analysis of the overlapping genes associated with BII. (A) The overlap between DEGs identified by the Deseq2 and WGCNA module genes consists of 386 genes, as depicted in the Venn diagram. (B) KEGG pathway analysis of these intersecting genes reveals various significant pathways, with different colours indicating distinct pathways and their associated genes. (C–E) The GO analysis encompasses biological processes, cellular components, and molecular functions. The y-axis indicates the GO terms, whereas the x-axis shows the proportion of genes associated with each term. The circle size reflects the number of genes, and their colour denotes the significance level, as indicated by the p value. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The KEGG analysis indicated that the CGs were enriched mainly in “cell adhesion molecules”, “axon guidance”, “tight junctions”, and “Staphylococcus aureus infection” (Figure 3B). GO analysis indicated that the CGs were enriched primarily in biological process (BP) terms such as “anatomical structure morphogenesis”, “tissue development”, and “biological adhesion” (Figure 3C). In terms of cellular component (CC) ontology, the common genes were associated with “intrinsic component of membrane”, “integral component of membrane”, and “plasma membrane part”, as illustrated in Figure 3D. With respect to molecular function (MF), “transporter activity” was the most significant term among the overlapping genes, as shown in Figure 3E.

The enrichment analysis revealed that the CGs associated with BII were predominantly associated with cell adhesion molecules and Staphylococcus aureus infection pathways, which corroborated some prior scholarly conjectures and discourses regarding the pathogenesis of BII.

Analysis of Enrichment in Breast Cancer Linked to Breast Implant Illness and Identification of Hub Genes via Protein‒Protein Interaction Network Assembly

The overlap of DEGs associated with BC and the key module genes related to BII encompassed 541 genes, as illustrated in Figure 4A. To investigate the connection between BC-related genes and the pathogenesis of BII, enrichment analysis was conducted on these genes. The KEGG analysis revealed that the 541 genes were predominantly enriched in the “PI3K–Akt signalling pathway”, “cell cycle” and “tight junction” terms, as demonstrated in Figure 4B. The GO analysis indicated that the genes were associated with “anatomical structure morphogenesis”, “tissue development” (BP), “extracellular region”, “extracellular region part” (CC), and “cell adhesion molecule binding” (MF), as shown in Figure 4C–E.

Figure 4 Functional enrichment analysis of genes associated with both breast cancer (BC) and breast implant illness (BII). (A) Venn diagram illustrating that 541 genes were identified from the intersection of genes identified in BC using Limma and in BII using WGCNA. (B) KEGG analysis of these 541 common genes. (C–E) The GO analysis revealed the biological process, cellular component, and molecular function of these 541 genes.

To further elucidate the potential interactions among proteins encoded by shared DEGs and pinpoint central genes, we constructed a PPI network of the DEGs using STRING (Supplementary Figure 2). To investigate genes that might be significant in the co-occurrence of BC and BII, the central genes for the overlapping DEGs were pinpointed with the aid of CytoHubba. Considering the multifaceted nature of biological networks, a variety of topological analysis algorithms were utilized in tandem to detect hub genes. The metrics MCC (Figure 5A), MNC (Figure 5B), Degree (Figure 5C), and EPC (Figure 5D) were applied to forecast and assess the top 20 prominent hub genes within the PPI networks. The intersection of these 20 genes identified by the four algorithms yielded 17 potential central genes: PROM1, KRT5, KRT18, MUC1, KRT8, CDKN2A, CD24, KRT14, EPCAM, CLDN7, KRT17, ERBB2, KIT, ALB, CLDN4, KRT19, and GATA3 (Figure 5E).

Figure 5 Illustration of the construction of a PPI network and the identification of key genes. (A) Identification of the 17 candidate hub genes using four different algorithms: DEGREE (B), EPC (C), MCC (D), and MNC (E). Abbreviation: PPI, protein-protein interaction network.

Identification of Candidate Hub Genes via Machine Learning

LASSO analysis was used to select candidate genes for the risk model, resulting in 16 genes being filtered based on the optimal λ value (Figure 6A). The RF algorithm further assessed the significance of these genes, ranking a total of 17 genes. From this ranking, we identified the top five most critical genes (Figure 6B). Additionally, a machine learning method that employed SVM offered insights into the top fourteen most significant genes (Figure 6C). Using genes identified through LASSO analysis, RF, and SVM, three genes, KRT14, KIT, and ALB, were identified for the final validation (Figure 6D).

Figure 6 Identification of key hub genes using LASSO analysis, the RF algorithm, and SVM. (A) LASSO regression analysis. (B) Application of the RF algorithm. (C) Machine learning strategy employing SVM. (D) Venn diagram highlighting the central genes pinpointed by LASSO, SVM-RFE, and RF. Abbreviations: LASSO, least absolute shrinkage and selection operator; SVM, support vector machine; RF, random forest.

Assessment of Diagnostic Value

A prognostic model was constructed utilizing the three identified hub genes (Figure 7A), and ROC curve was plotted to assess the diagnostic sensitivity and specificity of each gene along with the nomogram. We computed the AUC and 95% confidence interval for each item. The findings were as follows: KRT14 (AUC 0.846, 95% CI 0.798–0.894), KIT (AUC 0.903, 95% CI 0.869–0.937), ALB (AUC 0.806, 95% CI 0.755–0.856), and nomogram (AUC 0.992, 95% CI 0.981–1.000) (Figure 7B–E). All the identified genes demonstrated significant diagnostic potential for BC in the context of BII, and the constructed nomogram had the highest diagnostic value. We validated the model with the GSE31448 dataset (AUC 0.898, 95% CI 0.855–0.940), as demonstrated in Figure 7F. Both the genes and the nomogram demonstrated substantial diagnostic value for BC in conjunction with BII.

Figure 7 Creation of a nomogram and its diagnostic value assessment. (A) The nomogram designed for diagnosing BII in the context of breast cancer. (B-F) ROC curves for each candidate gene (KRT14, KIT, ALB), the nomogram, and its validation using the GSE31448 dataset. Abbreviation: ROC, receiver operating characteristic.

ScRNA-Seq Analysis of Key Genes

To investigate the differential expression of key genes across various cell types, we performed scRNA-seq analysis on BC samples and normal breast samples. After processing the data, we acquired gene expression profiles from 15,705 cells across three BC samples and 18,238 cells across three normal breast samples for further analysis. We utilized PCA to decrease the dimensionality of the data (Figure 8A and B). The single R package facilitated cell type annotation, revealing nine different cell types (Figure 8C). We determined the expression profiles of three key genes in different samples and various cell types: KRT14 and KIT were predominantly highly expressed in normal tissue samples, whereas ALB was predominantly highly expressed in tumor samples. In addition, KRT14 was expressed predominantly in epithelial cells; KIT was expressed mainly in common myeloid progenitor (CMP); and ALB was expressed primarily in induced pluripotent stem cells (iPS cells) (Figure 8D–H).

Figure 8 Overview of single-cell atlases in BC patients. (A) t-SNE clustering plot of 6 samples. (B) t-SNE clustering plot of BC and normal samples. (C) 9 cell types identified based on marker gene expression. T-SNE plot and Dot plot (H) highlighting the expression patterns of KRT14 (D), KIT (E) and ALB (F) for the 9 cell types and different samples (G). Abbreviation: BC, breast cancer.

Validation of Differentially Expressed Hub Genes in Clinical Samples and Cell Models

To further assess the credibility of the prognostic models related to the hub genes, we first measured the mRNA levels of three hub genes via quantitative RT‒PCR analysis in various BC cell lines alongside the normal breast epithelial cell line MCF-10A. KRT14 was predominantly highly expressed in MCF-10A cell lines, but lowly expressed in other breast cancer cell lines (Figure 9A). While KIT and ALB were predominantly highly expressed in MCF-10A and T47D, but lowly expressed in SUM-159PT, MCF-7 and UACC-812 (Figure 9B and C).

Figure 9 Evaluation of the diagnostic significance of key genes in breast cancer cell lines and paraffin-embedded tissues. Comparison of the transcription levels of KRT14 (A), KIT (B), and ALB (C) mRNAs in breast cancer cell lines to normal human epithelial breast cell line. (D) Immunohistochemical (IHC) staining with Diaminobenzidine (DAB) demonstrates strong positivity for the target antigen, indicating high KRT14 protein expression levels in breast cancer. (E) IHC staining with DAB demonstrates weak positivity for the target antigen, indicating low KRT14 protein expression levels in normal breast tissue. (F) IHC staining of KRT14 expression in breast cancer and normal breast tissue. (G) IHC staining with DAB demonstrates strong positivity for the target antigen, indicating high KRT14 protein expression levels in normal breast tissue with BII. (H) IHC staining with DAB demonstrates weak positivity for the target antigen, indicating low KRT14 protein expression levels in normal breast tissue with Non-BII. (I)IHC staining of KRT14 expression in normal breast tissue stratified by BII and Non-BII. * p<0.05; *** p<0.001; **** p<0.0001. Abbreviation: BII, breast implant illness.

Because KRT14 had the highest weight among the three genes in the random forest analysis, we performed immunohistochemical staining and investigated KRT14 expression in breast cancer tissues and normal breast tissues from 25 patients with breast cancer who underwent prosthetic implantation at Harbin Medical University Cancer Hospital. We randomly selected 25 patients with early-stage breast cancer for immunohistochemical analysis. Nine of these patients were BII while 16 were Non-BII. With the groups for pTNM stages I, II, and III consisting of 18 (72.0%), 6 (24.0%), and 1 (4.0%) patients, respectively. Among these patients, 10 (40.0%) were HER2-positive, and 15 (60.0%) were HER2-negative; additionally, 12 (48.0%) had Ki-67 levels above 20%, whereas 13 (52.0%) had levels below 20% (Table 2). KRT14 expression was primarily localized in the cytoplasm (Figure 9D and E,9G and H). Figure 9E and F shows that KRT14 is highly expressed in paired normal breast tissues from the same patient but lowly expressed in breast cancer tissues. The results were measured by assessing the positive area in comparison with the total area. We compared breast cancer carcinoma and normal breast tissues from 25 patients and found that KRT14 was significantly expressed in normal tissues (Figure 9F). Similarly, Figure 9G and H shows that KRT14 is highly expressed in normal breast tissues of BII patients but lowly expressed in Non-BII patients. We also compared KRT14 expression in normal breast tissue from 25 patients (Figure 9I). Compared with that in patients who did not have BII, the expression of KRT14 was also significantly greater in patients with BII, suggesting that KRT14 may be a relevant diagnostic marker for BII.

Table 2 Chi-Square Test Analysis of the Connection Between BII and Clinicopathological Features in BC of HMU Cancer Hospital Cohort

Comparison of the Immune Infiltration Microenvironments of Patients with BII and BC

Given that prior research has presented findings corroborating a link between breast implants and autoimmune conditions, it is plausible that BII might also rely on biofilm infection and the microenvironment surrounding the implant. The analysis of immune cell infiltration can provide deeper insights into the role of immunity in BII. To investigate the common pathogenic mechanisms and immune microenvironments in BC and BII, we employed the ssGSEA algorithm to thoroughly evaluate the levels of immune cell infiltration under both conditions.

We employed ssGSEA to assess the variations in the prevalence of 28 distinct immune cell populations between the healthy and breast cancer cohorts, as well as between the non-BII and BII subgroups (Figure 10A and B). In contrast to normal individuals, individuals with BC presented increased immune cell infiltration, including increased levels of activated CD4 T cells, activated CD8 T cells, activated dendritic cells, CD56dim natural killer cells, immature B cells, myeloid-derived suppressor cells (MDSCs), memory B cells, regulatory T cells, and type 2 T helper cells. In contrast, healthy controls had higher levels of CD56bright natural killer cells, effector memory CD8+ T cells, mast cells, neutrophils, and eosinophils. In the comparison between the non-BII and BII groups, the levels of eosinophils and CD56dim natural killer cells were significantly different. As shown in Figure 10C and D, we performed a correlation analysis of immune cells in both patients with BC and patients with BII, revealing significant associations among most of the immune cells in these groups. We performed distinct correlation analyses to explore the relationships between our identified hub genes KRT14, KIT, and ALB and the levels of immune cells. In patients with BC, a positive correlation was found between the expression of KRT14 and KIT and the number of plasmacytoid dendritic cells and the p-value is statistically significant (r = 0.40). In contrast, ALB expression was significantly negatively correlated with activated CD4+ T cells (r = −0.19) (Figure 10E).

Figure 10 Investigation of the immune landscape in BII complicating breast cancer (BC) and its link to pivotal genes using single-sample gene set enrichment analysis (ssGSEA). (A and B) Analysis of the levels of immune cell infiltration across 28 immune cell subsets in the BC group versus the normal group and the BII group versus the non-BII group. (C and D) Correlation between distinct immune cell populations in the BC group versus the normal group and the BII group versus the non-BII group. (E and F) Relationships between the infiltration of immune cells and two central genes in patients from both the BC group and the normal group and between the BII group and the non-BII group. * p<0.05; ** p<0.01;*** p<0.001; **** p<0.0001. Abbreviations: BC, breast cancer; BII, breast implant illness.

In individuals with BII, a significant inverse relationship was detected between the expression of ALB and immature B cells, with a correlation coefficient of −0.40. Conversely, the expression of KRT14 and KIT was strongly positively correlated with type 2 T helper cells, with a correlation significant of 0.54 (Figure 10F). Overall, these findings suggest that these pivotal genes contribute to bolstering the immune response as the disease progresses in individuals with BII.

Discussion

Breast cancer is one of the most frequently diagnosed malignancies affecting females globally, representing a significant health concern as a result of its high incidence and mortality rates. Over the last ten years, bioinformatics analysis has been extensively utilized to investigate prognostic indicators in breast cancer.^21,22 Conversely, BII, recognized by both patients and physicians, refers to the diverse array of systemic manifestations that occur subsequent to the insertion of breast prosthetics in the absence of a definitive clinical diagnosis for such symptoms.²³ BII is not classified as a disease by the World Health Organization, and there are currently no diagnostic procedures available to ascertain that the systemic symptoms are a result of silicone implants.¹¹ Additionally, clinical examinations often reveal no objective findings, and there have been no prior studies that have integrated the two diseases. In this research, we investigated the shared signalling pathways and mechanisms underlying BII and BC using an integrated approach that combined bioinformatics analysis with machine learning methodologies to construct a nomogram and evaluate its effectiveness in diagnosing BII among patients with breast cancer. A pivotal discovery was the identification of three principal candidate genes—KRT14, KIT, and ALB—which led to the creation of a nomogram designed for the diagnosis of BII in individuals with breast cancer.

Consistent functional enrichment analysis of the WGCNA and DEG datasets consistently highlighted the activation of pivotal signalling pathways, notably, cell adhesion molecules, axon guidance, tight junctions, and Staphylococcus aureus infections, in the context of BII. Our findings are corroborated by the increased expression of cell adhesion molecules and the incidence of Staphylococcus aureus infection in BII. Following breast implant surgery, proteins may adsorb on the implant surface, and the adsorption of these proteins is necessary for wound healing. However, if the thickness of the protein layer increases or if cell adhesion is enhanced by specific protein conformations, this may lead to excessive envelope formation.²⁴ In another study, 20 patients who underwent breast implant removal and total capsulectomy were analysed and found to have increased mast cell counts in the resected tissue of all subjects, indicating that cell adhesion molecules could be involved in the pathogenesis of BII. Colonization of the implant surface by bacteria and biofilm formation are critical steps in breast implant-related infections and complications. Cell adhesion molecules may be involved in the interaction between bacteria and implant surfaces, affecting biofilm formation and stability. In particular, Staphylococcus aureus infection can occur.¹⁷ The PI3K–Akt pathway is a key signalling pathway common to both BC and BII. Dysregulation of the PI3K–Akt signalling pathway has been linked to increased cell proliferation, metabolic changes, and resistance to treatment in breast cancer.²⁵ In addition, a study proposed the biofilm hypothesis for BII, which mentions the role of macrophages in the peri-implant tissue response, including possible biofilm formation and macrophage activation.²⁶ The activation of the PI3K–Akt signalling pathway may influence macrophage behaviour,²⁷ and these cells play a role in the response of tissues surrounding implants. In our study, we found a greater number of macrophages in patients in the BII group than in those in the non-BII group by immune infiltration analysis. More definitive conclusions still need to be justified in deeper studies.

The KRT14 gene encodes keratin 14, an intermediate filament protein that belongs to a member of the keratin family. KRT14 is expressed predominantly in epithelial cells and forms heterotetramers with keratin 5 (KRT5), which together constitute the cytoskeleton of cells,²⁸ and KRT14 is deregulated in a variety of cancer types and participates in a common protein‒protein interaction network, suggesting that it may play a role in cancer invasion and potentially in cancer invasion and metastatic processes. Elevated levels of KRT14 have been correlated with decreased progression-free survival (PFS) in patients with ovarian cancer, particularly among those in the early stages (stages I–II) of the disease. In addition, elevated KRT14 levels were linked to reduced PFS. Treatment with platinum- and paclitaxel-based chemotherapy is a negative prognostic factor for PFS. This finding implies that KRT14 levels serve as a distinct prognostic indicator for individuals with high-grade serous ovarian cancer.²⁹ In triple-negative breast cancer, upregulation of KRT14 is associated with the regulation of EZH2–H3K27me3, and this upregulation promotes peritoneal metastasis in TNBC. In addition, KRT14 is also a characteristic marker in basal-like breast cancers, which are often associated with BRCA1 mutations and exhibit increased aggressiveness and resistance to chemotherapy.³⁰ In our investigation, KRT14 exhibited notably elevated expression in healthy tissues compared with that in breast cancer tissues across both the training and validation cohorts, and the same findings were obtained through assessments of pathological tissues at our site. The biological function of KRT14 was contrary to its expression in cancer, which suggests that we need to evaluate it in a specific biological and pathological context. Whereas KRT14 expression was also low in BII according to correlative analysis of BII, it was also validated in pathological breast tissue from our medical centre, suggesting that KRT14 expression could be used as a relevant diagnostic marker for BII.

The KIT protein, which is encoded by the KIT gene, functions as a type III receptor tyrosine kinase and plays a role in cellular signalling; it is typically activated upon interaction with its ligand, stem cell factor.³¹ In gastrointestinal stromal tumours, approximately 85% of these malignancies feature mutations in the KIT gene, which is a key driver. These mutations lead to continuous activation of the KIT protein, thereby facilitating the proliferation of cancer cells.³² Owing to the pivotal role of the KIT gene in numerous malignancies, therapeutic agents that specifically inhibit the KIT gene have been the subject of investigation for cancer therapy.³³ However, tumour progression is linked to reduced levels of c-KIT expression, which is thought to contribute to the early phases of human breast cancer development. A prior investigation evaluated KIT expression levels in 60 samples of both breast cancer and healthy breast tissue through immunohistochemistry, suggesting that reduced KIT expression is correlated with the shift of breast epithelial cells to a malignant state. The prognostic and predictive significance of KIT has been examined in individuals with high-risk breast cancer. Specifically, KIT expression was examined in 236 patients using tissue microarrays, which revealed reduced expression in 12% of the patients with breast cancer.³⁴

The ALB gene produces albumin, the most plentiful protein circulating in the bloodstream, constituting approximately half of the total protein content in the plasma of healthy individuals.³⁵ Research indicates that reduced expression levels of the ALB gene might be linked to a poorer prognosis in numerous conditions, such as various types of cancer. As a negative acute-phase response protein, a reduction in the serum ALB gene is a threat parameter for long-term survival in a variety of clinical settings and a strong biomarker for poor prognosis in most diseases.³⁶ A study investigated the immunoregulatory mechanism of cordycepin in breast cancer and revealed that the ALB gene is a key target of cordycepin. ALB gene expression is suppressed in tumours, whereas cordycepin is able to increase ALB gene expression and thus exert antitumour effects.³⁷ The ALB gene has a potential role in cancer and is important in immune regulation, particularly in breast cancer. These findings suggest that the ALB gene may become a novel target for cancer therapy.

When the immune infiltration microenvironments of both patients with BC and patients with BII were examined, distinct patterns emerged. Different immune cells are activated in BC. Our study revealed abnormal upregulation of NK cells and B cells in breast cancer tissues. NK cells and B cells play complex roles in the immune microenvironment of breast cancer, and they can act as antitumour effector cells or promote tumour escape and metastasis in some cases. For example, NK cells play an important role in immune surveillance against tumour cell generation, and they are able to recognize and clear senescent, apoptotic cells as well as variant tumour cells.³⁸ Within the immunosuppressive milieu of the tumour, NK cells might experience dysfunction as a result of exposure to inhibitory molecules released by cancer cells, potentially resulting in tumour evasion.³⁹ However, tumour-dominant B cells can promote the lymph node metastasis of patients with breast cancer by secreting pathological antibodies targeting the tumour antigen HSPA4. These antibodies are able to bind the glycosylation site of HSPA4 and activate the downstream NF-κB pathway, thereby promoting the formation of the microenvironment before tumour metastasis and inducing the migration of tumour cells into the draining lymph nodes.⁴⁰ Our study revealed that eosinophils and CD56dim natural killer cells were significantly different between the two groups in the BII group. Eosinophils may play a role in antitumour function through direct cytotoxic effects or indirect immunomodulatory effects. For example, they are able to attract CD8+ T cells into cancer tissue by secreting chemokines and normalizing the tumour vasculature.⁴¹ Although there is no direct evidence linking eosinophils or NK cells to BII, the role of eosinophils and NK cells in the immune microenvironment may be related to chronic inflammation and immune responses in BII. Future studies may further investigate the link between these two factors.

In summary, our research indicates that BC and BII share common molecular pathways during development. Bioinformatics analysis revealed enrichment of several biological pathways, such as those involving cell adhesion molecules, axon guidance, and Staphylococcus aureus infection, indicating that disruptions in these pathways could play a role in the development of BII; conversely, the PI3K‒Akt signalling pathway, cell cycle, and tight junction pathways might be involved in the progression of BII among patients with breast cancer. Our study has several limitations, such as the lack of validation of additional gene sets associated with BII using our established prognostic model. Second, the relevant genes we screened were not validated at the cytological level. Nonetheless, our findings offer novel insights into the molecular underpinnings of BII–BC development and could aid in identifying therapeutic targets for this condition. Nevertheless, additional research is essential to confirm these results and investigate the functional impact of the identified genes and pathways in the pathogenesis of BII–BC.

Conclusion

Our research meticulously identified three potential key genes (KRT14, KIT, and ALB) and established a predictive tool for distinguishing BII from BC using an array of bioinformatics analysis and machine learning algorithms. Additionally, we have developed a predictive model for distinguishing BII from BC that can support clinical decision processes. In summary, our results might offer novel perspectives on the pathogenesis and diagnosis of BII with BC.