Construction of stomach adenocarcinoma prognostic signature based on anoikis-related lncRNAs and clinical significance

ABSTRACT

As a dominant type of gastric cancer, stomach adenocarcinoma (STAD) is characterized by high morbidity and mortality rates. Anoikis factors participate in tumor metastasis and invasion. This study was designed to identify prognostic risk factors in anoikis-related long non-coding RNAs (lncRNAs) for STAD. First, with STAD expression datasets and anoikis-related gene sets downloaded from public databases, anoikis-related prognostic lncRNA signatures (AC091057.1, ADAMTS9.AS1, AC090825.1, AC084880.3, EMX2OS, HHIP.AS1, AC016583.2, EDIL3.DT, DIRC1, LINC01614, and AC103702.2) were screened by Cox regression to establish a prognostic risk model. Kaplan–Meier and receiver operating characteristic curves were used to evaluate the survival status of patients and verify predictive accuracy of the model. Besides, risk score could be an independent prognostic factor to assess the prognosis of STAD patients. Nomograms of the prognostic model that combined clinical information and risk score could effectively predict survival of STAD patients, as validated by calibration curve. Gene ontology and Kyoto encyclopedia of genes and genomes enrichment analyses were performed for differentially expressed genes (DEGs) in high- and low-risk groups. These DEGs were related to neurotransmitter transmission, signal transmission, and endocytosis. Moreover, we analyzed immune status of different risk groups and found that STAD patients in low-risk group were more sensitive to immunotherapy. A prognostic risk assessment model for STAD using anoikis-related lncRNA genes was constructed here, which was proven to have high predictive accuracy and thus could offer a reference for prognostic evaluation and clinical treatment of STAD patients.

KEYWORDS:

• anoikis
• gastric adenocarcinoma
• lncRNA
• prognosis
• immunity

1. Introduction

Ranking fifth in the incidence rate among all cancers, gastric cancer (GC) is a malignancy of the digestive system, with stomach adenocarcinoma (STAD) being the most common subtype [1]. Radical surgery combined with chemotherapy at the early stage is the most efficacious treatment for STAD, although postoperative metastases and recurrence remain a problem for some STAD patients [2]. However, most patients with STAD are diagnosed at a rather advanced stage which is no longer the best time for surgery. Prognosis of patients with advanced STAD is far from desired with a median survival of only 8 months [3]. Conventional therapies yield limited efficacy in patients with advanced GC, and immunotherapy may emerge as another efficient treatment after surgery, chemotherapy, radiotherapy, and target therapy [3]. A prognostic model that enables efficient risk stratification of STAD patients will facilitate the prediction of the overall survival of patients and guide prognostic management.

As a specific form of programmed cell death, anoikis refers to apoptosis triggered by the detachment of normal epithelial cells or solid tumor cells without metastatic nature in situ into the blood [4]. When cells are detached from the extracellular matrix, anoikis prevents tumor cells isolated from the primary tumor from continuing to survive [5]. In this way, the anoikis factor is vital in the development and metastasis of STAD. Analysis of the prognosis of STAD based on the anoikis factor is conducive to the correct assessment of tumor development and metastasis risk. More than 200 nt in length, long non-coding RNAs (lncRNAs) are involved in many human cancers such as GC, lung cancer, liver cancer, pancreatic cancer, and colorectal cancer. Besides, lncRNAs are crucial in various biological processes such as cell proliferation, invasion, apoptosis, and stem cell maintenance [6–10]. Li et al. pointed out that lncRNAs may be potential targets for cancer diagnosis and therapy [10]. There is association between lncRNA expression and anoikis resistance and metastasis in cancer [11–13]. Investigating the effect of anoikis-related lncRNAs on STAD may help doctors better evaluate the prognosis of STAD patients and explore possible targets for disease treatment.

