Association of a novel seven-gene expression signature with the disease prognosis in colon cancer patients

Differential analysis

The mRNA, lncRNA, and miRNA expression profiles obtained were subjected to differential expression analysis. The screening thresholds for mRNA and miRNA were adjusted P value < 0.01 and |logFC| > 2, and for lncRNA was adjusted P value < 0.05 and |logFC| > 1. A total of 1476 (323 upregulated and 1153 downregulated) DEMs, 9 (5 upregulated and 4 downregulated) DELs, and 243 (113 upregulated and 130 downregulated) DEMis were identified. Top10 genes (in ascending order according to the adjusted P value) were selected from DEMs, DELs, and DEMis to draw the heatmaps, as shown in Figure 1.

GO enrichment analysis of DEMs

With a threshold of adjusted P value < 0.01, 13 BPs, 14 CCs, and 5 MFs were significantly enriched by the DEMs. The top 5 terms for the three categories are shown in Figure 2A. The seven most significant terms are shown in detail. According to the genes in term, the clustering of genes was carried out to construct the phylogenetic tree (Figure 2B). In terms of the involved BP terms, the DEMs were significantly enriched in chemical synaptic transmission, axon guidance, nervous system development, cell adhesion, and transport. From the perspective of CC terms, the encoding proteins were mainly located in an integral component of the plasma membrane, plasma membrane, and apical plasma membrane.

Pathway enrichment analysis

According to the KEGG pathway enrichment analysis of the DEMs obtained above by GSEA, a total of 67 significantly enriched pathways were identified (adjusted P value < 0.05), among which 16 were upregulated pathways (normalized enrichment score (NES) > 0) and 51 were downregulated pathways (NES < 0). KEGG ribosome biogenesis in eukaryotes had the largest NES in the significantly upregulated pathways (Figure 2C). In addition, DNA replication, RNA degradation, IL17 signaling pathway, p53 signaling pathway, and oxytocin signaling pathway were also significantly upregulated. Among the significantly downregulated pathways, drug metabolism cytochrome P450 had the largest NES (Figure 2D). Additionally, chemical carcinogenesis and mineral absorption pathways were also downregulated.

Co-expression analysis

For the screened DELs and DEMis, the correlation coefficients and statistical significance degree between them and DEMs were calculated, respectively, and the gene pairs with |r| > 0.6 and an adjusted P value < 0.01 were screened out. There were 193 nodes (199 mRNAs and 6 lncRNAs) and 229 edges in the co-expression network of mRNA-lncRNA, where lncRNA Pvt1 oncogene (PVT1) had the highest connectivity degree (Figure 3A). In the co-expression network of mRNA-miRNA, there were 184 nodes (392 mRNAs, 102 miRNAs) and 1,935 edges (Figure 3B).

miRNA prediction

A total of 200 lncRNA-miRNA pairs (3 lncRNAs and 141 miRNAs) were predicted, based on the lncRNA involved in lncRNA-mRNA co-expression pairs using starbase (Figure 4A). Additionally, the mRNA involved in lncRNA-mRNA co-expression pairs were subjected to mRNA-miRNA prediction, and 117,570 pairs of miRNA-mRNA (199 mRNAs and 2498 miRNAs) were obtained with Score = 1 and Position = CDS. These mRNA-miRNA pairs were intersected with the obtained mRNA-miRNA co-expression pairs, and 89 mRNA-miRNA pairs (36 mRNAs and 38 miRNAs) were further screened (Figure 4B).

LncRNA and miRNA pathway enrichment analyses

DKFZp779M0652 and PVT1 that had many target genes were subjected to functional enrichment. DKFZp779M0652 was associated with functions of cellular response to tumor necrosis factor and cellular response to interleukin-1. PVT1 was related to one-carbon metabolic process and reciprocal meiotic recombination.

The mRNAs in the miRNA-mRNA relation pairs obtained above were used as the target genes of miRNAs. Since there were few target genes and enrichment analysis could not be carried out, thus the functions of each target gene were inquired on the genecards to obtain the functions of all miRNA (data not shown).

ceRNA network construction

A total of 8,353 lncRNA-miRNA-mRNA pairs were obtained, involving 197 mRNAs, 135 miRNAs, and 3 lncRNAs. The constructed ceRNA network is shown in Figure 5A. After further screening of the network, only the aforementioned co-expressed mRNA-miRNA pairs were retained in the mRNAs-miRNA relation pairs. Finally, a total of 21 lncRNA-miRNA-mRNA pairs were obtained, including 7 miRNAs, 2 lncRNAs, and 10 mRNAs. The network is depicted in Figure 5B.

Establishment of mRNA prognostic risk model

