Introduction
Gastric cancer (GC) is a common malignancy located in the gastrointestinal tract with more than one million confirmed cases worldwide each year, thus making it the fifth most commonly diagnosed cancer and the seventh most prevalent cancer in the world.1 In recent years, the incidence of GC in people aged 40 years or younger has been on the rise in some countries or regions.2 To date, the exact mechanism of GC occurrence and development is still unclear.
Immunogenic cell death (ICD) is a form of regulated cell death (RCD) that is sufficient to activate adaptive immune responses in an immunologically active environment. ICD is a type of tumor cell death induced by certain chemotherapeutic drugs, lytic viruses, physical chemotherapy, photodynamic therapy, and radiation therapy.3 Through the function of immunogenic death mediating the body’s immune response, we could possibly analyze the relationship of ICD-related genes in GC and provide partial theoretical support for treating this form of cancer.
As a class of biomolecules with regulatory roles in many life activities, such as cell cycle and cell differentiation, lncRNA plays an important role in transcriptional regulation, post-transcriptional processing, and chromatin epigenetic control in both normal human operation and pathological conditions.4 At this stage, there are many studies on lncRNAs, but the link between ICD-related lncRNAs and GC is less studied. Therefore, we aimed to study the role and prognostic value of lncRNA in ICD in GC to fill the gap in this area. This study constructed prognostic signals for GC by screening lncRNAs associated with ICD, analyzed the associated immune escape and immunotherapy sensitivity, and screened potential drugs for the disease.
Materials and methods
Data acquisition and collation
We extracted RNA-seq transcriptomic data and corresponding clinical data from TCGA for GC patients, where all sample types belonged to adenocarcinoma (according to the World Health Organization’s staging model for GC). Through a literature search, we collected 20 genes associated with ICD to be added to this study. The tumor immune dysfunction and exclusion (TIDE) scores associated with the current selection of cases from the TIDE database were downloaded (Supplementary File 1).
Construction of immunogenic cell death-associated lncRNA signaling
We assessed the lncRNA associated with survival related to ICD by univariate Cox regression analysis. GC samples were randomly divided into training and test sets in a 1:1 ratio, and then the prognostic characteristics of the GC patients were constructed based on LASSO Cox regression (Table 1). The formula used for the risk scoring was:
Risk score=∑expression (lncRNA)×coef (lncRNA)
Table 1The base data for model formula building
| id | Coef# |
|---|
| AC005332.1 | −0.618092349510564 |
| AC116312.1 | 0.921584780441913 |
| LINC00705 | 0.453026000009769 |
| CEP250-AS1 | −0.816745292649202 |
| AC234775.2 | 0.607496025313933 |
| LINC01150 | 0.83819920086105 |
| FLJ16779 | 0.655190927969806 |
| UBL7-AS1 | −0.606773322017963 |
| AC010457.1 | 0.586370556932745 |
Clinicopathological staging and the construction of nomograms
Univariate and multifactorial Cox regression analyses were performed to investigate whether patient-related clinicopathological characteristics and risk scores were associated with the overall survival. The age, gender, grading, stage, TNM classification, and risk score of the GC patients were analyzed by a nomogram, and the total score was calculated to predict the probability of survival at one, three, and five years, respectively. The performance of the nomogram was tested using calibration curves (Supplementary File 1; clinically relevant data).
Risk difference analysis and enrichment analysis
Differential analysis was used to distinguish differentially expressed genes in the low-risk and high-risk groups. We analyzed the functions and channels of differentially expressed gene enrichment by using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.
Mutation frequency of related genes and tumor mutational burden
We collated the mutation data of the patients’ genes and then used the maftools package to display the mutation frequencies of the genes using waterfall plots. An optimal cut-off point was calculated using the survminer package to distinguish between the low and high tumor mutation load groups, and finally survival analysis was performed (Supplementary File 1; tumor mutational burden (TMB) data).
Immune function and prediction of immunotherapy sensitivity
By using the limma package,5 the differences in various immune-related functions between the low-risk and high-risk groups were revealed. The sensitivity between the low-risk and high-risk groups when receiving immunotherapy was also investigated (Supplementary File 1; TIDE analysis of the source data).
Screening for potential drugs to treat diseases
Screening for sensitivity to different drugs between the low-risk and high-risk groups was applied by the pRRophetic R package,6 which shed some light on the treatment of the disease.
Statistical analysis
All statistical analyses including Kaplan-Meier survival analysis and univariate and multivariate Cox analyses were performed based on R language version 4.1.2 (Now developed by the R Core Team in New Zealand). The p-value <0.05 was considered statistically significant in the different comparisons. The presence of “*”, “**”, and “***” in all images was defined as p-value <0.05, p-value <0.01, and p-value <0.001, respectively.
Results
Predicting the co-expression of lncRNA and ICD in GC patients
In this experiment, the RNA-seq transcriptome and related clinical data downloaded from TCGA database (a total of 407 stomach tissues, including 375 tumor tissues and 32 normal tissues) were used for the prediction of the gene co-expression with 20 ICD-related genes obtained through the literature search, and 18 genes (ADAR, ANXA1, APC, AXL, BAX, BLM, BTK, CALR, CASP8, CD47, COX10, CP, CXCL10, EPHB4, HMGB1, IDO1, STAT3, and TLR4) with a strong correlation were also obtained (Fig. 1).
Development and validation of an ICD-related lncRNA prognostic model
The GC tissue samples were assigned 1:1 as the training and test groups, and the prognostic models were constructed by Cox regression and Lasso regression. By judging the risk level of lncRNA by their median risk score, the differences in expression could be derived for 30 lncRNA (p < 0.05), of which 10 lncRNA were low-risks (HR < 1) and 20 lncRNA were high-risks (HR > 1) (Fig. 2a). A correlation heat map determined the correlation between the ICD-related genes and lncRNA involved in the construction (Fig. 2b).
Through survival analysis and visualization, we obtained the overall survival curves as well as survival curves for the training and test groups (Fig. 2c–e). The analysis of the progression-free survival (PFS) predicted significant differences in the PFS between patients in the low-risk and high-risk groups (p < 0.001), and the PFS in the low-risk group was more optimistic (Fig. 2f). The survival distribution and risk profile showed that more patients died with increasing risk scores in both groups (Fig. 3a–d). The low-risk and high-risk profiles of lncRNA in both groups is shown in Figure 3e and f.
Independent prognostic analysis and survival-related nomograms
The analysis of the various clinical trait factors predicted that age, stage, and the model we constructed (risk score) was significant (p < 0.05) in both the univariate prognostic and the multivariate prognostic analyses (Fig. 4a, b). By plotting the nomograms, we could predict the survival of the patients (Fig. 4c), and the nomograms were highly reliable, as verified by the calibration curve (Fig. 4d).
We validated the model by clinical grouping and could conclude that our model had some applicability in different clinical groupings (p < 0.05) (Fig. 5a, b). Furthermore, the consistency index clearly showed that the predictions of our constructed model had high consistency compared to the actual results (Fig. 5c). The accuracy of the experimentally constructed model was tested by plotting the ROC curves (the area under the ROC curve for one, three, and five years was 0.689, 0.730, and 0.707, respectively) (Fig. 5d). In addition, the model we constructed predicted the survival of patients with superior predictive performance compared to other clinical traits (Fig. 5e).
KEGG and GO analysis of differential genes
We analyzed the three aspects of the molecular functions (MF), biological processes (BP), and cellular components (CC) by GO analysis. We could see that the differential genes were enriched in the muscle system process (GO:0003012) for the BP, mainly in collagen containing (GO:0062023) for the CC and in the MF for actin binding (GO:0003779) (Fig. 6a, c, e). Moreover, we found that these differential genes were significantly enriched in the neuroactive ligand-receptor interaction (ID: hsa04080) pathway by the KEGG analysis (Fig. 6b, d, f).
Mutation frequency of related genes and tumor mutational burden
By using the maftools software package, a waterfall plot was used to visualize the mutation frequency of the gene, and the study found that the frequency of mutations was generally higher in the low-risk group, while the mutation frequency of the gene titin was higher in both groups (Fig. 7a, b). By analyzing the survival of the TMB, we found that the survival analysis was different between the low TMB group and the high TMB group (p < 0.05), and the high mutation load group was more superior, which could be related to the fact that cells with a high mutation load were more likely to be attacked by the immune system. A joint survival analysis of the TMB and risk scores also allowed us to conclude that they were significantly different (Fig. 7c, d).
Differences in immune function and sensitivity to immunotherapy
By exploring the differences between the two groups (p < 0.05), we obtained five immune-related functions that were different between the low-risk and high-risk groups; namely, MHC_class_I, Type_II_IFN_Response, CCR, APC_co_stimulation, and Parainflammation (p < 0.05) (Fig. 8a). The sensitivity analysis of immunotherapy showed that there was a difference between the low-risk and high-risk groups with the high-risk group being more sensitive to immunotherapy (Fig. 8b). The differential analysis of the immune checkpoints showed that the expression of the immune checkpoints differed between the low-risk and high-risk groups (Fig. 8c).
Drug sensitivity analysis and potential drug screening
By screening the drugs, we derived differences in sensitivity to different drugs between the high-risk and low-risk groups, thus expecting to discover potential drugs to treat the disease. For a demonstration, we selected three drugs: TGX221, XMD8-85, and CGP-60474 (Fig. 9; panels g to l are from the PubChem website). The sensitivity of all three drugs differed between the high-risk and low-risk groups (p < 0.05), and the inhibitory concentration 50 (IC50) was lower in the high-risk group, thereby indicating that the high-risk group was more sensitive to them.
Discussion
The current study identified nine lncRNA (AC005332.1, AC116312.1, LINC00705, CEP250-AS1, AC234775.2, LINC01150, FLJ16779, UBL7-AS1, and AC010457.1) to model the prognostic features of GC. In 2016, Liu et al. reported LINC00705 as one of the features to predict the recurrence of human breast cancer.7 Likewise, LINC01150 was reported by He et al.8 as one of the possible independent predictors of lung adenocarcinoma prognosis. Sun et al.9 reported that LINC01150 could be one of the lncRNA with an independent prognostic value in GC. Wang et al.10 reported that FLJ16779 was one of three lncRNAs associated with the prognosis of GC. Cao et al.11 reported that UBL7-AS1 was associated with cell cycle gene expression in a study based on glioblastoma samples. UBL7-AS1 also had prognostic value in squamous cell carcinoma of the cervix.12 In addition, UBL7-AS1 and AC010457.1 were associated with the prognosis of patients with GC.13,14 The drug TGX221 was tested in the treatment of prostate cancer in nude mice and the treatment of human renal cell carcinoma.15,16 Moreover, CGP-60474 is a potential drug for curing many cancers and could have a role in treating COVID-19 lung injury.17,18
As a malignant tumor with high incidence, GC could be caused by factors, such as H. pylori infection, Epstein-Barr virus infection, smoking, a high salt diet, and self-produced cytokines.19 Since GC usually has no apparent symptoms in the early stages, it would not be easy for patients to pay attention to it and be diagnosed. Patients diagnosed with GC are mostly advanced patients, which would make it more challenging to treat GC.20,21 Nowadays, the first choice of clinicians for GC treatment is still surgery22 through open laparoscopic surgery and robotic surgery.23 Chemotherapy and radiation therapy are also increasingly chosen, but these methods have disadvantages, such as technical limitations and more adverse effects.24 Additionally, molecular targeted therapy and perioperative treatment are not ideal for some patients with GC.25 Therefore immunotherapy is now desired for the treatment of cancer, and various modalities of immunotherapy are now being used in the clinic.26
Many types of regulated cell death, such as ferroptosis, necroptosis, and pyroptosis, have now been studied by many scholars for their role and prognosis in cancer. For example, Tang et al.27 investigated the role of ferroptosis in head and neck squamous cell carcinoma; Yang et al.28 analyzed the prognosis of necroptosis in hepatocellular carcinoma; Shao et al.29 investigated the role of prognostic analysis of pyroptosis in GC. However, reports on the prognosis of ICD are still scarce. ICD is a type of RCD that activates dead cells to express endogenous or exogenous antigens to subject them to an adaptive immune response.30 This has also been used as an immunotherapeutic approach to treat cancer in recent years through the efforts of numerous researchers and clinicians. Liu et al.31 treated melanoma with a nanomaterial-mediated immunotherapy approach. However, studies on the prognosis of ICD-related lncRNA in GC are still lacking.
In this study, we predicted ICD-related lncRNA by co-expression and distinguished the training and test groups by performing one random grouping of the results. The lncRNA in the training group was filtered for significant prognosis-related lncRNA using the univariate Cox regression method, and then the Cox model was constructed by Lasso regression analysis to derive a risk score formula for the training group. The test group was scored by the risk scoring formula obtained from the training group. We then used univariate and multifactorial independent prognostic analyses to conclude that the risk score could be used as an independent prognostic indicator (p < 0.001). There was a significant difference in survival between the high-risk and low-risk groups in the training group and the test group. The C-index and ROC curves showed that the model we constructed was reliable and highly accurate. The positions of differential gene enrichment were predicted by the GO and KEGG enrichment analysis. Simultaneously, the differential analysis of the immune checkpoints and the differential analysis of the immune function were used in the hope that this would shed light on the immunotherapy of GC. However, this study had some limitations. In particular, our analysis only included external data from TCGA database, thus lacking our own data. Additionally, we did not conduct the basic biological experiments to validate the results.
Conclusions
This study obtained prognostic and immunoreactive results of genes related to ICD in GC and new prognostic-related lncRNAs. The immune-related predictions and model could help predict the outcome of GC and could provide a reference for clinical practice.
Supporting information
Supplementary File 1
Data used in current cases.
(XLSX)
Abbreviations
- GC:
gastric cancer
- GO:
gene ontology
- ICD:
immunogenic cell death
- lncRNA:
long non-coding RNA
- KEGG:
Kyoto Encyclopedia of Genes and Genomes
- RCD:
regulatory cell death
- ROC curve:
receiver operating characteristic (ROC) curve
Declarations
Ethical Approval
Our data were sourced from open-source databases and did not require ethical approval.
Data sharing statement
All data used in the study are available in TIDE (http://tide.dfci.harvard.edu/) and TCGA (https://portal.gdc.cancer.gov/). The clinically relevant data, TMN data, gene expression data for gastric cancer, and TIDE analysis of source data can be found in the Supplementary File 1. The clinically relevant data were basic information of the cases selected for this study, including age, gender, and stage; The TMN data were values of the tumor mutation load of the selected cases. The gene expression data were the information of the gene expression detected in all samples; The TIDE analysis of the source data included the TIDE scores of the cases selected for this study.
Funding
This study was not supported by any funding agency from the public, commercial, or nonprofit sectors.
Conflict of interest
The authors declare that there was no conflict of interest in this study.
Authors’ contributions
YH, JQN, and LZS conceptualized this study, JQN and LZS performed the data compilation and data analysis. YH and LZS performed the literature review. JQN and LZS drafted the manuscript. All authors made a significant contribution to this study and approved the final manuscript.