However, to our knowledge, no studies have investigated the effect of anoikis-related lncRNAs on prognosis of STAD. With prognostic markers screened from anoikis-related lncRNAs in STAD, this study built a prognostic model of STAD to predict survival outcomes and draw nomograms with risk score to facilitate prognostic management of STAD patients. Besides, the immune status and response to immunotherapy in high- and low-risk groups were assessed to offer a reference for immune checkpoint blockade therapy and other immunotherapies for STAD.

2. 2. Materials and methods

2.1. Data extraction and differential analysis

From The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/), TCGA-STAD dataset (32 normal samples and 375 STAD samples) including RNA-seq expression data, SNV mutation data, and corresponding clinical data was downloaded. Anoikis-related mRNAs were acquired from Genecards (https://www.genecards.org/).

Differential analysis of mRNA expression retrieved from TCGA-STAD was performed using the edgeR package [14] to screen differentially expressed genes (DEGs) (FDR <0.05, logFC > 1). DEGs and anoikis-related mRNAs were intersected to obtain differentially expressed anoikis-related mRNAs.

2.2. Construction of an anoikis-related lncRNA-based STAD risk prognostic model

With clinical information collected from the dataset, samples that survived 30 days were screened and randomly assigned into training and validation sets at a ratio of 7:3. LncRNAs strongly correlated with differentially expressed anoikis-related mRNAs were selected as the study subjects (cor ≥0.5). Univariate Cox regression analysis of selected lncRNAs was carried out using the survival package [15]. The dimensionality of the model was reduced by LASSO analysis using the glmnet package [16] to avoid overfitting. Following that, multivariate Cox regression analysis of the screened lncRNAs was performed with the survival package to construct the final prognostic model.

2.3. Validation of the prognostic model for STAD

The STAD samples were assigned into high- and low-risk groups by the median risk score. Kaplan–Meier survival curves of two groups were drawn using survival and survminer packages (https://rdocumentation.org/packages/survminer/versions/0.4.9). The training set was then analyzed using the time ROC package [17] and the survival package. Receiver operating characteristic (ROC) curves for 1, 3, and 5 years were plotted using area under the curve (AUC) calculated. Finally, the validation set was utilized to assess predictive accuracy and performance of the model.

2.4. Independent prognostic factor analysis

Univariate and multivariate Cox regression analyses that included risk score of STAD samples and clinical information (age, gender, pathological stage, and other characteristics) were performed using the survival package for identifying independent prognostic factors. Forest plots were drawn using the forestplot package. The rms package [18] and survival package were used to draw nomograms based on risk score and clinicopathologic factors (age, gender, T, and N) to further confirm the clinical application of the model. Calibration curves were drawn with the rms package [18] to assess whether the predicted overall survival by nomogram was in accordance with the actual overall survival (1-year, 3-year, and 5-year).

2.5. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses

The edgeR package [14] was applied for the differential analysis of high- and low-risk groups to screen DEGs (FDR <0.05, logFC > 1). GO and KEGG enrichment analyses of DEGs in the two risk groups were performed separately using the clusterProfiler package [19].

2.6. Analyses of immune infiltrate and immune checkpoint

Tumor microenvironment (TME) of each STAD sample was assessed using the estimate package [20]. Immune score, stromal score, and ESTIMATE score of each sample were calculated, and the Wilcox test was performed. Violin plots were plotted for the two risk groups.

Then, single-sample gene set enrichment analysis (ssGSEA) using GSVA package [21,22] was conducted to compare immune-related cells and functions in the two risk groups, respectively. Correlations between risk score, immune-related cell score, and immune-related function score were calculated to reveal the relation between immune level and risk score.

The expression of PDCD1/CTLA4 in the two groups was analyzed to assess the effectiveness of immunotherapy for patients in different risk groups. Tumor Immune Dysfunction and Exclusion (TIDE) scores were calculated on TIDE (tide.dfci.harvard.edu) to predict anti-tumor immune escape level in STAD patients during immunotherapy.

3. Results

3.1. Differential analysis, anoikis-related signature gene screening, and lncRNA-based prognostic model construction

After subjecting samples in TCGA-STAD to gene differential expression analysis, a total of 4705 DEGs in STAD that met the threshold (FDR <0.05, logFC > 1) were screened. DEGs and anoikis-related mRNAs (Table S1) were intersected to obtain 185 anoikis-related differentially expressed mRNAs (AnoRGs) (Figure 1a).

Figure 1. Differential analysis and construction of a STAD prognostic model.

Note: (a) Upset plots between DEGs and anoikis-related mRNAs; (b-c) coefficient distribution plots generated for log (λ) sequences in LASSO analysis as well as change trajectories; (d) forest plots for the multivariate Cox regression model based on lncRNAs.

A total of 496 lncRNAs that were in relatively high correlation (cor ≥0.5) with AnoRGs were screened. Referring to expression of these 496 lncRNAs and clinical information, samples with survival time >30 days were screened out and randomly divided into a training set (Table S2) and a validation set (Table S3) at a 7:3 ratio. Subsequently, univariate Cox regression analysis screened out 17 candidate lncRNAs dramatically associated with survival (FDR <0.05), among which 11 important feature lncRNAs were eventually selected after the LASSO Cox analysis (Figure 1b-c). A multivariate Cox regression model was constructed based on feature lncRNAs yielded from LASSO analysis (Figure 1d). Below is the final prognostic risk model for STAD:

(1) $Riskscore=-0.1091\ast AC091057.1-0.0112\ast ADAMTS9.AS1+0.024\ast AC090825.1+0.0163\ast AC084880.3+0.0091\ast EMX2OS+0.0525\ast HHIP.AS1+0.035\ast AC016583.2+0.0891\ast EDIL3.DT+0.032\ast DIRC1+0.0552\ast LINC01614-0.0688\ast AC103702.2$(1)

3.2. Validation of the anoikis-related lncRNA risk prognostic model

Samples were split into high- and low-risk groups with the median risk score. From risk score distribution diagram (Figure 2a) and survival state diagram (Figure 2b), it can be seen that the mortality rate of patients gradually increased with the risk score. The heat map of gene expression (Figure 2c) showed a difference between high- and low-risk groups. In general, expression levels of AC084880.3, DIRC1, LINC01614, EMX2OS, AC016583.2, HHIP.AS1, ADAMTS9.AS1, AC090825.1, and EDIL3.DT were higher in high-risk group than in low-risk group, while AC091057.1 and AC103702.2 expression levels were lower in high-risk group. We plotted ROC curves with the training set to evaluate accuracy of the STAD prognostic model. The values of AUC for 1, 3, and 5 years were 0.62, 0.69, and 0.68, respectively (Figure 2d). Survival analysis displayed lower survival rate in high-risk group (p < 0.01) (Figure 2e).

Figure 2. Validation of the predictive performance of the STAD prognostic model based on anoikis-related lncRNAs.

Note: (a-b) Score distribution map and survival status distribution map of high- and low-risk groups; (c) heat map of 11 feature lncRNAs; (d-e) ROC curves and Kaplan–Meier survival curves of two groups from training set; (f-g) ROC curves and Kaplan–Meier survival curves of two groups from the validation set.

Predictive performance of the model was further validated by ROC and survival curves using the validation set. The AUC values were 0.65, 0.75, and 0.81 at 1, 3, and 5 years, respectively (Figure 2f). Survival curves of the validation set were similar to those of the training set, showing a remarkably lower survival rate in high-risk group (p < 0.01) (Figure 2g). It could be concluded that this established prognostic model had good predictive performance and high accuracy.

3.3. Analysis of independent prognostic factor and clinical significance

To further investigate whether the STAD prognostic model could be an independent prognostic factor, we performed univariate and multivariate Cox regression analyses combining with clinical information. The P values of risk score were less than 0.01, indicating that risk score could be an independent prognostic factor (Figure 3a-b). To predict survival in STAD patients, risk score was combined with other clinicopathological features to draw nomograms (Figure 3c). The results of calibration curves showed that predicted survival of the prognostic model fitted well with actual survival (Figure 3d-f). With high accuracy and fitting degree, the risk prognostic model constructed from these 11 signature lncRNAs can effectively assess the prognostic risk of STAD patients.

Figure 3. Analysis of the model as an independent prognostic factor.

Note: (a-b) Forest plot for univariate and multivariate Cox regression analyses of risk score combined with clinical information; (c) nomogram of risk score combined with clinical information; (d-f) calibration curves for 1, 3, and 5 years.

3.4. GO and KEGG enrichment analyses

Differential analysis was conducted in two groups to analyze the reasons for survival differences between high- and low-risk groups. After that, enrichment analyses of 1340 DEGs were completed using GO and KEGG to study their differences in biological functions as well as pathways. GO analysis uncovered that these genes were mainly enriched in functions such as muscle system processes, collagen-containing extracellular matrix, receptor–ligand activity, signaling receptor activator activity, and neuronal cell bodies (Figure 4a). KEGG results showed that these genes were mainly enriched in pathways such as PI3K-Akt signaling pathway, neuroactive ligand–receptor interaction, focal adhesions, cAMP signaling pathway, and calcium signaling pathway (Figure 4b). These enrichment results suggested that anoikis-related genes responsible for STAD risk score differences were mainly associated with neurotransmitter transmission, signal transmission, and endocytosis.

Figure 4. GO and KEGG enrichment analyses of DEGs between high- and low-risk groups.

Note: (a) GO enrichment analysis; (b) KEGG enrichment analysis.

3.5. Assessment of immune infiltration levels in different risk groups

We first scored stromal cell fraction and immune cell fraction of samples and performed the differential analysis of immune scores in high- and low-risk groups. It was found that stromal score, immune score, and ESTIMATE score were higher in high-risk group (p < 0.001) (Figure 5a). Results from the ssGSEA approach showed that macrophages, neutrophils, T helper cells, T follicular helper cells, and regulatory T cells (Treg) were significantly higher in high-risk group (p < 0.001) (Figure 5b). Differential analysis of immune function showed that CCR, inflammation-promoting, T cell co-inhibition, para-inflammation, and Type II IFN response were significantly higher in high-risk group than in low-risk group, which may correlate with tumor progression (p < 0.001) (Figure 5c). Based on the above-mentioned observations, prognostic risk score should be associated with the immune microenvironment.

Figure 5. Assessment of immune infiltration in the high- and low-risk groups.

Note: (a-c) Violin plots for differential immune score analysis, immune cell panel analysis, and differential immune function analysis in different risk groups, with blue representing the low-risk group and red representing the high-risk group.

3.6. Expression of immune checkpoint molecule and TIDE score results

Although antibodies against immune regulators CTLA4 and PD-L1/PD-1 are proven to be successful in a variety of cancers, only a subset of people present durable responses [23]. In this study, we further analyzed the response to immunotherapy in high- and low-risk groups by assessing expression of PDCD1 and CTLA4. The results found significantly higher expression of PDCD1 and CTLA4 in high-risk group (p < 0.01) (Figure 6a,b).

Figure 6. Expression of immune checkpoint molecule and TIDE scores for high- and low-risk groups.

Note: (a) Box plots of PDCD1 expression; (b) box plots of CTLA4 expression; (c) violin plots of TIDE scores.

According to the TIDE database, the TIDE scores of samples from different risk groups (Table S4) were calculated to predict the effect of immune checkpoint inhibition therapy. With lower TIDE scores in low-risk group and higher TIDE scores in high-risk group, patients in high-risk group were prone to immune escape (p < 0.001) (Figure 6C).

4. Discussion

Despite thorough research and substantial improvement in diagnosis and treatment of STAD, the morbidity and mortality rates remain high [24]. Several studies have demonstrated that anoikis exerts critical function in tumor invasion and metastasis [4]. Multiple mechanisms and TME contribute to anoikis resistance in cancer cells to counteract apoptosis mechanisms and maintain pro-survival signals to counteract allergic reactions after detachment from stromal surfaces. All these events are considered promising targets for antimetastatic pharmacological therapy. Therefore, the impact of anoikis-related lncRNAs on prognosis of STAD became the focus of this study.

Based on 11 feature lncRNAs related to the prognosis of STAD, namely, AC091057.1, ADAMTS9.AS1, AC090825.1, AC084880.3, EMX2OS, HHIP.AS1, AC016583.2, EDIL3.DT, DIRC1, LINC01614, and AC103702.2, screened from AnoRGs-related lncRNAs in this study, the prognosis model was constructed. Immune-related AC091057.1 has been used to construct the lung adenocarcinoma risk signature [25]. ADAMTS9.AS1 regulates colorectal cell migration and proliferation in vivo and in vitro and inhibits Wnt/β-catenin signaling pathway from participating in its mediated gene transcriptional regulation [26]. EMX2OS, an enhancer RNA from enhancer regions, bears some relation with poor survival outcomes in STAD [27]. The aberrant expression of HHIP-AS1, which is actively transcribed from a bidirectional promoter shared with the Sonic Hedgehog (SHH) regulator HHIP, is a hallmark of the SHH signaling pathway driving tumors [28]. Increased expression of DIRC1 in gastric cancer tissue may influence cancer cell proliferation, invasion, and migration through the AKT/mTOR signaling [29]. Since AC090825.1, AC084880.3, AC016583.2, EDIL3.DT, and AC103702.2 have been rarely reported in cancer studies, the molecular mechanism by which these lncRNAs are involved in cancer is unknown. In summary, most of these lncRNAs directly or indirectly affect tumor development and progression by regulating cell proliferation and migration or acting on signaling pathways.

Targeted immune checkpoint therapy is a novel treatment for GC, although the immune response is absent in some patients [30,31]. We assessed immune infiltration levels among STAD patients who did respond to immunotherapy. Immunosuppressive Treg cells can induce immune escape in TME, thereby promoting tumor metastasis and progression [32]. Significant differences in inflammation-promoting and part inflammation functions suggest that inflammation may greatly contribute to malignant tumor development and progression [33]. Enhanced enrichment of T cell co-inhibition function observed in high-risk group may inhibit T cell activation. Subsequently, we further investigated the response to immunotherapy in high- and low-risk groups. The difference in expression of PDCD1 and CTLA4 suggested that STAD tumor cells in high-risk group were more likely to evade attack by the host immune system, which was also confirmed by the TIDE score. Taken together, it can be seen that the prognostic risk of STAD constructed by anoikis-related factors can distinguish the immune level of patients and provide a reference for STAD patients considering immunotherapy in clinical practice.

Although the STAD prognostic model is constructed by standard procedures and had its performance proved by convincing evidence, unfortunately, there still exist some limitations. Firstly, the molecular function of some prognostic risk markers in STAD is unknown. Secondly, this study used only retrospective data to construct and validate the STAD prognostic risk model. Clinical observation should be included to further verify predictive ability of the prognostic model and evaluate its value in clinic.

Author contribution

Conceptualization: Lina Lu, Guofa Jiang

Data curation: Min Yu

Formal Analysis: Wei Huang

Acquisition: Wei Huang

Investigation: Hui Chen

Methodology: Hui Chen

Project administration: Guofa Jiang

Resources: Guofa Jiang

Software: Guofa Jiang

Supervision: Lina Lu

Validation: Gangxiu Li

Visualization: Gangxiu Li

Writing – original draft: Lina Lu, Min Yu

Writing – review & editing: All authors

Declaration of Conflicting Interests

The authors report no conflict of interest.

Ethics approval and consent to participate

Not applicable.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The datasets generated and analyzed during the current study are not publicly available but are available from the corresponding author on reasonable request.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19932820.2023.2220153.

Additional information

Funding

This study received no funding.