Identification of a novel MSI-related ceRNA network for predicting the prognosis and immunotherapy response of gastric cancer

Differentially expressed mRNAs (DEmRNAs), miRNA (DEmiRNAs) and lncRNAs (DElncRNAs) between MSIH and MSIL/S

Supplementary Table 1 shows the patient information collected in this study, and the proportion of patients with different clinical characteristics. DEmRNAs, DEmiRNAs and DElncRNAs between the MSIH and MSIL/S groups were analyzed based on the 343 samples with MSI information from TCGA-STAD. The screening conditions were set at |log2 (fold change) |>1 and FDR<0.05 for mRNA and lncRNA, and |log2 (fold change) |>0.5 and FDR<0.05 for miRNA. The results in Figure 1A–1C show that there were 900 DEmRNAs of which 44 and 856 were up and down regulated, respectively; 96 DEmiRNAs of which 29 and 67 were up and down regulated, respectively; and 60 DElncRNAs of which with 9 and 51 were up and down regulated, respectively.

Construction of an MSI-related lncRNA-miRNA-mRNA ceRNA network

Based on the screened DEmRNAs, DEmiRNAs and DElncRNAs, the “gdcCEAnalysis” function in the “GDCRNATools” package was used to construct the MSI-related ceRNA network. The cut-off P values for the hypergeometric and Pearson correlation tests were both set at 0.01 and the regulation similarity value was required to be higher than 0. The constructed ceRNA network consisting of 2,729 lncRNA-mRNA pairs, 13,334 edges and 895 nodes (808mRNA, 44 miRNAs and 43 lncRNAs), is shown in Supplementary Table 2. The table also contains information on LncRNA-miRNA-mRNA pairs, hypergeometric distribution, correlation test details, regulation similarity and sensitivity correlation. Enrichment analysis was conducted based on the 808 mRNAs. The findings in Figure 1D, 1E show that 44 KEGG pathways were enriched, including MAPK signaling pathway, focal adhesion and Wnt signaling pathway. On the other hand, GO analysis revealed that the most enriched biological process (BP) was the muscle system process, the most enriched cellular component (CC) was the collagen-containing extracellular matrix, and the most enriched molecular function (MF) was channel activity.

Screening of the hub MSI-related ceRNA with prognostic significance

We conducted survival analysis on 43 lncRNAs in the constructed MSI-related ceRNA network. Three lncRNAs with prognostic significance were identified, including MIR99AHG, HAGLR and LINC02381 (Figure 2A). Based on the sequences downloaded from NCBI, the lncLocator prediction results (Figure 2B) showed that only MIR99AHG was located in cytoplasm. Therefore, MIR99AHG was selected as the hub lncRNA and the MIR99AHG-related ceRNA network was shown in Figure 2C. This subnetwork contained 87 lncRNA-mRNA pairs, 1116 edges and 109 nodes (1 lncRNA, 21 miRNAs and 87 mRNAs). Additionally, Unicox analysis was conducted, targeting the nodes in this subnetwork (Supplementary Table 3). Only nodes with prognostic significance were retained for constructing the final hub ceRNA network (Figure 2D). This network is comprised of 26 lncRNA-mRNA pairs, 43 edges and 29 nodes (1 lncRNA, 2 miRNA and 26 mRNA). Of the two miRNAs, let-7f-5p was identified as a protective factor while mir-125a-5p was shown to be the risk factor. All the 26 mRNAs were risk-associated genes.

Establishment and validation of the prognosis prediction model based on the MSI-related ceRNA regulatory network

