Abstract
Background: Natural killer (NK) cells are involved in monitoring and eliminating cancers. The purpose of this study was to develop a NK cell-related genes (NKGs) in pancreatic cancer (PC) and establish a novel prognostic signature for PC patients.
Methods: Omic data were downloaded from The Cancer Genome Atlas Program (TCGA), Gene Expression Omnibus (GEO), International Cancer Genome Consortium (ICGC), and used to generate NKG-based molecular subtypes and construct a prognostic signature of PC. NKGs were downloaded from the ImmPort database. The differences in prognosis, immunotherapy response, and drug sensitivity among subtypes were compared. 12 programmed cell death (PCD) patterns were acquired from previous study. A decision tree and nomogram model were constructed for the prognostic prediction of PC.
Results: Thirty-two prognostic NKGs were identified in PC patients, and were used to generate three clusters with distinct characteristics. PCD patterns were more likely to occur at C1 or C3. Four prognostic DEGs, including MET, EMP1, MYEOV, and NGFR, were found among the clusters and applied to construct a risk signature in TCGA dataset, which was successfully validated in PACA-CA and GSE57495 cohorts. The four gene expressions were negatively correlated with methylation level. PC patients were divided into high and low risk groups, which exerts significantly different prognosis, clinicopathological features, immune infiltration, immunotherapy response and drug sensitivity. Age, N stage, and the risk signature were identified as independent factors of PC prognosis. Low group was more easily to happened on PCD. A decision tree and nomogram model were successfully built for the prognosis prediction of PC patients. ROC curves and DCA curves demonstrated the favorable and robust predictive capability of the nomogram model.
Conclusion: We characterized NKGs-derived molecular subtypes of PC patients, and established favorable prognostic models for the prediction of PC prognosis, which may serve as a potential tool for prognosis prediction and making personalized treatment in PC.
1 Introduction
Pancreatic cancer (PC) as a lethal malignancy shows a high mortality worldwide, causing over 331000 deaths per year globally (). Although advances in the treatment of PC, patients who received surgical resection have a five-year survival rate ranging from 10% to 25% (). PC was usually diagnosed at a late stage due to the impalpable symptoms at the early stage, and approximately 80%–85% of PC was unresectable or metastatic at the time of diagnosis (). Currently, chemotherapy is the main treatment for PC but remains an unsatisfactory prognosis, and more effective and precise therapies are required ().
Immunotherapy has been recently developed to help improve the prognosis of various cancer types, such as renal cell carcinoma (), non-small cell lung cancer (), hematologic malignancies (), and melanoma (). The principle of tumor immunotherapy is to fight against tumors through the activation of immune system, during which restarting and maintaining tumor-immune cycle plays a crucial role. Therapeutic targeting of immune checkpoints with immune checkpoint inhibitors has revolutionized cancer treatment (; ; ). It was reported that checkpoint blockade in combination with GVAX has the potential for clinical benefit for patients with PC (). T-cell immunity is associated with the exceptional outcome of the few long-term survivors. A study identified unique neoantigens as T-cell targets in PC patients, which might be used to guide the application of immunotherapies (). Pembrolizumab is a PD-1 inhibitor and has been approved for tumor patients with deficient mismatch repair or high microsatellite instability, including PC (). However, the efficacy was restricted to a rare population due to the complex, highly immunosuppressive tumor microenvironment of PC ().
The tumor immune microenvironment (TME) contains tumor cells, immune cells, cytokines, etc., and its heterogeneity can potentially impact the patient’s response to immunotherapy. Natural killer (NK) cells are a subset of innate immune cells and play a crucial role as effector cells against tumors. NK cell can directly kill malignant even at a relatively low ratio in the early presence of tumors () or promotes adaptive T-cell immunological responses to limit cancer cell aggressiveness (). The activation of NK cells is controlled by the integration of signals from cytokine receptors and a range of germline-encoded inhibitory and activating receptors (; ). Studies found that NK cell activity was significantly negatively correlated with the risk of malignancy (), and patients with a higher NK cell infiltration into cancers had better outcomes (; ; ). Cutting-edge immunotherapy targeting NK cells exerts great potential in the treatment of cancer and become an attractive alternative to T cell immunotherapies (; ). Accumulating evidence described the molecular characteristics of NK cells in various cancers (; ), but a comprehensive molecular characterization of NK cells in PC remains poorly understood.
In the present study, the PC patients were clustered on the basis of natural killer cell-related genes (NKGs), and further comparison of the clinicopathological, mutational, immunological and pathway characteristics among subtypes was conducted. In addition, we identified prognostic differentially expressed genes (DEGs) among subgroups and constructed a risk signature for prognosis prediction. The decision tree and nomogram model were built using clinicopathological features and the risk signature to assist in prognostic prediction and personalized treatment of patients with PC.
2 Materials and methods
2.1 Data collection and preprocessing
Transcriptome files and clinicopathological data of patients with PC were obtained from the Cancer Genome Atlas Program (TCGA) (https://tcga-data.nci.nih.gov/tcga/), Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), and the International Cancer Genome Consortium (ICGC) (https://www.icgc.org) databases. After removal of patients without complete clinical information and outcome status, as well as follow-up of fewer than 30 days, 176 PC patients from the TCGA pancreatic adenocarcinoma (TCGA-PAAD) cohort were retained as a training set. Ensembl was converted into gene symbol, and median value was kept when a genes had multiple gene symbols. The validation set contains 63 samples from the GSE57495 cohort and 215 patients of the PACA-CA cohort from the ICGC database. When multiple gene symbols appear or multiple probes appear for a gene, the median is taken as the gene expression value. A total of 134 human NKGs were downloaded from the ImmPort (https://www.immport.org/resource) database.
2.2 Consensus clustering
The prognostic NKGs were identified via univariate Cox regression analysis and were used to perform consensus clustering of PC patients. Consensus clustering analysis was conducted using the “ConsensusClusterPlus” R package to determine subgroups of PC patients based on the prognostic NKGs (). The best classification was determined using the partition around medoids (PAM) algorithm and 1-Pearson correlation distance, with 500 bootstraps.
2.3 Risk score
The DEGs among NKGs-derived clusters were screened out using “limma” package according to the false discovery rate (FDR) < 0.05 and |log2 [fold change (FC)]| > log2 (2) (). The univariate and the least absolute shrinkage and selection operator (LASSO) Cox regression analysis were adopted to identify and filter prognosis-related NKGs, respectively. Finally, by choosing the optimal penalty parameter lambda correlated with the minimum 10-fold cross-validation, multivariate Cox regression analysis was then implemented to establish the prognostic signature. The formula for the risk signature was as follows: risk score = . Where the βi represents the coefficient and Expi represents the normalized expression level of a gene. Two risk groups (high and low) were generated by a threshold of zero, and K–M analysis was conducted to compare overall survival (OS) differences between the high- and low-risk groups. The receiver operating characteristic (ROC) analysis was performed to estimate the predictive accuracy of the risk score.
2.4 Gene set enrichment analysis
GSEA was performed to analyze the differences in specific gene sets using the “GSVA” R package (). The hallmark gene sets from the Molecular Signatures Database (MSigDB), the inflammation-related gene sets (), and the angiogenesis-related gene set () were used to be analyzed. These pathways with the FDR <0.05 was considered to be significant. Functional enrichment analysis included Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) (biological process (BP), cellular component (CC), and molecular function (MF)) analysis was performed on DEGs in clusters using WebGestaltR package ().
2.5 Immune infiltration, chemotherapeutic sensitivity, and immunotherapy response predictions
The relative proportion of immune cells was calculated using the CIBERSORT algorithm (https://cibersort.stanford.edu/), which performs cell type enrichment analysis from gene expression data for 22 immune cells. The “ESTIMATE” R package was applied to estimate and extrapolate the fraction of stromal and immune cells in tumor samples (). The expression levels of the immune checkpoints were compared in different groups. To predict the chemosensitivity of osteosarcoma patients to several common anti-cancer drugs (methotrexate, paclitaxel, cisplatin, and doxorubicin), we adopted the “pRRophetic” R package to infer the half-maximal inhibitory concentration (IC50) values by constructing the ridge regression model based on Genomics of Drug Sensitivity in Cancer (GDSC) (www.cancerrxgene.org/) cell line expression spectrum and gene expression profiles ().
2.6 Establishment of a predictive nomogram
The decision tree model was applied to classify subgroups based on clinicopathogicial features and risk scores by using the “rpart” R package (https://cran.r-project.org/web/packages/rpart/index.html). The independent prognostic factors of OS for PC were identified by univariate and multivariate Cox regression analysis. A nomogram integrating the risk signature and independent prognostic clinicopathological factors was constructed in the TCGA cohort by the “rms” R package (https://cran.r-project.org/web/packages/rms/index.html). The calibration curves were utilized to evaluate the prediction accuracy between the predicted 1-, 2- and 3-year OS probabilities and the actual observations. The discriminate ability of the nomogram was assessed by time-dependent ROC curves. The decision curve analysis (DCA) was conducted to test the clinical utility of the nomogram using the “rmda” R package (https://cran.r-project.org/web/packages/rmda/index.html).
2.7 Mutation analysis
Tumor mutation burden (TMB) is was determined as the number of somatic indels and base substitutions per million bases in the coding region of the genome detected. Gene mutation data of PC patients were downloaded from the TCGA database and TMB was calculated using the “maftools” package () as previously described ().
2.8 Programmed cell death (PCD) analysis
12 PCD patterns were acquired from previous (). Altogether, 580 apoptosis genes, 52 pyroptosis genes, 87 ferroptosis genes, 367 autophagy genes, 14 cuproptosis genes, 9 parthanatos genes, 15 entotic cell death genes, 101 necroptosis genes, 8 netotic cell death genes, 7 alkaliptosis genes, 220 lysosome-dependent cell death genes, and 5 oxeiptosis genes were collected. Based on the expression data of above gene sets, ssGSEA analysis was conducted on tumor samples using the R package GSVA.
2.9 Statistical analysis
The R software (v3.6.3) was used for statistical analyses. Wilcoxon test compared differences between two groups. Survival differences were compared using K–M curves with a Log-rank test. The Cox proportional hazard model was performed to estimate the β regression coefficient, hazard ratios, p-value, and their corresponding 95% confidence interval for each of the selected risk predictors. a nomogram was constructed with the “rms” package in R. The C-index and calibration curve with the bootstrap method were used to evaluate the prediction performance of the nomogram. A p-value <0.05 was deemed to be a statistical significance.
3 Results
3.1 Molecular subtypes derived from natural killer cell-related genes
The flowchart is shown in Supplementary Figure S1. To obtain molecular subtypes of PC based on NKG, we first identified 32 NKGs that were significantly associated with the prognosis of PC (p < 0.05, Figure 1A). Notably, positive correlations among the expression of the 32 NKGs were observed in Figure 1B. Subsequently, consensus clustering of the 32 NKGs generated three stable clusters (C1, C2, and C3) in the TCGA-PAAD cohort (Figures 1C–E). Survival analysis demonstrated that the C3 cluster had a favorable prognosis whereas the C1 cluster had a poorer prognosis (Figure 1F). The individuals in the PACA-CA cohort were also divided into three clusters, which exerted similar prognosis characteristics as the clusters in the TCGA-PAAD cohort (Figure 1G). Among the 32 NKGs, the risk genes were generally overexpressed in the C1 cluster, and the protective genes were mainly elevated in the C3 clusters (Figure 1H).
FIGURE 1
3.2 Genomic landscapes among molecular subtypes
We compared defined three clusters with the molecular subtypes derived from a pan-cancer study and immune signatures (). As shown in Figure 2A, the C1 cluster presented with a higher TMB, aneuploidy, homologous recombination defects, and loss of heterozygosity. Meanwhile, a significantly higher proportion of immune signature-derived C3 subtype in our defined C3 subtype was observed (Figure 2B). The immune signature-derived C3 subtype was characterized by the overexpression of TH17 and Th1 genes, a low to moderate proliferation rate of tumor cells, and lower levels of aneuploidy and overall somatic copy number alterations. Meanwhile, the immune signature-derived C3 subtype showed a better prognosis than other subtypes, which is consistent with our defined C3 cluster showing the best prognosis, as shown in Figure 1F. The gene mutations in each cluster were compared and the top 20 genes with a lower p-value were illustrated in Figure 2C. Most mutations were present in KRAS, TP53, and SMAD4, accounting for 75.3%, 28.2%, and 19.7%, respectively. It was noticed that the C1 cluster with a poor prognosis had more gene mutations.
FIGURE 2
3.3 Pathway characteristics among molecular subtypes
GSEA was performed to elucidate the pathway features in each cluster by using the Hallmark candidate gene sets. As shown in Figure 3A, the C1 cluster was significantly enriched in 38 pathways in the TCGA cohort. Generally, the activated pathways mainly included cell cycle-related pathways, such as E2F_TARGETS, G2M_CHECKPOINT, MYC_TARGETS_V1, whereas the inhibited pathways primarily contained INFLAMMATORY_RESPONSE, COMPLEMENT, and INTERFERON_GAMMA_RESPONSE. Similar results were also observed in the PACA-CA cohort. In addition, we compared the pathway characteristics between clusters (Figures 3B–D). It revealed that PC patients with the 3 subtype had activated immune pathways, such as cell cycle-related pathways, indicating that the 32 NKGs might play vital roles in the regulation of cell cycle and TME.
FIGURE 3
3.4 Immune signatures between molecular subtypes and differences in immunotherapy/chemotherapy/PCD
Furthermore, we assessed the relative abundance of 22 immune cells in the TCGA-PAAD and PACA-CA cohorts using the CIBERSORT algorithm. As shown in Figures 4A, C, significant differences among three clusters were observed for several immune cell types, such as CD8+T cells and activated CD4+ memory T cells. Meanwhile, we observed a significantly higher immune score in the C3 cluster than in other clusters (Figures 4B, D), indicating that the C3 cluster had a higher immune infiltration. In addition, we investigated the 7 inflammation-related metagenes clusters in the three molecular subtypes. As a result, 6 of the 7 metagenes clusters were significantly differently expressed among subtypes, except for interferon (Figure 4E). Overall, the C1 cluster presented with a higher inflammation activity than other clusters. Meanwhile, we also observed a higher enrichment score of LCK and MHC-II, and STAT1 in the C1 cluster than the other two clusters in the PACA-CA cohort (Figure 4F). The ssGSEA analysis of 12 PCD patterns indicated that 10 PCD patterns had obviously differences among 3 subtypes, and in general, C1 or C3 subtype had higher ssGSEA scores (Figure 4G).
FIGURE 4
3.5 Immunotherapy response and drug sensitivity among clusters
Immunotherapy achieved favorable therapeutic effects in various cancers and immune checkpoint genes (ICG) play vital roles in these processes. Therefore, we evaluated the expression of ICGs among clusters and found an elevated expression of PD-1, PD-L1, and CTLA4 in the C3 cluster, as shown in Figure 5A. Meanwhile, we assessed the capability of clusters in predicting immunotherapy response using the T cell inflamed GEP score and observed a higher score in the C3 cluster than in other clusters (Figure 5B). INF-γ is a cytokine that plays a key role in immune regulation and anticancer immunity (), therefore, we calculated the ssGSEA score of the GOBP_RESPONSE_TO_INTERFERON_GAMMA gene set and found a significantly higher score of INF-γ response in the C3 cluster (Figure 5C). In addition, we also observed a higher CYT score in the C3 cluster than in other clusters (Figure 5D), which was used to reflect cytotoxic effects. Moreover, our data showed that the C1 cluster was more sensitive to cisplatin, gemcitabine, and erlotinib.
FIGURE 5
3.6 Establishment of a risk signature
A total of 294 DEGs among clusters were identified, as shown in Figures 6A–C. Enrichment analysis on the DEGs was performed and the results showed that the C3 cluster contained DEGs that were significantly associated with immune-related pathways (Figure 6D). Univariate COX analysis showed that 122 of the 293 DEGs were significantly associated with the prognosis of PC (p < 0.01), including 84 risk genes and 38 protective genes (Figure 7A). Subsequently, lasso COX regression was adopted to compress the gene number and found 9 candidate genes when lambda = 0.0666 (Figures 7B, C). Finally, four genes were identified after stepwise multivariate regression analysis on the 9 candidate genes and were used to construct a prognosis model (Figure 7D), RiskScore = +0.306*MET+0.299*EMP1-0.225*NGFR+0.182*MYEOV. The four gene expressions were negatively correlated with methylation level (Supplementary Figure S2). The risk score was calculated for each patient in the TCGA cohort and was used to divided the patient into the high and low group (Figure 8A). ROC analysis demonstrated a favorable predictive capability in forecasting the 1-, 3-, and 5-year survival rates (Figure 8B). Survival analysis showed a significantly difference in prognosis between the high and low groups (Figure 8C). In addition, we evaluated the robustness of the prognosis model in the PACA-CA and GSE57495 cohorts, which had similar results as the TCGA cohort (Figures 8D–G).
FIGURE 6
FIGURE 7
FIGURE 8
3.7 Differences in clinicopathological features and clusters between the high and low groups
The correlations between risk score and clinicopathological characteristics were analyzed in the TCGA and PACA-CA cohorts, and the results found significant associations between risk score and grade, but not stage, age, and gender (Figures 9A, C). Meanwhile, the risk score was significantly different among the three clusters, which manifested by a higher risk score in the C1 cluster and a lower risk score in the C3 cluster (Figures 9B, D). In addition, K-M curves showed that the risk score exhibited a favorable capability in the prognostic prediction of PC in sub-populations with specific clinicopathological features (Figures 9E, F).
FIGURE 9
3.8 Immune infiltration and pathway characteristics in different risk groups
As shown in Figure 10A, we observed a significantly difference in the relative abundance of four immune cells, including naive B cells, CD8 T cells, monocytes, and M0 macrophages, between the high and low groups in the TCGA cohort. The correlations between risk score and immune cells were illustrated in Figure 10B. In addition, a higher immune score was observed in the low group than the high group, indicating a higher immune infiltration in the low group (Figure 10C). The ssGSEA scores on each pathway were calculated for individuals and were compared between two risk groups. The results demonstrated that the High group was significantly associated with cell cycle-related pathways (Figures 10D, E).
FIGURE 10
3.9 Immunotherapy response, chemotherapy sensitivity and PCD between two risk groups
As shown in Figure 11A, we observed a significantly higher T cell inflamed GEP score in the Low group as compared with those in the High group. Our data also revealed a significantly higher response to IFN-γ and cytolytic activity in the Low group, when compared with the high group (Figures 11B, C). In addition, we found elevated expression of PD-1 and CTLA4, but not PD-L1, in the low group (Figure 11D), suggesting potential differences in immunotherapy response between the two risk groups. Chemotherapy sensitivity in different risk groups was analyzed and found that the patients in the high group were more likely to be sensitive to gemcitabine, cisplatin, and erlotinib, as shown in Figure 11E. In addition, four of 12 PCD patterns had increased ssGSEA score in low group, while 3 PCD had higher ssGSEA score in high group (Figure 11F, left). Furthermore, we analyzed the correlation between RiskScore, four model genes and 12 PCD patterns, and there were different degrees of correlation with each other (Figure 11F, right).
FIGURE 11
3.10 Improvement of the prognostic model
As shown in Figure 12A, a decision tree was constructed based on the risk score and clinicopathological features and generated four groups (Lowest, Low, Mediate, High) using three parameters (risk score, N stage, age). Survival analysis demonstrated significant differences in prognosis among the four groups (Figure 12B, p < 0.001). The correlations between the decision tree-derived groups and risk groups were illustrated in Figures 12C, D. Univariate regression analysis showed that T stage, N stage, age, and risk score was associated with the prognosis of PC, and three of them (N stage, age, and risk score) were identified as independent risk factors via multivariate regression analysis (Figures 12E, F). Therefore, a nomogram was built using the three factors (Figure 12G). It was observed that the predicted values were close to the observed values in terms of the 1-, 2, and 3-year OS (Figure 12H), indicating that the nomogram had good prediction performance. In addition, a decision curve was used to evaluate the reliability of the model, and it was observed that the risk signature and nomogram model had a higher standardized net benefit as compared with other clinicopathological features (Figure 12I).
FIGURE 12
4 Discussion
Tumor immunotherapy has brought hope for cancer treatment, and more and more studies have shown that innate immune cells, including NK cells, have unique advantages in anti-tumor immunotherapy. However, most of the current research focuses on adaptive immune cells, and the role of innate immune cells has not been paid enough attention. Studies have shown that the abundance of tumor infiltrating NK cells is closely related to the prognosis of patients with various solid tumors (; ; ). The prognostic model based on NKG has the potential ability to predict prognosis and immunotherapy response (). Meanwhile, a novel human NK cell-based immunotherapy was developed and showed efficacy in human metastatic PC models (). Inspired by these findings, we attempted to investigate the molecular subtypes of PC based on prognosis-related NKGs using transcriptomic data in this study. Distinct differences in prognosis, immunotherapy response, and drug sensitivity among subtypes were observed, indicating the crucial role of NK cells in the progression and treatment of PC. Functional enrichment analysis showed that NKGs involved in activated immune pathways, such as cell cycle-related pathways, indicating that the those NKGs might play vital roles in the regulation of cell cycle and TME. Furthermore, we developed a novel prognostic prediction signature based on DEGs that were found among NKGs-derived molecular subtypes of PC, which exerts a favorable capability of prognostic prediction.
Herein, we identified 32 prognosis-related NKGs in PC, including 12 protective genes and 20 risk genes, and the expression of most of these genes was significantly correlated. A number of studies had proposed potential roles of these NKGs in PC. For instance, the tumor necrosis factor ligand superfamily member 10 (TNFSF10), also known as TRAIL, encodes a cytokine that belongs to the tumor necrosis factor (TNF) ligand family, it preferentially induces apoptosis in transformed and tumor cells and was proposed as a prognostic indicator of PC (; ). As a well-known driver gene, KRAS frequently mutated in PC patients (), our data revealed that KRAS was the most mutated gene in the TCGA-PAAD cohort. KRAS gene mutations has been reported to be involved in the invasion and metastasis of tumor cells, as well as chemoresistance (; ). It was found that TMB was associated with the sensitivity of immunotherapy response and was more effective than ICG expression in screening patients suitable for immunotherapy (). This finding may result from the enrichment of immune cells due to the elevated production of “non-self” neoantigen under high TMB (). In addition, it was observed that the phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit beta isoform (PIK3CB) was involved in metastasis of PC cells (). Therefore, further investigation on these prognostic NKGs and their mutations might provide clues for the development of novel treatment of PC.
Three stable clusters with distinct differences in prognosis were generated based on the prognostic NKGs, and GSEA results found significant differences in cell cycle pathways and immunity-related pathways among clusters. Therefore, the inferior prognosis of patients in the C1 cluster may be partly attributed to the disturbance of cell cycle regulation, which is closely related to tumor proliferation and progression (). Meanwhile, these data indicated that the prognostic NKGs used for molecular typing play important roles in the cell cycle process and tumor immune microenvironment. For example, Rac1 plays an important role in regulating cell function, and its activation affects cell morphology (), cell cycle and gene expression (), survival and apoptosis (). Tyrosine kinase FYN was reported to be associated with mediating mitogenic signals and involved in regulating cell cycle and proliferation (). Besides, we observed significant differences in immune cells infiltration among NKG-derived clusters. The C1 cluster was characterized as so-called “cold tumor” since it presented with a lower immune cell infiltration. The tumor-infiltrating immune cells participated in tumor development and influence prognosis (), and anti-tumor activity of “cold tumor” is decreased because low immune cell infiltration could increase tumor cell escape from immune surveillance and contribute to tumor progression (). These finding may partly contribute to the significant reduction in survival of the C1 and C2 clusters. Meanwhile, a lower stromal score was observed in the C1 and C2 clusters, which was suggested to be associated with a poor OS of osteosarcoma ().
Since GSEA revealed significant inhibition of inflammatory response among clusters, we further evaluated the relationships between NKG-derived clusters and inflammatory activities by analyzing inflammatory-related metagenes. Notably, significant differences in hematopoietic cell kinase (HCK), IgG, MHC-II, src-family kinases p56 (LCK), MHC-I, and were observed among clusters. HCK plays a pivotal role in innate immunity and was overexpressed in various cancers. It could regulate the phagocytosis of neutrophils and macrophages (), as well as immune cell infiltration in the TME (). LCK is critical for proximal T-cell antigen receptor (TCR) signal transduction and is involved in the earliest steps of TCR-mediated T-cell activation (). MHC-I and MHC-II are two pivotal molecules presenting with the function of antigenpresentation, and their loss of expression would make tumor cells escape T-cell killing (). Therefore, a lower level of these inflammatory-related metagenes may partly account for the immunosuppressive microenvironment in the C1 and C2 clusters.
Discrepancy between inflammatory activities and immune cell infiltration among clusters prompted us to explore the immunotherapy response. It was suggested that ICG expression partly contribute to the success of immune checkpoint blockade therapy. Herein, we revealed significant differences in ICG expression among clusters, indicating potential differences in the response to immunotherapy among clusters. In addition, a T cell-inflamed gene expression profile (GEP) was found to be effective in predicting response to anti-PD-1-directed therapy (). Our data showed that the C3 cluster had a significantly higher T cell-inflamed GEP score, indicating that PC patients in the C3 cluster might be more sensitive to anti-PD-1 therapy. Cytokine IFN-γ plays a key role in anticancer immunity and immune regulation, and the C3 cluster presented with a higher elevated expression of the gene set that responds to IFN-γ. Moreover, the cytolytic activity score (CYT) has been considered as a useful tool to evaluate anti-tumor immunity. It has been revealed that high CYT was associated with better prognosis of colorectal cancer, which could be explained by increased immunity and cytolytic activity of T cells and M1 macrophages (). In this study, elevated cytolytic activity was observed in the C3 cluster. Moreover, our data also revealed the differences in chemotherapeutic drug sensitivity among clusters. The clusters derived from NKG have significant differences in immunotherapy and chemotherapy responses, which has potential value to guide individualized treatment strategies.
On the basis of NKG-derived clusters, we established a novel prognostic signature using the DEGs found among clusters. This prognostic signature has satisfactory prognostic performance and has shown good predictive power in immunotherapy response and chemotherapeutic drug sensitivity. Despite the promising findings obtained, several limitations in this study should be acknowledged. First, due to the high heterogeneity of the tumor immune microenvironment, the prognosis-predicting ability of NKG-derived molecular subtypes and subsequent prognostic models was limited. Second, analysis of NK cell characteristics based on single cell sequencing will help to further understand its role in PC. Finally, further study is required to investigate the underlying mechanism of the genes in the risk signature and PC patients’ outcomes.
5 Conclusion
In conclusion, we established three molecular clusters of PC using 32 prognosis-related NKGs and revealed differences in clinicopathological and genomic features, pathways, immunotherapy response, and drug sensitivity among clusters. Furthermore, a prognostic signature with robust prognosis-predicting ability was built and validated. The NKG-derived molecular clusters and prognostic signature might serve as a useful tool for assisting in the decision of individualized treatment and the selection of suitable individuals for chemotherapy.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
YL: Writing; YL, QJ, and MF collecting data; PZ analyzed the data; MZ supervised and submitted the paper.
Funding
The present study was supported by the National Natural Science Foundation of China (81972002); Natural Science Foundation of Shandong Province (ZR2022MC174).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1100020/full#supplementary-material
SUPPLEMENTARY FIGURE S1Working flow chart
SUPPLEMENTARY FIGURE S2The correlation of expression and Methylation level in 4 genes
Abbreviation
NK, Natural killer; NKGs, NK cell-related genes; PC, pancreatic cancer; TME, tumor microenvironment; TCGA, The Cancer Genome Atlas Program; GEO, Gene Expression Omnibus; PAM, partition around medoids; DEGs, differently expressed genes; FDR, false discovery rate; FC, fold change; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; MSigDB, Molecular Signatures Database; IC50, half-maximal inhibitory concentration; DCA, decision curve analysis; TMB, tumor mutation burden; ICG, immune checkpoint genes; TNFSF10, tumor necrosis factor ligand superfamily member 10; TNF, tumor necrosis factor; PIK3CB, phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit beta isoform; HCK, hematopoietic cell kinase; LCK, src-family kinases p56; TCR, T-cell antigen receptor; GEP, gene expression profile; CYT, cytolytic activity score.
References
1
AlvesP. M.de ArrudaJ. A. A.ArantesD. A. C.CostaS. F. S.SouzaL. L.PontesH. A. R.et al (2019). Evaluation of tumor-infiltrating lymphocytes in osteosarcomas of the jaws: A multicenter study. Virchows Arch.474 (2), 201–207. 10.1007/s00428-018-2499-6
2
AyersM.LuncefordJ.NebozhynM.MurphyE.LobodaA.KaufmanD. R.et al (2017). IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J. Clin. Invest.127 (8), 2930–2940. 10.1172/JCI91190
3
BalachandranV. P.ŁukszaM.ZhaoJ. N.MakarovV.MoralJ. A.RemarkR.et al (2017). Identification of unique neoantigen qualities in long-term survivors of pancreatic cancer. Nature551 (7681), 512–516. 10.1038/nature24462
4
BarnesT. A.AmirE. (2017). HYPE or HOPE: The prognostic value of infiltrating immune cells in cancer. Br. J. Cancer117 (4), 451–460. 10.1038/bjc.2017.220
5
BonaventuraP.ShekarianT.AlcazerV.Valladeau-GuilemondJ.Valsesia-WittmannS.AmigorenaS.et al (2019). Cold tumors: A therapeutic challenge for immunotherapy. Front. Immunol.10, 168. 10.3389/fimmu.2019.00168
6
BoyiadzisM. M.KirkwoodJ. M.MarshallJ. L.PritchardC. C.AzadN. S.GulleyJ. L. (2018). Significance and implications of FDA approval of pembrolizumab for biomarker-defined disease. J. Immunother. Cancer6 (1), 35. 10.1186/s40425-018-0342-x
7
BuscailL.BournetB.CordelierP. (2020). Role of oncogenic KRAS in the diagnosis, prognosis and treatment of pancreatic cancer. Nat. Rev. Gastroenterol. Hepatol.17 (3), 153–168. 10.1038/s41575-019-0245-4
8
ChalmersZ. R.ConnellyC. F.FabrizioD.GayL.AliS. M.EnnisR.et al (2017). Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med.9 (1), 34. 10.1186/s13073-017-0424-2
9
ChengF.LiQ.WangJ.WangL.LiW.ZengF. (2022). HCK is a potential prognostic biomarker that correlates with immune cell infiltration in acute myeloid leukemia. Dis. Markers2022, 3199589. 10.1155/2022/3199589
10
ChoY. H.KimM. S.ChungH. S.HwangE. C. (2017). Novel immunotherapy in metastatic renal cell carcinoma. Investig. Clin. Urol.58 (4), 220–227. 10.4111/icu.2017.58.4.220
11
ChoucairK.MorandS.StanberyL.EdelmanG.DworkinL.NemunaitisJ. (2020). Tmb: A promising immune-response biomarker, and potential spearhead in advancing targeted therapy trials. Cancer Gene Ther.27 (12), 841–853. 10.1038/s41417-020-0174-y
12
CocaS.Perez-PiquerasJ.MartinezD.ColmenarejoA.SaezM. A.VallejoC.et al (1997). The prognostic significance of intratumoral natural killer cells in patients with colorectal carcinoma. Cancer79 (12), 2320–2328. 10.1002/(sici)1097-0142(19970615)79:12<2320::aid-cncr5>3.0.co;2-p
13
ConsidineB.HurwitzM. E. (2019). Current status and future directions of immunotherapy in renal cell carcinoma. Curr. Oncol. Rep.21 (4), 34. 10.1007/s11912-019-0779-1
14
CursonsJ.Souza-Fonseca-GuimaraesF.ForoutanM.AndersonA.HollandeF.Hediyeh-ZadehS.et al (2019). A gene signature predicting natural killer cell infiltration and improved survival in melanoma patients. Cancer Immunol. Res.7 (7), 1162–1174. 10.1158/2326-6066.CIR-18-0500
15
Etienne-MannevilleS.HallA. (2002). Rho GTPases in cell biology. Nature420 (6916), 629–635. 10.1038/nature01148
16
GarridoF.AptsiauriN. (2019). Cancer immune escape: MHC expression in primary tumours versus metastases. Immunology158 (4), 255–266. 10.1111/imm.13114
17
GeeleherP.CoxN.HuangR. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One9 (9), e107468. 10.1371/journal.pone.0107468
18
GuillereyC.HuntingtonN. D.SmythM. J. (2016). Targeting natural killer cells in cancer immunotherapy. Nat. Immunol.17 (9), 1025–1036. 10.1038/ni.3518
19
HänzelmannS.CasteloR.GuinneyJ. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma.14, 7. 10.1186/1471-2105-14-7
20
HellmannM. D.CiuleanuT. E.PluzanskiA.LeeJ. S.OttersonG. A.Audigier-ValetteC.et al (2018). Nivolumab plus ipilimumab in lung cancer with a high tumor mutational burden. N. Engl. J. Med.378 (22), 2093–2104. 10.1056/NEJMoa1801946
21
HuntingtonN. D.VosshenrichC. A.Di SantoJ. P. (2007). Developmental pathways that generate natural-killer-cell diversity in mice and humans. Nat. Rev. Immunol.7 (9), 703–714. 10.1038/nri2154
22
ImaiK.MatsuyamaS.MiyakeS.SugaK.NakachiK. (2000). Natural cytotoxic activity of peripheral-blood lymphocytes and cancer incidence: An 11-year follow-up study of a general population. Lancet356 (9244), 1795–1799. 10.1016/S0140-6736(00)03231-1
23
IshigamiS.NatsugoeS.TokudaK.NakajoA.XiangmingC.IwashigeH.et al (2000). Clinical impact of intratumoral natural killer cell and dendritic cell infiltration in gastric cancer. Cancer Lett.159 (1), 103–108. 10.1016/s0304-3835(00)00542-5
24
KomatsubaraK. M.CarvajalR. D. (2017). Immunotherapy for the treatment of uveal melanoma: Current status and emerging therapies. Curr. Oncol. Rep.19 (7), 45. 10.1007/s11912-017-0606-5
25
LanierL. L. (2008). Evolutionary struggles between NK cells and viruses. Nat. Rev. Immunol.8 (4), 259–268. 10.1038/nri2276
26
LeD. T.LutzE.UramJ. N.SugarE. A.OnnersB.SoltS.et al (2013). Evaluation of ipilimumab in combination with allogeneic pancreatic tumor cells transfected with a GM-CSF gene in previously treated pancreatic cancer. J. Immunother.36 (7), 382–389. 10.1097/CJI.0b013e31829fb7a2
27
LiangJ.OyangL.RaoS.HanY.LuoX.YiP.et al (2021). Rac1, A potential target for tumor therapy. Front. Oncol.11, 674426. 10.3389/fonc.2021.674426
28
LiuQ.ChengR.KongX.WangZ.FangY.WangJ. (2020). Molecular and clinical characterization of PD-1 in breast cancer using large-scale transcriptome data. Front. Immunol.11, 558757. 10.3389/fimmu.2020.558757
29
López-SotoA.GonzalezS.SmythM. J.GalluzziL. (2017). Control of metastasis by NK cells. Cancer Cell.32 (2), 135–154. 10.1016/j.ccell.2017.06.009
30
MasieroM.SimoesF. C.HanH. D.SnellC.PeterkinT.BridgesE.et al (2013). A core human primary tumor angiogenesis signature identifies the endothelial orphan receptor ELTD1 as a key regulator of angiogenesis. Cancer Cell.24 (2), 229–241. 10.1016/j.ccr.2013.06.004
31
MayakondaA.LinD. C.AssenovY.PlassC.KoefflerH. P. (2018). Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res.28 (11), 1747–1756. 10.1101/gr.239244.118
32
MengN.GlorieuxC.ZhangY.LiangL.ZengP.LuW.et al (2019). Oncogenic K-ras induces mitochondrial OPA3 expression to promote energy metabolism in pancreatic cancer cells. Cancers (Basel)12 (1), 65. 10.3390/cancers12010065
33
MizrahiJ. D.SuranaR.ValleJ. W.ShroffR. T. (2020). Pancreatic cancer. Lancet395 (10242), 2008–2020. 10.1016/S0140-6736(20)30974-0
34
MorettaL.BottinoC.PendeD.CastriconiR.MingariM. C.MorettaA. (2006). Surface NK receptors and their ligands on tumor cells. Semin. Immunol.18 (3), 151–158. 10.1016/j.smim.2006.03.002
35
MuellerS.EngleitnerT.MareschR.ZukowskaM.LangeS.KaltenbacherT.et al (2018). Evolutionary routes and KRAS dosage define pancreatic cancer phenotypes. Nature554 (7690), 62–68. 10.1038/nature25459
36
NarayananS.KawaguchiT.YanL.PengX.QiQ.TakabeK. (2018). Cytolytic activity score to assess anticancer immunity in colorectal cancer. Ann. Surg. Oncol.25 (8), 2323–2331. 10.1245/s10434-018-6506-6
37
NelsonM. H.PaulosC. M. (2015). Novel immunotherapies for hematologic malignancies. Immunol. Rev.263 (1), 90–105. 10.1111/imr.12245
38
OkashaH.ElkholyS.El-SayedR.WifiM. N.El-NadyM.El-NabawiW.et al (2017). Real time endoscopic ultrasound elastography and strain ratio in the diagnosis of solid pancreatic lesions. World J. Gastroenterol.23 (32), 5962–5968. 10.3748/wjg.v23.i32.5962
39
O'ReillyE. M.OhD. Y.DhaniN.RenoufD. J.LeeM. A.SunW.et al (2019). Durvalumab with or without tremelimumab for patients with metastatic pancreatic ductal adenocarcinoma: A phase 2 randomized clinical trial. JAMA Oncol.5 (10), 1431–1438. 10.1001/jamaoncol.2019.1588
40
PulluriB.KumarA.ShaheenM.JeterJ.SundararajanS. (2017). Tumor microenvironment changes leading to resistance of immune checkpoint inhibitors in metastatic melanoma and strategies to overcome resistance. Pharmacol. Res.123, 95–102. 10.1016/j.phrs.2017.07.006
41
QuJ.ZhengB.OhuchidaK.FengH.ChongS. J. F.ZhangX.et al (2021). PIK3CB is involved in metastasis through the regulation of cell adhesion to collagen I in pancreatic cancer. J. Adv. Res.33, 127–140. 10.1016/j.jare.2021.02.002
42
RawlaP.SunkaraT.GaduputiV. (2019). Epidemiology of pancreatic cancer: Global trends, etiology and risk factors. World J. Oncol.10 (1), 10–27. 10.14740/wjon1166
43
RibasA.WolchokJ. D. (2018). Cancer immunotherapy using checkpoint blockade. Science359 (6382), 1350–1355. 10.1126/science.aar4060
44
RitchieM. E.PhipsonB.WuD.HuY.LawC. W.ShiW.et al (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res.43 (7), e47. 10.1093/nar/gkv007
45
RoseweirA. K.PowellA. G. M. T.HorstmanS. L.InthagardJ.ParkJ. H.McMillanD. C.et al (2019). Src family kinases, HCK and FGR, associate with local inflammation and tumour progression in colorectal cancer. Cell. Signal56, 15–22. 10.1016/j.cellsig.2019.01.007
46
SalmondR. J.FilbyA.QureshiI.CasertaS.ZamoyskaR. (2009). T-cell receptor proximal signaling via the Src-family kinases, Lck and Fyn, influences T-cell activation, differentiation, and tolerance. Immunol. Rev.228 (1), 9–22. 10.1111/j.1600-065X.2008.00745.x
47
SchumacherT. N.SchreiberR. D. (2015). Neoantigens in cancer immunotherapy. Science348 (6230), 69–74. 10.1126/science.aaa4971
48
SiegelR. L.MillerK. D.Goding SauerA.FedewaS. A.ButterlyL. F.AndersonJ. C.et al (2020). Colorectal cancer statistics, 2020. CA Cancer J. Clin.70 (3), 145–164. 10.3322/caac.21601
49
Souza-Fonseca-GuimaraesF.CursonsJ.HuntingtonN. D. (2019). The emergence of natural killer cells as a major target in cancer immunotherapy. Trends Immunol.40 (2), 142–158. 10.1016/j.it.2018.12.003
50
SunY.SedgwickA. J.PalarasahY.MangiolaS.BarrowA. D. (2021). A transcriptional signature of PDGF-DD activated natural killer cells predicts more favorable prognosis in low-grade glioma. Front. Immunol.12, 668391. 10.3389/fimmu.2021.668391
51
SunY.SedgwickA. J.KhanM. A. A. K.PalarasahY.MangiolaS.BarrowA. D. (2021). A transcriptional signature of IL-2 expanded natural killer cells predicts more favorable prognosis in bladder cancer. Front. Immunol.12, 724107. 10.3389/fimmu.2021.724107
52
TangD.LiuH.ZhaoY.QianD.LuoS.PatzE. F.Jret al (2020). Genetic variants of BIRC3 and NRG1 in the NLRP3 inflammasome pathway are associated with non-small cell lung cancer survival. Am. J. Cancer Res.10 (8), 2582–2595.
53
TengK. Y.MansourA. G.ZhuZ.LiZ.TianL.MaS.et al (2022). Off-the-Shelf prostate stem cell antigen-directed chimeric antigen receptor natural killer cell therapy to treat pancreatic cancer. Gastroenterology162 (4), 1319–1333. 10.1053/j.gastro.2021.12.281
54
ThorssonV.GibbsD. L.BrownS. D.WolfD.BortoneD. S.Ou YangT. H.et al (2018). The immune landscape of cancer. Immunity48 (4), 812–830.e14. 10.1016/j.immuni.2018.03.023
55
VillegasF. R.CocaS.VillarrubiaV. G.JimenezR.ChillonM. J.JarenoJ.et al (2002). Prognostic significance of tumor infiltrating natural killer cells subset CD57 in patients with squamous cell lung cancer. Lung Cancer35 (1), 23–28. 10.1016/s0169-5002(01)00292-6
56
WangX.DouX.RenX.RongZ.SunL.DengY.et al (2021). A ductal-cell-related risk model integrating single-cell and bulk sequencing data predicts the prognosis of patients with pancreatic adenocarcinoma. Front. Genet.12, 763636. 10.3389/fgene.2021.763636
57
WangX.NiM.HanD. (2022). Identification of a novel risk model: A five-gene prognostic signature for pancreatic cancer. Evid. Based Complement. Altern. Med.2022, 3660110. 10.1155/2022/3660110
58
WatersA. M.DerC. J. (2018). Kras: The critical driver and therapeutic target for pancreatic cancer. Cold Spring Harb. Perspect. Med.8 (9), a031435. 10.1101/cshperspect.a031435
59
WilkersonM. D.HayesD. N. (2010). ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinformatics26 (12), 1572–1573. 10.1093/bioinformatics/btq170
60
YoshidaT.ZhangY.Rivera RosadoL. A.ChenJ.KhanT.MoonS. Y.et al (2010). Blockade of Rac1 activity induces G1 cell cycle arrest or apoptosis in breast cancer cells through downregulation of cyclin D1, survivin, and X-linked inhibitor of apoptosis protein. Mol. Cancer Ther.9 (6), 1657–1668. 10.1158/1535-7163.MCT-09-0906
61
YoshiharaK.ShahmoradgoliM.MartinezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun.4, 2612. 10.1038/ncomms3612
62
YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics16 (5), 284–287. 10.1089/omi.2011.0118
63
ZhangX.ZengY.QuQ.ZhuJ.LiuZ.NingW.et al (2017). PD-L1 induced by IFN-γ from tumor-associated macrophages via the JAK/STAT3 and PI3K/AKT signaling pathways promoted progression of lung cancer. Int. J. Clin. Oncol.22 (6), 1026–1033. 10.1007/s10147-017-1161-7
64
ZhengJ.LiH.XuD.ZhuH. (2017). Upregulation of tyrosine kinase FYN in human thyroid carcinoma: Role in modulating tumor cell proliferation, invasion, and migration. Cancer Biother Radiopharm.32 (9), 320–326. 10.1089/cbr.2017.2218
65
ZouY.XieJ.ZhengS.LiuW.TangY.TianW.et al (2022). Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int. J. Surg.107, 106936. 10.1016/j.ijsu.2022.106936
Summary
Keywords
natural killer cells, pancreatic cancer, consensus clustering, nomogram, methylation, programmed cell death, prognosis
Citation
Lan Y, Jia Q, Feng M, Zhao P and Zhu M (2023) A novel natural killer cell-related signatures to predict prognosis and chemotherapy response of pancreatic cancer patients. Front. Genet. 14:1100020. doi: 10.3389/fgene.2023.1100020
Received
16 November 2022
Accepted
13 March 2023
Published
23 March 2023
Volume
14 - 2023
Edited by
Xiang Xue, University of New Mexico, United States
Reviewed by
Mingxin Yu, China Medical University, China
Tiansheng Cao, Southern Medical University, China
Updates
Copyright
© 2023 Lan, Jia, Feng, Zhao and Zhu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Min Zhu, zhumin@fybjsd.org.cn
This article was submitted to RNA, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.