Expression characterization of the guanylate-binding protein gene family in breast cancer and its association with the immune microenvironment

Data collection

We obtained BRCA clinical data and RNA-Seq data in FPKM from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). This was achieved by removing samples without survival time or status, ensuring that the survival time of each patient was longer than 30 days but shorter than 10 years. Ultimately, 994 BRCA samples and 114 paraneoplastic control samples were included. Here, we focused on protein-coding genes due to their well-defined biological functions, stronger interpretability and comprehensive annotation information. Additionally, GSE20685 data was obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. In total, 327 BRCA samples were analyzed after the probes were converted to symbol using the annotation file while samples without overall survival (OS) data or clinical follow-up were eliminated. GSE20685 dataset and TCGA-BRCA dataset served as independent validation set and training set, respectively. Furthermore, seven GBPs were collected from previous research literature for further investigation (Ning et al., 2023).

Differences in the mutation features of GBPs in BRCA samples

Based on the TCGA dataset, variations in the expressions of GBP genes between BRCA and paraneoplastic samples were analyzed. Next, ssGSEA (Hanzelmann, Castelo & Guinney, 2013) method in the R package “GSVA” was employed to calculate the GBP score for each sample in the TCGA_BRCA cohort, with the candidate GBPs as the background gene set. The BRCA samples were assigned into low- and high-GBP scores by the median value. The GBP gene mutations for each sample in the TCGA cohort were analyzed as genomic mutations are strongly linked to the progression of disease (Zahn & Travis, 2015). MuTect 2 software (Beroukhim et al., 2010) was utilized to process the mutation dataset of BRCA samples and paraneoplastic samples downloaded from the TCGA database and to map gene mutations in the low- and high-GBP subgroups.

Construction of weighted gene co-expression network

Gene modules closely linked to the GBP score were classified by WGCNA. To more effectively detect strong correlations between modules, the pickSoftThreshold function of “WGCNA” was employed to decide the optimal soft threshold power (β = 7) (Langfelder & Horvath, 2008). Under the threshold of minModuleSize = 60, hierarchical cluster analysis was performed to identify gene modules. According to the first principal components (PCs) of module expression, the “Heatmap” package (Saunders, Liang & Li, 2007), several module-trait genes were selected. Then, the correlation between module-trait genes and clinical trait for diagnosis was analyzed to evaluate the association between the modules and GBP scores. For modules showing the strongest module-trait relationships, the genes contained in the modules were tested for further analysis.

Enrichment analysis

Subsequently, the “limma” package was employed to screen DEGs between tumor samples and paraneoplastic normal samples in the TCGA_BRCA cohort under |log₂ (Fold Change)| > log₂(2) with an adjusted p < 0.05 as the screening criteria. Then, genes (p < 0.05) in the intersection between the selected modular genes and DEGs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional analyses using the R package “clusterProfiler” (Yu et al., 2012; Xia et al., 2024). Top functions enriched in biological process (BP) term in GO analysis and the top enriched KEGG pathways were visualized into bar graphs.

Risk modeling and validation

The intersecting genes were subjected to univariate Cox proportional risk regression in the R package “survival” (Therneau & Lumley, 2015) to select genes closely (p < 0.05) linked to the prognostic outcomes of BRCA patients in the TCGA_BRCA cohort (Zhang et al., 2025). The final genes for developing an effective risk model were refined applying 10-fold cross-validation and LASSO Cox regression analysis with the “glmnet” package (Friedman, Hastie & Tibshirani, 2010). The RiskScore model was developed by Multivariate stepwise regression analysis to screen important genes that showed a correlation coefficient independently linked to the prognosis of BRCA in the TCGA_BRCA dataset. The formula of the risk model was: RiskScore = Σβi × Expi, where βi represents the multivariate regression Cox analysis coefficient of each gene, and Expi represents the expression of each gene. Following zscore normalization, the optimal threshold value of Riskscore was used to divide the patients in the TCGA_BRCA dataset into low- and high-risk groups. The R package “survminer” (Ozhan, Tombaz & Konu, 2021) was then used to conduct survival analysis between high- and low-risk groups, and Kaplan–Meier (KM) survival curves were plotted for prognostic analysis. Furthermore, the diagnostic accuracy of the Riskscore in the TCGA-BRCA training set was assessed according to receiver operating characteristic analysis (ROC) using the R package “timeROC” (Blanche, 2015) and area under ROC curve (AUC) of 1-, 3- and 5-year. Next, the performance of the Riskscore was similarly validated in GSE20685.

Functional characterization and immune infiltration analysis of BRCA patients in different risk groups

