Long non-coding RNA signature for predicting gastric cancer survival based on genomic instability

Data download and processing

We downloaded somatic mutation data (MuTect2 Variant mutation and Masking data) of 437 gastric cancer samples, transcriptome data of 407 samples (including STAD-Counts and STAD-FPKM), and clinical data of 348 gastric cancer patients from The Cancer Genome Atlas (TCGA) database [22]. We also obtained data of 300 gastric cancer patients with clinical data and transcriptional information from the GEO database (data set GSE62254) [23]. For subsequent analysis, we used the COUNT data type for differential analysis, and the FPKM data type was converted to the TPM (transcripts per kilobase of exon model per million mapped reads) data type [24]. All expression data were transformed by log2 for further analysis [25]. After excluding patients without clinical data or with 0 survival time, we obtained a final cohort of 348 gastric cancer patients with both clinical information and expression data. We randomly divided these patients into a training group and a test group, as well as a combined dataset (Table 1). PCR belongs to the second half (validation part), using samples from patients of The Affiliated Huaian No. 1 People’s Hospital of Nanjing Medical University.

LncRNAs associated with genomic instability were obtained

We utilized the R-package “maftools” to analyze mutation data of gastric cancer samples, calculating the total number of mutations in each sample. Patients were then sorted based on their mutation load, defining the top 25% as genome-stable (GS-like) patients and the bottom 25% as genome-unstable (GU-like) patients. Differences in total mutation load and the expression of UB1LN4, a gene related to genomic instability, were examined between the two groups. The lncRNA expression data for the two groups were obtained by matching transcriptome data, and the “DEseq2” package was used for differential analysis of lncRNAs between the two groups. Genomic instability-related lncRNAs (GlnLncRNAs) were identified (|logFC| > 2 and p < 0.05), and genetic variations were visualized through heatmaps and volcano plots, with the five most up-regulated and down-regulated lncRNAs labeled. Spearman correlation analysis was performed between the 101 differential lncRNAs and coding proteins, with the top 10 mRNAs most associated with each lncRNA being selected. The lncRNA-mRNA network was then mapped to illustrate their relationships. Additionally, functional and pathway enrichment analysis of differential genes between the genome-stable and genome-unstable groups was conducted through the “clusterProfiler” package.

Building a prognostic model

We finally obtained 13 prognosis-related lncRNAs by the univariate Cox analysis. In the training group, we then performed Lasso regression analysis on these 13 genes to construct a prognostic model. We then summed up the risk score by multiplying the amount of gene expression in each patient’s model by the corresponding correlation coefficient in the training set, validation set, and the entire data set. After obtaining a risk score for each patient, we divided the patients in the three datasets into low-risk and high-risk groups based on the median risk value of the patients in the training group. The risk score was calculated with the formula as follows: $\text{Riskscore}={\sum }_{i=0}^n\text{coefi}\ast \text{X}i$. The “coefi” and “X_i” represent the coefficient and expression level of each prognostic lncRNAs.

Evaluation of the model

We assessed the relationship between risk score and prognosis in the training group, test group, and entire dataset, as well as the accuracy of predicting 5-year survival using the area under the ROC curve. Differences in the number of mutations and the expression of UB1LN4, a gene associated with genomic instability, were examined between the high-risk group or the genomically stable group and the genomically unstable group. We analyzed whether the model had the same value across different populations, such as gender, age, and early and late-stage tumors. Univariate and multivariate Cox analyses were used to evaluate whether the model could be used as an independent prognostic indicator for gastric cancer patients in these three datasets. We created a dynamic nomogram representing patient A in the dataset to evaluate the 1-, 3-, and 5-year mortality of patients with gastric cancer. Additionally, calibration curves were plotted to evaluate the accuracy of the model in predicting outcomes at 1, 3, and 5 years.

Correlation between the model and immune cells and checkpoint

