DEGs between TMB_high and TMB_low PDAC patients
Originally, the transcriptome data of 186 patients were downloaded from TCGA-PAAD. A total of 150 PDAC patients were included after 32 patients with other pancreatic neoplasms were excluded. Among these PDAC patients, 136 patients had TMB data and hence were enrolled in the following analysis. In addition, one patient whose TMB score deviated extremely from that of other patients was also removed from the present study. The patients were further divided into two groups (TMB_high and TMB_low) based on the median value of TMB across the whole cohort (Fig. 1A). First, we compared the baseline characteristics between the two groups. The proportion of tumors in the pancreatic head was larger in the TMB_low group (P = 0.007), while tumor purity was increased in the TMB_high group (P < 0.0001) (Table 1). Second, we investigated whether the TMB score was associated with the prognosis of PDAC patients. The results showed no significant difference between the two groups, although the median OS was slightly prolonged in patients with lower TMB scores (Fig. 1B; P = 0.14). However, TMB was to some extent determined by tumor purity. Given the obvious difference in baseline tumor purity we observed between the two groups, we compared OS between the TMB_high and TMB_low groups after adjusting for tumor purity. Interestingly, the median OS was significantly increased in the TMB_low group, which had adjusted tumor purity, compared with the TMB_high group (Fig. 1C; P = 0.03). Third, we identified TRGs and their potential functions. With strict screening criteria, a total of 718 DEGs were identified (LogFC > 2 and FDR < 0.05; Fig. 1D).
Table 1
The baseline characteristics of PDAC patients included in this study.
| TMB_low (n = 67) | TMB_high (n = 68) |
Age (year)# | 66.22 ± 1.33 | 63.82 ± 1.40 | p = 0.22 |
Gender (female %) | 32 (49.2%) | 27 (40.3%) | p = 0.38 |
Location (head %) | 59 (88.1%) | 46 (67.7%) | p = 0.007* |
T3T4 (%) | 56 (83.6%) | 61 (89.7%) | p = 0.32 |
N1 (%) | 55 (80.9%) | 44 (66.7%) | p = 0.07 |
AJCC pathologic tumor stage (a2~) | 56 (83.6%) | 46 (68.7%) | p = 0.07 |
KRAS | 61 (91%) | 64 (94.1%) | p = 0.53 |
Tumor purity# | 0.27 ± 0.02 | 0.47 ± 0.02 | p < 0.0001* |
#Data was presented as mean ± deviance. |
*Statistically significant. |
The KEGG analysis revealed that the TRGs were enriched in some classic oncologic signaling pathways, such as MAPK and HIF-1 signaling. In addition, we also observed that these TRGs may be involved in some pancreatic physiopathologies, such as pancreatic secretion and diabetes. TMB is a common biomarker for predicting the response to immunotherapy in multiple cancers. Here, we showed that some TRGs may regulate the remodeling of the immune microenvironment in PDAC. For example, TRGs may affect Th17 cell differentiation, leukocyte transendothelial migration, antigen presentation and processing and IgA production (Fig. 1E). The GO analysis also demonstrated that the TRGs were involved in neutrophil-mediated immunity, negative regulation of immune system process and MHC class II protein complex (Supplementary Fig. 1A). Furthermore, we performed GSEA to identify the upregulation or downregulation of certain gene sets in groups with different TMB scores, and the top five enriched gene sets in two gene lists (KEGG and GO) were visualized. According to the GSEA results, metabolic remodeling might be a potential downstream factor for TMB variation, since terms such as drug metabolism, cytochrome P450, glycine, serine and threonine metabolism, pentose and glucuronate interconversions and hormone regulation were enriched in PDAC patients with decreased TMB scores (Fig. 1F and Supplementary Fig. 1B).
Next, we compared the most frequent somatic mutations between the two groups. Overall, four driver genes, TP53, KRAS, SMAD4 and CDKN2A, were similar between the TMB_high and TMB_low cohorts in terms of their mutation frequencies. The ranking of other genes showed slight changes, as shown in Fig. 2A-B. For example, the mutation frequency of DAMTS15 was ranked 10th in the TMB_high group (3%), but it dropped out of the top 20 most frequently mutated genes. In addition, the co-occurrence and mutual exclusion between mutated genes were significantly different between the TMB_high and TMB_low groups. In TMB_high samples, the co-occurrence between mutated genes was extremely common, while mutual exclusion was observed for only the KRAS-KMT2C, KRAS-GNAS, KRAS-COL6A2, KRAS-ATM, TP53-GNAS, TP53-ARID1A and TP53-ATM pairs (Fig. 2C). In contrast, in TMB-low samples, less co-occurrence and mutual exclusion were observed, and a common trend of mutual exclusion widely existed in this cohort (Fig. 2D). In addition, we visualized the mutational landscape of the two groups in Supplementary Fig. 2.
The correlation between TRG expression and the prognosis of PDAC patients
We conducted univariate Cox regression to identify survival-related TRGs in PDAC patients. The expression of 9.3% (67/718) of TRGs was associated with OS, where 34 genes were favorable survival factors, while the other 33 genes were unfavorable survival factors (Supplementary table 2). Given the capability of these TRGs to predict the OS of PDAC, we then constructed a prognostic model based on their expression levels using lasso regression (Supplementary Fig. 3). Fifty-one genes were removed after lasso regression to avoid the overfitting phenomenon. Finally, 16 genes were retained for subsequent model construction. The coefficient of each gene in the model is provided in Supplementary table 3. The patients were divided into high- and low-risk groups based on their risk score calculated by the model (Fig. 3A). Their survival time and status varied along with an increased risk score (Fig. 3B). OS was significantly prolonged in the low-risk group (Fig. 3C). ROC curves were calculated to evaluate the accuracy of the model. The area under the curve (AUC) of this model was 0.849, which demonstrated good accuracy in predicting the OS of PDAC patients (Fig. 3D). Then, we created a validation cohort consisting of 655 PDAC patients from five independent datasets (Supplementary table 4). Using the same genes, coefficients and cutoff values, we divided the patients into high- and low-risk groups (N = 350 and 305, respectively; Supplementary Fig. 4A-B). The survival analysis showed that our model could accurately distinguish patients with dismal prognosis from the whole cohort (P = 0.03) (Supplementary Fig. 4C-E).
To further decipher the role of the genes used in the prognostic model, we investigated the differential expression of these genes between tumor and normal samples using bulk sequencing data from the TCGA and Genotype-Tissue Expression (GTEx) databases. Among them, four genes were upregulated in tumor samples and associated with dismal prognosis (Supplementary Fig. 5). Given the cell heterogeneity in tumor tissues, the differential expression of specific genes may not be caused by tumor cells themselves. Hence, we confirmed the differential expression of genes in a pancreatic cancer cell line and a normal pancreatic ductal cell line. The relative mRNA levels of these genes in the cell lines showed a similar trend as the bulk sequencing results, except no difference was found in terms of the mRNA expression of MMP28 between the Capan-1 cell line and the HPDE cell line (Supplementary Fig. 5).
The TMB score is associated with the remodeling of the immune microenvironment in PDAC
Although TMB is regarded as an effective biomarker for predicting the response to immunotherapy in patients with solid tumors, the effectiveness of TMB in some immunologically cold tumors, such as PDAC, remains controversial. In this context, we analyzed the association between the TMB score and immune cell infiltration using multiple algorithms. First, we used ssGSEA to calculate the activity of 29 immune signatures and further analyzed their correlation with TMB (Fig. 4A). The results showed that TMB was negatively associated with many anticancer signatures, such as CD8 + T cells and cytolytic activity. However, TMB was also negatively correlated with some immunoinhibitory factors, such as Treg cells and APC co-inhibition. Of note, different algorithms for the estimation of immune infiltration may yield conflicting conclusions. For example, while TMB was negatively associated with the fraction of CD8 + T cells using ssGSEA, Timer and CIBERSORT (Supplementary table 5), we found a positive correlation between CD8 + T cells and the TMB score using the EPIC algorithm (Fig. 4B). Some results also seemed to be complex and conflicting. For instance, more cancer-associated fibroblasts, which are normally seen as protumoral factors, were infiltrated in the TMB_low group. However, several anticancer factors, such as NK cells and cytotoxic scores, were also enriched in the TMB_low group (Fig. 4B). Negative correlations were observed between TMB and stromal (r = -0.34, P < 0.001) and immune scores (r = -0.29, P < 0.001) (Fig. 4C). Overall, the TMB score was associated with multiple components in the tumor microenvironment; however, whether it
is an effective biomarker reflecting anticancer immunity remains obscure in PDAC in view of the complex relationship between TMB and various immune signatures. Under this circumstance, we performed a more precise PDAC classification and focused on single gene-level regulation that mediated the influence of TMB on PDAC development in the following analysis.
A PDAC subgroup featuring TMB low MSI high was associated with prolonged OS
Previous studies have indicated that cancers with high MSI respond very well to immune checkpoint inhibitors[15, 16]. We hence explored the correlation between survival-related TRGs and MSI. Five genes were found to be negatively associated with MSI (Fig. 5A). Meanwhile, all these genes were upregulated in tumor samples with lower TMB. Among the 5 genes, PDX1 is a classical negative regulator of PDAC initiation, as shown in a PDX1-deleted PDAC animal model, but the roles of the other genes in PDAC remain unclear. Then, we confirmed their differential expression between tumor and normal tissues and survival relevance in the validation cohort. HHEX was identified as a gene of interest because it was downregulated in tumor tissues and was associated with prolonged OS (Fig. 5A).
To obtain more information on the influence of TMB and MSI on PDAC survival, we divided the patients into four groups based on the median TMB and MSI scores and then compared the OS among the subtypes. We found that the TMBlowMSIhigh group had the longest OS with marginal statistical significance (log-rank test P = 0.05; Gehan-Breslow-Wilcoxon test P = 0.007; Fig. 5B). Interestingly, the TMBhighMSIlow subtype had the shortest OS. Hence, we analyzed the DEGs between the TMBlowMSIhigh and TMBhighMSIlow groups. A total of 14 genes were deemed to be differentially expressed, and only one of them (ANKRD55) was associated with patient OS (Fig. 5C). In this context, we raised the possibility that ANKRD55 mediated the survival benefits of the TMBlowMSIhigh subtype. Therefore, we further studied whether ANKRD55 was differentially expressed between tumor and normal tissues. Immunohistochemical staining demonstrated that ANKRD55 was universally downregulated or even not expressed in PDAC samples. In contrast, it had medium expression across normal pancreatic tissues (Fig. 5D-E). To further validate our findings, we detected the expression level of ANKRD55 in patients from our center. Overall, ANKRD55 was highly expressed in stromal and normal ductal structures but rarely expressed in malignant ductal structures, which is consistent with the HPA results.
Next, we sought to determine whether ANKRD55 affected the survival of PDAC patients through immune regulation. Interestingly, ANKRD55 expression was positively correlated with CD8 + T cell infiltration not only in PDAC (r = 0.70, P < 0.001) but also in most other tumors (Fig. 6A-B). In addition, its expression was negatively associated with myeloid-derived suppressor cell (MDSC) infiltration (r = -0.65, P < 0.001; Fig. 6C). Then, we confirmed this association in the GEO cohort (Fig. 6C). Therefore, it is plausible that ANKRD55 inhibited PDAC development through CD8 + T cell enrichment and MDSC exclusion. We further explored the relationships between ANKRD55 expression and immune checkpoints, DNA repair-related genes and DNA transmethylase in the pan-cancer profile. The results showed that ANKRD55 expression was positively associated with most immune checkpoints in cancers, suggesting that although ANKRD55 predicted more intratumorally infiltrated CD8 + T cells, cancer cells may still evade immune system-mediated killing through immune checkpoint overexpression (Supplementary Fig. 6A). Among the four DNA transmethylases, only DNMT1 and DNMT2 were associated with ANKRD55 expression in pancreatic cancer (Supplementary Fig. 6B). Additionally, MLH1 was positively associated with ANKRD55 expression, while EPCAM was negatively associated with ANKRD55 expression in pancreatic cancer (Supplementary Fig. 6C).