A total of 277 samples (68 deceased and 209 living) with clinical information of overall survival time were selected. 197 mRNAs from ceRNA network were used as candidate mRNAs. The expression level of mRNA in the sample set was screened, and univariable Cox regression analysis was conducted based on the clinical prognosis information of the sample. The regression threshold of P value was set at 0.05 to screen the prognostic mRNAs, as shown in Table 1. The risk score models of prognostic mRNAs are shown in Table 2. The drawn survival curves are shown in Figure 6A. All models were significantly correlated, and the ROC curves were drawn based on the risk scores of all the models and survival information (Figure 6B). AUC values were calculated (Table 3), and the time periods were 1, 3, and 5 years, respectively. The largest AUC of 5-year prognosis scoring model was used as the final risk score model: Risk score = 0.58343 * exp (SGCG) + (-0.11860) * exp (CLDN23) + (-0.09726) * exp (SLC4A4) + 0.18416 * exp (CCDC78) + 0.13586 * exp (SLC17A7) + 0.40269 * exp (OTOP3) + (-0.23459) * exp (SMPDL3A). The parameters of mRNAs (SGCG, CLDN23, SLC4A4, CCDC78, SLC17A7, OTOP3, and SMPDL3A) associated with the optimal prognosis are shown in Table 4.

Cox regression analysis of prognostic clinical factors

The clinical factors of gender, age at initial pathologic diagnosis anatomic neoplasm subdivision, and pathologic stage were selected for univariable Cox regression analysis combining with the survival information, and the influence factor of pathologic stage (P < 0.0001) was selected. The samples without pathologic stage were removed from 277 samples in the previous step, and the remaining 268 samples were used for subsequent analysis. Pathologic stage contained the following types: “Stage I”, “Stage IA”, “Stage II”, “Stage IIA”, “Stage IIB”, “Stage IIC”, “Stage III”, “Stage IIIA”, “Stage IIIB”, “Stage IIIC”, “Stage IV”, “Stage IVA”, and “Stage IVB”. These types were converted into stage “1”, “2”, “3”, and “4”, and were denoted as “stage n” which was added into the model above for multivariable Cox regression analysis (p = 0.00026) to obtain the new scoring model as follows: Risk score = 0.51396 * exp (SGCG) + (-0.16881) * exp (CLDN23) + (-0.08028) * exp (SLC4A4) + 0.11820 * exp (CCDC78) + 0.09786 * exp (SLC17A7) + 0.44172 * exp (OTOP3) + (-0.22250) * exp (SMPDL3A) + 0.66818 * stage_n. The parameters of scoring model are shown in Table 5. The survival curve and ROC curve of survival rate at 1, 3, and 5 years are shown in Figure 7.

Model validation using independent datasets

Risk score values of all the samples in two independent validation datasets were calculated based on the above risk score model. All samples were divided into high risk group and low risk group according to the median score of risk score. Kaplan-Meier survival curve analysis showed that the developed prognostic signature could well distinguish high-risk and low-risk groups; the survival time of patients in low-risk group was significantly longer than that of patients in high-risk group (P = 0.00083 and P < 0.0001) (Figure 8A and 8B). The AUCs of this prognostic signature in predicting 1/3/5 survival rate of patients in GSE17538 were 0.755, 0.799 and 0.758, respectively (Figure 8C).

Public data processing

RNA-seq expression profile, miRNA expression data, and clinical phenotype information of all the samples of TCGA colon adenocarcinoma were downloaded from the University of California Santa Cruz (UCSC) Xene database [16] on May 27^th, 2019. RNA-seq data from 329 samples, including 288 tumor tissues and 41 normal tissues and miRNA expression data from 261 samples, including 253 tumor tissues and 8 normal tissues were downloaded. Besides, clinical information of 551 patients, including age, gender, tumor histological grading, survival time and survival status, were also downloaded.

LncRNA/mRNA re-annotation

GENCODE database is a scientific project in genome research which was used to identify and map all protein-coding genes within the ENCODE regions [24]. The gtf gene annotation file (Release 26, GRCh38.p10) provided by GENCODE database was downloaded, and the downloaded RNA-seq expression data were re-annotated with mRNA and lncRNA to obtain the mRNA expression profile and lncRNA expression profile, respectively. The gene with annotation information of "protein coding" was retained as mRNA, and with annotation information of "antisense", "sense intronic", "lincRNA", "sense overlapping", "processed transcript", "3prime overlapping ncRNA", or "non-coding" was retained as lncRNA. The clinical phenotype information corresponding to the gene expression profile of the samples were screened.

Differential analysis

The log2(count+1) data were transformed into raw count data and normalized into distribution of same mean value and equal variance for each sample using betaqn function in R package. The differences between tumor and normal samples in the three expression profiles were analyzed using the Bayesian method in limma package (version 3.40.0) [12]. Significance test was performed using the paired t-test to obtain the P values of all the genes. The P values were then adjusted for multiple test using Benjamini & Hochberg (BH) method to obtain the adjusted P values. Log2 fold change (FC) and adjusted P values were used to select the differentially expressed mRNAs (DEMs), and the final DEMs, DELs, and DEMis were screened out.