The DEGs between the two risk groups were filtered by the “limma” package under the filtering criterion |log₂ (Fold Change)|> log₂(2) and an adjusted p < 0.05. The DEGs were then subjected to functional annotation through GO and KEGG pathway enrichment analyses using the DAVID online tool (https://davidbioinformatics.nih.gov/). Briefly, the DEGs were loaded into DAVID and the organism was selected as “Homo sapiens” to conduct GO (in three terms) and KEGG functional enrichment analysis. The gene-correlated pathways were examined by KEGG enrichment analysis. Only entries with p < 0.05 were selected.

The ssGSEA function of the R package “GSVA” was employed to calculate the scores of 28 types of immune cells (Charoentong et al., 2017) for the two risk subgroups, so as to analyze the correlation between Riskscore and immune function in different BRCA patients. The immune infiltration for the TCGA-BRCA dataset samples was assessed by calculating StromalScore, ImmuneScoreh and ESTIMATEScore with the R package “estimate” (Yoshihara et al., 2013). Then, we collected 29 gene signatures that represent the primary functional components of the tumor and immune, stromal, and other cell populations (Bagaev et al., 2021). The correlations between gene signatures and Riskscore were analyzed using ssGSEA (Liu & Huang, 2023).

Cell culture and plasmid transfection

From the American Type Culture Collection (ATCC, Manassas, VA, USA), we ordered human BRCA cells (MDA-MB-231 and AU565) and human normal mammary epithelial cells (MCF-10A). DMEM medium (11965092, Gibco, Waltham, MA, USA) and RPMI 1640 (11875093, Gibco, USA) medium were used to culture MDA-MB-231 cells and AU565 cells, respectively. All the cultures were added with 1% penicillin streptomycin (15140122, Gibco, USA) and 10% heat-inactivated fetal bovine serum (FBS, 10099141, Gibco, USA). Cells were cultured at 37 °C and 5% CO₂ in an incubator.

Following the instructions, DACT2 overexpression plasmid (oe-DACT2) or control plasmid (oe-NC) purchased from GenePharma (Shanghai, China) was transfected into MDA-MB-231 and AU565 cells in logarithmic growth phase (2 × 10⁴ cells/well) using Lipo3000 Liposome Transfection Reagent (L3000-001, Thermo Fisher Scientific, Waltham, MA, USA). After transfection for 48 h, the cells were collected for further study.

RNA extraction and quantitative real-time PCR

Total RNA was separated from MDA-MB-231, MCF-10A, and AU565 cells using the RNA Extraction Kit (TRIzol, Invitrogen, Waltham, MA, USA), following the protocols. The purity and concentration of the total RNA was quantified, and the cDNA templates were generated using the HiScript II kit (R233-01, Vazyme, Nanjing, China). The qRT-PCR was conducted utilizing specific primers and the KAPA SYBR^® FAST kit (Sigma Aldrich, St Louis, MO, USA). Data were processed using the 2^−ΔΔCT method, with GAPDH as an internal control (Zhang et al., 2023). See Table 1 for the primer sequences used in the study.

Cell viability

Following the instructions, CCK-8 assay (Dojindo, Rockville, MD, USA) was applied to assess the effect of DACT2 on the viability of MDA-MB-231 and AU565 cells. Briefly, the cells were cultured for 0 (7,000 cells/well), 24 (5,000 cells/well), 48 (3,000 cells/well) and 72 (2,000 cells/well) h in 96-well microtiter plates. After washing the cells with PBS twice, 100 µL of fresh culture media and 10 µL of CCK-8 solution were added into each well for cell incubation at 37 °C with 5% CO₂ for 3 h. The absorbance at 450 nm was read by SPECTROstar^® Nano (BMG LABTECH GmbH, Ortenberg, Germany) (Tang et al., 2024).

Cell migration and invasion assays

The impact of DACT2 expression on the migration and invasion capacities of MDA-MB-231 and AU565 cells was also investigated. Wound healing assay was conducted to measure the cell migratory ability. Briefly, 6-well plates were added with transfected cells (5 × 10⁵/mL) and two mL of cell suspension for cell incubation at 37 °C with 5% CO₂. A 10 µL plastic pipette tip was utilized to scratch a uniform wound on the monolayer once the cell density reached approximately 80%. Following PBS wash, the monolayers were cultured in non-FBS medium. The wound edge between the two edges of the migrating cell sheet were measured and photographs at 0 and 48 h. Each experiment was carried out in triplicate. The upper transwell chamber coated with 10% Matrigel (Corning, Corning, NY, USA) contained 1 × 10⁵ cells. After incubating the cells for 24 h, those still in the upper chamber were removed. After that, the cells in the lower chamber was fixed by 4% paraformaldehyde and colored by 0.1% crystal violet solution. Finally, the invaded or migrated cells were quantified from six distinct fields of view under a microscope.

Statistical analysis

All the statistical analyses were conducted using R software version 3.6.0 (R Foundation, Vienna, Austria) and GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA). Differences between continuous variables in the two groups were calculated by Wilcoxon rank-sum test. Correlations were computed by Spearman method, and the log-rank test was applied to compare survival differences between two groups of patients. In cellular assays, multiple group comparisons were performed using Student’s t-test or two-way ANOVA. For data that did not conform to normal distribution, non-parametric tests such as Mann–Whitney U test or the Kruskal–Wallis test were used for statistical analysis. A p < 0.05 stood for statistically significant difference.