The Lasso Cox regression algorithm was applied to establish the prognosis prediction model, based on the 26 mRNAs in the final hub ceRNA network (Supplementary Figure 2). The constructed model comprised of 7 genes and the risk score was calculated as follows: RiskScore=IL1RL1*0.042604+SPAG16*0.069518+FAM110B*0.050744+ANKRD6*0.017904+ACSS3*0.068896+CORO2B*0.034964+TNFAIP8L3*0.036276. Patients were stratified into the high- and low-risk groups based on the predicted score. In addition, KM curves were plotted, and survival difference was analyzed. The results revealed that high-risk patients had worse OS (Figure 3A). Moreover, ROC curve analysis (Figure 3B) showed that the AUC of the prediction model was 0.76. Results from PCA (Figure 3E) showed that the distribution of these 2 groups was diverse, and could not be separated linearly. Moreover, the scatter plots (Figure 3E) showed that, with increased risk scores, patients had a shorter lifespan and a higher possibility of death. Finally, a heatmap was used to show the expression patterns of the 7 genes used in the model. As shown in Figure 3F, all the genes were risk factors for GC patients. Additionally, the precited risk score was associated with some clinical characteristics, therefore correlation tests were conducted. The findings revealed that a high risk was associated with MSIL/S, a higher grade and a younger age (Supplementary Figure 3A).

To test the generalization ability of the established model, the above tests were repeated in an independent GEO dataset (n=689). The overall results were consistent with those obtained from the TCGA database, further validating our model. Figure 4A shows that high-risk patients had worse OS and Figure 4B shows that the AUC was 0.7. The PCA results (Figure 4C) showed that the distribution of the 2 groups was different and the scatter plot (Figure 4E) confirmed that the mortality rate increased with a higher risk score. Finally, the heatmap showed that the 7 genes were all risk factors, and a higher risk was associated with a younger age (Supplementary Figure 3B).

To further test whether the predicted risk score was an independent prognosis index for GC patients, this study integrated the risk score with other clinical factors into Unicox and Multicox analyses. Results for the TCGA database are shown in Figure 5A, 5B. Unicox analysis revealed that age, stage, T, M, N, and the risk score, were risk factors for GC patients (Figure 5A). On the other hand, Multicox analysis confirmed that age, M and risk score, were independent prognostic factors for GC (Figure 5B). Results for the GEO database were shown in Figure 5C, 5D. Unicox analysis showed that age and risk score were risk factors for GC patients (Figure 5C). Moreover, Multicox analysis identified these 2 factors (age and risk score) as independent prognostic factors for GC (Figure 5D). Furthermore, Differential expression gene (DEG) analysis was conducted in the TCGA database between the high- and low-risk groups. The results showed that 913 genes were up regulated in the high-risk group while 30 genes were downregulated. KEGG and GO enrichment analyses were then conducted on these DEmRNAs. A total of 52 pathways (Figure 5E) were enriched in the KEGG analysis and the most interesting included the PI3K-Akt signaling pathway, focal adhesion and vascular smooth muscle contraction. After comparing these pathways with the 43 enriched pathways in Figure 1D (representing the MSI-related ceRNA regulatory network), 28 pathways were shown to be conserved (Supplementary Table 4). This demonstrated that the constructed model generally represented the MSI-related ceRNA regulatory network. Finally, GO analysis showed that the most enriched BP, CC, and MF were extracellular matrix organization, collagen-containing extracellular matrix and the extracellular matrix structural constituent, respectively (Figure 5F).

Correlation between risk score and IPS

315 IPS scores for TCGA-STAD patients were downloaded from the TCIA database and integrated into patient data. The IPS score can be classified into 2 types, where one type represents the possibility of response to the PD1/PDL1/PDL2 blocker, while the other represents the possibility of response to the CTLA4 blocker. Additionally, patients can be divided into 3 groups based on the IPS score: 4-6, 7-8 and 9-10. With respect to the PD1/PDL1/PDL2 blocker, these 3 groups of patients accounted for 18%, 66% and 16%, respectively, while for the CTLA4 blocker, these 3 categories of patients accounted for 7%, 57% and 36%, respectively. Differences in risk scores among the 3 groups were compared and the results are shown in Figure 6A, 6B. In both therapies, the risk score of the 9-10 group was significantly lower than that of the 7-8 category, indicating that a lower risk score was associated with a better ICI response. These results demonstrated that our calculated signature could not only predict OS but also forecast the ICI response to a certain degree. To further explore this interesting observation, we assessed the correlation between the expression of the hub nodes (MIR99AHG and 7 genes used in the model) and 3 immune checkpoint genes (PD-L1(CD274), PD1(PDCD1), CTLA4). The findings showed that MIR99AHG, ANKRD6, ACSS3, CORO2B, FAM110B, and SPAG16 were negatively associated with PD-L1 (Figure 6C). Additionally, TNFAIP8L3 was positively correlated with PD1, while IL1RL1 was positively correlated with CTLA4. However, SPAG16 was negatively correlated with CTLA4.

