Gene coexpression network approach to develop an immune prognostic model for pancreatic adenocarcinoma

Abstract

Background: Pancreatic adenocarcinoma (PAAD) is a nonimmunogenic tumor, and very little is known about the relationship between the host immune response and patient survival. We aimed to develop an immune prognostic model (IPM) and analyze its relevance to the tumor immune profiles of patients with PAAD.

Methods: We investigated differentially expressed genes between tumor and normal tissues in the TCGA PAAD cohort. Immune-related genes were screened from highly variably expressed genes with weighted gene correlation network analysis (WGCNA) to construct an IPM. Then, the influence of IPM on the PAAD immune profile was comprehensively analyzed.

Results: A total of 4902 genes highly variably expressed among primary tumors were used to construct a weighted gene coexpression network. One hundred seventy-five hub genes in the immune-related module were used for machine learning. Then, we established an IPM with four core genes (FCGR2B, IL10RA, and HLA-DRA) to evaluate the prognosis. The risk score predicted by IPM was an independent prognostic factor and had a high predictive value for the prognosis of patients with PAAD. Moreover, we found that the patients in the low-risk group had higher cytolytic activity and lower innate anti-PD-1 resistance (IPRES) signatures than patients in the high-risk group.

Conclusions: Unlike the traditional methods that use immune-related genes listed in public databases to screen prognostic genes, we constructed an IPM through WGCNA to predict the prognosis of PAAD patients. In addition, an IPM prediction of low risk indicated enhanced immune activity and a decreased anti-PD-1 therapeutic response.

Keywords: Immune profile; Immune prognostic model; Immunotherapy; Pancreatic cancer; WGCNA.

Fig. 1

Construction of the gene coexpression network for PAAD. a Checking the scale-free topology when β = 12. b The consensus gene dendrogram and corresponding module colors are shown. The vertical axis represents the gene expression value, and the horizontal axis represents the genes. Each vertical line in the dendrogram relates to a gene, and each branch indicates highly coexpressed genes as a module (one color). Twelve modules were detected and merged into 10 main modules. c Module-trait relationships. Each row represents a ME, the two columns represent the immune score and tumor purity, and each cell contains the corresponding correlation and P value. The matrix is color-coded by correlation according to the color legend. d Scatterplot of gene significance (y-axis) vs. module membership (x-axis) in the most significant module (pink module, see panel c)

Fig. 2

Analysis of enriched GO terms for the hub genes in the pink module. The analysis of enriched GO terms was performed using the function “enrichGO” in the “clusterProfiler” package. Biological process (BP, panel a); cellular component (CC, panel b); molecular function (MF, panel c). The y-axis represents the number of genes associated with the GO term. The intensity and color of dots are indicated on the right side of the heatmap and are represented by their corresponding adjusted P values

Fig. 3

Analysis of the prognostic value of the IPM. The risk scores (a) and OS (b) of each patient. The patients were ranked by risk score. The dot plot shows the survival status of each patient. Red: deceased; pink: alive. Kaplan-Meier survival curves showing the OS times of patients stratified into low/high-risk groups (c). P values were obtained from the log-rank test. Forest plot of the multivariate Cox regression proportional hazards regression analysis of OS in TCGA PAAD cohort (d). CI, confidence interval; HR, hazard ratio

Fig. 4

Associations between low/high risks and immune profiles. The activated CD8⁺ T cell fraction (a) and cytolytic activity (b) in the low-risk PAAD group were significantly increased compared with the high-risk group. The IPRES signatures of the low-risk PAAD group were significantly decreased compared with the high-risk group (c). P values were calculated with the Wilcoxon test; the box shows the upper and lower quartiles (*P < 0.05, **P < 0.01, and ***P < 0.001). Heatmap showing scores for IPRES signatures in TCGA PAAD cohort (d)