Calcium-related genes predict prognosis, immune characteristics, and response to immunotherapy in LUAD patients.

Introduction
Lung cancer (LC) is the most commonly diagnosed type of cancer and is characterized by its complexity and global impact [1]. According to the Global Cancer Statistics 2024 (GLOBOCAN 2024), LC currently ranks first in the world in terms of prevalence and mortality [2]. Lung adenocarcinoma (LUAD) is the main subtype of LC, accounting for about 35%−40% of LC cases [3, 4]. Despite early surgical interventions and late systemic whole-body treatments having ameliorated LUAD symptoms and prolonged overall survival (OS) of patients to some extent, the overall prognosis of LUAD remains less than satisfactory [5–7]. This challenge is not only due to the limitations of existing treatments but also intricately related to the complex mechanism of tumor pathogenesis [8]. Furthermore, the high heterogeneity of tumors further enhances the difficulty of disease evaluation and treatment [9]. Therefore, innovating reliable prognostic biomarkers is necessary to provide strong support for clinical research [5]. In recent years, various models have been created to predict the survival of patients with LUAD based on different signatures such as lipid metabolism [10], iron metabolism [11], ferroptosis [12], and calbindin S100 [13]. However, predictive models based on calcium-related genes (CRGs) have not been created.
As an essential metal element in the human body, calcium plays a key regulatory role in biological processes such as cell differentiation, aging, and signaling [14, 15]. With the deepening of the study of calcium-related signaling mechanisms in recent years, researchers have discovered that the dynamic changes of calcium ions are crucial for the occurrence and development of tumor diseases. These changes drive and influence a variety of key cellular processes, including cell proliferation, invasion, death, epithelial-mesenchymal transition, and the development of therapeutic resistance [16]. CRGs can not only effectively predict the prognosis of endometrial cancer, but also function as potential therapeutic targets [17]. Moreover, calcium plays an indispensable part in immune modulation and immunotherapy [18]. These findings suggested that calcium is crucial for tumor formation, tumor immune microenvironment (TIME) modulation, and disease prognosis. A thorough investigation into the prognostic value of CRGs holds great implications in exploring the mechanisms underlying tumors and developing corresponding therapeutic strategies.
In this work, we created a risk-scoring model based on CRGs for the first time to evaluate the prognosis of patients with LUAD. At the same time, we systematically analyzed the characteristics of the TIME and gene mutation landscape of different risk groups and predicted the efficacy of immunotherapy. Moreover, we identified several anticancer drugs with potential therapeutic value. Overall, this work provides essential theoretical foundations for the prognosis evaluation of LUAD patients, the analysis of the TIME, and the future development of precision medicine.

Materials and methods

