Construction of a Competitive Endogenous RNA Network for Pancreatic Adenocarcinoma Based on Weighted Gene Co-expression Network Analysis and a Prognosis Model

Abstract

Pancreatic adenocarcinoma (PAAD) is a pancreatic disease with considerable mortality worldwide. Because of a lack of obvious symptoms at the early stage, most PAAD patients are diagnosed at the terminal stage and prognosis is usually poor. In this study, we firstly obtained RNA sequencing data of 181 patients with PAAD from The Cancer Genome Atlas (TCGA) database to identify early diagnostic biomarkers for PAAD. Survival-related mRNAs were identified using a weighted gene co-expression network analysis (WGCNA), and then a linear prognostic model of seven long non-coding RNAs (lncRNAs) was established using univariate and multivariate Cox proportional hazards regression analyses, which is verified using a time-dependent receiver operating characteristic (ROC) curve analysis. Finally, according to the survival analysis, we constructed a survival-related competing endogenous RNA (ceRNA) network. Our results showed that: (1) The upregulated genes related to cell cycle-related pathway (including homologous recombination, DNA replication and mismatch repair) in PAAD can increase the proliferation ability of cancer cells; (2) The 7-lncRNA signature can predict the overall survival (OS) of PAAD patients; and (3) The key mRNAs and lncRNAs are involved in mutual regulation in the ceRNA network.

Keywords: competing endogenous RNA network; pancreatic adenocarcinoma; prognostic signature; the cancer genome atlas; weighted gene co-expression network analysis.

Figure 1

Description of DEmRNAs, DEmiRNAs, and DElncRNAs. (A) Volcano plot showing DEmRNAs (green), DEmiRNAs (yellow), and DElncRNAs (red), The X axis represents logFC, the Y axis represents log10 (P value). Genes with |logFC| ≥ 1 and P ≤ 0.05 were defined as DEGs. Heatmap showing the normalized expression of (B) DEmRNAs, (C) DEmiRNAs and (D) DElncRNAs. T indicates cancer tissue, N indicates paracancerous tissue.

Figure 2

Weighted gene co-expression network analysis (WGCNA) network. (A) Cluster dendrogram of DEmRNAs. Each branch represents a single gene. Height indicates the Euclidean distance. Each color indicates a single module. (B) Hierarchical clustering tree of the TCGA-PAAD samples. Dendrogram tips are labeled with the TCGA-PAAD unique name. Height indicates the Euclidean distance. (C) Heatmap showing the Pearson correlation between modules and patients' clinical information. The numbers represent correlation coefficients and P-value. (D) Correlation analysis showed that genes in green module were significantly correlated with OS (cor = 0.76) (E) KEGG enrichment analysis of genes in the green module, ranked by P-value. The size of each dot represents the number of genes enriched in the pathway.

Figure 3

Four key genes in the green module. (A) Expression of four genes by grade in PAAD patients (due to the lack of expression data for G4, only the expression for G1–G3 is shown). P-value, 95% CI, and other information are shown at the top of the figure. (B) Relationships between the four genes and the cancer status among PAAD patients. (C) Survival analysis of the four genes. Patients were divided into high-expression and low-expression groups according to the median gene expression.

Figure 4

Kaplan–Meier survival curves of the training set, test set and total set. Survival curves of (A) training set (B), test set and (C) total set. Patients were divided into high- and low-risk groups according to the median prognostic score. log-rank, Wilcoxon, Fleming-Harrington, Tharone-Ware, and Peto-Peto was used to compare the differences of survival curves between two groups and get consistent results (P < 0.0001). The P-value of plots based on log-rank test.

Figure 5

LncRNA prognostic signature model and least absolute shrinkage and selection operator (Lasso) regression results. (A) Patients were divided into high- and low-risk groups according to the median prognostic score. The prognostic score, OS and expression levels of 7 lncRNAs in the two groups are shown. (B) Regression coefficient diagram using Lasso regression. (C) ROC curves of the 7-lncRNA prognostic signature for predicting 1-, 3-, and 5-year survival. (D) Diagnostic efficiency, based on the AUC, of 7-lncRNA signature models constructed using different machine learning models (random forest, support-vector machine [SVM] and XGBoost). (E) Boxplot showing the expression levels of 7 lncRNAs in the high- and low-risk groups. **p < 0.01, ****p < 0.0001.

Figure 6

CeRNA network related to PAAD overall survival (OS). (A) Kaplan–Meier curves of 3 mRNAs (TOB1, GOLGA8B and ANLN),3 miRNAs (203a-3p, 424-5p, and 135b-5p) and 3 lncRNAs (AP004608.1, AC005674.2, and AL365277.1) were selected for display. (B) CeRNA network related to PAAD OS.