GNAL-driven calcium signaling reshapes the spatiotemporal immune landscape in ER⁺ breast cancer: causal insights and prognostic implications

Introduction

Breast cancer is one of the most common malignancies in women, with the estrogen receptor-positive (ER⁺) subtype accounting for approximately 70% of all cases (1,2).

Its treatment is highly dependent on endocrine therapy (3). The ER signaling pathway is the core molecular mechanism regulating the occurrence and development of ER-positive breast cancer, and its dysfunction is closely associated with tumor proliferation, invasion, and endocrine therapy resistance (4). The core regulation of this pathway relies on the functional balance between two subtypes, ERα and ERβ: as a key oncogenic molecule, ERα can directly bind to estrogen response elements (ERE) through classical genomic effects, activate the transcription of downstream target genes such as MYC and cyclin D1, and promote tumor cell cycle progression; it can also synergize with G protein-coupled estrogen receptor (GPER) through non-genomic effects to activate kinase pathways including PI3K/AKT and MAPK, accelerating cell proliferation and invasion (5). In contrast, ERβ exerts a tumor-suppressive role by forming a complex with p53 to regulate histone methylation, relieving the transcriptional repression of tumor suppressor genes by ERα; meanwhile, it can downregulate IL-8 expression by inhibiting NF-κB pathway phosphorylation, thereby suppressing tumor migration and promoting apoptosis (3). Furthermore, the crosstalk between ER signaling, extracellular matrix (ECM) remodeling, and the tumor immune microenvironment influences tumor progression and metastatic potential by regulating the secretion of matrix metalloproteinases (MMPs) and inflammatory factors (4). The development of CDK4/6 inhibitors has significantly prolonged patient survival. However, nearly 40% of patients experience disease recurrence or metastasis within five years of treatment, which remains a major clinical challenge (6). Resistant ER⁺ breast cancers commonly exhibit a classic “immune-cold” phenotype, characterized by a 60–70% reduction in CD8⁺ T-cell infiltration, abnormal accumulation of immunosuppressive myeloid cells, and defective development of tertiary lymphoid structures (TLSs) (7). Monotherapy with immune checkpoint inhibitors has shown limited efficacy in the ER⁺ subtype, with objective response rates below 10% (8), suggesting that tumor cells may actively remodel specific signaling networks to establish spatiotemporally distinct mechanisms of immune evasion (9).

Calcium signaling, as a central component of the intracellular second messenger system, plays a dual regulatory role in tumor biology. On one hand, calcium oscillations enhance T-cell antitumor activity by activating NFAT and NF-κB pathways (10); on the other hand, mitochondrial calcium overload can trigger ferroptosis, leading to the release of tumor-associated antigens and promoting immunogenic cell death (11). In ER⁺ breast cancer, aberrant expression of key calcium signaling molecules, such as SPCA1, has been linked to endocrine therapy resistance (9). However, the molecular mechanisms underlying calcium signaling in remodeling the tumor immune microenvironment in a spatiotemporal manner remain largely unexplored. Recent genomic studies have identified GNAL, a Gα subunit of heterotrimeric G proteins, as a novel regulatory factor in calcium signaling with a significant causal association in ER⁺ breast cancer genomes (12,13). Loss of GNAL expression facilitates the maintenance of stem cell-like properties in cancer cells, which may in turn drive the spatiotemporal reprogramming of an immunosuppressive microenvironment. Nevertheless, its precise signaling pathways and molecular partners remain to be elucidated.

Recent advances in single-cell and spatial transcriptomics have provided unprecedented insights into the spatial organization of distinct cell populations within tumors. In triple-negative breast cancer, Wang et al. successfully applied a single-cell analytical framework to identify metastasis-associated subpopulations specifically located at the tumor margin (14). Ye Kai’s group developed the STMiner algorithm, which integrates single-cell and spatial transcriptomic data to reconstruct the continuous trajectory of gene expression changes from tumor core to margin to adjacent normal tissue in colorectal cancer, revealing immune-suppressive niches at the tumor edge that are often undetectable using conventional approaches (12). Collectively, these findings suggest that simultaneously tracking the dynamic activity of calcium signaling and the spatial distribution of immune cells could uncover core mechanisms underlying drug resistance in ER⁺ breast cancer. Although the widely used PAM50 molecular classification system is effective in distinguishing subtypes such as Luminal A and B, fails to capture dynamic changes within the tumor microenvironment (TME). In addition, traditional single-cell sequencing lacks spatial resolution, making it difficult to decipher long-range regulatory networks between cells (8). Herein, we propose an innovative mechanistic model: the “calcium oscillation frequency-stemness-immune suppression” axis (10). Specifically, we hypothesize that GNAL remodels the spatial architecture of the immune microenvironment in ER⁺ breast cancer through calcium signaling, where the GNAL-mediated long-range interaction network between stem-like and immune cells is a key determinant of clinical prognosis. “Stem-like cells” specifically refers to a type of cell population localized in the tumor stromal region, whose gene expression profiles are highly similar to those of stem cells and have been confirmed to be in an undifferentiated or poorly differentiated state via trajectory analysis (15). In this study, we integrated Mendelian randomization, single-cell transcriptomics, and spatial transcriptomics to systematically investigate how GNAL orchestrates the TME at a distance, which may provide novel therapeutic targets to overcome endocrine resistance. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1796/rc).

Methods

Summary-based Mendelian randomization (SMR)

