Four drug metabolism-related subgroups of pancreatic adenocarcinoma in prognosis, immune infiltration, and gene mutation

2.1 Data acquisition and screening

A total of 181 PAAD case samples and clinical data were downloaded from the TCGA dataset (http://xena.ucsc.edu) [18], and 177 PAAD samples included in this study (all tumors samples) were obtained after excluding the adjacent samples (TCGA.H6.A45N.11A, TCGA.H6.8124.11A, TCGA.YB.A89D.11A, TCGA.HV.A5A3.11A). For the data quality control (QC), the count value of the PAAD gene expression matrix was used for log2(count + 1) standardization, and then 60,489 transcripts were merged with the same transcripts through Python (the mean method for processing). After removing the genes whose expression level was 0, 31,186 Ensembl_IDs were obtained. Gencode (ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/) was used for gene annotation of the transcriptome. In addition, the conversion between Ensembl_IDs and Symbol ID was performed on Python. The “normalizeBetweenArrays” function in the R language limma 3.9.19 package was used to normalize the data [19]. Finally, 298 drug metabolism-related genes were downloaded from PMID 33202946. Besides, through matching with the TCGA data set (PAAD) and eliminating genes with missing values less than 50%, the drug metabolism-related gene expression matrix was set for subsequent clustering studies.

3.1 Construction of protein–protein-interaction (PPI) network and pathway enrichment analysis

In order to explore the interaction of the expressed proteins of drug metabolism-related genes, PPI prediction analysis was performed on the prognostic-related drug metabolism proteins obtained above. STRING (v9.1) (https://www.string-db.org) was applied in the construction of the PPI network, and a threshold (score threshold = 0.4) was set to draw the online protein interaction network. In addition, the PPI score (score threshold = 0.4) was adopted to show the PPI interaction relationship more clearly. The Cytoscape (v3.8.2) tool software (https://github.com/cytoscape/cytoscape.js/tree/v3.8) was used to draw the network diagram.

The R language ClusterProfiler (v3.16.1) package (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) was used for pathway enrichment analysis [21]. Moreover, through gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, pathway enrichment analysis was performed on prognostic-related drug metabolism genes, and the false discovery rate (FDR) method was used to correct the P value. Finally, the ClusterProfiler (v3.16.1) package was used to draw the pathway enrichment map.

3.2 Unsupervised clustering analysis

To further determine the regulation mechanism of different drug metabolism-related genes in patients with PAAD, unsupervised clustering analysis was performed by the R language NbClust package (v 3.0), which was used in all TCGA samples (177 cases) (https://cran.r-project.org/web/packages/NbClust/index.html) to determine the best classification group [22]. Then, clustering was performed on the k-mean method provided by the Nbclust package. Finally, principal components analysis (PCA, R language dplyr package v 1.0.5, https://cran.r-project.org/web/packages/dplyr/, prcomp function) and linkage clustering analysis (R language pheatmap package v 1.0 .12, https://cran.r-project.org/web/packages/pheatmap/) were used to explore the differences between different subgroups.

3.3 Gene set variation analysis (GSVA) of different subgroups

In order to explore the differences in the pathway scoring status of the four different subgroups, we used GO (http://software.broadinstitute.org/gsea/msigdb) as a reference set and conducted the subgroup of 177 samples through the ‘gsva’ method provided by the R language GSVA (v 3.10.3, http://www.bioconductor.org/packages/release/bioc/html/GSVA.html). First, linkage cluster analysis was used to explore whether different subgroups have associated clustering patterns on GSVA scores and whether there was a linkage between subgroups. If linkage existed, subgroups were considered to be combined and participated in the analysis of the pathway score. Limma package (v 3.44.3) of R was used to analyze the different pathways among subgroups.

3.4 Prognostic survival analysis and clinical information association analysis among subgroups

The R language Survival (v3.2-10) package (Version3.10.3, http://www.bioconductor.org/packages/release/bioc/html/survival.html) was used for prognostic risk difference analysis based on the survival matrix of the subgroups and subtype groupings. In addition, the Kaplan–Meier (K–M) survival prognostic curve (P value < 0.05) was used to explore the difference between different subgroups.

Then, TCGA clinical data were used to count the number of WHO grades (G), sex ratios, age distributions, and responses to chemotherapy (complete response, CR; partial response, PR; progressive response, PD; stable disease, SD) included in each subtype. In addition, the Chi-square test in the MASS (v 7.3-53.1) package of R language was used to calculate whether there were differences in different clinical characteristics of each subgroup. Because each subgroup contained different number of people, the above clinical characteristics were counted as “percentage = count/total number of people” as the final index included in the study. Finally, ggplot2 was applied in the drawing of bar graphs.

3.5 Analysis of the abundance of immune infiltration in subgroups

In order to explore the differences in the abundance of immune infiltration between subgroups, 177 TCGA samples were included in this part of the study. Immune infiltration analysis was performed by the CIBERSORT (R 4.0.2) package [23], which used a leukocyte signature matrix (LM22). Besides, a variety of immune cell types in all PAAD samples were analyzed by the deconvolution method [24]. After calculating the proportion of immune infiltrating cells, the Mann–Whitney U-test was used to compare and analyze the differences in the proportion of immune infiltration of patients in different subgroups.

Furthermore, we searched for markers of different M2 macrophages subtypes (M2a, M2b, M2c, M2d) from the published literature [25,26]. The number of markers for M2a, M2b, M2c, and M2d are 6 (CD163, CD200R1, CD206, TGM2, IL1R2, CD209), 2 (CD163, CD86), 3 (CD163, TLR1, TLR8), and 2 (VEGF, IL10), respectively. Next, based on the expression matrix of all genes, the ssGSEA (single sample enrichment analysis) algorithm was used to estimate the relative infiltration level of each M2 macrophages subtype. Then, the t-test was used to compare whether there was a significant difference between the M2 macrophages subtypes of each two clusters. P < 0.05 was considered as a statistically significant threshold.

3.6 Analysis of gene mutations in subgroups

In order to explore the differences in gene mutations between different subgroups, PAAD genome Maf files (somatic mutations) and corresponding clinical signature files in the TCGA database were adopted. At first, the genomic Maf files of different subgroups were separated. Then, the maftools package of R language was used to perform gene mutation analysis [27], and the mutation waterfall chart of different subtypes was drawn. Finally, the mutation information of different subgroups was counted. Furthermore, we selected the top 2 genes with mutation frequency in each subgroup and analyzed the KEGG pathways of these genes that were significantly different between mutant and nonmutant samples.

4.1 Eighty-three drug metabolism genes were related to the prognosis of PAAD

A total of 177 cases of TCGA expression matrix were included in this study (Table S1). The TCGA expression matrix included 31,186 genes and the clinical phenotypes, namely, WHO grade, gender, age, response to chemotherapy (SD, PD, CR, and PR), TNM classification (TNM), and tumor stage ((Table S2).

After the intersection of the 298 reported drug metabolism genes with TCGA (Table S3), a total of 270 genes were used as the drug metabolism-related gene expression matrix for subsequent clustering studies (Table S4). The results of univariate Cox regression analysis between the survival data of 177 samples in TCGA and 270 drug metabolism genes suggested that 83 genes had a significant correlation with the prognosis of PAAD (P < 0.05) (Table S5).

4.2 PPI network construction and pathway enrichment analysis results

The PPI network of 83 PAAD prognosis-related drug metabolism proteins was constructed using Cytoscape (Figure 1a). A total of 248 interactions relationship pairs were obtained with an interaction score >0.4 (Table S6). It can be seen that in the PPI network the two relationship pairs with the highest degree of connection were ABCC8 (ATP binding cassette subfamily Cc member 8) – KCNJ11 (potassium inwardly rectifying channel subfamily j member 11) and PPARG (peroxisome proliferator activated receptor gamma) – RXRA (retinoid × receptor alpha).

Figure 1

Protein–protein interaction network (PPI) and enrichment analysis of 83 drug metabolism genes: (a) PPI network; (b) enrichment results of Kyoto encyclopedia of genes and genomes (KEGG); and (c) enrichment results of gene ontology (GO).

Besides, the enrichment analysis of KEGG and GO pathways was performed on the 83 genes, and a total of eight KEGG and eight GO pathways were enriched. The eight enriched KEGG pathways were as follows: drug metabolism-cytochrome P450, metabolism of xenobiotic by cytochrome P450, chemical carcinogenesis, retinol metabolism, drug metabolism-other enzymes, bile secretion, ABC transporters, and tyrosine metabolism (Figure 1b). Moreover, the GO pathway included xenobiotic metabolic process, response to xenobiotic stimulus, cellular response to xenobiotic stimulus, hormone metabolism, cellular hormone metabolic process, steroid metabolic process, retinoic acid metabolic process, and drug catabolic process (Figure 1c).

4.3 PAAD subgroups were obtained from unsupervised clustering analysis

Based on the 83 drug metabolism-related genes, patients were divided into different subgroups through unsupervised cluster analysis. As shown in Figure 2a, as the number of clusters increased, the total within sum of square (WSS) gets smaller. When the WSS decreased slowly, the effect of further increasing the number of clusters was considered not to be better, and this point was the optimal number of clusters. The final result showed that the samples can be divided into four subgroups (Figure 2a). Moreover, Cluster 1, Cluster 2, Cluster 3, and Cluster 4 contained 48, 55, 8 and 66 people, respectively (Table S7), and the linkage clustering heat map revealed that although there was a certain correlation between samples in different groups, the correlation between the same subgroup was greater, and the correlation was also reflected in different clinical characteristics (Figure 2b). Furthermore, it can be seen from PCA that the four groups were separated well (Figure 2c, Dim1 = 28.1%).

Figure 2

Unsupervised clustering of 83 drug metabolism genes: (a) the best type was determined; (b) linkage association clustering analysis (Clusters1–4); and (c) principal components analysis (PCA) of subgroup clustering.

4.4 Prognostic K–M survival analysis and clinical information association analysis results

For the above four subgroups (Clusters 1–4), the survival data of TCGA were used for K–M survival analysis. The results suggested that there were significant differences among the four subgroups (Figure 3a, P = 0.0017), confirming that the obtained subgroups were real. In addition, the results of clinical information association analysis showed significant differences in the four subgroups (Figure 3b–f). These differences were reflected in the WHO classification (G1, G2, and G3), gender ratio, survival status, and response to chemotherapy (CR, PR, PD, and SD) (Table 1, p < 0.001). In the Cluster 1 subgroup, there were more individuals with the age of <65 years, WHO G2 grade, and PR, but the survival rate was relatively lower. Similar to Cluster 1, in Cluster 2, the WHO classification results showed that G2 > G3 > G1 and the survival rate was relatively low. The difference was that more individuals in this subgroup were >65 years old and PD. Moreover, the age of all patients in the clusters was less than 65. All patients in Cluster 3 alive were alive with CR. After WHO classification, the G1 level had a clear advantage. In Cluster 4, the proportion of SD was slightly higher than that in other subgroups.

Figure 3

Prognostic analysis and clinical phenotype of the subgroup: (a) Kaplan–Meier survival analysis; (b) the distribution of WHO classifications in different subgroups; (c–f) distribution of chemotherapy response (c), gender (d), age (e), and tumor stage (f) in the four subgroups (because 53% of the samples in Figure 4c lacked data on chemotherapy response, each group of samples cannot constitute 100% ratio).

4.5 GSVA analysis of different subgroups

The correlation between the significantly enriched GO pathways was explored through linkage association cluster analysis (Figure 4a). From the perspective of subgroups, we can clearly see that the four subgroups formed interlocking modules, indicating significant differences in GO pathways between the four groups (Figure 4b). In addition, multigroup square difference analysis was used to explore significantly different GO pathways among the four subgroups, and the results indicated 17 significantly different pathways including cell cycle, metabolism and synthesis of organic matter, response to stimulants such as drugs, and transmembrane transport (Figure 4c).

Figure 4

Gene set variation analysis: (a) the linkage-heatmap describes the existence of linkage correlation modules between the GSVA scores of the GO pathway; (b) the linkage-heatmap describes the correlation modules that are linked between the GSVA scores within different subtypes of TCGA; (c) significant differences in GO pathways (P value <0.05) in the heatmap can see obvious stratification in the four subgroups.

4.6 Abundance analysis for immune infiltration among subgroups

The CIBERSORT algorithm was used to calculate the 22 immune infiltration abundances among the above subgroups, and more than 50% of immune cells with no immune abundance were excluded. Finally, 12 immune cells were included in the study (Table S8). Furthermore, the 12 immune cells were subjected to the Whitney U-test (Table S9). The results indicated that the immune score of immune cells varies among different subgroups (Figure 5a). Compared with the other three clusters, M2 macrophages had the highest immune score in Cluster 3, while M0 macrophages and M1 macrophages had the lowest immune score in Cluster 3 (Figure 5a). Further analysis of different subtypes of M2 macrophages (M2a, M2b, M2c, M2d) showed that the infiltration levels of the four subtypes were basically consistent among different subgroups (Figure 5b). In the same subgroup, M2b macrophages had the highest infiltration level, followed by M2a macrophages, then M2c macrophages, and finally, M2d macrophages.

Figure 5

Analysis of the differences in immune infiltration among subgroups: (a) immune scores of 12 kinds of immune cells in different subgroups; and (b) infiltration levels of M2a, M2b, M2c, and M2d in subgroups.

4.7 Gene mutation analysis for subgroups

Four sets of different mutation waterfall charts were drawn to count the mutation information of different subgroups (Figure 6, Table 2). The results showed that there were differences in the mutant genes and the mutation frequency among the 4 subgroups (Table 3). Besides, KRAS and TP53 with the highest mutation frequency in Clusters 1, 2, and 4 did not appear in the Top 10 mutation frequencies of the Cluster 3 subgroup.

Figure 6

Somatic mutation analysis of Cluster 1 (a), Cluster 2 (b), Cluster 3 (c), and Cluster 4 (d).

To further analyze whether the genes with missense mutation might have functional changes and further affect the pathway, we selected the genes with high missense mutation frequency (KRAS, TP53, RGPD3) in different subgroups for further analysis. The pathways with differences between mutant and nonmutant KRAS samples in Cluster 1 included taurine and hypotaurine metabolism, regulation of autophagy, pancreatic cancer, etc. (Figure 7a). KRAS mutations in Cluster 2 involved changes in pathways such as taurine and hypotaurine metabolism, immunodeficiency, antigen processing and presentation (Figure 7b). For Cluster 4, there were 27 patients with mutations in the KRAS (Figure 6). There were differences in multiple pathways between patients in the KRAS gene mutation group and patients in the nonmutation group, including immunodeficiency, natural killer cell-mediated cytotoxicity, JAK stat signaling pathway, etc. (Figure 7c). TP53 mutations in Cluster 1 patients may involve changes in pathways including pancreatic cancer and TCA cycle (Figure 7d). In Cluster 2, patients with TP53 mutation and nonmutation had changes in multiple pathways, including phenylalanine metabolism, antigen processing and presentation, etc. (Figure 7e). For TP53 mutations in Cluster 4, there were some changed pathways between mutant and nonmutant patients, including base excision repair, homologous recombination, cell cycle, etc. (Figure 7f). We only detected RGPD3 mutations in Cluster 3 patients. According to whether RGPD3 was mutated or not, the patients in Cluster 3 were divided into the mutation group and nonmutation group. It was found that there were differences between various pathways, including ubiquitin-mediated proteolysis, taurine and hypotaurine metabolism, MTOR signaling pathway, etc. (Figure 7g).

Figure 7

Significant differences in KEGG pathways between samples with genetic mutations and nonmutated samples. Significant differences in KEGG pathways between KRAS-mutated and -nonmutated samples in Clusters 1 (a), 2 (b), and 3 (c). Significant differences in KEGG pathways between TP53-mutated and -nonmutated samples in Clusters 1 (d), 2 (e), and 3 (f). Significant differences in KEGG pathways between RGPD3 mutated and nonmutated samples in Clusters 3 (g).