Function and pathway enrichment analyses of DEMs

Using the enrichment analysis tool DAVID [25] (version 6.8, https://david.ncifcrf.gov/), Gene Ontology (GO) (biological process (BP), cellular component (CC) and molecular function (MF)) analysis was performed, and the results were visualized using GOplot (version 1.0.2) [26]. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted based on the KEGG database [27] using the gene set enrichment analysis (GSEA, version 3.0) [28].

Co-expression analysis

For DELs and DEMis, the Pearson correlation coefficients between them and with the DEMs were calculated, and correlation tests were performed using the corr.test method in R package psych [29] (ci = F, adjust = "BH"). Multiple test was performed using the BH method. The co-expression pairs were obtained according to the correlation coefficient and significance degree, and the co-expression network was constructed using Cytoscape software (version 3.7.1) [30].

Online prediction of miRNA

Using starbase (version 3.0, http://starbase.sysu.edu.cn/) [31], the lncRNAs involved in the obtained lncRNA-mRNA co-expression pairs were subjected to lncRNA-miRNA prediction, thus acquiring the lncRNA-miRNA relation pairs. Additionally, the mRNA-miRNA relation pairs were also predicted based on the mRNAs involved in the obtained lncRNA-mRNA co-expression pairs, using the online tool mirwalk (version 3.0) [32]. Moreover, the intersection of these mRNA-miRNA relation pairs with the above mRNA-miRNA co-expression pairs was screened to obtain the final mRNA-miRNA pairs.

LncRNA and miRNA pathway enrichment analyses

The mRNAs that had co-expression relations with lncRNAs were considered as the target genes of lncRNAs. The lncRNAs with the maximum target genes were subjected to pathway and function enrichment analyses using DAVID [25], and for the lncRNAs with fewer target genes, their functions were searched in genecards (https://www.genecards.org/) [33]. Furthermore, the mRNAs in the mRNA-miRNA pairs were used as the target genes of miRNAs, and the related functions of each miRNA were analyzed in the same way as described above.

ceRNA relation integration and network construction

From the lncRNA-mRNA pairs with co-expression relations, we selected the lncRNA-mRNA pairs regulated by the same miRNA, which were then integrated with miRNAs to construct the ceRNA network. Cytoscape software (version 3.7.1) [30] was used for network construction.

Establishment of mRNA prognostic risk model

The mRNAs in the ceRNA network were used as candidate mRNAs, and the univariable Cox regression analysis in R survival package (version 2.44-1.1) [34] was used to analyze the regression coefficient and P value of each candidate mRNA in relation to survival time and status. The mRNAs with P value < 0.05 were initially considered as the mRNAs related to prognosis.

To further screen the prognostic mRNAs, the calculation model of risk score was defined as follows: Risk score = β gene1*expr (gene 1) + β gene2 * expr (gene 2) + ... + β genen * expr (gene n), Where, β is the prognostic correlation coefficient beta estimated by Cox analysis which equals to log (Hazard Ratio), and expr represents the expression value of corresponding gene.

Thereafter, based on the P values of mRNAs in Cox regression analysis ranking in ascending order, the mRNA with the minimum P value was used as the starting point, followed by adding the other mRNAs (Table 2). The samples were divided into high-risk and low-risk groups according to the median value of risk core. Survival analysis for high-risk and low-risk groups was performed by log-rank test after the addition of a certain mRNA. Then, the receiver operation characteristic (ROC) curves of 1, 3, and 5 years were calculated using the R survival ROC package (version 1.0.3) [35]. The risk score model with the highest Area Under the Curve (AUC) was taken as the best scoring model.

Cox regression analysis of prognostic clinical factors

Important clinical factors (gender, age at initial pathologic diagnosis, anatomic neoplasm subdivision, and pathologic stage) were selected from the phenotypic information corresponding to the samples. Thereafter, an univariable Cox regression analysis was performed by combining the survival information, and regression coefficient, and statistical P value between each clinical factor, and survival time and state were calculated. The influencing factors with threshold of P < 0.05 were further subjected to multivariable Cox regression analysis to obtain the final risk score model.

Model validation using independent datasets

The prognostic signatures developed were further validated in two independent datasets. The first dataset was the RNA-seq expression profile of colon cancer from the UCSC Xene database. This dataset included HTSeq-FPKM data of 512 samples, clinical phenotype information of 570 samples and survival data of 546 samples. The samples in the training dataset or samples without corresponding survival information or tumor stage information were excluded. At last, 164 samples were included for the validation.

The second dataset was GSE17538 downloaded from Gene Expression Ominibus (GEO), which included 213 colon cancer samples with complete information of survival time and survival status [36, 37].

Risk score was calculated for each individual, and the samples in validation datasets were divided into high-risk group and low-risk group according to the median score of risk score. Kaplan-Meier survival curve analysis was then performed to calculate the difference in the survival prognosis time between the samples of the high and low risk groups.