Identification and Validation of an Endoplasmic Reticulum Stress-Related lncRNA Signature for Colon Adenocarcinoma Patients

Introduction

Colon adenocarcinoma (COAD) is a malignant tumor and is the third-most lethal tumor of the digestive system. The 5-year survival rate of COAD is approximately 50%.^1,2 With the continuous development of diagnostic and treatment methods in recent years, the survival rate of COAD has improved, but the mechanism of COAD progression remains unclear. Because the detection rate is still not perfect, some COAD biomarkers are being used or tested in clinical trials.^3–5 There is a need to explore additional potentially informative biomarkers to improve COAD prognoses.

Endoplasmic reticulum (ER) is the organelle responsible for the synthesis, folding, modification, and control of secreted proteins and membrane proteins in eukaryotic cells.^6–9 Many disturbances can impair a cell’s ability to process protein folding and then lead to ER stress.^10,11 When an unfolded protein response cannot maintain the balance of ER stress, it can induce cell apoptosis and affect the growth, survival, differentiation, and maintenance of protein homeostasis in tumor cells.^12–16 Therefore, the utilization of ER stress-related genes may be an effective strategy for the treatment of COAD. Long non-coding RNAs (lncRNAs) have been generally defined as non-coding RNAs > 200 nucleotides in length with incomplete functional open reading frames.^17,18 They participate in the regulation of COAD-related cell signaling pathways.^19,20 Several studies indicate that lncRNAs are related to ER stress, and recent studies suggest that lncRNA DANCR protects cardiomyocytes from ER stress injury by sponge-adsorbed microRNA-6324.²¹ LncRNA MEG3 promotes ER stress through the MEG3/miR-103a-3p/PDHB pathway.²² We therefore hypothesized that the identification of ER stress-related lncRNAs may be useful for the prediction of COAD prognoses, which provide the basis for individualized treatment.

In this study, COAD RNA sequencing data downloaded from the Cancer Genome Atlas (TCGA) database were used as a training set, then data downloaded from the Gene Expression Omnibus (GEO) database were used as a validation set. A risk prognosis model including eight ER stress-related lncRNAs was established via Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses. COAD patients were then divided into low-risk and high-risk groups based on the median risk score, and survival analysis, tumor mutation burden (TMB) and copy number variation (CNV) analyses, and the gene set enrichment analysis (GSEA) analysis were conducted in the two groups. The predictive capacity of the model was verified in the GEO validation cohort. Lastly, lncRNA expression in the model was validated via RT-PCR using clinical tissue samples.

Materials and Methods

Data Sources

