Establishment of Differentially-Expressed ARGs.
The flowchart of our study process is displayed in Fig. 1a. Initially, the mRNA expression profiles of four primary pancreatic ductal adenocarcinoma tissues (the primary group) and five liver metastasis tissues (the metastasis group) were extracted from GSE19279, setting p. adjusted < 0.05 and |log2FC| > 1 as the cut-off, identifying 1113 differentially-Expressed mRNAs (DEmRNAs). Then, we downloaded 243 human autophagy-related genes from the HADb. Finally, twenty-two differentially-Expressed ARGs were obtained by taking the intersection of DEmRNAs with the 243 ARGs (Fig. 1b), containing 4 downregulated and 18 upregulated genes (Fig. 1c). DEARGs associated with survival prognosis were further selected.
Construction Of The 5-args Prognostic Model
To confirm the relationship between DEARGs and OS outcomes, the Kaplan-Meier method and univariate proportional hazards Cox regression analysis were enforced and identified 7 DEARGs that were notably associated with patient prognosis (Fig. 2a-h). Next, to identify stable DEARGs, LASSO Cox regression analysis was used in the TCGA-PAAD cohort. Five autophagy-related prognostic genes were built when the model arrived at the minimum lambda value (Figs. 2i-j). Five DEARGs expressions in TCGA-STAD, TCGA-CHOL, TCGA-COAD, and TCGA-LIHC were verified (Supplementary Fig. 1). Then, using real time-PCR to exam the mRNA expression of CASP1, CHMP2B, MYC,RHEB and HDAC6 in pancreatic normal cell line (HPDE6-C7), five pancreatic cancer cell lines (Capan-1, CFPAC-1, PANC-1, MIA PaCa-2, COLO357) and five pancreatic cancer and matched adjacent normal tissues. The reult showed that the mRNA expression of CASP1, CHMP2B, MYC and RHEB was high, whereas that of HDAC6 was low in cancer cell lines and tissues (Supplement Fig. 2). The results of the two verification methods are consistent with that in GEO database.Eventually, five independent prognostic ARGs (CASP1, RHEB, CHMP2B, MYC, HDAC6) were screened as hub genes and were involved in structuring the autophagy-related prognostic model. The Lasso Cox regression analysis coefficient was applied as a weighting divisor, and the risk score was calculated using the following equation: 0.1919*RHEB expression quantity + 0.4924*CHMP2B expression quantity + 0.0096*MYC expression quantity + 0.1333*CASP1 expression quantity + (-0.6840)*HDAC6 expression quantity. Then, a risk score was computed for every patient and selected the median as the cut-off value. It divided the patients into low- and high-risk groups and observed the difference in survival between them to determine whether the model could effectively predict pancreatic cancer patients. By computing the risk score, clinicians could assess the overall survival probability predicted for each patient preliminarily.
Prognostic Ability Of The 5-args Prognostic Model And Application In Clinicopathological Stratification
Ranking the risk scores of the low-risk and high-risk groups, the distribution diagrams of the risk scores, survival statuses, and heatmaps of 5 ARGs were then plotted (Figs. 3a-c). The results revealed that the higher the risk scores, the more probably the patient was to die. Expressions of MYC, CHMP2B, RHEB, and CASP1 were lower, and that of HDAC6 was higher in the low-risk group. Kaplan-Meier survival curves indicated that the OS was significantly higher in the low-risk group than in the high-risk group (Fig. 3d). The risk scores value for the prediction (AUCs) for 1-year, 3-year, and 5-year OS were 0.68, 0.69, and 0.70, separately (Fig. 3e). Besides, the stratification based on risk score could efficiently distinguish pancreatic cancer patients with favorable or unfavorable prognoses. At different ages (≤ 65-year-old vs. > 65-year-old), gender (female vs. male), stage (I vs. II), and T staging (T 1 + 2 vs. T 3 + 4), patients with low risk scores had significantly favorable prognoses (Additional File 2). Therefore, the 5-ARGs prognostic model correlated with the survival prognosis of pancreatic cancer patients, and low risk score indicated a favorable prognosis.
Prognostic Value Of The 5-args Model And Clinicopathological Features
To evaluate prediction efficiency between clinicopathological features and the risk score, the univariate Cox regression analysis was applied to recognize prognosis-related factors. Factors in the univariate analysis with p < 0.05 were considered as significant factors and further involved in the multivariate Cox regression analysis. After above analysis, the forest plots showed that only the risk score (HR = 2.385, 95% CI = 1.549–3.672, p < 0.001) and stage (HR = 2.278, 95% CI = 1.238–4.226, p = 0.008) were the signficant independent risk factors (Figs. 4a). ROC curves were plotted for the pathological characteristics and risk score, displayed in Fig. 4b (The computation of T and M staging were neglected due to too few data). The AUCs of the pathological indexes were much lower than those of the risk score in pancreatic cancer patients (AUC = 0.692, 1 year). Decision curve analysis (DCA) was worked and further suggested that compared with other indicators, the risk score was a better prognostic signature in clinical decision-making at 1-, 3-, and 5-year (Additional File 3). Moreover, calibration curves were applied to assess the prediction veracity of the 5-ARGs prognostic model (Fig. 4c). The calibration curves revealed that at 1-, 2-, and 3-year, the prognostic model’s predicted consequences were closed to the criteria curve, which suggested that the 5-ARGs model had benefit prognostic ability.
Functional Enrichment Analysis And Screening Small-molecule Drugs
Further, GO and KEGG pathway analyses of five hub genes were carried out. The results indicated that the top three enriched BPs (Biological Processes) were macroautophagy, autophagy, and process utilizing autophagic mechanism; the top three enriched CCs (Cellular Components) were late endosome, multivesicular body, and ESCRT III complex; the top three enriched MFs (Molecular Function) were activating transcription factor binding, cysteine-type endopeptidase regulator activity involved in apoptotic process and peptidase activator activity involved in apoptotic process. At the same time, the top three enriched KEGG pathways were necroptosis, cellular senescence, and thyroid hormone signaling pathway (Fig. 5a). The results indicated that hub genes might participate in carcinogenesis and metastasis by regulating autophagy.
To identify potential drugs for pancreatic cancer, five hub genes (one down-regulated gene and four up-regulated genes) were uploaded to the L1000FWD website and matched small-molecule chemotherapies. The top ten predicted signficant small-molecule compounds and their similarity scores are listed in Table 1. The top three negative fractions drugs were apicidin, simvastatin, and PD-0325901, which were expected to inverted regulate the expression of five hub genes. Their 3D and 2D constructions were exteriorized via the PubChem website (Additional File 4). These potential small-molecule compounds might invert gene expression by inducting autophagy, guiding the advance of targeted therapies for pancreatic cancer. The therapeutic value of these candidate small-molecule compounds in pancreatic cancer treatment needs to be further verified.
Table 1 The screened drugs for PAAD treatment
drug
|
similarity score
|
p-value
|
q-value
|
Z-score
|
combined score
|
MOA
|
Predicted MOA
|
apicidin
|
-0.6
|
6.90E-05
|
2.74E-01
|
1.59
|
-6.6
|
Unknown
|
HDAC inhibitor
|
simvastatin
|
-0.6
|
5.36E-05
|
2.74E-01
|
1.63
|
-6.98
|
HMGCR inhibitor
|
glucocorticoid receptor
|
PD-0325901
|
-0.6
|
4.55E-05
|
2.74E-01
|
1.76
|
-7.63
|
MEK inhibitor
|
MEK inhibitor
|
auranofin
|
-0.6
|
5.11E-05
|
2.74E-01
|
1.81
|
-7.78
|
NFkB panthway inhibitor
|
calcium channel
|
panobinostat
|
-0.4
|
3.69E-03
|
2.76E-01
|
1.55
|
-3.77
|
HDAC inhibitor
|
HDAC inhibitor
|
PP-2
|
-0.4
|
2.83E-03
|
2.74E-01
|
1.62
|
-4.13
|
scr inhibitor
|
cyclooxyenase inhibitor
|
BMS−387032
|
-0.4
|
3.12E-03
|
2.74E-01
|
1.6
|
-4.01
|
CDK inhibitor, MCL1 inhibitor
|
CDK inhibitor
|
LY-294002
|
-0.4
|
2.94E-03
|
2.74E-01
|
1.65
|
-4.19
|
mTOR inhibitor, PI3K inhibitor
|
PI3K inhibitor
|
olvanil
|
-0.4
|
3.11E-03
|
2.74E-01
|
1.68
|
-4.2
|
TRPV agonist
|
adrenergic receptor
|
THM-I-94
|
-0.4
|
3.49E-03
|
2.74E-01
|
1.64
|
-4.04
|
HDAC inhibitor
|
HDAC inhibitor
|
Build of a TF/miRNA/mRNA Regulatory Network in Pancreatic Cancer Metastasis and Correlation Analysis of five Hub Genes and MAP1LC3B
Next, to investigate the regulatory theory of five hub genes in pancreatic cancer metastasis, a hub genes-related TF/miRNA/mRNA regulatory network was established. Interacting the miRNA-DEARG forecast from the Starbase and TF-DEARG predicted from the chEA3, a network consisting of ten TFs, five miRNAs, and five targeted mRNAs was obtained(Fig. 5b). Particularlytia, most of the hub genes were regulated and controlled by two TFs (FOS and MYC).
MAP1LC3B was recognized as a star molecule in autophagy [26]. Further, to realize the mechanism of hub genes in autophagy, the correlation of MAP1LC3B and five hub genes were investigated via Pearson correlation analysis (Fig. 5c). The result showed the expression levels of four hub genes, including RHEB, CHMP2B, MYC, and CASP1, had significant associations with MAP1LC3B, which further indicated that hub genes might participate in the metastasis process of pancreatic cancer through autophagy.
Immune Cell Infiltration Analysis and Co-expression Analysis of five Hub Genes and Metastasis Related Immune Cells
Next, we investigated how five hub genes affected the prognoses of pancreatic cancer patients. Autophagy has been reported to be related to immune cell infiltration [5, 27–28], which has previously been identified as one of the mechanisms of cancer invasion and metastasis [29, 30]. We explored the association between the expressions of five hub genes and the proportions of stromal and immune components. The immune score, stromal score, ESTIMATE score, and tumor purity in the primary and metastasis groups were calculated using the Sangerbox online analysis tool. The results displayed the scores of the metastasis group were observably lower in stromal scores (Fig. 6a), immune scores (Fig. 6b), ESTIMATE scores (Fig. 6c), and higher in tumor purity(Fig. 6d) than the primary group. The results suggested that the variation in five hub genes expression was related to the tumor microenvironment.
We examined diversities in the abundance of 22 immune cell subsets between the primary and metastasis groups in GSE19279 using CIBERSORT (Fig. 6f). T cell follicular helpers and CD8 T cells were the most abundant subpopulation of 22 immune cell forms in the primary group, and CD8 T cells and Treg cells were the amplest immune cell subpopulations in the metastasis group (Fig. 6i). The proportions of the immune cells from primary and metastasis samples represented signally intergroup bias and individual variations via principal component analysis (Fig. 6e). There were observable differences in two immune cell subsets between the two groups — the expression levels of plasma cells were much higher in the primary group, and the levels of Treg cells were much higher in the metastasis group. The two types of immune cells with significantly different infiltration were further included in the co-expression analysis.
Then, the potential communicative mechanisms of immune cells and five hub genes were investigated by Pearson correlation analysis. Firstly, we probed underlying relationships between twenty-two immune cell subgroups (Fig. 6g). The result demonstrated that the infiltrating centage of different immune cell subsets displayed moderate to strong associations. Next, the correlation analysis was carried out between the significantly different immune cell subsets and five hub genes (Fig. 6h). The results indicated that the expression level of CASP1, HDAC6, and RHEB had prominent correlations with the levels of infiltrating Plasma cells, and the expression levels of RHEB, MYC, and CHMP2B had noticeable correlations with the levels of infiltrating Treg cells in metastasis samples. These findings strongly suggested that hub genes have particular regulatory functions in immune cell infiltration in pancreatic cancer, notably Treg cells.
Summary of significant results show the following:
- Five hub genes-based risk model was correlated with the pancreatic patients’ prognosis and had good prediction performance.
- Hub genes were enriched in cellular components, biological functions, and related pathways related to autophagy, identifying small molecule compounds with potential therapeutic effects through autophagy.
- Hub genes were significantly correlated with MAP1LC3B and Treg cells.