Multiple Perspectives Reveal the Role of DNA Damage Repair Genes in the Molecular Classification and Prognosis of Pancreatic Adenocarcinoma

Abstract

Pancreatic adenocarcinoma (PAAD) is a highly heterogeneous and immunosuppressive cancer. This study investigated the diversity of DNA damage repair (DDR) and immune microenvironment in PAAD by transcriptomic and genomic analysis. Patients with PAAD were divided into two DDR-based subtypes with distinct prognosis and molecular characteristics. The differential expression genes were mostly enriched in DDR and immune-related pathways. In order to distinguish high- and low-risk groups clinically, a DDR- and immune-based 5-gene prognostic signature (termed DPRS) was established. Patients in the high-risk group had inferior prognosis, a low level of immune checkpoint gene expression and low sensitivity to DDR-associated inhibitors. Furthermore, single-cell sequencing was used to observe the performance of the DDR-based signature in a high dimension, and immunohistochemistry was used to verify the relationship between the genes we identified and the prognosis of patients with PAAD. In conclusion, the DDR heterogeneity of PAAD was demonstrated, and a novel DDR- and immune-based risk-scoring model was constructed, which indicated the feasibility of DPRS in predicting prognosis and drug response in PAAD patients.

Keywords: DNA damage repair; immune; molecular classification; pancreatic adenocarcinoma; prognosis.

Figure 1

Workflow of the current study.

Figure 2

Consensus clustering for DNA damage repair (DDR) related genes in the TCGA-PAAD cohort (n = 168). (A) The consensus matrix of the DDR-related genes. (B) Cumulative distribution function (CDF) plot with consensus values ranging from 0 to 1. (C) The Nbclust plot represented the chosen optimal cluster number (k = 2) for DDR genes. (D) Heatmap of DDR-related genes in DDR subtypes. (E) Overall survival (OS) of patients in DDR-subtype1 and subtype2.

Figure 3

Somatic mutation and immune features between subgroups. Landscape of mutation profiles in DDR-subtype1 (A) and DDR-subtype2 (B). Mutation information of each gene in each sample is shown in the waterfall plot. Top panel shows individual tumor mutation burden. (C) Patients with DDR-subtype1 showed higher mutation frequencies of KRAS and SMAD4; (D) Patients with DDR-subtype2 showed higher mutation frequencies of TP53 and CDKN2A. Immune profile alterations (E) and human leukocyte antigen (F) between the DDR-subtypes. * represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001, ns represents no significant difference.

Figure 4

Identification of differential expression genes (DEGs) and related biological processes of DDR subtypes. Volcano plot (A) and heatmap (B) of DEGs in PAAD based on data from TCGA. (C) The top 20 hallmark gene sets of Gene Set Enrichment Analysis of DEGs.

Figure 5

Generation of DDR- and immune-based signature. (A) Venn-diagram-described DDR-based DEGs (n = 1081) were intersected with the immune genes (n = 2483) and altogether 191 overlapping genes were screened. (B) Forest plot of univariate Cox regression analysis showed that 13 DEGs were associated with prognosis. (C) LASSO coefficient profiles. (D) Selection of the tuning parameter (lambda) in the LASSO model by 10-fold cross-validation based on minimum criteria for OS. (E) Heatmap of relationship between the 5-gene signature and immune infiltration. (F) Heatmap of 5-gene signature by unsupervised clustering. The DDR subtype, risk group and risk score as gene annotations were correlated.

Figure 6

Construction and validation of DDR- and immune-based risk score (DPRS) model. (A–D) Construction and validation of TCGA training set. (A) The OS of training set. (B) Distribution of DPRS and OS of training set. Time-dependent ROC curves (C) and calibration curves (D) validation at 1, 2, and 3 years of prognostic value in TCGA cohort. (E–H) Construction and validation of GEO external validation set (GSE85916). (E) The OS of validation set. (F) Distribution of DPRS and OS of validation set. Time-dependent ROC curves (G) and calibration curves (H) validation at 1, 2, and 3 years of prognostic value in GEO cohort.

Figure 7

Correlation of risk score with TMB and clinical characteristics (A) Comparison of TMB between two DDR subtypes. (B) The OS of patients in the high- and low-risk groups. (C) The OS of patients in the risk groups combined with the TMB groups. (D) Univariate Cox analysis and (E) multivariate Cox analysis of clinical characteristics. (F) Nomogram predicting OS for PAAD patients. Time-dependent ROC curves (G) and calibration curves (H) validation at 1, 2, and 3 years of prognostic value in nomogram. *** represents p < 0.001.

Figure 8

The boxplots and scatter plots of CTLA4 (A) and PDCD1 (B) in the low-and high-risk groups. (C) Boxplots showing estimated IC50 values for Olaparib, Niraparib, Berzosertib, AZD6738 and MK8776 in the TCGA-PAAD dataset (D) The uniform manifold approximation and projection (UMAP) plot demonstrates main cell types in PAAD. (E) The distribution of each type and DDR-based score expression in PAAD. (F) DDR-based scores in different cells are various (p < 0.05). (G) Staining scores of MET and ERAP2 were graded at four levels and two groups. Scale bar, 100 μm. (H) The OS of patients in the high and low expression of MET or ERAP2.