GBP gene expression and mutation in BRCA

In the TCGA cohort, analysis on the expressions of GBPs showed that the levels of GBP2, GBP4 and GBP6 were lower in BRCA samples than in paraneoplastic samples (p < 0.0001), while the expressions of GBP3 and GBP5 were notably higher in BRCA samples than paraneoplastic samples (p < 0.05, Fig. 1A). The GBP score for each sample in the TCGA_BRCA dataset was computed by the ssGSEA method, and the patients were divided by the median value into high and low GBP score groups. Prognostic study demonstrated that BRCA patients with high GBP scores had longer survival time, whereas those with low GBP values had a worse OS (p = 0.014, Fig. 1B). Next, analysis on the GBP gene mutations between BRCA and paraneoplastic samples in TCGA cohort showed that only GBP4 (2%) and GBP5 (1%) genes were mutant in paraneoplastic tissues (Fig. 1C). In contrast, BRCA tissues had numerous mutations in GBP3 and GBP5 genes in the same sample. Additionally, we identified mutations in GBP7 (1%) and GBP3 (1%) genes, along with the previously mentioned GBP4 (1%) and GBP5 (1%) (Fig. 1D).

WGCNA identified GBP-related gene modules

Gene expression modules linked to GBPs in the TCGA-BRCA cohort were identified by the R package “WGCNA”. The soft threshold power of 7 was selected to develop a topological network to ensure the scale-free topology of the network (Fig. 2A). Ultimately, we obtained 16 co-expression modules, with 60 genes in each module at least (Fig. 2B). Notably, the number of genes in the magenta, grey, and turquoise modules was comparatively high (Fig. 2C). Genes in the grey module could not be merged into other modules. In order to select clinically significant modules, we calculated the correlation of each module with GBP score and plotted a heat map of module-shape correlation. Light cyan and magenta modules were closely positively connected with the GBP score (cor = 0.73, p = 1.42e−165; cor = 0.77, p = 3.97e−192, Fig. 2D). Moreover, a strong association between MM and GS was detected in the lightcyan module (cor = 0.72, p = 3.6e−14) and magenta module (cor = 0.9, p < 1e−200) (Figs. 2E–2F).

Differential gene analysis and functional characterization

Subsequently, we identified DEGs between BRCA samples and paraneoplastic samples in the TCGA dataset and finally screened 4,215 DEGs (Fig. 3A). Subsequently, we found 393 present in the intersection between the lightcyan and magenta module genes and the DEGs (Fig. 3B). GO and KEGG was performed on the common genes to investigate their regulatory functions in BRCA pathogenesis. GO-BP analysis showed significant enrichment of the intersecting genes in immune-related functions (Fig. 3C), including leukocyte migration, T-cell activation, positive regulation of cytokine production. The KEGG analysis revealed that the intersecting genes were considerably enriched in the viral protein interaction with cytokines and cytokine receptors and cytokine-cytokine receptor interaction pathways (Fig. 3D).

Establishment of a risk model and verification