We downloaded seven types of algorithms from the TIMER database (http://timer.cistrome.org/) to assess immune infiltration in each patient. We analyzed the expression of immune cells in the high and low-risk groups, identified cells with differential expression (p < 0.05), and developed a heatmap. We also studied the expression of immune checkpoint genes between the high and low-risk groups, extracting genes with different expressions (p < 0.05) and creating a boxplot. Additionally, we used the DREIMT database (http://www.dreimt.org/) to further analyze the correlation between the model and immune cells. We input 80 down-regulated genes and the first 199 up-regulated genes in ascending order of p-value into the website for verification.

Model to predict potential drug candidates for gastric cancer

We utilized the “pRRophetic” package and the expression matrix of gastric cancer patients to predict the minimum drug inhibition concentration (IC50) of drugs in gastric cancer patients. Based on the model, we identified drugs that could become candidates for the treatment of gastric cancer. Additionally, we used the DREIMT database (http://www.dreimt.org/) to predict potential drugs related to the risk for gastric cancer by inputting 80 differential genes in the high and low-risk groups and the first 199 up-regulated genes in ascending order of p-value into the website.

Explore the pathways associated with the model

We downloaded the Hallmarks data set from the Gene Set Enrichment Analysis database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). We then used the Gene Set Variation Analysis (GSVA) algorithm to calculate the enrichment fraction of a special gene set in the Hallmarks data set, based on expression data from the high-risk group of patients with gastric cancer. We analyzed the differences between the high and low-risk groups using LIMMA (p< 0.05), and obtained 42 significantly enriched pathways, which were displayed in the form of a histogram.

External dataset validation and model comparison

We attempted to further assess the accuracy of the signature by exploring data sets associated with gastric cancer from the GEO and ICGC databases. However, we could not find any data set that included all of the lncRNAs in the signature. Therefore, we chose the GSE62254 data set with clinical data to explore the role of LINC00163 in gastric cancer. Additionally, we compared the 5-year accuracy under the curve of our signature with other models.

PCR to verify the expression of LINC00163 in GC

PCR experiments were conducted to compare gene expression levels between 4 pairs of GC and Para-GC tissues from The Affiliated Huaian No. 1 People’s Hospital. Total RNA was isolated from the tissues using the TRIzol reagent. cDNA was synthesized using the PrimeScript™ RT reagent kit. Real-time PCR was performed using the SYBR Green Real-time PCR Master Mix on an Applied Biosystems QuantStudio 7 Flex Real-time PCR System. The relative gene expression levels were calculated using the 2^−ΔΔCt method, with GAPDH as the internal reference gene. The primer sequence of LINC00163 is Forward: GAGCAAGCTCTAGCTCTCGG; Reverse: AGAGCTTTGGGAAGACACCG. GAPDH Forward: AATGGGCAGCCGTTAGGAAA; Reverse: GCCCAATACGACCAAATCAGAG. All experiments were performed in triplicate, and the data were analyzed using GraphPad Prism software.

Statistical analysis

Chip-squared and Rank sum test (Wilcoxtest) test was used to assess the difference in categorical data between different groups as well as datasets. The version of R software was 4.0.5. The statistical significance was defined with the two-tailed p<0.05.

The flow chart

Our research process is summarized in Figure 1.

Figure 1. The flow chart.

Data download and processing

We randomly divided the 348 patients with gastric cancer into the training group and the test group, as well as the combined dataset. The specific information is shown in Table 1, and there were no significant differences in clinical data between the training group and the test group (p > 0.05).

LncRNAs associated with genomic instability were obtained

We obtained 101 genomic instability-related lncRNAs (GlnLncRNAs) with | logFC | > 2 and p< 0.05. We created a heatmap and volcano plot to visualize the differential expression of the GlnLncRNAs, and labeled the top five most significantly up-regulated and down-regulated lncRNAs. The top five down-regulated lncRNAs were MR143HG, PGM5P4-AS1, AC003090.1, ZFHX4-AS1, and the up-regulated lncRNAs were HOXA11-AS, RP11-297P16.4, CQ0ORF91, ELFN1-AS1, RP4-694A7.2. Additionally, we found that the somatic mutation count and the expression of UBQLN4 were different in the GS-like and GU-like groups, and were up-regulated in the GU-like group (p< 0.05) (Figure 2C–2D). We also identified that the lncRNAs and mRNAs formed a network and interacted with each other in gastric cancer (Figure 2E). Furthermore, we conducted functional enrichment analysis of the differentially expressed genes between the GS-like and GU-like groups, and found that they were mainly associated with receptor activity, G protein-coupled peptide receptor activity, cAMP signaling pathway, cell adhesion molecules, calcium signaling pathway, and cGMP-PKG signaling pathway, which are involved in the tumorigenesis and development of gastric cancer (Figure 2F–2G).

Figure 2. (A) Expression of differential genes in GS-like group and GU-like group. The darker the red, the higher the expression level, and the darker the blue, the lower the expression level. (B) The display of differential genes in the volcano map: the top five most up-regulated lncRNAs and the top five most down-regulated lncRNAs were labeled respectively. The blue indicated down-regulated lnRNA, and the red indicated up-regulated lncRNAs. (C) The number of somatic mutations in GS-like group and Gu-like group was different (*** p <0.001). (D) The expression of UBQLN4 gene was different between the GS-like group and the Gu-like group (*** p <0.001). (E) The differentially expressed lncRNAs and the top 10 most related mRNAs formed a co-expression network, with yellow representing lncRNAs and green representing mRNAs. (F) GO enrichment analysis in GS-like group and GU-like group. The top five most significantly enriched functions were extracted from BP, CC and MF (p<0.05). (G) KEGG enrichment.

Construction of the prognostic model

After performing univariable Cox regression, we identified 13 lncRNAs that were significantly associated with prognosis (Figure 3A). Subsequently, we performed Lasso regression analysis on these 13 lncRNAs in the training group, with an optimal Lambda value of 0.01149079 (Figure 3B, 3C). Ultimately, we identified seven GlnLncRNAs, including PTENP1-AS, LINC00163, RP11-169F17.1, C8ORF87, RP11-389G6.3, LINCO1210, and RP11-115H13.1 (Table 2). Using the weights of these seven lncRNAs in the training set, validation set, and combined dataset, we calculated the risk score for each patient. Patients were then divided into high-risk and low-risk groups based on the median score in each dataset. The risk score calculation was as follows: risk score = PTENP1-AS * 0.00635 + LINC00163 * 0.51230 + RP11-169F17.1 * 0.02056 + C8ORF87 * 0.13023 + RP11-389G6.3 * 0.36579 + LINCO1210 * (-0.15732) + RP11-115H13.1 * 0.05948.

Figure 3. (A) After univariate Cox regression, 11 prognostic related ginlncRNAs were obtained, of which 10 were deleterious genes and one was protective gene. (B, C) Lasso regression results showed that when the best lambda value was 0.01149079, the following curve tended to be stable, and 7 lncRNAs were selected for the model construction.

Evaluation of the model

We performed a survival analysis of patients in the high-risk and low-risk groups. Outcomes were significantly worse in the high-risk group than in the low-risk group in the training set (Figure 4A, p=0.009), significantly worse in the validation set (Figure 4B, p<0.001), and similar results were obtained in the merge set (Figure 4C, p <0.001). The ROC curve was used to verify the accuracy of the survival analysis. The area under the curve (AUC) was 0.741 in the training set (Figure 4D), 0.798 in the validation set (Figure 4E), and 0.747 in the merge set (Figure 4F). As pictures shows (Figure 5A, 5F, 5K), LINC00163 was up-regulated in low-risk group, while the others were up-regulated in high-risk group. With the rising of risk score, the probability of death is up. We then analyzed the number of somatic mutations between the high-risk and low-risk groups, and between the Gu-like and GS-like groups. The results (Figure 5B, 5D, 5G, 5I, 5L, 5N) showed that in the three data sets, the number of somatic mutations was higher in the high-risk group or the Gu-like group, and the difference was statistically significant (p<0.05). The expression analysis of UBQLN4 also showed (Figure 5C, 5E, 5H, 5J, 5M, 5O) that the expression of UBQLN4 was mainly upregulated in the high-risk group compared with the low-risk group or the GU-like group comparing with the GS-like group (p<0.05). After univariate and multivariate Cox regressions of risk score and patient clinical characteristics (gender, age, stage) in three datasets, we found that risk score and stage were independent prognostic factors for patients (Figure 6A–6F, p<0.05). The heatmap (Figure 7A) showed the relationship between 7 lncRNAs and stage, age, sex, risk score, GS-like, and Gu-like, finding that there is a difference in GS-like and GU-like group as well as age group (p<0.001). Then, all patients were divided into two groups according to gender, age (greater than 65 years old, less than or equal to 65 years old), and stage respectively. The results showed that high-risk patients had a poor prognosis in groups by sex, age, and stage (Figure 7B–7G, p<0.05), which suggested that the signature had a similar result in different types of gastric patients. Dynamic Nomogram (Figure 8A) was used to predict the 1,3,5year overall survival of patients with gastric cancer and the calibration curve shows that the signature has a high prediction accuracy (Figure 8B–8D).

Figure 4. (A–C) Prognostic differences between the high-risk and low-risk groups were investigated in the training set, validation set, and the entire dataset. High risk was found to indicate poor prognosis in all three datasets (p<0.001). (D–F) In the training set, validation set and the entire set, the AUC values obtained by the prognostic model to predict the accuracy of the 5-year survival rate of patients were 0.741, 0.798 and 0.747, respectively.

Figure 5. (A) In the training set, we can see the distribution of patient survival and death as the risk score increased, as well as the heat map of lncRNAs expression in the high-risk and low-risk groups. (B, C) In the training set, the number of somatic mutations and the expression of UBQLN4 gene were different between the GS-like group and the Gu-like group (** p<0.01, *** p<0.001). (D, E) In the training set, the number of somatic mutations and the expression of UBQLN4 gene were different between the high-risk group and the low-risk group (** p<0.01, *** p<0.001). (F) In the validation set, we can see the distribution of patient survival and death as the risk score increased, as well as the heat map of lncRNAs expression in the high-risk and low-risk groups. (G, H) In the validation set, the number of somatic mutations and the expression of UBQLN4 gene were different between the GS-like group and the Gu-like group (** p <0.01, *** p<0.001). (I, J) In the validation set, the number of somatic mutations and the expression of UBQLN4 gene were different between the high-risk group and the low-risk group (** p <0.01, *** p <0.001). (K) In the entire set, we can see the distribution of patient survival and death as the risk score increased, as well as the heat map of lncRNAs expression in the high-risk and low-risk groups. (L, M) In the entire set, the number of somatic mutations and the expression of UBQLN4 gene were different between the GS-like group and the Gu-like group (** p<0.01, *** p <0.001). (N, O) In the entire set, the number of somatic mutations and the expression of UBQLN4 gene were different between the high-risk group and the low-risk group (** p <0.01, ***p <0.001).

Figure 6. (A) Univariate Cox regression for risk-score, gender, age, and tumor stage in the training set found that risk score and stage were independent prognostic factors (p<0.05). (B) Multivariate Cox regression for risk-score, gender, age, and tumor stage in the training set found that risk score and stage were independent prognostic factors (p<0.05). (C) Univariate Cox regression for risk-score, gender, age, and tumor stage in the validation set found that risk-score and stage were independent prognostic factors (p<0.05). (D) Multivariate Cox regression for risk-score, gender, age, and tumor stage in the validation set found that risk score, age, and stage were independent prognostic factors (p<0.05). (E) Univariate Cox regression for risk-score, gender, age, and tumor stage in the entire set found that risk-score, age, and stage were independent prognostic factors (p <0.05). (F) Multivariate Cox regression for risk-score, gender, age, and tumor stage in the entire set found that risk score, age, and stage were independent prognostic factors (p <0.05).

Figure 7. (A) Heat maps of the relationship between ginlncRNAs in high-risk group and low-risk group, GU-like group and GS-like group, tumor stage, age, and sex (* p<0.01, ** p<0.01, *** p<0.001). (B) The risk values of the constructed model were associated with poorer prognosis in male gastric cancer patients (p<0.05). (C) The risk value of the established model was associated with poor prognosis in patients with gastric cancer aged <=65 years (p<0.05). (D) The risk value of the constructed model was associated with poor prognosis in both early gastric cancer patients (p<0.05). (E) In female patients (p<0.05); (F) in aged >65 patients (p<0.05); (G) in late stage patients (p<0.05).

Figure 8. (A) A Nomogram for predicting prognosis based on risk value, age, tumor stage, and sex shows the clinical characteristics of the patient (TCGA-CG-571) and the predicted 1 -, 3 - and 5-year mortality of 0.177, 0.472, and 0.605, respectively. (B–D) The 1 -, 3 - and 5-year calibration curve constructed by risk value was used to evaluate the accuracy of prognosis prediction based on risk score, and the results were all satisfactory.

Correlation between the model and immune cells

The immune-related heatmap displayed the immune landscape based on high- and low-risk groups, where the infiltration level of various immune cells was obtained using different algorithms such as TIMER, CIBERSORT, QUANTISEQ. As shown in Figure 9A, there was a noticeable difference in the expression of immune cells between the two groups, with a trend of immune cells being up-regulated in the high-risk group compared to the low-risk group, such as macrophages, T cells CD4, and T cells CD8. To further verify the accuracy of this finding, we predicted the potential immune cells using the DREIMT database, which showed a similar result indicating that macrophages were highly associated with the risk groups (Figure 9B, 9C, p<0.05). We also analyzed the relationship between the high-risk and low-risk groups and immune checkpoint genes, where we found that the expression of immune checkpoint genes was mainly up-regulated in the high-risk group (Figure 9D, p<0.05).

Figure 9. (A) Heat maps of immune cells that differ between the high-risk and low-risk groups. (B, C) Macrophage was the most different immune cells predicted by Dreimt database among the high-risk and low-risk groups (p<0.05). (D) Immune checkpoint genes that were differentially expressed in the high-risk and low-risk groups were selected and presented in a boxplot. (* p< 0.05, **p < 0.01, * * * p< 0.001).

Potential drug candidates for gastric cancer

Our analysis showed that certain drugs, including Sunitinib and Shikonin, had lower IC50 concentrations in the high-risk group (Figure 10A). These findings suggest that patients in the high-risk group may be more sensitive to these drugs. In addition, using the DREIMT database, we predicted and identified candidate drugs for future treatment strategies (Figure 10B–10D). Our top 5 candidate drug types were adrenergic receptor antagonists, cyclooxygenase inhibitors, dopamine receptor antagonists, estrogen receptor agonists, and serotonin receptor antagonists. These results could inform personalized treatment approaches in the future.

Figure 10. (A) Chemotherapy drugs with differential IC50 expression in the high-risk and low-risk groups predicted by the Prophetic package, and drugs with low IC50 expression in the high risk group were screened out (* * p < 0.01, ***p < 0.001). (B) The drugs predicted by DREMIT website are different in the high and low risk groups, in which the best candidates in the first and second quadrants mean the candidates with more reliable results, while the remaining drugs mean better candidates. (C) Classification of predicted chemotherapeutic drug validation levels. (D) Drug classification of predicted chemotherapeutic agents.

Exploration of the pathways associated with the model

We conducted a GSVA enrichment analysis and identified 42 significantly enriched pathways, as shown in Figure 11. This analysis allowed us to explore differences in signaling pathways between the high-risk and low-risk groups. Among the pathways associated with cancer, we found that the EMT signaling pathway, P13-AKT-TOR signaling pathway, and WNT-BETA-CATEIN signaling pathway were highly enriched. These pathways are known to play important roles in tumorigenesis and development in gastric cancer.

Figure 11. Through GSVA analysis, significant enrichment pathways with differences in the high-risk and low-risk groups were obtained (p <0.05). Pathways enriched in the high-risk group were marked in red, those enriched in the low-risk group were marked in blue, and those enriched in -2

External dataset validation and model comparison

We investigated the role of LINC00163 in the GSE62254 dataset and found that it was associated with prognosis, as shown in Figure 12A. Specifically, up-regulated LINC00163 expression was associated with a good prognosis in gastric cancer. This finding is consistent with experimental studies, which have shown that up-regulated LINC00163 can suppress gastric carcinoma development [26]. In our signature, the coefficient of LINC00163 was -0.15732, indicating that high expression of LINC00163 may result in a low-risk score and a good prognosis. As shown in Figure 12B–12D, we found that LINC00163 expression was similar across different gender, age, and stage groups.

Figure 12. (A) KM survival curve of LINC00163 in gastric cancer:High expression of LINC00163 was associated with poor prognosis (p<0.001). (B–D) Expression of LINC00163 in different stages, genders and ages of gastric cancer. (E) Our model was compared with other models by plotting a 5-year ROC curve associated with prognosis.

We also compared our signature to other models, such as ChengLncRNAs and SongLncRNAs, and found that our model had a better prognostic value, with an area under the curve of 0.747 for 5-year prognosis (Figure 12E). This suggests that our signature may be a useful tool for diagnosing prognosis in gastric cancer patients [27, 28].

PCR to verify the expression of LINC00163 in GC

PCR was performed to verify the expression of LINC00163. The results showed that the expression of LINC00163 was significantly up-regulated in 4 adjacent normal tissues compared with the GC tissues. (Figure 13, **p<0.01).

Figure 13. PCR showed that LINC00163 was significantly up-regulated in 4 adjacent normal tissues compared with the GC tissues. **p<0.01.