Analysis of immune cell infiltration between different risk groups

This study simultaneously investigated immune cell infiltration and chemokine gene expression among different risk groups in the TCGA and GEO databases. The results revealed that there were 8 immune cells with consistent infiltration in the 2 databases (Figure 7A, 7B). T cells CD4 memory resting, Mast cells resting, B cells naïve and Monocytes were infiltrated in the high-risk group, while T cells CD4 memory activated, Mast cells activated, Neutrophils and Macrophages M1 were mainly infiltrated in the low-risk group. In addition, the stromal score and the estimate score in the high-risk group were markedly higher, demonstrating that the patients had lower tumor purity. And the risk score was positively correlated with the stromal score and the immune score (Supplementary Figure 4). Additionally, 17 chemokines were consistently expressed in the 2 databases (Figure 7C, 7D). A total of 14 chemokines (XCL1, CCL2, CCL11, CCL13, CCL16, CCL17, CCL19, CCL21, CCL22, CCL23, CXCL12, CXCL13, CXCL14 and CX3CL1) were highly expressed in the high-risk group, while three others (CCL20, CXCL1 and CXCL2) were mainly expressed in the low-risk group. The correlation of the above immune cells of interest with chemokines was further studied. Eight co-expression patterns were identified among 5 immune cells and 6 chemokines (Co-expression patterns were denoted as “chemokine ~ immune cell”), as shown in Figure 7E, 7F. They included CX3CL1 ~ T cells CD4 memory activated, CXCL14 ~ B cells naïve, CXCL12 ~ mast cells resting, CXCL1 ~ neutrophils, CXCL2 ~ neutrophils, CXCL1 ~ mast cells activated, CXCL2 ~ mast cells activated and CCL19 ~ B cells naïve. This study also investigated the association between the hub lncRNA MIR99AHG and the 6 identified chemokines. MIR99AHG was observed to be positively related to CXCL14, CXCL12 and CCL19, but negatively correlated to CX3CL1, CXCL1 and CXCL2 (Figure 7G). Finally, survival analyses were performed on the five identified immune cells in the 2 databases. The results showed that Mast cells resting and T cells CD4 memory activated were significantly correlated with the OS of patients in the GEO database (Figure 7H). Moreover, Mast cells resting was identified as a risk factor while T cells CD4 memory activated had the opposite effect.

Validation of the expression and regulatory relationship of the ceRNA network

To further validate our findings, we quantified the expression levels of MIR99AHG and the crucial genes in the prognostic model (IL1RL1, SPAG16, FAM110B, ANKRD6, ACSS3, CORO2B and TNFAIP8L3) in 15 GC samples and their adjacent nontumor tissues. The PCR results revealed that MIR99AHG, IL1RL1, SPAG16, ANKRD6 and ACSS3 were upregulated in GC tissues, while CORO2B and TNFAIP8L3 were downregulated (Figure 8). There was no significant difference in the expression of FAM110B between GC samples and the corresponding nontumor tissues. Due to the limited number of high-MSI samples, additional MSIH samples are needed to verify the expression levels of the above genes in different MSI types.

What’s more, we confirmed that MIR99AHG was overexpressed in GC cell lines compared with normal gastric epithelial cell (Supplementary Figure 5A). Knockdown or overexpression of MIR99AHG could negatively regulate the expression level of let-7f-5p and miR-125a-5p, and positively regulate the expression level of IL1RL1, SPAG16, FAM110B, ANKRD6, ACSS3, CORO2B and TNFAIP8L3 (Supplementary Figure 5B, 5C). The above results confirmed the regulatory relationship among the ceRNA network molecules. In addition, we analyzed the effect of MIR99AHG knockdown and overexpression on PD-L1 mRNA and protein expression, and also confirmed that MIR99AHG could negatively modulate PD-L1 expression (Supplementary Figure 5D–5F). The specific mechanism needs to be further elaborated.

Research design and data preprocessing