COAD RNA sequencing and clinical characteristics data were obtained from two databases. A training cohort dataset data from 434 patients was downloaded from the TCGA project (https://portal.gdc. cancer.gov/repository). The external validation cohort dataset containing data from 556 patients was derived from the GEO databases (GSE39582)(https://www.ncbi.nlm.nih.gov/geo/). RNA sequencing transcriptome data of HTseq-counts were transformed into TPM for subsequent analysis. Patients with unknown overall survival (OS) or incomplete information were excluded from both datasets.

Screening of ER Stress Gene-Associated lncRNAs

ER stress-related genes were obtained from the GSEA/MSigDB platform (https://www.gsea-msigdb.org/gsea/login.jsp). Three gene sets associated with ER stress were included²³ (GO_PERK_MEDIATED_UNFOLDED_PROTEIN_RESPON

ESPONSE, GO_IRE1_MEDIATED_UNFOLDED_PROTEIN_RESPON

ESPONSE, and REACTOME_ATF6_ATF6_ALPHA_ACTIVATES_CHAPE

CHAPERONE_GENES). The limma package was used to normalize the presentation data. Pearson’ s correlations were calculated to assess relationships between ER stress-related genes and lncRNAs (|R| >0.2).

Establishment of a Risk Signature

Univariate Cox regression analysis was used to identify ER stress-related lncRNAs associated with prognosis (p < 0.05). The lncRNAs identified were then used in LASSO analysis and multivariate Cox regression analysis to generate a set of prognostic characteristics associated with ER stress. Based on Akaike Information Criteria values, the best ER stress-related lncRNA prognostic signature was selected for further analysis. Each patient’s risk score was calculated via the following formula, where Coef(i) and x(i) respectively represent estimates of each regression coefficient and expression value of ER stress-related lncRNA: risk score = ∑Coef(i)×x(i). COAD patients in the TCGA and GEO databases were divided into high-risk and low-risk groups based on the median risk score.

TMB and CNV Analyses

Mutation data were downloaded from the cbioportal database, and the “maftools” R package was used to visualize somatic mutation data in the mutation annotation format. The TMB value was then calculated by adding the number of somatic mutations to the length of the exon. Finally, differences in TMB values between the high-risk group and the low-risk group were analyzed.

The GISTIC algorithm was used to analyze masked CNV data derived from different transcriptome clusters. CNV data associated with Y chromosome aberrations in germ cells were excluded from this analysis. Differences in CNV values between the high-risk group and the low-risk group were analyzed.

Survival Analysis

To construct survival curves for patients at high and low risk, Kaplan-Meier curves and Log rank tests were used to assess associations between risk model scores in high-risk and low-risk groups and patient survival. Because COADs are closely associated with genetic mutations, survival in high-risk and low-risk groups was also analyzed with respect to BRAF, KRAS, NRAS, MSI, and MSS mutational status. By using the “survivorROC” package in R language, time-dependent receiver operating characteristic (ROC) curves were generated to evaluate the predictive accuracy of prognostic features.

Patients and Tissue Samples

Altogether 6 COAD tissues and paired normal tissues were collected from patients who provided informed consent who had undergone radical surgical resection of COAD between September 2016 and January 2017 at the Third Xiangya Hospital of Central South University, China. All dissected samples were immediately stored in liquid nitrogen until preparation for total RNA extraction. The collection of tumor samples from patients with COAD was approved by the Ethics Committee of the Third Xiangya Hospital of Central South University, China (application number: 2017-S259).

Reverse Transcription and Quantitative Real-Time PCR

Total RNA was isolated from tissue using RNAiso Plus (TRIzol) (Takara, Dalian, China, Code No:9109) and reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (Takara, Dalian, China, Code No: RR047A) in accordance with the manufacturer’s instruction. The reaction was performed in a Light Cycler@ 480 II Sequence Detection System (Roche, Basel, Switzerland). Relative gene expression was calculated using the 2^-ΔCt method with glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as an internal control. All primers were synthesized by Sangon Biotech (Shanghai, China). The primers sequences were: GAPDH FP: 5’-TCAAGATCATTGCTCCTCCTGA-3’, RP: 5’-ACATCTGCTGGAAGGTGGACA-3’; MIR210HG FP: 5’-FGGTTCTGGCTTGCTGACAC-3’, RP: 5’-CAACTCGGCTTGGTTATTTC-3’; ALMS1-IT1 FP: 5’-AGGATACCTGGGATGGCTGT-3’, RP: 5’-ACAGCCATCCCAGGTATCCT-3’; LINC00638 FP: 5’-ACAATTCGACCCGTAACAGC-3’, RP: 5’-TGCTCGATATTCCCATGTCA-3’; ZEB1-AS1 FP: 5’-CCGTGGGCACTGCTGAA-3’, RP: 5’-CTGCTGGCA AGCGGAAC-3’; FLJ21408 FP: 5’-GCCCGCTCTCGATCCTTAAT-3’, RP:

5’-ACTACATTGTGCCAGCGGAA-3’; LINC00648 FP: 5’-CGCATTTTCACCCTCTTCGC-3’, RP:5’-CCTGAGAAGGCATAGCGACT-3’.

Statistical Analysis

R software (version 3.6.3) and GraphPad Prism v7.0 (GraphPad Software Inc.) were used as analysis tools. Expression levels of single lncRNAs in tumor tissues and normal tissues were compared via the paired t-test. P < 0.05 was deemed to indicate statistical significance.

Results

Identification of an ER Stress-Related lncRNA Signature

A total of 990 patients were included in the study, 434 derived from the TCGA, and 556 derived from the GEO. The clinicopathological characteristics of COAD patients are detailed in Table 1. A flow chart representing the entire study is shown in Figure 1. Pearson’s correlational analysis yielded 675 ER stress-related lncRNAs. In univariate Cox analysis, 109 ER stress-related lncRNAs exhibited prognostic value in COAD patients (p < 0.05). Eight ER stress-related lncRNAs were subsequently identified via LASSO regression. In multivariate Cox regression analysis, 8 lncRNAs were independent prognostic factors. Lastly, an 8-ER stress-related lncRNA signature was generated (Figure 2A–C).

Table 1 The Clinicopathological Characteristics of COAD Patients

Figure 1 Flowchart of data preparation, processing, analysis, and validation.

Figure 2 The construction of risk model and identification of ER stress-related lncRNAs. (A) The forest plot of univariate COX regression analysis among eight lncRNAs. (B) Partial likelihood deviation was plotted relative to the logarithm of lambda in 10-fold cross-validation. (C) The trajectory graph of each variable. (D and E) The expression heat map of eight lncRNAs and clinical features in TCGA (D) and GEO (E) cohorts. (F and G) The histogram of the expression of eight lncRNAs in the TCGA (F) and GEO (G) cohorts.

The median risk score was used as the cut-off for low-risk and high-risk COAD patients. Expression of lncRNAs in the training and validation cohorts is shown in Figures 2D and E. Heatmaps were generated representing gene expression profiles in the high-risk and low-risk groups in the TCGA-training cohort (Figure 2F) and the GEO-validation cohort (Figure 2G). In the TCGA-training cohort, the high-risk and low-risk groups differed with respect to TNM stage and pathological stage. Specifically, there were more advanced malignancies in the high-risk group (Table 2). Similarly, in the GEO-validation cohort, TNM stage and pathological stage were significantly higher in the high-risk group than in the low-risk group. In the TCGA cohort, there were no significant differences in tumor site or BRAF mutation between the high-risk and low-risk groups, but there were significant differences between the high-risk and low-risk groups in the GEO cohort (Table 3). ER stress-related lncRNAs were more likely to be expressed in the high-risk groups in both the training cohort and the validation cohort. Online database GEPIA analysis indicated that the differential gene expression identified was related to the prognosis, whereby high expression indicated a worse prognosis (Figures 3A–G).

Table 2 Demographic and Molecular Comparison Between HRisk and LRisk in TCGA

Table 3 Demographic and Molecular Comparison Between HRisk and LRisk in GSE38582

Figure 3 The prognostic analysis of a single lncRNA in GEPIA. (A–G) The survival curves of 7 lncRNAs.

Prognostic Effects of the lncRNA Signature in the Training and Validation Cohorts

Risk score was significantly associated with OS in patients with COAD. Kaplan-Meier curves indicated that prognoses were worse in high-risk group (Figure 4A). In ROC curve analysis, the areas under the curve (AUCs) corresponding to 1-year, 3-year, and 5-year survival were 0.7, 0.73, and 0.81. In the validation cohort, there were differences in OS between the high-risk and low-risk groups, The AUCs corresponding to 1-year, 3-years, and 5-years survival were 0.71, 0.61, and 0.6 (Figure 4B). The distribution patterns of risk scores and survival status are shown in Figures 4C and D. As risk score increased, the mortality rate gradually increased.

Figure 4 Analysis of prognostic survival of patients with COAD by risk model. (A and B) Kaplan-Meier curve and ROC curve drawn by dividing high- and low- risk groups based on the risk model in TCGA (A) and GEO (B) cohorts. (C and D) The risk score value of each sample, the survival status ranked from low to high-risk score in TCGA (C) and GEO (D) cohorts. (E) Nomogram to predict the 1-, 3-, 5-year OS for eight lncRNAs. (F) The calibration curves of the TCGA dataset. (G) The forest plot of multivariate COX regression analysis includes the risk model and clinical features.

A nomogram of the eight genes was constructed, predicting 1, 3, and 5-year survival (Figure 4E). A calibration curve indicated that the predicted and actual 5-year survival values were in good agreement, but the predicted and actual 1- and 3-year survival values were somewhat discrepant (Figure 4F). In multivariate Cox regression analysis performed based on the risk scores in the training cohort and validation cohort, risk score was a prognostic indicator (training group: p=1.77e-03, hazard ratio: 2.169; Validation group: p=5.57e-03, hazard ratio: 1.613) (Figure 4G). Because the prognostic model, tumor position stage, and BRAF/KRAS mutation types were independent indicators of OS in the two separate datasets, and NRAS mutation and MSI/MSS were also independent prognostic indicators in the TCGA training cohort, stratified analyses were conducted to investigate whether the prognostic model effectively predicted OS in COAD patients with different tumor position stages, BRAF/KRAS/NRAS mutation types, and MSI/MSS types. The results are shown in Figures 5A–I.

Figure 5 Analysis of prognostic survival of patients with COAD by risk model. (A and B) Kaplan-Meier curve based on the tumor location in TCGA (A) and GEO (B) cohort. (C–F) Kaplan-Meier curve based on BRAF and KRAS status in TCGA (C and D) and GEO (E and F) cohorts. NRAS status in TCGA (G). MSI status in TCGA (H). MSS status in TCGA (I).

COAD Tumor Mutation Profiles

Missense mutations were the most frequent type in the training cohort and the validation cohort. With regard to mutation types, the frequency of single nucleotide polymorphisms was higher than that of deletion and insertion mutations, and C > T was the major type of single nucleotide variant class. In the COAD samples, the top 10 mutated genes in the high-risk group were APC (71%), TP53 (60%), KRAS (44%), PIK3CA (32%), CSMD3 (22%), RYR2 (21%), ZFHX4 (21%), FBXW7 (19%), HTDIN (18%), BRAF (16%), in the low-risk group, they were APC (71%), TP53 (52%), KRAS (42%), PIK3CA (33%), RYR2 (21%), ZFHX4 (21%), SPTA1 (18%), DCSH2 (17%), DMD (17%), CSMD3 (16%) (Figures 6A and B). Details of mutations in the COAD samples in the high-risk and low-risk groups are represented as waterfall diagrams in Figures 6C and D.

Figure 6 The summary information of somatic mutation in colon cancer samples from the TCGA database. (A and B) Somatic mutation in high- and low-risk groups, including frequency of different mutation classification, frequency of mutational variant types, frequency of SNV class, TMB level of each colon cancer sample, and top 10 mutated genes in colon cancer samples. SNP, single nucleotide polymorphism; INS, insertion; DEL, deletion; SNV, single nucleotide variants. (C and D) The waterfall plot of the mutation profile of a single gene in each COAD sample.

Copy Number Profiles in Between Low-Risk and High-Risk Groups

Tthe GISTIC algorithm was used to analyze copy number data derived from 426 COAD patients downloaded from the database. Frequently changing areas were identified in the high-risk score and low-risk score groups’ genomes (Figures 7A and B). Many areas with significant copy gain and loss were detected. The copy number-altered genomes in the high-risk group were greater than those in the low-risk group because the copy number-lost genomes in the high-risk group were greater than those in the low-risk group (Figure 7C, p<0.05). These results suggested that lncRNA copy deletion may be associated with the occurrence and development of COAD.

Figure 7 Continued.

Figure 7 Gene copy number variations in COAD. (A and B) The amplification and deletion of gene fragments in COAD patients among two groups. (C) The comparison of copy number of the altered genome, lost genome, and gained genome between high- and low-risk groups.

Enrichment Analysis

To investigate the biological functions potentially associated with the lncRNA signature, GSEA was conducted in the high-risk and low-risk groups in the TCGA and GEO cohorts. Several tumor process pathways were enriched in the high-risk group, such as degradation of the extracellular matrix (ECM), ECM proteoglycans, ECM organization, and MET activation of PTK2 signaling (Figures 8A and B). The scores of these pathways in the samples are represented as heat maps in Figures 8C and D. Gene ontology and Kyoto Encyclopedia of Genes and Genomes analyses revealed enrichment mainly of genes involved in cell-ECM organization and extracellular structure organization (Figures 8E–H). Generally, ER stress triggers the activation of an adaptive response known as the Unfolded Protein Response. Studies indicate that the mechanism of the unfolded protein response in cell migration and invasion is closely related to ECM remodeling and cell adhesion modification.²⁴ Concordantly, the relatively highly expressed genes in the high-risk group were mainly involved in ECM degradation.

Figure 8 Continued.

Figure 8 Enrichment analysis. (A and B) GSEA analysis reveals pathways in the regulation of high-expressed genes in the high-risk group. (C and D) The heat maps of the scores of enrichment pathways in each sample. (E–H) GO analysis and KEGG analysis.

Confirmation of Expression of the Eight Genes

After generating the model based on tumor sample data analysis, differences between the eight lncRNAs identified were analyzed in normal samples and tumor samples. Expression of LINC00648 and ALMS1-IT1 were lower in tumor samples than in normal samples, whereas expression of FLJ21408, LINC00638, and MIR210HG were higher (Figures 9A–E).

Figure 9 Gene expression profiles in clinical samples. (A) FLJ21408 expression profiles in clinical samples. (B) LINC00638 expression profiles in clinical samples. (C) MIR210HG expression profiles in clinical samples. (D) LINC00648 expression profiles in clinical samples. (E) ALMS1-IT1 expression profiles in clinical samples. In clinical samples, the target gene expression was normalized to GAPDH (ΔCt), data are presented as 2−ΔCt.

Discussion

Despite advances in surgery and chemotherapy in recent years, the prognoses of patients with advanced and metastatic COAD remains unsatisfactory.^25,26 In addition, due to molecular heterogeneity, patients with the same TNM stage or similar risk factors may have different clinical outcomes.^27,28 Therefore we aimed to establish a risk-scoring model for the early diagnosis of COAD.

LASSO regression and Cox regression were utilized to identify eight ER stress-related lncRNAs, and a new signature that can be used as an independent prediction factor was constructed. By applying this risk prediction signature to the patients, the outcomes of patients with high-risk and low-risk scores can be predicted. Patients in the high-risk group had significantly worse outcomes than those in the low-risk group. These results suggested that the prognostic value of the ER stress-related lncRNA signature was reliable for survival prediction in COAD patients. In the analysis of TMB and CNVs there were differences between patients with high and low scores. The high-risk group had higher CNVs. Genomic CNVs contribute to various types of human dysfunction or tumorigenicity in many ways, suggesting that high-risk groups are at higher risk of tumors.

GSEA revealed potential signal pathways that the 8 ER stress-related lncRNAs were associated with. ECM-related processes were upregulated in the high-risk group. Previous studies indicate that in the course of tumor occurrence and development ECM undergoes a remodeling process.^29,30 The reconstructed ECM can create a favorable microenvironment for tumor growth, leading to tumor cell proliferation, invasion, and metastasis.^31,32 The relevant characteristics of ECM remodeling can be used as important indicators for tumor clinical staging, early diagnosis, and prognosis.

Eight ER stress-related lncRNAs (LINC01354, MIR210HG, HOTAIR, ALMS1-IT1, and ZEB1-AS1) have been reported to be associated with cancer. LINC01354 promotes the proliferation and invasion of lung cancer cells.³³ It is also closely related to the development of colorectal cancer. Studies have shown that LINC01354 promotes the proliferation and metastasis of colorectal cancer,³⁴ and promotes osteosarcoma cell invasion via upregulation of integrin β1.³⁵ MIR210HG is highly expressed in a variety of cancers and it can promote tumor growth in triple-negative breast cancer and non-small cell lung cancer,^36,37 as well as the proliferation and invasion of cervical cancer cells.³⁸ HOTAIR acts as an oncogene in human cancer and some studies indicate that HOTAIR enhances liver cancer cell proliferation via the promotion of epithelial-to-mesenchymal transition.³⁹ A recent study indicated an association between a genomic variant within HOTAIR and the risk of psoriasis.⁴⁰ ALMS1-IT1 reportedly has prognostic value, and some studies indicate that high ALMS1-IT1 expression is associated with a poor prognosis in head and neck squamous cell carcinoma.⁴¹ Subsequent studies have suggested that ZEB1-AS1 is highly expressed in glioma and bladder cancer, among others,^42,43 and that it is upregulated in hepatocellular carcinoma.⁴⁴ ZEB1-AS1 regulates ZEB1 expression, promoting tumor cell proliferation and invasion in hepatocellular carcinoma, leading to a poor prognosis. There are few reports on LINC00638, FLJ21408, and LINC00648, thus our future studies will focus on those lncRNAs.

The advantage of the current study is that we generated a prognostic indicator based on an 8-gene signature that predicted survival with relatively high AUCs in both a training cohort and a validation cohort. Nonetheless, the study had some limitations. Neither in vitro nor in vivo mechanism investigations were conducted to confirm the reliability of the signature analysis. Therefore, a large number of experiments will be conducted in the future to demonstrate the mechanistic association between these genes and the progression of COAD.

Conclusions

After a series of bioinformatics analyses and verifications, we constructed a prognostic model consisting of 8 ER stress-related lncRNAs with strong predictive value. To the best of our knowledge, this is the first report of the construction of an ER stress-related lncRNA prognostic model for malignancies. This study provides novel research targets for studying the pathogenesis and progression of COAD.