Data acquisition and processing
RNA sequencing data, somatic mutation data, demographics, and clinical characteristics of LUAD patients were sourced from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database. 59 normal samples (adjacent non-cancerous tissues) and 541 LUAD tumor samples were collected, and these 600 samples were used as the training set. A validation set GSE31210 (n = 226) was collected from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database and CRGs were from GeneCards (https://www.genecards.org). A total of 537 CRGs were obtained by screening and downloading genes with a Relevance Score greater than or equal to 8. Calcium was the keyword for screening (Table S1).

Identification and functional analysis of differentially expressed CRGs (DE-CRGs) in LUAD
The R package edgeR [19] was applied in the differential analysis of normal and tumor groups of LUAD (|logFC|>1, FDR < 0.05). The differentially expressed genes (DEGs) and CRGs were intersected to yield DE-CRGs. Subsequently, GO enrichment analysis was performed on these CRGs to probe into their functions.

Building and testing of a prognostic model
In the training set, the survival package [20] was utilized to develop the univariate Cox regression analysis on DE-CRGs to identify survival-associated genes in LUAD. The p-value of 0.01 was the threshold for significance. Subsequently, LASSO Cox regression analysis was performed to further narrow down the range of CRGs associated with survival. Finally, multivariate Cox regression analysis was utilized to screen characteristic genes and construct a calcium-related risk-scoring model. The calculation formula for the risk-scoring model is as follows:
Model=∑Coefficientgene*Expressionvalue (gene).
Patients in the training and validation sets were subdivided into high-risk (HR) and low-risk (LR) groups with the median riskscore as the cut-off point. Subsequently, the survival information of patients in each group was recorded, and the K-M survival curve was plotted. Additionally, a time-dependent receiver operating characteristic (ROC) curve was established to verify the effectiveness of the risk-scoring model in the training set and the validation set. Finally, we used an expression heatmap of the model gene in the training set to preliminarily investigate the expression differences between different risk groups.

Functional enrichment analysis
The edgeR package was employed to screen the DEGs in the HR and LR groups, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was carried out with the software gene set enrichment analysis (GSEA).

Clinical independence assessment of the nomogram
Combined with riskscores and clinicopathological information (including clinical stage, tumor grade, and lymph node status), the differences in riskscore among patients with different clinicopathological features were assessed using the Wilcoxon test. Moreover, clinical indicators and riskscores were included in univariate Cox regression analysis and multivariate Cox regression analysis to test the independence of the risk-scoring model. The rms package [21] was used to generate a nomogram to predict the 1-, 3-, and 5-year survival of LUAD patients. The corresponding calibration curves were graphed to verify the prediction effect of the nomogram.

Assessment of TIME and prediction of immunotherapy response
Based on the sequencing data in the training set, the stromal score, immune score, and ESTIMATE score were assessed using the ESTIMATE algorithm [22] to evaluate the expression of immune cells and stromal cells in the TIME. Then, the enrichment of immune cells and immune function in the HR group and LR group was analyzed by using the single sample GSEA (ssGSEA) algorithm and the CIBERSORT algorithm [23], respectively. At the same time, the expression of immune checkpoints and HLA genes in the two different groups was counted. The immunophenoscore (IPS) and tumor immune dysfunction and exclusion (TIDE) scores of LUAD were downloaded from the Atlas of Cancer Immunome (TCIA, https://tcia.at) online platform and the TIDE website (http://tide.dfci.harvard.edu/), respectively. The differences between the IPS and the TIDE scores were compared in combination with the information from the HR and LR groups of the sample.

Somatic mutation analysis
Mutation data for LUAD derived from the TCGA database. Moreover, the maftools package [24] was applied in the mutation analysis of all samples and the comparison of tumor mutational burden (TMB) between the two risk groups. Volcano maps were graphed to visualize the mutations of the top 30 mutated genes in the HR and LR groups and the genes in the risk-scoring model.

Drug sensitivity analysis
The CellMiner database (https://discover.nci.nih.gov/cellminer/) was applied in the screening of antitumor drugs whose sensitivity was greatly correlated with prognostic genes and whose IC50 was significantly different in the HR and LR groups.

Results

Identification of DE-CRGs in LUAD
To screen for DE-CRGs in LUAD, we compared the mRNA expression profiles of normal and tumor samples in the training set. A total of 5461 DEGs were obtained, of which 196 genes were CRGs (including 116 up-regulated CRGs and 80 down-regulated CRGs) (Fig. 1A). In GO enrichment analysis, these CRGs were mainly in transmembrane transporter activity (calcium ion transmembrane transporter activity, metal ion transmembrane transporter activity) and ion channel activity (ion channel regulator activity, calcium channel activity, calcium channel regulator activity) (Fig. 1B).

Prognostic model development and validation of CRGs
To probe into the prognostic value of CRGs, we developed a relevant risk-scoring model in the training set and validated it in the validation set. First, in the training set, using univariate Cox regression analysis, we identified 22 CRGs associated with LUAD survival (Table S2). To reduce the complexity of the model, we launched a LASSO Cox regression analysis (Fig. 2A-B). Finally, we identified eight CRGs as prognostic biomarkers for LUAD through multivariate Cox regression analysis and created corresponding risk-scoring models. The formula is provided here:

To validate the reliability and predictive power of the model, riskscores were calculated in both cohorts, and samples were clustered into HR and LR groups based on the median riskscore (Fig. 3A-B). The K-M curve manifested that the survival outcome was better in the LR group (Fig. 3C-F). The ROC curve displayed that the model had excellent sensitivity and specificity for predicting OS in LUAD patients, with an AUC value greater than 0.69 in both cohorts (Fig. 3G-H). The heat map displayed that GRIA1, BTK, and CLCA1 were highly expressed in the LR group, and PDGFB, S100P, TRPA1, F2RL1, and FBN2 were highly expressed in the HR group (Fig. 3I).

GSEA results
Using GSEA software, we conducted the KEGG pathway enrichment analysis in the HR group. The results manifested that the HR group was greatly enriched in Adherens junction, Cell cycle, Axon guidance, and P53 signaling pathway, as well as pathways related to Escherichia coli infection and various types of cancer (Fig. 4).

Assessment of the clinical independence of riskscore and formulation and validation of the nomogram
Multivariate and univariate Cox regression analyses were processed by combining risk-scoring models and other clinical factors. The model functioned as an independent prognostic indicator (p < 0.001) (Fig. 5A-B). Violin plot analysis revealed that the riskscores of patients with LUAD with different clinical features were highly consistent with disease severity, with higher riskscores associated with advanced stages (stage III-IV), tumor enlargement or invasion (T3-T4), and more lymph node metastases (N1-N3) (Fig. 5C). Based on riskscores and clinical characteristics, we graphed a comprehensive nomogram. The calibration curves displayed an excellent consistency of the survival rates predicted by the nomogram with actual survival rates at one, three, and five years (Fig. 5D-G), further validating the effectiveness and reliability of both the model and the nomogram.

Immunological characteristics of the risk-scoring model and evaluation of immunotherapy effect
The comparison of immune-related scores in the training set uncovered that the LR group had higher immune scores and ESTIMATE scores, implying that their immune cells and functions were more abundant (Fig. 6A). ssGSEA analysis further demonstrated significant enhancement of multiple immune cells and functions (e.g., Mast cells, B cells, and Tfh) in the LR group (Fig. 6B). In CIBERSORT analysis, the LR group was found to be significantly enriched with immune cells such as B cells naive and T cells CD4 memory resting, while the HR group was enriched with Tregs and Macrophages M0 (Fig. 6C). Comparisons between immune checkpoints and HLA genes revealed that the majority of these checkpoints (such as CD28, CD80, CTLA4, TNFSF15, etc.), along with all HLA genes, exhibited high expression in the LR group, whereas some checkpoint genes (such as CD200, CD276, VTCN1, etc.) displayed higher expression in the HR group (Fig. 6D-E). Moreover, the LR group had lower TIDE scores and higher IPS scores. They had greater sensitivity to immunotherapy and a higher likelihood of benefiting from treatment (Fig. 6F-G).

Somatic mutation analysis
By visualizing the mutational data in a training set comprising 483 tumor samples, we discovered that a considerable number of missense mutations existed in LUAD samples, resulting in the conversion of corresponding codons into different amino acid codes. Furthermore, the higher frequency of nucleotide base exchanges (such as cytosine (C) to adenine (A)), contributed to the increased number of SNPs. Furthermore, the waterfall plot revealed that mutations in the genes TTN and MUC16 occurred in over 40% of the tumor samples (Fig. 7A). Furthermore, the analysis of the mutation landscape of the two groups revealed that they had different TMBs and that the mutation landscape of the HR group was more complex than that of the low-level group (Fig. 7B). Correspondingly, the high-mutation genes TP53 and TTN exhibited mutation rates of 55% and 48% in the HR group, respectively, and 42% and 40% in the LR group, respectively (Fig. 7C-D). Moreover, a total of 165 samples (86.67%) had mutations in these eight genes in the tumor samples, among which the mutations in the FBN2 gene were the most common in the samples, with a mutation frequency of 50%, followed by TRPA1 (29%), BTK (9%), CLCA1 (6%), GRIA1 (5%), F2RL1 (2%), PDGFB (1%), and S100P (1%) (Fig. 7E). These mutations may disrupt DNA or RNA’s normal structure and function in a particular way. Illuminating these mechanisms is crucial for the development of precision medicine and targeted therapies.

Drug sensitivity analysis
We evaluated the potential correlation between model genes and IC50 of drugs based on the CellMiner database. The expression of PDGFB was significantly positively correlated with the IC50 of Sapitinib, Erlotinib, and Afatinib. S100P was significantly positively correlated with the IC50 of Barasertib (Fig. 8A). In the PDGFB high-expression group, the IC50 of Sapitinib, Erlotinib, and Afatinib was significantly higher (Fig. 8B). Taken together, the expression levels of PDGFB and S100P may be key regulators of tumor cell sensitivity to these drugs. High expression of PDGFB and S100P is significantly associated with cell resistance to these drugs.

Discussion
LUAD is one of the subtypes with the highest mortality rate in LC and poses a daunting threat to human health [25]. Despite tremendous advancements in the diagnosis and treatment of (LUAD) in recent years, biological markers for early detection are still limited to fully meet clinical demands. Therefore, innovating novel and effective biomarkers is essential to improve the prognosis and treatment of LUAD. Calcium, as an essential metal element, plays a critical role in cell proliferation, differentiation, and signal transduction in various biological processes. Serum calcium levels serve as a risk factor for bone metastasis in bladder cancer (BCa) patients and can be combined with hemoglobin levels to predict the likelihood of bone metastases in BCa patients [26]. Furthermore, in a prostate cancer study, calcium levels are found to correlate with adverse prognosis [27]. Though the current literature on calcium’s role in the prognosis of LUAD is scarce, evidence suggests that serum calcium-binding protein S100A10 is closely correlated with the disease stage and lymph node metastasis of LUAD. It activates the Akt-mTOR signaling pathway to reinforce tumor cell proliferation and invasion, indicating its potential as a tumor biomarker [28, 29]. Based on this, the present study systematically dissected the potential value of CRGs in the prognosis of LUAD using bioinformatics methods, thereby designating a new perspective for further understanding its biological mechanisms.
In this work, we identified 8 CRGs (GRIA1, BTK, CLCA1, PDGFB, S100P, TRPA1, F2RL1, FBN2) that were involved in the prognosis of LUAD through regression analysis, and successfully created a risk-scoring prognostic model. The ROC curve verified the reliability of the model. Additionally, the analysis demonstrated that the independent prognostic value of the model was significant and highly correlated with the clinical characteristics of LUAD patients. LUAD patients with HR scores exhibited a preference for disease status in an advanced stage, higher tumor burden, and broader lymph node metastasis, further suggesting a correlation between HR scores and worse prognosis. Furthermore, we observed that GRIA1, BTK, and CLCA1 were protective genes identified in the model, with elevated expression in the LR group. Conversely, PDGFB, S100P, TRPA1, F2RL1, and FBN2 were risk genes identified in the model, with heightened expression in the HR group. Notably, this work uncovered that all eight model genes exhibited varying degrees of gene mutations. Therefore, it is speculated that the abnormal expression or variation of these genes may drive diverse tumor characteristics, ultimately affecting disease progression. GRIA1 is a subunit encoding AMPA-type glutamate receptors, whose mutation can lead to neurological developmental disorders [30]. Despite the limited specific studies on GRIA1’s role in LUAD currently available, numerous research findings have suggested that GRIA1 can serve as an important prognostic biomarker for LUAD, which is consistent with the conclusion of this study [31–33]. As a member of the Tec kinase family, BTK is implicated in the regulation of multiple immune pathways [34]. At present, BTK inhibitors have made tremendous progress in the treatment of B-cell malignancies [35]. In non-small cell lung cancer (NSCLC), overexpression of BTK has been found to curb tumor cell proliferation, invasion, and migration [36]. Furthermore, in LUAD, BTK is considered a critical prognostic factor and is closely related to the reprogramming of the TIME [37]. These findings suggest a potential role of BTK in immune regulation and tumor progression in LC. CLCA1 is a member of the calcium-activated chloride channel regulator family, regulates epithelial cell chlorination currents, and is implicated in the pathogenesis of respiratory and gastrointestinal diseases associated with mucus hypersecretion, including chronic obstructive pulmonary disease, asthma, pneumonia, etc [38]. In hepatocellular carcinoma, CLCA1 is thought to have a tumor-suppressive effect [39]. However, its function in LUAD is unclear and requires further research. PDGFB has been revealed to play a key part in tumor angiogenesis [40]. Research in LUAD has uncovered that PDGFB expression is implicated in fibrosis and epithelial-mesenchymal transition processes, which exert influence on the invasion and metastasis of LUAD [41]. S100P is a member of the calcium-binding protein family S100, associated with malignant phenotypes [42]. S100P is upregulated in LUAD and plays a role both intracellularly and extracellularly, where it can interact with the advanced glycation end product receptor (RAGE) to facilitate tumor development [43]. Recent evidence supports that the ion channel TRPA1 is associated with LUAD, which reinforces tumor progression and metastasis through direct interaction with FGFR2 [44]. F2RL1 is an F2R-like trypsin receptor, and its overexpression in LUAD tissues and cells not only facilitates cancer cell proliferation and stem cell spheroid formation but also boosts the phosphorylation of epidermal growth factor receptor (EGFR) [45] by upregulating VEGFA expression. FBN2 is an extracellular matrix glycoprotein that exhibits heightened expression in LC, with great correlations to TNM stage and lymph node metastasis status, serving as a marker of adverse patient prognosis [46, 47]. Currently, although the function of FBN2 in LUAD remains unclear, research suggests that it can be a key prognostic biomarker for LUAD [48]. Collectively, the abnormal expression or variation of these model genes identified in this work may drive the progression and metastasis of LUAD. Moreover, KEGG pathway enrichment analysis demonstrated that the HR group was considerably enriched in signaling pathways such as Adherens junction and P53 signaling pathway, which may synergistically drive tumor progression and metastasis. The enrichment of adhesion pathways may reflect changes in cell-cell adhesion, which is a key initial step in the distant metastasis of tumors [49]. Abnormalities in the P53 signaling pathway may further boost tumor development and progression by disrupting genomic stability [50]. Comprehensive analysis manifested that CRGs play an instrumental part in the occurrence and development of LUAD, and may impact the prognosis of patients by modulating specific signaling pathways, such as adhesion junction and P53 signaling pathway. These findings shed fresh insights into the relationship between CRGs and tumor signaling networks for further research.
Calcium signaling, which is at the core of numerous cellular processes that can lead to cancer treatment-induced immune responses, cancer growth, and apoptosis, has been integrated into current approaches to cancer immunotherapy [51]. However, the correlation between calcium and immune cell infiltration in LUAD remains unknown. Therefore, we herein integrated ESTIMATE, ssGSEA, and CIBERSORT algorithms to systematically analyze the immune-related features of different risk groups, revealing great disparities between the LR group and the HR group in terms of their TIME and responsiveness to immunotherapies. The LR group exhibited enhancement in multiple immune cells and immune function (e.g., Mast cells, B cells, and Tfh), while the HR group had an increase in specific immunosuppressive cells (e.g., Tregs). In previous research, mast cell abundance is related to a good prognosis in patients with LUAD [52]. The cooperation of B cells and Tfh can reinforce anti-tumor CD8 T cell responses, which in turn facilitate anti-tumor immunity [53]. The enrichment of Tregs reinforces the formation of an immunosuppressive microenvironment, which in turn drives the disease progression of LUAD [54]. Furthermore, we discovered that in the LR group, a substantial elevation in the expression of immune checkpoint genes was prevalent, which not only corresponds to their complex TIME but may also reflect their higher potential for anti-tumor immunity. In particular, elevated expression of immune checkpoint genes (e.g., CD28, CD80, etc.) suggested that LR groups may be more likely to benefit from immune checkpoint blockade therapy. The presence of CD28 in LUAD is impactful for the functionality of PD1+ CD8+ T cells, manifested by the multifunctionality of PD1 + CD28 + T cells at tumor sites and the lower expression of inhibitory molecules (such as TIGIT, TIM-3, and LAG-3), all of which can maintain the antitumor activity of T cells [55]. Similarly, the upregulation of CD80 may interact with CTLA4, which in turn facilitates tumor progression in LUAD [56]. Further analysis of IPS and TIDE scores demonstrated that the LR group had higher IPS scores and lower TIDE scores, indicating that the LR group was considerably more sensitive to immunotherapy than the HR group and had greater clinical treatment potential. Taken together, we hypothesized that LR patients may achieve more obvious efficacy when receiving immunotherapy, while HR patients may require a combination therapy strategy to overcome their immunosuppressive status.
We also predicted some drug candidates, such as Erlotinib, Afatinib, and Barasertib. Erlotinib is an inhibitor targeting EGFR and has been approved in many countries for the second and third-line treatment of advanced or locally metastatic NSCLC [57]. Afatinib is an irreversible ErbB family inhibitor that is effective in patients with NSCLC with uncommon mutations [58]. Barasertib (AZD1152) has been validated to curb the growth of small-cell lung cancer cell lines in vivo and in vitro [59].
In general, this work systematically investigated the potential value of CRGs in the prognosis of LUAD and constructed a riskscore prognostic model based on eight CRGs. We discovered that the model score was closely associated with the clinical features of LUAD, with HR patients more likely to have advanced-stage tumors, heavier tumor burdens, and lymph node metastasis, suggesting that the model can effectively predict prognosis. The analysis of TIME suggested that patients in the LR group had a stronger immune response, for whom immunotherapy may be a better option. Furthermore, we predicted the clinical potential of some of our drug candidates. In conclusion, this work brings insights into the role of the TIME of LUAD and CRGs in tumor prognosis and treatment. Nevertheless, this study relied primarily on bioinformatic analyses but lacked experimental validation, limiting its widespread applicability. Further research is necessary to illuminate the relationship between CRGs and the TIME, particularly the mechanisms by which calcium signals modulate immune cell infiltration. In the future, further experimental studies and clinical validations will aid in confirming the specific roles of these CRGs in LUAD, particularly their potential in immune evasion and tumor immunotherapies. Furthermore, the study of calcium signaling and immunotherapy is expected to offer more accurate guidance for the individualized treatment of LUAD patients.

Supplementary Information