The intersecting genes were subjected to univariate Cox regression to identify GBP-linked DEGs with the most significant influence on the prognostic ourcomes of BRCA in TCGA-BRCA while removing redundant confounding genes. The screened genes were then further compressed based on LASSO analysis (Figs. 4A–4B). Multivariate stepwise regression analysis determined four distinctive genes (PSME2, DACT2, PIGR and STX11) independently linked to the prognosis of TCGA-BRCA patients (Fig. 4C). Then, the prognostic outcomes of TCGA-BRCA patients was predicted by the risk model: RiskScore = (−0.28*PSME2) + (−0.157*DACT2) + (−0.066*PIGR) + (−0.237*STX11). Additionally, based on the optimal critical value of RiskScore, TCGA-BRCA patients were classified into low- and high-risk groups. The OS of low-risk TCGA-BRCA patients was more unfavorable than that the high-risk patients, as shown by the KM curves (p < 0.0001, Fig. 4D). Utilizing the “timeROC” R package, ROC analysis for 1-, 3-, and 5- year prognostic prediction was performed to evaluate the effectiveness of RiskScore in prognostic characterization. The AUC results verified the classification accuracy of the RiskScore, with an AUC value of 0.73, 0.66 and 0.69 for 1-, 3-, and 5- year prognostic prediction, respectively (Fig. 4E). Analysis on the expressions of the identified genes in TCGA-BRCA patients indicated that the high-risk BRCA group had lower expressions of PSME2, DACT2, PIGR and STX11 than the low-risk group (Fig. 4F). We divided TCGA-BRCA patients by the median expression of the four genes into low and high expression groups, and examined the correlation between the expression of the four genes and patient prognosis. It was found that BRCA patients in the low expression groups of PSME2 (p = 0.00064), DACT2 (p = 0.00015), PIGR (p = 0.00093) and STX11 (p = 0.0001) had poorer prognosis and shorter survival time (Fig. 4G).

The robustness of the RiskScore in the GSE20685 dataset was tested applying the RiskScore and equivalent coefficients used in the analysis on the training set. The validation set showed similar results to those in the training set that high-risk BRCA patients had a worse prognostic outcome than BRCA patients with a low risk (p < 0.015, Fig. 4H). The AUC value of 1-, 2-, 3-, 4-, and 5- year prognostic prediction in GSE20685 validation set reached 0.69, 0.6, 0.58, 0.58 and 0.57, respectively (Fig. 4I). All the four prognostic genes were low-expressed in the high-risk category in the GSE20685 validation cohort (Fig. 4J), suggesting the accuracy of these genes as prognostic predictors for BRCA.

TCGA-BRCA enrichment analysis on the DEGs between the two risk groups

Next, 813 DEGs were filtered between low- and high-risk subgroups. The GO-BP enrichment analysis revealed that these DEGs were mainly implicated in inflammatory response, signal transduction, inflammatory response immune response, adaptive immune response, and adaptive immune response pathways (Fig. 5A). In CC terms, these genes mainly localized at the extracellular space, cell surface, external side of the plasma membrane, and some other structures (Fig. 5B). In MF terms, these genes were closely involved in signaling receptor activity and cytokine activity (Fig. 5C). The DEGs among different risk subgroups of BRCA mainly influenced cytokine activity, hematopoietic cell lineages, cell adhesion molecules, transmembrane signaling receptor activity, viral protein-cytokine, cytokine-cytokine receptor interactions, and cytokine receptor interactions (Fig. 5D).

Relationship between RiskScore and immune characteristics

Immune cell infiltration in TCGA-BRCA patients was analyzed to examine the differences in the immunological milieu of the patients across risk groupings. Analysis of the infiltration abundance of 28 immune cells in BRCA patients revealed that high-risk BRCA patients had higher infiltration of activated CD8 T cells, activated CD4 T cells, activated B cells than the low-risk group, while the infiltration of mast cells, macrophages, and effector memory CD8 T cells was lower in the high-risk group (Fig. 6A), indicating that high-risk BRCA patients had lower immunoreactivity. The ESTIMATE showed a negative correlation between RiskScore and StromalScore, ImmuneScore, and ESTIMATEScore (Fig. 6B). This suggested that immune cell infiltration may be lowered in patients in the BRCA high-risk group. Next, applying the ssGSEA method, it was observed that the majority of the 29 gene signatures were closely linked to PSME2, DACT2, PIGR and STX11. Furthermore, the RiskScore had a negative correlation with anti-tumor immune infiltrate but was irrelevant to tumor proliferation (Fig. 6C).

In vitro cell-based model to validate the expressions and potential biological functions of the key signature genes

According to the results of qPCR, MDA-MB-231 and AU565 cells had remarkably higher PSME2 levels than in MCF-10A cells, while both AU565 and MDA-MB-231 cells had significantly lower levels of DACT2 and STX11. However, only AU565 cells had notably downregulated expression of PIGR (Fig. 7A). Research showed that DACT2 could cause cellular G1/S phase inhibition and suppress BRCA cell proliferation, demonstrating a strong research basis and biological feasibility (Li et al., 2017). For this reason, we selected DACT2 for further verification. The CCK-8 assay revealed that AU565 and MDA-MB-231 cell viability was considerably downregulated by DACT2 overexpression (Figs. 7B–7C). The wound healing and transwell assays demonstrated that DACT2 overexpression markedly suppressed MDA-MB-231 and AU565 cell migration and metastasis (Figs. 7D–7G).