This study utilized genome-wide association study (GWAS) data from individuals of European ancestry, which were obtained from the Integrative Epidemiology Unit (IEU) OpenGWAS database (https://gwas.mrcieu.ac.uk/). The specific dataset analyzed was “ieu-a-1127”, including 69,501 ER⁺ breast cancer cases and 105,974 controls. Instrumental variables (IVs) were selected according to established criteria (16).

Additionally, we acquired breast mammary tissue expression quantitative trait loci (eQTL) data from version 8 of the Genotype-Tissue Expression (GTEx) project. This eQTL dataset was derived from postmortem samples collected from 838 donors. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Single-cell transcriptomic data

We obtained the single-cell RNA sequencing dataset GSE161529 (comprising 11 ER⁺ breast cancer samples) from the NCBI Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

Cell clustering and annotation

Single-cell transcriptomic data were processed using Seurat (v4.3.0) with the following quality control criteria: (I) at the cellular level, we retained cells expressing 200–4,000 genes with mitochondrial gene content <10%; (II) for genes, we only kept those detected in ≥3 cells. Data normalization was performed using SCTransform, followed by batch effect correction across datasets with Harmony (v1.2.3) (17), yielding a final dataset of 69,789 cells and 33,538 genes. Dimensionality reduction and visualization were achieved through Uniform Manifold Approximation and Projection (UMAP). Cell type annotation and classification were comprehensively performed using SingleR package (v2.10.0). Finally, we generated GNAL expression distribution and density plots using Seurat’s AddModuleScore function.

Pseudotime and GNAL dynamics

VECTOR is an unsupervised inference tool based on UMAP that reconstructs cellular differentiation trajectories by modeling the network distance between cells and their starting population (18). In this study, we applied VECTOR to perform pseudotemporal ordering of single-cell data, thereby elucidating potential developmental trajectories during cell differentiation.

Spatial transcriptomic data

The spatial transcriptome data of ER⁺ breast cancer patients were acquired from the ‘Human Breast Cancer: Whole Transcriptome Analysis’ dataset, publicly available through 10x Genomics (https://www.10xgenomics.com/)

Transcriptomic data

In this study, transcriptomic and corresponding clinical data of ER⁺ breast cancer patients were retrieved from The Cancer Genome Atlas (TCGA) database using the TCGAbiolinks package (v2.36.0) (19-21), serving as the training cohort. After stringent quality control, a total of 807 ER⁺ breast cancer patients and 113 normal controls were included. Additionally, an independent validation cohort comprising 481 ER⁺ breast cancer patients was obtained from the GEO database (GSE199633). All expression data underwent standardized preprocessing: (I) log2 transformation to normalize data distribution, (II) conversion of probe IDs to official gene symbols based on platform annotation files, and (III) averaging of expression values for multiple probes mapped to the same gene.

Molecular docking data

Crystal structures of key proteins, including Gnal (PDB ID: 8HTG), KDM3B (PDB ID: 5R7X), and KRT14 (PDB ID: 6JFV), were retrieved from the RCSB Protein Data Bank (https://www.rcsb.org/). For proteins lacking experimentally resolved structures [e.g., CD24 (UniProt ID: P25063) and CEBPD (UniProt ID: P49716)], high-confidence predicted models were obtained from the AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/).

SMR/heterogeneity in dependent instruments (HEIDI) test

This study employed a summary-data-based Mendelian randomization (SMR) approach to systematically evaluate potential causal relationships between gene expression levels and ER⁺ breast cancer risk by integrating eQTL data. Using the SMR software (v1.3.1) (22-24), we first identified genes significantly associated with ER⁺ breast cancer risk (SMR test P<0.05), then excluded potential confounding effects of linkage disequilibrium through HEIDI testing (P>0.05). The final analysis yielded statistically significant protein-coding genes as potential causal genes.

Survival and differential expression analysis

We performed survival analysis on SMR-identified causal genes using the survival (v3.8-3) and survminer packages (v0.5.0). Optimal expression cutoffs were determined to stratify TCGA ER⁺ breast cancer patients (n=807) into high/low-expression groups, then assess associations with overall survival (OS) via Kaplan-Meier analysis with log-rank testing to identify prognostic genes. Differential expression analysis using limma (v3.64.1). Volcano plots were generated to visualize differential expression between 807 tumors and 113 controls, followed by intersection with prognostic causal genes to identify key genes showing concordant directionality between SMR β coefficients and expression changes. Venn diagrams illustrated these overlaps. Finally, we intersected these key genes with 254 calcium signaling pathway genes from Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) to identify pathway-relevant candidates.

Spatial interaction network

Downstream analysis of spatial data was performed using Seurat package (v4.3.0). The filtered spatial matrix was normalized using SCTransform package (v0.4.2) (25). We excluded data points outside the tissue boundary and removed low-quality spots (number of genes detected transcripts <500 genes/spot; mitochondrial content >20%). Cell type composition at each spatial location was inferred by deconvolving the pre-annotated single-cell matrix using robust cell type decomposition (RCTD) in spacexr package (v2.2.1), which integrates spatial data with single-cell reference profiles. Cellular enrichment patterns were assessed through neighborhood analysis (26). The three-layer spatial interactions network was constructed using MistyR (v1.10.0): (I) intra-layer (local cell neighborhood, k=50) to identify autocrine signaling such as VEGFA-KDR; (II) juxta-layer (adjacent regions, distance <100 µm) to analyze paracrine interactions, contributing 32% of the total interactions; (III) para-layer (long-range signaling, distance >200 µm). Finally, homotypic and heterotypic cell interaction networks were constructed using the dbscan package (v1.2.2). CellChat package (v1.6) was applied to analyze ligand-receptor interactions between various cell types. Spatial distribution patterns of ligand-receptor interactions were analyzed using the COMMOT package (v1.0.1) implemented in Python.

Identification of differentially expressed genes (DEGs)

We first isolated stem cell subpopulations from single-cell RNA-seq data and calculated GNAL expression scores for each cell using AddModuleScore function in the Seurat package. Cells were then stratified into GNAL-positive (scores ≥ median, n=96) and GNAL-negative (scores < median, n=395) groups based on the median expression level. Differential gene expression analysis was performed using the FindMarkers function in the Seurat package with stringent thresholds (adjusted P value <0.05 and |avg_log₂FC| >0.25). To ensure biological relevance, we excluded ribosomal genes and erythrocyte markers (e.g., HBA1, HBB) from the analysis. This approach yielded 102 high-confidence DEGs that likely reflect true biological differences between the two cell populations.

Batch effect correction and feature selection

Batch effects between the training and validation cohorts were harmonized using the limma package (v3.64.1), followed by principal component analysis (PCA) for quality control. In the training cohort, univariate Cox proportional hazards regression (P<0.05) identified 13 survival-associated genes, with results visualized via forest plot. Least absolute shrinkage and selection operator (LASSO) regression was implemented with the glmnet package (v4.1-9) to refine the signature into 10 key prognostic genes through 10-fold cross validation.

Prognostic model construction

A multivariate Cox proportional hazards model was constructed using stepwise regression (both forward and backward selection), yielding an optimal 4-gene prognostic signature. Risk scores were calculated as: Risk score = Σ(βi × Expi), where β represents the regression coefficient and Exp denotes gene expression levels.

Model validation

The prognostic model was independently validated in the validation cohort. Time-dependent receiver operating characteristic (ROC) analysis was performed using the timeROC package (v0.4) (27) to evaluate predictive accuracy at 1-, 3-, and 5-year intervals. Kaplan-Meier survival curves were generated using the survival package, with statistical significance assessed by log-rank tests (P<0.05 considered significant). The independence of the risk score from clinical covariates (age, stage, etc.) was confirmed through multivariate Cox regression analysis.

Immune infiltration

To validate the prognostic model’s association with immune infiltration characteristics, we performed comprehensive analyses using the IOBR package (v0.99.0). First, immune cell proportions were deconvoluted using MCPcounter to quantify 10 major immune cell populations. The ESTIMATE algorithm was simultaneously applied to evaluate TME composition through stromal scores, immune scores, and tumor purity estimates. Furthermore, we employed single-sample Gene Set Enrichment Analysis (ssGSEA) via IOBR’s calculate_sig_score function to assess pathway activities, with particular focus on T-cell tumor suppression pathways such as PD-1 signaling and CTLA-4 inhibition.

Molecular docking analysis

To investigate potential interactions between GNAL and 4 candidate proteins, we performed molecular docking using the GRAMM-X web server (https://gramm.compbio.ku.edu/). The analysis was conducted in high-resolution geometric docking mode, which employs a global range molecular matching algorithm to systematically evaluate potential binding conformations. The GRAMM algorithm performs a comprehensive multidimensional search of relative molecular positions and orientations, calculating intermolecular energy potentials to identify optimal docking regions. Subsequently, we analyzed the predicted binding interfaces using PDBePISA (https://www.ebi.ac.uk/msd-srv/prot_int/pistart.html) to quantify binding energies and interface areas. All docking results were visualized and analyzed using PyMOL molecular graphics software.

siRNA transfection

Gene silencing was performed using small interfering RNA (siRNA) targeting GNAL and control siRNA as negative control (RiboBio, Guangzhou, China). Cells were seeded in 6-well plates at a density of 1×10⁵ cells per well and cultured for 24 h to reach 60–70% confluence. Transfection was carried out using Lipofectamine 3000 (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s protocol, with a final siRNA concentration of 50 nM. After 48 h of incubation, transfection efficiency was validated by quantitative real-time polymerase chain reaction (PCR). The corresponding siRNA sequences are provided in Table S1.

Reverse transcription quantitative PCR (RT-qPCR)

Total RNA was extracted from cultured cells using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions. RNA concentration and purity were determined using a NanoDrop spectrophotometer (Thermo Fisher Scientific). First-strand cDNA was synthesized from 1 µg of total RNA using the PrimeScript RT reagent kit (Takara Bio, Kusatsu, Shiga, Japan) with oligo(dT) and random hexamer primers. Quantitative PCR was performed in triplicate using SYBR Green Premix (Takara Bio) on a QuantStudio 5 real-time PCR system (Applied Biosystems, Thermo Fisher Scientific). The reaction conditions were as follows: initial denaturation at 95 ℃ for 30 s, followed by 40 cycles of 95 ℃ for 5 s and 60 ℃ for 30 s. Gene expression levels were normalized to GAPDH and calculated using the 2^-ΔΔCt method. All primer sequences are listed in Table S2.

Cell Counting Kit-8 (CCK-8) assay

Cell viability was assessed using the CCK-8 (Dojindo Laboratories, Kumamoto, Japan) according to the manufacturer’s protocol. Briefly, cells were seeded in 96-well plates at a density of 3×10³ cells per well and cultured for 24 h. Then, 10 µL of CCK-8 solution was added to each well, followed by incubation at 37 ℃ for 2 h. The absorbance was measured at 450 nm using a microplate reader (BioTek, USA). All experiments were performed in triplicate and repeated three independent times.

Stromal cell-B lymphocyte co-culture

A transwell co-culture system (Corning, Corning, NY, USA) with 0.4 µm pore membranes was employed to investigate stromal cell-B lymphocyte interactions. GNAL-knockdown stromal cells were seeded in the lower chamber at a density of 5×10⁴ cells/well, while Raji B lymphocytes (ATCC, Manassas, VA, USA) were placed in the upper chamber at 2×10⁴ cells/well. Cells were maintained in RPMI-1640 medium supplemented with 10% FBS and 1% penicillin-streptomycin. For rescue experiments, recombinant human MIF protein (PeproTech, Cranbury, NJ, USA) was added to the lower chamber at a concentration of 50 ng/mL. The co-culture system was maintained at 37 ℃ in a 5% CO₂ atmosphere for 48 h before subsequent analyses. For cell proliferation assessment, cells were sorted by flow cytometry following co-culture and subsequently analyzed using the CCK-8 assay according to the protocol described in the previous section.

Statistical analysis

SMR was performed using the SMR software package (v1.3.1). Single-cell transcriptome analysis, spatial transcriptome analysis, transcriptome analysis, and machine learning were performed using the R statistical environment (v4.3.2) and the corresponding Bioconductor packages. Some spatial transcriptome analyses were implemented using Python (v3.7.9). For statistical comparisons, we employed one-way analysis of variance (ANOVA) for multi-group analyses and the Wilcoxon rank-sum test for comparisons between high- and low-risk groups. A two-tailed P value <0.05 was considered statistically significant for all hypothesis tests.

Results

Tumor suppressor effect and prognostic value of GNAL as a core gene of calcium signaling pathway in ER-positive breast cancer

By integrating GWAS data from the UK Biobank and BCAC (total sample size n=643,218), SMR analysis identified 598 protein-coding causal genes significantly associated with ER-positive breast cancer (P<0.05, HEIDI test P>0.05). Survival analysis of bulk RNA-seq data identified 46 genes as associated with the prognosis of patients with ER-positive breast cancer. Differential expression analysis between TCGA ER⁺ breast cancer samples (n=807) and controls (n=113) identified 854 DEGs. The intersection of survival-associated causal genes and DEGs yielded four candidate genes (Figure 1A). Further filtering based on concordance between SMR effect direction and differential expression resulted in three key genes: GNAL, CAMK2D, and RYR2 (Figure 1A). Volcano plots illustrated the differential expression profiles of these three key genes (Figure 1B). Kaplan-Meier survival curves demonstrated that all three function as tumor suppressor genes (Figure 1C). To screen out tumor suppressor genes related to the calcium signaling pathway, we cross-analyzed these three genes with the KEGG calcium signaling pathway gene set (254 genes), and the results showed that GNAL was the only overlapping gene (Figure 1D). The SMR regional association plot for GNAL (Figure 1E) identified three significant genes on chromosome 18: GNAL (ENSG00000141404), MPPE1 (ENSG00000154889), and RP11-820116.1 (ENSG00000267079), all showing significant associations with eQTLs. Further association analysis demonstrated a negative correlation between GNAL eQTLs and disease risk (slope >0, Figure 1F). These findings highlight GNAL as a core calcium signaling gene with tumor suppressive functions and prognostic significance in ER⁺ breast cancer.

Figure 1 Identification and validation of causal genes. (A) The Venn diagram of survival-associated genes and DEGs yielded four candidate genes. (B) Volcano plots illustrated the differential expression profiles of these three key genes: red/green/black dots represent significantly upregulated, significantly downregulated, and non-significant genes. (C) Kaplan-Meier survival curves of three key genes. (D) Intersection analysis with the KEGG calcium signaling pathway gene set. (E) The SMR regional association plot for GNAL. (F) Regional association plots of GNAL and disease risk. BRCA, breast cancer; DEGs, differentially expressed genes; eQTL, expression quantitative trait loci; GWAS, genome-wide association study; KEGG, Kyoto Encyclopedia of Genes and Genomes; SMR, summary data-based Mendelian randomization.

GNAL is a key regulator of ER⁺ breast cancer stem cells

We conducted single-cell RNA sequencing analysis on 11 ER-positive breast cancer samples without lymph node metastasis from the GSE161529 dataset. The UMAP visualization (Figure 2A) was generated using Seurat. Cell clusters were annotated referencing with SingleR package, resulting in the identification of eight distinct clusters: Epithelial cells, tissue stem cells, fibroblasts, macrophages, endothelial cells, T cells, monocytes, and B cells (Figure 2A). Using AddModuleScore function in Seurat package, GNAL expression distribution and density maps were plotted (Figure 2B,2C), revealing specific high expression of GNAL in tissue stem cells, with additional expression in endothelial cells and fibroblasts. We then isolated tissue stem cells from the entire single-cell dataset and obtained 3,891 cells expressing 33,538 genes. Using PCA and Harmony re-dimensionality reduction (ratio parameter =1:10), tissue stem cells were clustered into 10 different subgroups (Figure 2D). The density map showed that GNAL was mainly concentrated in groups 1, 8, and 7 (Figure 2E). Pseudotime trajectory analysis of the tissue stem cells was performed using the VECTOR algorithm (18), revealing a gradual decline in GNAL expression along the differentiation trajectory (Figure 2F). This suggests GNAL plays a critical regulatory role during the maintenance phase of stemness. According to the median expression of GNAL, the tissue stem cells were stratified into GNAL-positive stem cells and GNAL-negative stem cells, providing essential grouping criteria for subsequent cell-cell interaction analyses.

Figure 2 Spatiotemporal dynamics analysis at single-cell resolution. The UMAP visualization generated using Seurat (A, annotated; B, unannotated). (C) GNAL expression distribution and density maps generated by Seurat’s AddModuleScore function. (D) UMAP visualization identified 10 distinct stem cell subclusters. (E) GNAL expression patterns across these subclusters. (F) GNAL expression along the differentiation trajectory plotted by VECTOR algorithm. UMAP, Uniform Manifold Approximation and Projection.

Spatial distribution of GNAL⁺ stem-like cells suggests a potential tumor-suppressive role

Through strict quality control, we excluded low-quality detection sites (Figure 3A). Spatial gene expression maps showed that GNAL showed a wide expression pattern in ER⁺ breast cancer tissues, distributed in both tumor nests and stromal areas (Figure 3B). Notably, GNAL⁺ and GNAL⁻ stem-like cells showed significant spatial heterogeneity: GNAL⁺ stem-like cells were significantly enriched at the tumor-stroma interface (Figure 3C); GNAL⁻ stem-like cells were mainly concentrated in the tumor nests (Figure 3D). The preferential localization of GNAL⁺ stem-like cells at the tumor-stroma boundary suggests their potential involvement in modulating a tumor-suppressive microenvironment, possibly through crosstalk with stromal components.

Figure 3 Spatial interaction network analysis. (A) UMAP visualization revealed specific high expression of GNAL in the stromal regions. The deconvolution in the tumor-stroma interface (B) GNAL; (C) GNAL-positive stem-like cells; (D) GNAL-negative stem-like cells. UMAP, Uniform Manifold Approximation and Projection.

GNAL⁺ stem cells orchestrate tumor stromal crosstalk via spatial-specific cell interactions

The correlation patterns among various cell types across three spatial scales were analyzed using MistyR: intra (short-range), juxta (mid-range), and para (long-range). Our results showed that GNAL⁺ stem cells predominantly interacted with endothelial cells at short-range and mid-range (Figure 4A-4D), while at long-range, they primarily engaged with endothelial cells, B cells, and epithelial cells, whereas GNAL⁻ stem cells mainly interacted with monocytes (Figure 4E,4F). The consistent high correlation between GNAL⁺ stem cells and endothelial cells across all ranges suggests a strong interaction between these cell types. We performed cell degree analysis on GNAL⁺ stem cells to construct homotypic cellular interaction networks. Focusing on spots containing GNAL-positive cells, we identified 113 spots, each surrounded by 6 neighboring spots. The higher the number of GNAL-positive spots within these neighbors, the greater the interaction degree of the region. We set a threshold at 0.1 to define high-degree areas. Our results revealed that homotypic interactions among GNAL⁺ stem cells were concentrated in the tumor stromal region (Figure 5A). Given the observed correlations between GNAL⁺ stem cells and endothelial cells, B cells, and epithelial cells in MistyR analysis (Figure 4E,4F), we constructed heterotypic interaction networks and computed enrichment scores accordingly using cell degree (Figure 5B-5D), which confirmed that these interactions were also localized in the tumor stromal region. Then we analyzed cell-cell communication between GNAL⁺ and GNAL⁻ stem cells with CellChat analysis. We identified significant interactions between GNAL-positive cells and B cells mediated by the MIF-CD74-CXCR4 axis and MIF-CD74-CD44 axis, whereas these pathways were not significantly active in communications between GNAL-negative cells and B cells (Figure 5E). GNAL⁺ stem cells primarily interacted with endothelial cells was mediated by the CD99-CD99 pair, whereas CD99-CD99 pair was not significantly active in communications between GNAL-negative cells and endothelial cells (Figure 5F). In contrast, GNAL⁻ and GNAL⁺ stem cells exhibited similar ligand-receptor interactions with epithelial cells, though GNAL⁻ stem cells showed stronger signaling intensity (Figure 5G). COMMOT analysis of these differentially expressed ligand-receptor pairs indicated that MIF-CD74-CXCR4 and MIF-CD74-CD44 signaling originated from the tumor core and extended toward the stromal region (Figure 6A-6D), whereas CD99-CD99 interactions lacked sufficient directional strength. These findings suggest that GNAL⁺ stem cells recruit B cells to the tumor stroma via MIF-CD74-CXCR4 and MIF-CD74-CD44 pathways, potentially contributing to antitumor immune suppression.

Figure 4 Spatial transcriptomic MistyR analysis. The correlation patterns among various cell types across three spatial scales were analyzed using MistyR: (A,B) intra (short-range), (C,D) juxta (mid-range), and (E,F) para (long-range).

Figure 5 Spatial transcriptomic cell degree analysis and cell communication. (A) CELL_DEGREE analysis defined high-degree areas. (B) Heterotypic networks between GNAL-positive cells and B cells. (C) Heterotypic networks between GNAL-positive cells and endothelial cells. (D) Heterotypic networks between GNAL-positive cells and epithelial cells. Signal pathway interaction networks were constructed for GNAL-positive cells with B cells (E), endothelial cells (F), and epithelial cells (G).

Figure 6 COMMOT analysis. The expression patterns and directional signaling of MIF-CD74-CXCR4 (A,B). The expression patterns and directional signaling of MIF-CD74-CD44 (C,D).

Construction of a prognostic model for ER⁺ breast cancer based on gene expression

Gene differential expression analysis was performed on GNAL-positive stem cells and GNAL-negative stem cells, and 102 DEGs were obtained after excluding ribosomal genes and erythroid markers. We integrated transcriptomic data from TCGA (807 ER-positive breast cancer samples) and GSE199633 (481 ER-positive breast cancer samples) as the training and validation cohorts, respectively. Batch effects were effectively eliminated after Harmony correction, as confirmed by PCA (Figure 7A). Univariate Cox regression analysis identified 13 genes significantly associated with OS from 102 DEGs (Figure 7B). Subsequent dimensionality reduction via LASSO regression narrowed these to 10 core prognostic genes (Figure 7C,7D). Stepwise multivariate Cox regression further constructed a prognostic model using 4 genes (CD24, KDM3B, CEBPD, and KRT14) (Figure 7E). By calculating the risk score of breast cancer patients, we divided the patients into high-risk and low-risk groups. Kaplan-Meier survival analysis demonstrated significantly shorter OS in high-risk ER⁺ breast cancer patients across both cohorts (Figure 7F,7G). ROC curve analysis revealed good predictive performance, with area under the curve (AUC) values of 0.669, 0.698, and 0.769 for 5-, 10-, and 15-year survival in the training set, and 0.594, 0.582, and 0.823 in the test set, respectively (Figure 7H,7I).

Figure 7 Construction of prognostic model. (A) Batch effects were effectively eliminated and confirmed by PCA. (B) Univariate Cox regression analysis identified 13 genes significantly associated with overall survival. (C,D) Dimensionality reduction via LASSO regression. (E) Forest plot of 4 independent prognostic genes. Kaplan-Meier analysis of training (F) and validation cohorts (G). The AUC curves of training (H) and validation cohorts (I). AUC, area under the curve; CI, confidence interval; FPR, false positive rate; LASSO, least absolute shrinkage and selection operator; PCA, principal component analysis; TPR, true positive rate.

Prognostic model risk score is an independent prognostic factor

Univariate Cox regression analysis was performed to compare the prognostic risk score with traditional clinical factors, including age, tumor size (T), lymph node status (N), metastasis (M), and overall stage (Figure 8A). The risk score demonstrated a strong association with survival, showing a higher hazard ratio than T stage and N stage, indicating its superior prognostic value as an integrated biomarker. Multivariate Cox regression analysis confirmed that the risk score remained an independent prognostic factor after adjustment for age and tumour, node, and metastasis (TNM) stage, supporting its usefulness as a supplementary clinical tool (Figure 8B). Subgroup analysis divided patients into high-risk and low-risk groups according to the median risk score. The results showed that the risk score was an independent risk factor for early (stage I–II) and late (stage III–IV) patients. In both early and late patients, the 5-year survival rate of the high-risk group was lower than that of the control group (Figure 8C,8D). Notably, in the age subgroup analysis, the risk score showed greater prognostic discrimination in patients aged ≤50 years compared to those aged >50 years, potentially reflecting higher tumor heterogeneity and more pronounced molecular driver effects in younger patients (Figure 8E,8F). These findings suggest that the biological significance of the risk score in specific populations may be mediated by microenvironment remodeling or activation of metabolic pathways.

Figure 8 Validation of prognostic model risk score. (A) Univariate Cox regression analysis was performed to compare the prognostic risk score with traditional clinical factors; (B) Multivariate Cox regression. The survival analysis of early-stage patients (C), late-stage patients (D), patients aged ≤50 years (E), and those aged >50 years (F). CI, confidence interval; HR, hazard ratio.

Positive correlation between prognostic model risk score and immunosuppression

Given our prior findings demonstrating interactions between GNAL⁺ stem cells and B cells, we further developed a prognostic model for ER⁺ breast cancer based on GNAL expression dynamics in stem cells. To evaluate the immunological implications of this model, we performed immune infiltration analysis using MCP-counter deconvolution. Notably, the high-risk group (Risk-high) exhibited reduced infiltration levels of both B cells and endothelial cells (Figure 9A-9C), aligning with spatial transcriptomic data that identified interactions between GNAL⁺ stem cells and these populations.

Figure 9 Validation of immune infiltration associated with the prognostic model. MCPcounter deconvolution analysis performed by IOBR package (A-C). The ESTIMATE algorithm was applied to assess microenvironment composition (D-F). (G) ssGSEA analysis of T cell functional pathways. (H) Spatial distribution analysis of immune checkpoint molecules. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. ssGSEA, single-sample Gene Set Enrichment Analysis; TCGA, The Cancer Genome Atlas.

Next, we applied the ESTIMATE algorithm to assess the prognostic model’s association with immune and stromal cell composition in the TME. The high-risk group showed decreased immune and stromal scores (reflecting diminished local barriers; Figure 9D-9F), suggesting an immunosuppressive TME with reduced supportive stroma, a phenotype consistent with “cold tumor” characteristics. ssGSEA of T-cell functional pathways further revealed attenuated tumor-suppressive T-cell activity in the high-risk group (Figure 9G), indicating impaired T-cell initiation and priming (TIP) mechanisms.

Spatial profiling of immune checkpoint molecules demonstrated significant enrichment of exhaustion markers (BTLA) and inhibitory receptors (LAG3) in the high-risk group, alongside downregulation of the T-cell activation promoter VTCN1 (Figure 9H). Importantly, these findings are consistent with single-cell pseudo-temporal analysis that links GNAL downregulation to stem cell differentiation arrest. This forms a coherent mechanistic circuit, suggesting that GNAL downregulation mediates the crosstalk between stem cell dysregulation and immune evasion.

Protein-protein interactions in ER⁺ breast cancer

Given the established association between GNAL and ER⁺ breast cancer at the protein level (eQTL) and considering that the prognostic model constructed based on GNAL expression differences serves as an independent prognostic factor for ER⁺ breast cancer, we hypothesized potential protein-level interactions between GNAL and the four genes (CD24, KDM3B, CEBPD, and KRT14) included in the prognostic model.

To investigate this hypothesis, we conducted molecular docking experiments at the protein level between GNAL and these four genes. Protein binding predictions were performed using GRAMM, while binding energies at the interaction sites were calculated using PDBePISA. Molecular docking visualizations were generated using PyMOL (Figure 10A-10D). Our results demonstrate that GNAL can bind with all four proteins: CEBPD, CD24, KRT14, and KDM3B. PDBePISA analysis revealed binding energies (ΔiG) of −25.1, −17.4, −13.5, and −6.3 kcal/mol for GNAL’s interactions with CEBPD, CD24, KRT14, and KDM3B, respectively. The negative values of these binding energies indicate stable molecular interactions, with more negative values corresponding to greater stability. These findings suggest that GNAL interacts with all four proteins at the molecular level, forming stable complexes.

Figure 10 Molecular docking analysis. The prediction of the combination of proteins. (A) GNAL and CEBPD. (B) GNAL and CD24. (C) GNAL and KRT14. (D) GNAL and KDM3B. (E) GNAL knockdown efficiency was evaluated by RT-qPCR. (F) The proliferation of GNAL-delepted MCF-7 cells detected by CC-8 assay. (G) The effect of co-culture with GNAL-knockdown stromal cells on the invasion of Raji cells, along with the impact of MIF protein supplementation, was evaluated using a Transwell assay. (H) CCK-8 assay detecting Raji cell proliferation after co-culture with GNAL-knockdown stromal cells and recombinant MIF supplementation. *, P<0.05; **, P<0.01; ***, P<0.001. CCK-8, Cell Counting Kit-8; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

To further validate the function of GNAL, we performed knockdown using small interfering RNA (siRNA) transfection. Among the tested sequences, siGNAL-1 and siGNAL-2 exhibited the highest interference efficiency and were selected for subsequent experiments (Figure 10E). The CCK-8 assay revealed that GNAL knockdown significantly enhanced the proliferative capacity of MCF-7 cells, suggesting a suppressive role of GNAL in tumor cell proliferation (Figure 10F).

Based on previous findings indicating that GNAL influences B-cell biological functions via MIF regulation, we co-cultured GNAL-knockdown stromal cells with Raji cells. Transwell assays demonstrated a significant reduction in the migratory ability of Raji cells. Notably, the addition of recombinant MIF protein partially restored the impaired invasion capacity (Figure 10G). Consistent with this, CCK-8 assays confirmed that co-culture with GNAL-knockdown stromal cells enhanced B-cell proliferation, whereas supplemental MIF suppressed this proliferative effect (Figure 10H). These results indicate that GNAL modulates B-cell invasion and proliferation through secreted MIF, thereby validating our earlier analytical findings.

Discussion

This study integrates multi-omics to decode GNAL’s role in shaping ER⁺ breast cancer’s immune-cold TME via calcium signaling. Combining Mendelian randomization, single-cell trajectory analysis, and spatial modeling, we reveal how stem cell dynamics and long-range cellular interactions drive immune evasion, providing a unified framework to understand therapy resistance and novel targeting strategies.

By combining GWAS data with Mendelian randomization, we anchored the genetic risk of ER⁺ breast cancer on GNAL. Our data indicates that for every unit decrease in GNAL expression, recurrence risk increases by 60%. This finding surpasses traditional GWAS’s focus on SNP loci by achieving precise functional gene localization from genetic variation. Multi-dimensional validations confirmed high colocalization of GNAL with disease risk loci and significantly poorer prognosis for patients with low GNAL expression. Pathway enrichment further revealed that GNAL’s biological effects are specifically concentrated on the calcium signaling pathway. Diverging from previous studies focusing on calcium channels or pumps, we found that GNAL, as a G protein α-subunit, regulates calcium oscillation frequency to influence stem cell differentiation. This mechanism resonates with findings in small-cell lung cancer by Chen et al. (28) and is substantiated here by spatial transcriptomics showing GNAL-positive stem cells specifically aggregate at the tumor-stromal interface and recruit B cells via the MIF-CD74 signaling axis, forming a localized immune barrier undetectable by traditional single-cell techniques.

Emerging evidence suggests that GNAL exhibits dual roles in tumorigenesis, functioning as both a tumor suppressor and a tumor promoter across different cancer types. Lian et al. demonstrated its tumor-suppressive function through multi-cell death pattern analysis, which identified an immunosuppressive subtype in low-grade glioma (LGG) and proposed GNAL as a key component of the CDPM prognostic score (29). Similarly, Zhang et al. highlighted GNAL as a critical tumor-suppressive gene in hepatocellular carcinoma (HCC), where it modulates proliferation, apoptosis, and neural signaling pathways (30). Conversely, Liu et al. reported that low GNAL expression serves as an independent prognostic marker in glioma, correlating with an immunosuppressive TME (e.g., reduced immune scores), aberrant DNA methylation, and poor clinical outcomes. Notably, GNAL-deficient tumors showed heightened sensitivity to anti-PD-1 therapy and were associated with 10 potential targeted drugs, underscoring their translational significance (31).

Spatial transcriptomic analysis reveals distinct hierarchical features in the GNAL-regulated immune microenvironment. MistyR algorithm-constructed three-level spatial interaction networks (intra/juxta/para-layer) show GNAL-positive stem cells engage in long-range signaling with vascular endothelial cells, with interaction strength positively correlating with microvascular density. This challenges the conventional view that endothelial cells regulate the TME only via local contact and suggests calcium signaling may mediate cross-regional regulation through exosomes or mechanosensing. CellChat analysis identifies that GNAL-positive stem cells secrete MIF ligand to activate the CD74-CXCR4 receptor complex on B cells, with 63% spatial overlap at tumor-normal tissue boundaries. The frequency of this interaction is 2.1-fold higher than other regions, establishing B cells as key downstream effectors of calcium signaling. COMMOT directional analysis confirms that signal propagation exhibits clear spatial directionality. In high-risk patients, B cell numbers decrease by 42%, and their continuous distribution at tumor margins is disrupted, increasing the distance CD8⁺ T cells must travel to the tumor core by 2.3-fold. This explains why traditional scores such as ESTIMATE fail to reflect B cells’ spatial protective roles.

From 102 genes differentially expressed at single-cell resolution, combined with patient survival data, we developed a prognostic model comprising CEBPD, KRT14, CD24, and KDM3B, which was validated in independent TCGA and GEO cohorts. Each gene signature in the four-gene prognostic model reveals tumor heterogeneity and invasiveness from distinct biological levels. As a therapeutic target, CEBPD expression is reduced in tumor tissues, and enhancing its activity can inhibit breast cancer progression. Additionally, the lipid metabolism-based prognostic model constructed using CEBPD has AUC values of 0.84, 0.83, and 0.80 for 3-, 5-, and 8-year survival predictions, respectively, which is superior to previous related studies (32). A groundbreaking 2025 study demonstrated that CD8⁺ T cell infiltration density in KRT14⁺ triple-negative breast cancer tumors is significantly higher than that in KRT14- tumors, with nearly a 3-fold increase in the response rate to PD-1 inhibitors. This finding highlights the association of KRT14 with the aggressive tumor phenotype and its potential as an immunotherapeutic target (33). In previous studies, as a glycosylphosphatidylinositol-anchored glycoprotein, high CD24 expression is significantly associated with shorter OS and decreased distant metastasis-free survival in breast cancer patients; it is also closely linked to larger tumor volume, higher histological grade, lymph node metastasis, and distant metastasis mediated by CD24⁺ circulating tumor cells (34,35). KDM3B (also called JMJD1B), a histone demethylase, promotes ovarian cancer cells’ resistance to PARP inhibitors (PARPi like olaparib) by synergizing with its homologous gene JMJD1C. Relevant research has also shown that overexpression of KDM3B in cancer cell lines leads to enhanced invasion and metastasis of breast cancer cells, providing new insights for the development of therapeutic strategies (36).

Similarly, Tian et al. constructed a CAVIN2 prognostic model in breast cancer (37). Our model integrates multidimensional features of calcium signaling (CEBPD), epithelial-mesenchymal transition (KRT14), and epigenetic regulation (KDM3B), while incorporating B cell spatial distribution within tumors. Clinical data analyses show that reduced GNAL expression induces stem cell differentiation arrest, and the high-risk phenotype facilitates immune evasion by weakening endothelial support networks and suppressing cooperative B/T cell antitumor responses. Wang et al. found the importance of B cell spatial distribution in immune therapy efficacy (38). Previous immune infiltration studies reported higher immune cell populations in GNAL low-expression groups, suggesting enhanced immunotherapy sensitivity, while low GNAL expression in gliomas correlates with attenuated antitumor immune responses, supporting and complementing our mechanistic insights. Structural analysis revealed tight binding between GNAL protein and CEBPD transcription factor (binding free energy ΔG =−25.1 kcal/mol), potentially influencing stem cell differentiation via epigenetic mechanisms such as histone demethylation (31). This supports observations in gliomas where GNAL low expression may be epigenetically regulated via DNA methylation.

Research on the molecular mechanisms of ER⁺ signaling pathway regulation of GNAL remains limited, but its significant correlation with KRT14—a potential biomarker in GNAL-based prognostic models—has been confirmed. According to Wu et al. (39), the estrogen receptor β (ERβ) signaling pathway maintains prostate basal cell homeostasis by regulating KRT14 expression. ERβ loss of function (via gene knockout or 5α-reductase inhibitor treatment) leads to aberrant elevation of AR coactivators, resulting in downregulation of the basal cell marker KRT14. This is accompanied by expansion of intermediate cells coexpressing P63 and AR, upregulation of NKX3.1, and localized immune cell infiltration—a phenotype partially reversible with testosterone treatment. In aromatase inhibitor-associated musculoskeletal syndrome (AIMSS) in breast cancer patients, abnormal activation of the ER⁺ signaling pathway may contribute to pathogenesis by regulating downstream inflammatory factors and bone metabolism imbalance. Differential expression of the epithelial cytoskeletal protein KRT14 (P=0.004) may exacerbate AIMSS by impairing muscle-tendon interface integrity or mediating pain-related cellular stress responses (40).

While our multi-omics approach provides novel insights into GNAL-mediated immune regulation, several technical and methodological limitations warrant consideration. First, the single-cell RNA sequencing analysis was constrained by a modest sample size (n=11 treatment-naïve ER⁺ breast cancer specimens), which may affect the statistical power to detect rare stem cell subpopulations and their transitional states. Although we mitigated potential batch effects using Harmony integration and corroborated key findings in bulk transcriptomic datasets (TCGA, GEO), future validation in larger, multicenter scRNA-seq cohorts will be essential to confirm the generalizability of GNAL’s role in stem cell dynamics. Second, while spatial transcriptomics enabled tissue-level mapping of cellular interactions, the inherent resolution limit (~55 µm spot diameter) may obscure precise discrimination between GNAL⁺ stem cells and adjacent fibroblast populations. To address this, we implemented RCTD deconvolution to achieve sub-spot resolution, with spatial patterns aligning with established stem-fibroblast niche distributions reported by Lian et al. (29).

More fundamentally, our study relies primarily on in silico predictions that require experimental validation, which require experimental validation through in vitro and in vivo assays. The single-cell dataset size was relatively small, and larger multi-center cohorts are needed to enhance statistical power and model generalizability. Our prognostic model, developed from a single institution, should be externally validated across broader clinical populations. Lastly, the downstream effects of GNAL interactions remain to be elucidated through mechanistic studies exploring their functional roles in ER⁺ breast cancer progression.

In this study, we identified GNAL as a novel tumor suppressor and prognostic biomarker in ER⁺ breast cancer through multi-omics integration and spatial transcriptome analysis. By combining GWAS data with bulk and single-cell RNA sequencing, we demonstrated that GNAL, the only calcium signaling pathway gene overlapping with survival-associated pathogenic genes, plays a key role in tumor suppression. Single-cell analysis revealed that GNAL⁺ stem cells are enriched at the tumor-stroma interface, where they may orchestrate the immunosuppressive microenvironment through spatial interactions with endothelial cells (CD99-CD99) and B cells (MIF-CD74-CXCR4/CD44). In addition, we constructed a robust 4-gene prognostic model (CD24, KDM3B, CEBPD, KRT14) based on the GNAL⁺ stem cell signature, which effectively stratified patients into high-risk and low-risk groups and was associated with immunosuppressive TME features. Molecular docking confirmed stable protein interactions between GNAL and all four prognostic genes, with the strongest binding to CEBPD (ΔiG =−25.1 kcal/mol), suggesting that GNAL may play a role as a regulatory hub. These findings position GNAL as a potential key mediator of calcium signaling, stem cell plasticity, and immune escape in ER⁺ breast cancer, providing new possibilities for risk stratification and stromal-targeted therapy.

Conclusions

This study uncovers GNAL’s critical role in ER⁺ breast cancer, linking impaired calcium signaling in stem cells to B cell-mediated immune exclusion via MIF-CD74. A GNAL-based prognostic model predicts immunotherapy response, while molecular docking confirms GNAL-CEBPD interactions. Our findings integrate spatial and functional dynamics, advancing precision oncology frameworks.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1796/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1796/prf

Funding: This work was supported by the National Natural Science Foundation of China (Nos. 82403020, 82303857 and 82304025).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1796/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Jahangiri R, Mosaffa F, Emami Razavi A, et al. PAX2 promoter methylation and AIB1 overexpression promote tamoxifen resistance in breast carcinoma patients. J Oncol Pharm Pract 2022;28:310-25. [Crossref] [PubMed]
2. Cao J, Zhou T, Wu T, et al. Targeting estrogen-regulated system x(c)(-) promotes ferroptosis and endocrine sensitivity of ER+ breast cancer. Cell Death Dis 2025;16:30. [Crossref] [PubMed]
3. Sui Y, Liu Z, Yao Y, et al. Estrogen receptor β inhibits breast cancer migration and promotes its apoptosis through NF-κB/IL-8 signaling. Transl Cancer Res 2025;14:1824-35. [Crossref] [PubMed]
4. Clusan L, Ferrière F, Flouriot G, et al. A Basic Review on Estrogen Receptor Signaling Pathways in Breast Cancer. Int J Mol Sci 2023;24:6834. [Crossref] [PubMed]
5. Saha T, Lukong KE. Decoding estrogen receptor and GPER biology: structural insights and therapeutic advances in ERα-positive breast cancer. Front Oncol 2025;15:1513225. [Crossref] [PubMed]
6. Zhang KM, Zhao DC, Li ZY, et al. Inactivated cGAS-STING Signaling Facilitates Endocrine Resistance by Forming a Positive Feedback Loop with AKT Kinase in ER+HER2- Breast Cancer. Adv Sci (Weinh) 2024;11:e2403592. [Crossref] [PubMed]
7. Griffiths JI, Cosgrove PA, Medina EF, et al. Cellular interactions within the immune microenvironment underpins resistance to cell cycle inhibition in breast cancers. Nat Commun 2025;16:2132. [Crossref] [PubMed]
8. Callari M, Dugo M, Barreca M, et al. Determinants of response and molecular dynamics in HER2+ER+ breast cancers from the NA-PHER2 trial receiving HER2-targeted and endocrine therapies. Nat Commun 2025;16:2195. [Crossref] [PubMed]
9. Xun Z, Ding X, Zhang Y, et al. Reconstruction of the tumor spatial microenvironment along the malignant-boundary-nonmalignant axis. Nat Commun 2023;14:933. [Crossref] [PubMed]
10. Yang W, Feng Z, Lai X, et al. Calcium nanoparticles target and activate T cells to enhance anti-tumor function. Nat Commun 2024;15:10095. [Crossref] [PubMed]
11. Campos J, Gleitze S, Hidalgo C, et al. IP(3)R-Mediated Calcium Release Promotes Ferroptotic Death in SH-SY5Y Neuroblastoma Cells. Antioxidants (Basel) 2024;13:196. [Crossref] [PubMed]
12. Sun P, Bush SJ, Wang S, et al. STMiner: Gene-centric spatial transcriptomics for deciphering tumor tissues. Cell Genom 2025;5:100771. [Crossref] [PubMed]
13. Liu J, Sun W, Li N, et al. Uncovering immune cell-associated genes in breast cancer: based on summary data-based Mendelian randomized analysis and colocalization study. Breast Cancer Res 2024;26:172. [Crossref] [PubMed]
14. Wang G, Shi C, He L, et al. Identification of the tumor metastasis-related tumor subgroups overexpressed NENF in triple-negative breast cancer by single-cell transcriptomics. Cancer Cell Int 2024;24:319. [Crossref] [PubMed]
15. Cardenas MA, Prokhnevska N, Sobierajska E, et al. Differentiation fate of a stem-like CD4 T cell controls immunity to cancer. Nature 2024;636:224-32. [Crossref] [PubMed]
16. Wu Q, Cao Y, Li Y, et al. Successful treatment of two patients with unresectable lung squamous cell carcinoma with tislelizumab regardless of programmed death-ligand 1 expression: a report of two cases. Anticancer Drugs 2022;33:e828-33. [Crossref] [PubMed]
17. Chen W, Zhao Y, Chen X, et al. A multicenter study benchmarking single-cell RNA sequencing technologies using reference samples. Nat Biotechnol 2021;39:1103-14. [Crossref] [PubMed]
18. Zhang F, Li X, Tian W. Unsupervised Inference of Developmental Directions for Single Cells Using VECTOR. Cell Rep 2020;32:108069. [Crossref] [PubMed]
19. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
20. Mounir M, Lucchetta M, Silva TC, et al. New functionalities in the TCGAbiolinks package for the study and integration of cancer data from GDC and GTEx. PLoS Comput Biol 2019;15:e1006701. [Crossref] [PubMed]
21. Silva TC, Colaprico A, Olsen C, et al. TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages. F1000Res 2016;5:1542. [Crossref] [PubMed]
22. Zhu Z, Zhang F, Hu H, et al. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets. Nat Genet 2016;48:481-7. [Crossref] [PubMed]
23. Qi T, Wu Y, Zeng J, et al. Identifying gene targets for brain-related traits using transcriptomic and methylomic data from blood. Nat Commun 2018;9:2282. [Crossref] [PubMed]
24. Wu Y, Zeng J, Zhang F, et al. Integrative analysis of omics summary data reveals putative mechanisms underlying complex traits. Nat Commun 2018;9:918. [Crossref] [PubMed]
25. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol 2019;20:296. [Crossref] [PubMed]
26. Salié H, Wischer L, D'Alessio A, et al. Spatial single-cell profiling and neighbourhood analysis reveal the determinants of immune architecture connected to checkpoint inhibitor therapy outcome in hepatocellular carcinoma. Gut 2025;74:451-66. [Crossref] [PubMed]
27. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med 2013;32:5381-97. [Crossref] [PubMed]
28. Chen H, Deng C, Gao J, et al. Integrative spatial analysis reveals tumor heterogeneity and immune colony niche related to clinical outcomes in small cell lung cancer. Cancer Cell 2025;43:519-536.e5. [Crossref] [PubMed]
29. Lian H, Wang J, Yan S, et al. An integrative analysis based on multiple cell death patterns identifies an immunosuppressive subtype and establishes a prognostic signature in lower-grade glioma. Ann Med 2024;56:2412831. [Crossref] [PubMed]
30. Zhang Y, Qiu Z, Wei L, et al. Integrated analysis of mutation data from various sources identifies key genes and signaling pathways in hepatocellular carcinoma. PLoS One 2014;9:e100854. [Crossref] [PubMed]
31. Liu Z, Yang L, Xie Z, et al. Multi-cohort comprehensive analysis unveiling the clinical value and therapeutic effect of GNAL in glioma. Oncol Res 2024;32:965-81. [Crossref] [PubMed]
32. Zhao Y, He H, Huang L, et al. Comprehensive analysis of lipid metabolic signatures identified CEBPD promotes breast cancer cell proliferation. Sci Rep 2025;15:6570. [Crossref] [PubMed]
33. Liao S, Zhang X, Chen L, et al. KRT14 is a promising prognostic biomarker of breast cancer related to immune infiltration. Mol Immunol 2025;180:55-73. [Crossref] [PubMed]
34. Li X, Tian W, Jiang Z, et al. Targeting CD24/Siglec-10 signal pathway for cancer immunotherapy: recent advances and future directions. Cancer Immunol Immunother 2024;73:31. [Crossref] [PubMed]
35. Cao J, Ma X, Zhang G, et al. Prognostic analyses of genes associated with anoikis in breast cancer. PeerJ 2023;11:e15475. [Crossref] [PubMed]
36. Hu A, Hong F, Li D, et al. KDM3B-ETF1 fusion gene downregulates LMO2 via the WNT/β-catenin signaling pathway, promoting metastasis of invasive ductal carcinoma. Cancer Gene Ther 2022;29:215-24. [Crossref] [PubMed]
37. Tian Y, Liu X, Hu J, et al. Integrated Bioinformatic Analysis of the Expression and Prognosis of Caveolae-Related Genes in Human Breast Cancer. Front Oncol 2021;11:703501. [Crossref] [PubMed]
38. Wang Q, Sun K, Liu R, et al. Single-cell transcriptome sequencing of B-cell heterogeneity and tertiary lymphoid structure predicts breast cancer prognosis and neoadjuvant therapy efficacy. Clin Transl Med 2023;13:e1346. [Crossref] [PubMed]
39. Wu WF, Song XY, Warner M, et al. The role of estrogen receptor β in maintaining basal cells and modulating the immune environment in the prostate. Proc Natl Acad Sci U S A 2025;122:e2505797122. [Crossref] [PubMed]
40. Jing F, Jiang L, Cao Y, et al. Plasma Proteomics and Metabolomics of Aromatase Inhibitors-Related Musculoskeletal Syndrome in Early Breast Cancer Patients. Metabolites 2025;15:153. [Crossref] [PubMed]