The study was designed according to the flow diagram (Supplementary Figure 1).

Three GEO microarray datasets and the TCGA-STAD dataset were downloaded. Specifically, expression data from 689 samples and their corresponding clinical information were extracted from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), including 200 samples from GSE15459, 56 from GSE34942 and 433 samples from GSE84437. The “Combat” function in the “sva” R package was used to remove the batch effect. Additionally, the mRNA, lncRNA and miRNA expression data from 348 patients and their clinical characteristics, were obtained from the TCGA GDC pipeline (https://portal.gdc.cancer.gov/). Moreover, immunophenoscore (IPS) and MSI data for patients in the TCGA database, were obtained from the Cancer Immunome Atlas (TCIA, https://tcia.at/home).

Screening for differentially expressed mRNAs (DEmRNAs), miRNA (DEmiRNAs) and lncRNAs (DElncRNAs)

Referring to the revised Bethesda guidelines for colorectal cancer, samples were divided into two groups: samples displaying instability at two or more of the five recommended loci (BAT-25, BAT-26, NR-21, NR-24, NR-27) were interpreted as MSIH (19.2%), while tumors with only one or no locus altered were identified as MSIL/S (80.8%). The “gdcDEAnalysis” function in the “GDCRNATools” package was applied to screen for differentially expressed genes between the MSIH and MSIL/S groups [40]. The screening conditions were set at |log2 (fold change) |>1 and FDR<0.05 for mRNA and lncRNA, and |log2 (fold change) |>0.5 and FDR<0.05 for miRNA. A volcano plot was used to display the DEmRNAs, DEmiRNAs, and DElncRNAs.

Construction of the lncRNA-miRNA-mRNA network

The “gdcCEAnalysis” function in the “GDCRNATools” package [16] was used to construct the ceRNA regulatory network based on the genes screened in step 2 (above). This function screens credible lncRNA-mRNA pairs based on the following criteria [41]: (1) a significant number of miRNAs must be shared between the lncRNA-mRNA pair (hypergeometric test, p<0.01); (2) a positive correlation must exist between this lncRNA-mRNA pair (Pearson correlation, Cor>0 and p<0.01); and (3) the common miRNAs should play similar roles in regulating the expression of the lncRNA and mRNA (regulation similarity>0) [42]. Notably, the database used to search for possible lncRNA -miRNA and mRNA-miRNA pairs can be set in the “gdcCEAnalysis” function. To explore as many pairs as possible, the parameter was set as “miRcode”, which constructed the largest network. Finally, the “clusterProfiler” package was used to conduct enrichment analysis targeting the mRNAs in the construed network, to reveal the potential functions and mechanisms of this ceRNA regulatory network.

Screening for the hub ceRNA network with prognostic significance

The constructed ceRNA network from the previous step comprised of 2729 lncRNA-mRNA pairs and 895 nodes (808 mRNAs, 44miRNAs and 43 lncRNAs). Survival analysis was conducted on the 43 lncRNAs to identify the hub lncRNAs with prognostic significance. According to literature, lncRNAs exert their function based on their cellular position and sponging lncRNAs are mainly located in the cytoplasm [43]. Therefore, location prediction was conducted following survival analysis. In addition, nucleotide sequences for the lncRNAs of interest were downloaded from the NCBI nucleotide website (https://www.ncbi.nlm.nih.gov/nuccore/), then uploaded to the lncLocator database (http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/) [44] to predict the possible positions of the lncRNAs. The result from this database revealed 5 possible locations and the corresponding possibility. Notably, MIR99AHG was targeted as the hub lncRNA, as it was located in the cytoplasm and was a significant risk factor for the overall survival of GC patients. Finally, we screened for the MIR99AHG-related mRNAs and miRNAs, then performed Unicox analysis on the related genes. Only genes with prognostic significance were retained and the final hub ceRNA regulatory network was constructed.

Construction and validation of a prognosis prediction model

Based on the mRNAs in the final hub ceRNA regulatory network, the “glmnet” package was used to construct the prognosis prediction model. Thereafter, patients were stratified into the high- and low-risk groups based on the median value of the predicted risk score. In addition, survival differences between these 2 groups were analyzed using the Kaplan-Meier curve method. The efficiency of the model was tested through the ROC curve which was plotted using the “ROCR” package. Moreover, PCA analysis was conducted to test whether the predicted risk score could successfully stratify patients based on risk. MultiCox analysis was also conducted to show whether the predicted risk score was an independent prognostic factor in GC patients. To test the generalization ability of the established model, all the above analyses were conducted in the GEO datasets which have been detailed in the “Data Acquisition and Preparation” section. Furthermore, correlation analysis was conducted to check for consistencies between the risk score and existing clinical characteristics. Finally, differentially expressed genes between the 2 risk groups were screened using FDR≤0.05 and | logFCfilter | ≥1, and GO and KEGG enrichment analyses were conducted based on these differentially expressed genes.

Correlation analysis between the predicted risk score and ICI response predictor

The IPS scores for 415 patients were downloaded from the TCIA database, which uses expression data and a machine learning-based algorithm to predict ICI (PD1/PDL1/PDL2 blocker or CTLA4 blocker) response in 20 types of solid tumors (range from 0-10, higher score means higher possibility to response to ICI) [10]. The downloaded TCGA-STAD IPS scores ranged from 4-10. Therefore, we divided patients into 3 different groups (4-6, 7-8 and 9-10) and the Wilcoxon rank test was used to show differences in the risk scores between the 3 categories. In addition, the study assessed the correlation between the hub genes in the final prognostic model and the PD-L1/PD1/ CTLA4.

Analysis of immune cell infiltration between different risk groups

Basic data files for 22 types of immune cells and the R script for infiltration scoring were obtained from the CIBERSORT website (https://cibersort.stanford.edu/) [45]. The immune-cell infiltration score was calculated and compared between different risk groups. The immune, stromal, and the estimate score were calculated for each sample using the “estimate” package. And then we used “ggplot2” package to reveal differences in the immune score and stromal score between the high- and low-risk groups and the relationship between risk score, immune score, and the stromal score [46]. Notably, chemokines are a type of small cytokines with chemotactic functions. They guide cell migration and localization between tissues and play an important role in embryogenesis, tissue development and immune response [47]. Studies have shown that chemokines play an important role in immune cell infiltration in the tumor microenvironment (TME) [48]. Therefore, this study obtained the list of 40 human chemokines, and analyzed differences in their expression between different risk groups. Pearson correlation analysis was used to find the potential regulatory chemokine and immune cell pairs (denoted as “chemokine ~ immune cell”) based on the differentially expressed immune cells and chemokines between the 2 risk groups. Furthermore, the correlation between MIR99AHG and these chemokines was analyzed. Finally, patients were divided into two groups according to the median infiltration score of the immune cells of interest, and differences in survival between the two categories were analyzed using the KM curve.

RNA isolation and qRT-PCR

Gastric cancer specimen and adjacent nontumor tissues were obtained from 15 individuals with GC who had undergone gastric cancer surgery between January 2021 and January 2022 in the Xijing Hospital of Digestive Diseases. The TRIzol reagent (Invitrogen, Waltham, MA, USA) was used to isolate total RNA from both GC and the adjacent nontumor tissues, according to the manufacturer’s protocol. Thereafter, 2 μg of the extracted total RNA was used for reverse transcription using the PrimeScript™ RT Master Mix kit (Takara). qRT-PCR analyses for the expression of the MIR99AHG, IL1RL1, SPAG16, FAM110B, ANKRD6, ACSS3, CORO2B and TNFAIP8L3 genes were conducted on the CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA), including with 45 cycles of 95° C for 5 s and 60° C for 30 s. The primers are listed in Supplementary Table 5. The 2^−ΔΔCt method was employed to analyze the relative changes in gene expression.

Statistics

In this study, the Wilcoxon rank test was applied to compare indexes between the high and low risk groups. Pearson analysis was used to test for the correlation between variables of interest and the KM curve was employed for survival analysis. All the statistical analyses were conducted in R (Version 4.0.2 (2020-06-22)), and p<0.05 was considered to be statistically significant, unless otherwise stated. A Github page (https://github.com/itkwffc/MSI-STAD.git) were created and all scripts and supported data were uploaded.