Identification of EMTlnc in PAAD Patients
Firstly, through the gene annotation file from the TCGA database, we identified 14805 lncRNAs in the TCGA dataset and 1440 lncRNAs in the ICGC dataset for the next analysis. Then we conducted GSVA analysis on EMT path-related genes in TCGA PAAD samples and obtained GSVA variation values. After applying WGCNA on the database of PAAD patients in TCGA, coexpression network by WGCNA analysis shown that the EMT related lncRNA model in PAAD samples were grouped into 48 models, contained MElightcyan and MEturquoise model with abs(cor)>0.5 and p.value <0.05, which were both included into further analysis in this research. The lncRNAs in those two models were defined as EMTlnc.
Combined with the prognostic information, univariate Cox regression was then implemented to screen EMT-related prognostic lncRNAs from the EMTlnc in both the TCGA and ICGC datasets (p < 0.05), respectively. Finally, we found that 33 EMTlnc were significantly correlated with the OS of PAAD patients in both two datasets through taking intersection. The work flow was shown in Figure 1A and the WGCNA in PAAD samples was shown in Figure 1B, the correlation between membership in the lightcyan module and membership in the EMT pathway of the eigengenes in lightcyan module by Pearson correlation analysis, and correlation between membership in the turquoise module and membership in the EMT pathway of the eigengenes in turquoise module by Pearson correlation analysis was shown in Figure 1C, Figure 1D respectively. The results of univariate Cox analysis of the 33 EMTlnc were shown in Table 1.
Table 1 | The thirty-three EMT-related prognostic lncRNAs
EMT related
lncRNAs
|
TCGA
|
ICGC
|
HR
|
HR.95L
|
HR.95H
|
p-value
|
HR
|
HR.95L
|
HR.95H
|
p-value
|
AC017002.1
|
2.9985
|
1.0379
|
8.6629
|
4.25E-02
|
1.4835
|
1.0744
|
2.0483
|
1.66E-02
|
AC093850.2
|
1.2967
|
1.0605
|
1.5855
|
1.13E-02
|
1.1653
|
1.0152
|
1.3376
|
2.96E-02
|
LINC00152
|
1.9583
|
1.3243
|
2.8959
|
7.60E-04
|
1.9065
|
1.2703
|
2.8615
|
1.84E-03
|
LINC01116
|
1.7473
|
1.1670
|
2.6161
|
6.73E-03
|
1.5378
|
1.1805
|
2.0033
|
1.42E-03
|
MIR4435-1HG
|
2.3464
|
1.4175
|
3.8839
|
9.10E-04
|
2.3326
|
1.4546
|
3.7404
|
4.39E-04
|
RP11-274H2.3
|
5.8172
|
1.4789
|
22.8821
|
1.17E-02
|
1.3529
|
1.0352
|
1.7681
|
2.69E-02
|
RP11-400N13.3
|
1.6652
|
1.3155
|
2.1079
|
2.24E-05
|
1.2304
|
1.0762
|
1.4068
|
2.41E-03
|
RP11-417E7.1
|
1.5019
|
1.0033
|
2.2482
|
4.82E-02
|
1.3728
|
1.1670
|
1.6149
|
1.31E-04
|
RP11-554I8.2
|
1.6514
|
1.2114
|
2.2512
|
1.51E-03
|
1.1728
|
1.0457
|
1.3153
|
6.46E-03
|
UCA1
|
1.4171
|
1.2289
|
1.6341
|
1.62E-06
|
1.3475
|
1.1798
|
1.5391
|
1.10E-05
|
AC009506.1
|
0.2309
|
0.0902
|
0.5911
|
2.24E-03
|
0.6498
|
0.4526
|
0.9330
|
1.95E-02
|
AC096772.6
|
0.4648
|
0.2784
|
0.7762
|
3.40E-03
|
0.5803
|
0.3373
|
0.9985
|
4.93E-02
|
DANCR
|
0.3418
|
0.2100
|
0.5562
|
1.56E-05
|
0.6598
|
0.4665
|
0.9331
|
1.87E-02
|
FLJ37035
|
0.0148
|
0.0007
|
0.3109
|
6.69E-03
|
0.7340
|
0.5878
|
0.9165
|
6.35E-03
|
GS1-358P8.4
|
0.3759
|
0.2246
|
0.6289
|
1.95E-04
|
0.3935
|
0.2300
|
0.6734
|
6.67E-04
|
HNF1A-AS1
|
0.5779
|
0.4157
|
0.8034
|
1.11E-03
|
0.7402
|
0.5912
|
0.9267
|
8.70E-03
|
LINC00261
|
0.6925
|
0.5461
|
0.8780
|
2.41E-03
|
0.8859
|
0.8096
|
0.9693
|
8.33E-03
|
LINC01128
|
0.2883
|
0.1461
|
0.5691
|
3.37E-04
|
0.3458
|
0.1789
|
0.6681
|
1.58E-03
|
PART1
|
0.0495
|
0.0026
|
0.9468
|
4.59E-02
|
0.8745
|
0.7691
|
0.9944
|
4.09E-02
|
PP7080
|
0.6783
|
0.4953
|
0.9288
|
1.55E-02
|
0.7559
|
0.5729
|
0.9974
|
4.79E-02
|
PRKAG2-AS1
|
0.6626
|
0.4526
|
0.9702
|
3.44E-02
|
0.7983
|
0.6602
|
0.9652
|
2.00E-02
|
RP11-16P6.1
|
0.0868
|
0.0184
|
0.4096
|
2.02E-03
|
0.6257
|
0.3945
|
0.9926
|
4.64E-02
|
RP11-226L15.5
|
0.4160
|
0.2031
|
0.8521
|
1.65E-02
|
0.5727
|
0.3337
|
0.9830
|
4.31E-02
|
RP11-244O19.1
|
0.1963
|
0.0790
|
0.4874
|
4.50E-04
|
0.7268
|
0.5435
|
0.9720
|
3.14E-02
|
RP11-384L8.1
|
0.5399
|
0.3382
|
0.8617
|
9.78E-03
|
0.6399
|
0.5037
|
0.8129
|
2.55E-04
|
RP11-700H6.1
|
0.0668
|
0.0087
|
0.5136
|
9.31E-03
|
0.8289
|
0.7284
|
0.9433
|
4.45E-03
|
RP1-193H18.2
|
0.3373
|
0.1930
|
0.5895
|
1.36E-04
|
0.7516
|
0.6052
|
0.9334
|
9.77E-03
|
RP11-968O1.5
|
0.5038
|
0.2815
|
0.9017
|
2.10E-02
|
0.7650
|
0.5933
|
0.9863
|
3.88E-02
|
RP5-1033H22.2
|
0.6204
|
0.3941
|
0.9766
|
3.92E-02
|
0.8018
|
0.6884
|
0.9338
|
4.50E-03
|
RP5-1085F17.3
|
0.3766
|
0.2402
|
0.5904
|
2.07E-05
|
0.5024
|
0.2875
|
0.8779
|
1.56E-02
|
RP5-894A10.2
|
0.5887
|
0.3842
|
0.9020
|
1.49E-02
|
0.6040
|
0.4265
|
0.8554
|
4.52E-03
|
SLC25A25-AS1
|
0.6439
|
0.4690
|
0.8839
|
6.45E-03
|
0.7420
|
0.5867
|
0.9383
|
1.27E-02
|
XXbac-B135H6.15
|
0.4622
|
0.2895
|
0.7377
|
1.22E-03
|
0.6998
|
0.5027
|
0.9740
|
3.44E-02
|
Color shaded lncRNAs were risky lncRNAs and others were protective lncRNAs.
Construction of the EMT-LPS in the TCGA Dataset and Validation of the EMT-LPS in the ICGC Dataset
To build the EMT-LPS for forecasting the OS of PAAD patients, we performed a LASSO Cox analysis on the basis of the 33 EMT-related prognostic lncRNAs in the TCGA cohort and it generated the EMT-LPS which contains 11 EMTlnc and coefficient of each (Figures 2A,B). The EMT-LPS involved 11 lncRNAs and, for each patient in the TCGA dataset, a risk score was calculated based on the coefficient for each lncRNA (Figure 2C). Patients in the TCGA cohort were divided into low and high-risk subgroups based on the median value of risk scores. Kaplan-Meier survival curves depicted that PAAD patients with higher risk scores had worse clinical outcomes (lower OS rates and a shorter OS time) (Figure 2D). Risk score and survival status distributions are plotted in Figure 2F. And the ROC curves demonstrated that EMT-LPS harbored a promising ability to predict OS in the TCGA cohort (1-year AUC = 0.81, 2-year AUC = 0.86, 3-year AUC = 0.9; Figure 2H).
To validate the prognostic ability of EMT-LPS, we calculated risk scores for patients in the ICGC cohort using the same formula. PAAD patients in the ICGC dataset were assigned to low- and high-risk groups based on the median risk score. The results were consistent with the findings in the TCGA dataset: PAAD patients with higher risk scores had lower OS rates and a shorter OS time in the ICGC dataset (Figure 2E). Risk score and survival status distributions are shown in Figure 2G and it showed that patients with higher risk scores had shorter overall survival time and dead status. The ROC analysis also indicated that EMT-LPS had a strong prognostic value for PAAD patients in the ICGC dataset (1-year AUC = 0.8, 2-year AUC = 0.85, 3-year AUC = 0.89; Figure 2I). These results showed that the EMT-LPS had a robust and stable OS-predictive ability.
Prognostic Analysis of the Eleven EMTlnc
Eleven EMTlnc were included in the EMT-LPS and univariate Cox regression analysis was used to evaluate their prognostic roles. The forest plot shows that PP7080, xxbac-B135H6.15, RP5-1085F17.3, DANCR, AC096772.6, LINC01128 are protective factors with HR (Hazard ratio) < 1, while LINC01116, UCA1, RP11-400N13.3 are risk factors with HR < 1 in EMT patients (Figure 3A). The heatmap (Figure 3B) shows that RP11-55418.2, RP11-400N13.3, UCA1, LINC00152 and LINC01116 expression increased with increasing risk score, whereas the expression of the AC096772.6, xxbac-B135H6.15, LINC01128, RP5-1085F17.3, PP7080 and DANCR decreased with increasing risk score. Their expression levels were also related to the clinicopathological features of PAAD, including N_stage, T_stage, gender, age and WHO grade (Figure 3B). The Kaplan-Meier survival curves confirmed that higher expression of AC096772.6, xxbac-B135H6.15, LINC01128, RP5-1085F17.3, PP7080 and lower expression of RP11-55418.2, RP11-400N13.3, UCA1, LINC00152, DANCR and LINC01116 were associated with better OS in the TCGA dataset (Figures 3C-M). In the ICGC dataset, the heatmap (Supplementary Figure S2C) also showed that RP11-55418.2, RP11-400N13.3, UCA1, LINC00152 and LINC01116 expression increased with increasing risk score, whereas the expression of the AC096772.6, xxbac-B135H6.15, LINC01128, RP5-1085F17.3, PP7080 and DANCR decreased with increasing risk score. Their expression levels were related to the clinicopathological features of PAAD, such as N_stage, T_stage, gender, and age.
Stratification Analysis of the EMT-LPS
We attempted to identify whether clinicopathological features were associated with the risk score. In TCGA dataset, the results revealed that PAAD patients with WHO grade III + Ⅳ had higher risk scores, while the risk score was not associated with age, gender, N_stage and T_stage (Figures 4A-E). To better assess the prognostic ability of the EMT-LPS, we performed a stratification analysis to confirm whether it retains its ability to predict OS in various subgroups. In contrast with patients with lower risk, higher risk PAAD patients had worse OS in age ≥ 60 (Figures 4F, G). Likewise, we confirmed that EMT-LPS retained its ability to predict OS for patients by female or male (Figures 4H,I) , patients with grade Ⅰ+Ⅱ or grade Ⅲ+Ⅳ, and patients with TNM stage N0, N1 or T3+T4(Figures 4J-N). In ICGC dataset, the results showed that PAAD patients with N1_stage, T3_stage and T4_stage had higher risk scores (Supplementary Figure S2D-G) . These data indicated that it could be a potential predictor for PAAD patients.
Principal Component Analysis
Based on the expression value of the 200 EMT-related genes, principal component analysis (PCA) was performed to assess the differences between the low- and high-risk subgroups (Supplementary Figure S2A,B). The results showed that the low and high-risk patients in both the TCGA and ICGC datasets were distributed in distinct directions. These results may suggest that differential EMT statuses exist in different risk subgroups.
Pathway and Process Enrichment Analysis and Gene Set Enrichment Analysis (GSEA)
For investigating the potential biological process and pathway involving in the molecular heterogeneity between the low- and high-risk subgroups, we identified 710 differential expression genes (DEGs) [|log2 (fold change)| > 1 and p < 0.05] between the low- and high-risk subgroups in the TCGA cohort. These DEGs were primarily enriched in these terms: chemical synaptic transmission, Neuronal System, regulation of membrane potential, behavior, plasma membrane bounded cell projection morphogenesis, GABAergic synapse(Figures 5A). Gene set enrichment analysis revealed that two tumor hallmarks were enriched in the high-risk subgroup, the interferon alpha response and the interferon gamma response (Figure 5B,C) . These results may give us some insights into the cellular biological effects related to the EMT-LPS.
EMT-LPS Was an Independent Prognostic Factor for PAAD Patients
We used univariate and multivariate Cox analyses to assess whether the EMT-LPS was an independent prognostic factor for patients with PAAD. Based on the data of PAAD patients in the TCGA dataset, univariate Cox analysis indicated that EMT-LPS was remarkably associated with OS [Hazard Ratio (HR): 3.832, 95% CI: 2.435-6.030, p < 0.001; Figure 6A] and multivariate Cox analysis further showed that EMT-LPS was an independent predictor of OS (HR: 3.573, 95% CI: 2.248-5.681, p < 0.001; Figure 6A). The conclusion was validated in the ICGC dataset, which confirmed that EMT-LPS was an independent predictor of OS for PAAD patients in the ICGC validation dataset (univariate: HR: 2.667, 95% CI: 1.923-3.697, p < 0.001; multivariate: HR: 2.787 95% CI: 1.945-3.94, p < 0.001; Figure 6B). These results indicated that our EMT-LPS, as an independent prognostic indicator, might be useful for clinical prognosis evaluation.
Construction and Validation of the EMT-LPS Based Nomogram
To create a clinically applicable quantitative tool to predict the OS of PAAD patients, we established a nomogram using the risk status (based on EMT-LPS), gender, age, T_stage, N_stage and WHO grade in the TCGA dataset and it was also tested in the ICGC dataset (Figure 6C). Calibration plots showed that the observed vs. predicted rates of 1-, 2- and 3-year OS showed perfect concordance in the TCGA (Figures 7A-C) and ICGC cohorts (Supplementary Figure S1B-D). Then time-dependent ROC curves were used to assess the prognostic predictive ability of the nomogram and other predictors (risk score, gender, age, T_stage, N_stage and WHO grade) in the TCGA (Figures 7D-F) and (Supplementary Figure S1E-G) and the results revealed that, compared with the other predictors, the nomogram had excellent accuracy regarding 1-, 2- and 3- year OS (AUC = 0.79, 0.83, and 0.86 in the TCGA, and 0.8, 0.84, and 0.89 in the ICGC , respectively). These data indicated that the nomogram has a robust and stable ability to predictive the OS for PAAD patients.
Construction of the ceRNA Network and Functional Enrichment Analysis
To further understand how the EMTlnc regulate mRNA expression by sponging miRNAs in PAAD, we constructed a ceRNA network based on the EMTlnc. Two of twelve lncRNAs were extracted from the miRcode database and thirteen pairs of interaction between the two lncRNAs and thirty-four miRNAs were identified. Then we used three databases (miRTarBase, miRDB, and TargetScan) to search target mRNAs based on the thirty-four miRNAs and totally 1539 mRNAs were identified in all the three databases. What’s more, these target mRNA were intersected with DEGs to botain differential expressed target mRNA. Ultimately, two lncRNAs, twelve miRNAs and thirteen mRNAs were included in our ceRNA network (Figure 8A). Furthermore, the 1539 target mRNAs were used to implemented functional enrichment analysis in the Metascape online tool and we found that these genes were enriched in vasculature development, pathway in cancer, regulation of cellular response to stress, Wnt signaling pathway, tissue morphogenesis, insulin signaling, regulation of cellular protein localization, response to growth factor, nagative regulation of cell differentiation (Figures 8B-D). These data may provide us some clues for finding the potential functions of these EMTlnc in PAADs.