Development and Validation of an ADME-Related Gene Signature for Survival, Treatment Outcome and Immune Cell Infiltration in Head and Neck Squamous Cell Carcinoma

ADME genes are a set of genes which are involved in drug absorption, distribution, metabolism, and excretion (ADME). However, prognostic value and function of ADME genes in head and neck squamous cell carcinoma (HNSCC) remain largely unclear. In this study, we established an ADME-related prognostic model through the least absolute shrinkage and selection operator (LASSO) analysis in the Cancer Genome Atla (TCGA) training cohort and its robustness was validated by TCGA internal validation cohort and a Gene Expression Omnibus (GEO) external cohort. The 14-gene signature stratified patients into high- or low-risk groups. Patients with high-risk scores exhibited significantly poorer overall survival (OS) and disease-free survival (DFS) than those with low-risk scores. Receiver operating characteristic (ROC) curve analysis was used to confirm the signature’s predictive efficacy for OS and DFS. Furthermore, gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analyses showed that immune-related functions and pathways were enriched, such as lymphocyte activation, leukocyte cell-cell adhesion and T-helper cell differentiation. The Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) and other analyses revealed that immune cell (especially B cell and T cell) infiltration levels were significantly higher in the low-risk group. Moreover, patients with low-risk scores were significantly associated with immunotherapy and chemotherapy treatment benefit. In conclusion, we constructed a novel ADME-related prognostic and therapeutic biomarker associated with immune cell infiltration of HNSCC patients.

Introduction

Head and neck squamous cell carcinoma (HNSCC), which originates from the squamous epithelium of the mouth; nasal cavity and paranasal sinuses; oropharynx; larynx and hypopharynx, is one of the most prevalent malignant cancers with more than 800,000 new cases each year (1, 2). Despite multidisciplinary treatment of locally advanced HNSCC by surgery combined with adjuvant chemoradiation or platinum-based concurrent chemoradiation, less than 50% of patients can be cured (3, 4). Tumor-node-metastasis (TNM) staging system helps guide therapeutic decisions and predict clinical outcomes. However, this anatomic staging system does not adequately reflect tumor complexity and heterogeneity in an individual patient. Therefore, identification of effective and sensitive biomarkers for the diagnosis, therapeutic response and prognosis of HNSCC is urgently needed.

Genes involved in drug absorption, distribution, metabolism and excretion are defined as ADME genes (5–7). According to the PharmaADME consortium, ADME genes consist of 32 core genes and 266 extended genes (http://www.pharmaadme.org). These genes are divided into several groups based on their roles in the pharmacokinetic process, such as phase I and II drug-metabolizing enzymes, drug transporters and modifiers (8–10). Growing evidence indicates that ADME gene polymorphism contributes to interindividual variability in carcinogenesis and drug response (11, 12). Consistently, ADME gene expression regulated at the transcriptional, translational and epigenetic levels also shows variability in the population (13–15).

Recent studies have revealed that ADME genes can act as biomarkers to predict therapy response, adverse drug reactions, drug resistance and survival outcomes based on their critical roles in the pharmacokinetic process. Previous research has demonstrated that CYP1B1 and ABCB1 could predict the clinical response to taxane therapy in breast cancer (16). Suthandiram et al. showed that ABCB1 was associated with an increased probability of methotrexate adverse events in patients with haematological malignancies (17). Zhang et al. reported that UGT1A1 expression correlated with 5-fluorouracil resistance in esophageal squamous cell carcinoma (18). Hu et al. identified a set of core ADME genes that could predict overall survival (OS) in many cancer types (19). However, the prognostic prediction values and biological functions of ADME genes in HNSCC have not yet been systematically evaluated.

In this study, we used gene expression data and clinical information from The Cancer Genome Atlas (TCGA) (n = 494) to identify differentially expressed ADME genes in HNSCC. We then defined a 14-gene signature for predicting survival outcomes by the least absolute shrinkage and selection operator (LASSO) analysis in TCGA training cohort (n = 247) and validated it through the internal TCGA validation cohort (n = 247) and an external Gene Expression Omnibus (GEO) cohort (n = 88). Finally, we performed functional enrichment, immune cell infiltration, immunotherapy response, mutation and chemotherapy response analyses between the high-risk group and the low-risk group. Our results indicated that the ADME-related signature was associated with immune cell infiltration and could predict prognosis and therapeutic response in HNSCC.

Materials and Methods

HNSCC Dataset Collection

We obtained the gene expression data and clinical annotation of the TCGA-HNSC and TCGA-SKCM datasets from UCSC Xena (https://xena.ucsc.edu/). Patients with no survival information or RNA expression data were removed from further analysis. An independent dataset (GSE102349, n = 113, 88 with survival outcome) from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) was analyzed to validate the prognostic value of the gene signature. Data and clinical information of IMvigor 210 cohort were obtained from the IMvigor210CoreBiologies R package. The probe signal intensity data were log2 transformed and quantile normalized. We transformed the ENSEMBL Gene ID to the Gene Symbol ID. Genes were excluded if they did not express in more than half of the samples. We conducted TCGA and GEO dataset analysis using R software (https://www.r-project.org/).

Extraction of Differentially Expressed ADME Genes

The 298 ADME genes that are currently defined by the PharmaADME Consortium (http://www.pharmaadme.org) consist of 32 core genes and 266 extended genes (8, 9). Their full names are shown in Table S1. Differentially expressed genes between tumor tissues and adjacent normal tissues in the TCGA-HNSC cohort was identified by the ‘limma’ R package. Genes with P value < 0.05 were identified as dysregulated genes. Differentially expressed ADME genes were extracted through VENNY (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

Identification of a Prognostic ADME-related Gene Signature

We first estimated the prognostic value of the overlapping dysregulated ADME genes through univariate Cox regression analysis. Next a protein–protein interaction (PPI) network for survival-related genes was established through the STRING database (20). We then chose these dysregulated prognostic genes to develop possible risk score through LASSO Cox regression analysis with the ‘glmnet’ R package. Penalty parameter lambda (λ) of the model was chosen according to 10-fold cross validation. The risk score of each patient was calculated according to the normalized expression of the candidate genes (Expi) and their corresponding regression coefficients (Coei). The formula of the risk score was constructed as follows:

$$\text{Risk score}={\displaystyle {\sum }_{i=1}^N\left(Expi \times Coei\right)}$$

HNSCC patients were stratified into high-risk and low-risk groups according to the median cutoff value. Afterwards Kaplan–Meier and ROC curve analyses were conducted to evaluate the prognostic performance of the novel gene signature with the ‘survminer’, ‘survival’ and ‘survivalROC’ R packages. Finally, univariate and multivariate Cox regression analyses were carried out to estimate prognostic independence of the ADME-related risk score and other clinical parameters in patients with HSNCC.

Functional Enrichment Analysis

The differentially expressed genes between the high-risk and low-risk groups were defined by the threshold of |fold change| ≥ 1.5 and P < 0.05. These genes were selected to perform gene ontology (GO) analysis and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis using the ‘clusterProfiler’ R package. The P-values of GO terms or KEGG pathways were corrected by the Benjamini–Hochberg method. Geneset variation analysis (GSVA) was carried out using the ‘GSVA’ package.

Estimation of Immune Cell Infiltration

The extent of immune cell infiltration of the tumor microenvironment (TME) in each HNSCC patient was estimated by the Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm (21). The relative abundance of 22 tumor-infiltrating immune cell types (naive B cell, memory B cell, CD8+ T cell, naive CD4+ T cell, resting/activated memory CD4+ T cell, regulatory T cells (Treg), follicular helper T cell, gamma delta T cell, resting/activated dendritic cell, M0/M1/M2 macrophage, eosinophil, neutrophil, resting/activated mast cell, resting/activated NK cell, monocyte and plasma cell) in HNSCC samples was evaluated using the ‘cibersort’ R package. Moreover, the relationship between the risk score and immune cell infiltration was calculated using Spearman correlation analysis, and the results were displayed through the ‘ggstatsplot’ R package.

The stromal score, immune score and ESTIMATE score of each sample were computed through ‘ESTIMATE’ R package (22). The xcell score which reflects the integrative immune activity in the tumor microenvironment was calculated using ‘xcell’ R package (23). Moreover, the other immune response-related scores (CYT score, MHC score, CD8 T effector score) were calculated using the related gene sets according to previous studies (24–26).

Drug Sensitivity Prediction

Drug sensitivity was estimated through ‘pRRophetic’ R package (27). The ridge regression was performed to calculate the 50% of maximum inhibitory concentration (IC50) based on the Genomics of Drug Sensitivity in Cancer (GDSC) database.

Statistical Analysis

A Student’s t-test was used to compare the continuous variables and a chi-square test was used to compare the categorical variables. The Kaplan–Meier method was conducted to evaluate OS and DFS, through the log-rank test. Univariate and multivariate Cox proportional hazards regression analysis was carried out to identify independent prognostic factors. A P value <0.05 was considered as significant. The above statistical analyses were performed using R software (Version 3.6.3).

Results

Identification of Dysregulated and Prognostic ADME Genes in HNSCC

The flow chart of this study is shown in Figure 1. We determined the gene expression profiles by comparing HNSCC samples and normal samples in the TCGA cohort and found 18353 dysregulated genes (P < 0.05). The heatmap and volcano plot of differentially expressed mRNAs were presented in Figure 2A and Figure S1. The Venn diagram showed that more than half of the ADME genes (160/298, 53.6%) were differentially expressed in HNSCC (Figure 2B). We performed the univariate COX analysis and identified 19 prognostic ADME genes. The Kaplan–Meier survival curves for the 19 ADME genes in HNSCC are shown in Figures 2C–H and Figure S2 (all P < 0.05). Moreover, we constructed the PPI network of 19 prognostic ADME genes (Figure 2I), and the correlation between these genes is presented in Figure 2J.

FIGURE 1

Figure 1 Flow chart of data collection, analysis and validation.

FIGURE 2

Figure 2 Identification of prognostic genes involved in drug absorption, distribution, metabolism, and excretion (ADME) in the Cancer Genome Atla (TCGA) cohort. (A) Heatmap showing the gene expression profiles of tumor (T) and normal (N) tissues. High, high expression; Low, low expression. (B) Venn diagram to show differentially expressed ADME genes. (C–H) Kaplan–Meier curves of overall survival (OS) for 6 representative prognostic ADME genes in HNSCC, all P <0.01. (I) A PPI network indicating the interactions among the 19 prognostic ADME genes. (J) Spearman correlation analysis of the 19 ADME genes.

Construction of a Prognostic ADME Gene Signature in the TCGA Training Cohort

To further investigate the prognostic roles of ADME genes, 494 patients of head and neck cancer from the TCGA project were divided into the training cohort (n = 247) and the internal validation cohort (n=247). We conducted LASSO penalized Cox regression analysis using the expression profile of the 19 ADME genes mentioned above (Figure 3A) in the TCGA training cohort. We then selected 14 genes (4core ADME genes and 10 extended ADME genes) according to the optimal λ value and the coefficients of candidate genes were derived from the LASSO algorithm (Figure 3B). The risk score of each patient was calculated as follows:

$$\begin{array}{l}\text{Risk score}=0.112\times \text{ UGT}1\text{A }-0.129\times \text{ SULT}1\text{B}1+0.223\times \text{ SOD}1+0.062\times \text{ SLCO}1\text{B}3+0.617\times \\ \text{SLCO}1\text{B}1+0.138\times \text{ PPARG }-1.328\times \text{ CYP}8\text{B}1-0.042\times \text{ CYP}7\text{B}1+0.052\times \text{ CYP}4\text{F}12-0.898\times \\ \text{CYP}2\text{D}6-0.009\times \text{ ALDH}2-0.076\times \text{ ABCB}8-0.278\times \text{ ABCB}4-0.507\times \text{ ABCB}1.\end{array}$$

FIGURE 3

Figure 3 Construction of the 14-gene signature in the TCGA train cohort. (A) Forest plot of the univariate Cox regression analysis between expression of 19 ADME genes and OS. (B) The coefficients of 14 ADME genes. (C) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the training cohort. Log-rank test, P <0.0001. (D) Time-dependent ROC curve of risk score in train cohort on OS. (E) Kaplan–Meier curves for the disease-free survival (DFS) of high-risk and low-risk groups in the training cohort. Log-rank test, P=0.0034. (F) Time-dependent ROC curve of risk score in train cohort on DFS.

The HNSCC patients in the TCGA training cohort were stratified into a high-risk group or a low-risk group based on the median cut-off value (0.774). The Kaplan–Meier survival curves indicated that patients in the high-risk group exhibited significantly poorer OS (P < 0.0001) and DFS (P = 0.0034) than those in the low-risk group (Figures 3C, E). Time-dependent ROC curves were conducted to evaluate the predictive efficacy of the novel signature for survival. The area the curve (AUC) for OS reached 0.666 at 1 year, 0.727 at 3 years and 0.702 at 5 years (Figure 3D). AUC for 1-, 3- and 5-year DFS were 0.698, 0.648 and 0.625, respectively (Figure 3F).

Validation of the 14-Gene Signature in 2 Cohorts

The robustness of the ADME gene signature was tested in the internal validation TCGA cohort (n = 247) and the external validation GEO cohort (GSE102349, n = 88) using the same risk score formula. The validation cohorts were also categorized into high- or low-risk groups by the same cut-off value from the training cohort (0.774). Consistent with the results of the TCGA training cohort, patients in the high-risk group had a shorter overall survival and disease-free survival than those in the low-risk group (Figures 4A, C, E). Moreover, ROC curves were applied and the AUC for survival were shown in Figures 4B, D, F, respectively.

FIGURE 4

Figure 4 Validation of the ADME gene signature in the internal and external cohort. (A) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the internal cohort. Log-rank test, P <0.0001. (B) Time-dependent ROC curve of risk score in internal cohort on OS. (C) Kaplan–Meier curves for the DFS of high-risk and low-risk groups in the internal cohort. Log-rank test, P=0.0001. (D) Time-dependent ROC curve of risk score in internal cohort on DFS. (E) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the GSE102349 cohort. Log-rank test, P=0.0015. (F) Time-dependent ROC curve of risk score in GSE102349 cohort on OS.

Independent Prognostic Value of the 14-Gene Signature

The clinical information of the cohort from the TCGA database is shown in Table 1. Univariate Cox regression analysis indicated that the risk score, gender, lymphovascular invasion and perineural invasion were all significantly associated with OS (Table 2). After adjustment for other confounding factors, the risk score was still correlated with poor overall survival in the multivariate Cox regression analysis (HR = 2.84, 95% CI = 1.92-4.21, P < 0.001, Table 2). Chi-square test showed no statistical significance concerning lymphovascular invasion (P = 0.7781) and perineural invasion (P = 0.1995) between different risk groups (Figures 5A, D).

TABLE 1

Table 1 Clinical characteristics of the patients stratified by low and high expression of ADME risk score in TCGA-HNSC cohort.

TABLE 2

Table 2 Univariate and multivariable Cox regression analysis of prognostic factors in 494 patients with HNSCC.

FIGURE 5

Figure 5 Association of risk score with lymphovascular invasion and perineural invasion in TCGA cohort. (A) Rate of patients with lymphovascular invasion (Yes with No) in high and low risk groups in the TCGA-HNSC cohort. (B) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the patients with lymphovascular invasion TCGA cohort. Log-rank test, P <0.0001. (C) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the patients without lymphovascular invasion TCGA cohort. Log-rank test, P =0.0011. (D) Rate of patients with perineural invasion (Yes with No) in high and low risk groups in the TCGA-HNSC cohort. (E) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the patients with perineural invasion TCGA cohort. Log-rank test, P <0.0001. (F) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the patients without perineural invasion TCGA cohort. Log-rank test, P <0.0001.

Stratification analysis was performed to further evaluate whether the ADME-related score was independent of lymphovascular invasion and perineural invasion. Patients in the TCGA dataset were stratified into the lymphovascular invasion group and the non-lymphovascular invasion group. The Kaplan–Meier survival curves showed a significant difference in OS between the high-risk and low-risk groups in the lymphovascular invasion group (P <0.0001, Figure 5B) as was in the non-lymphovascular invasion (P = 0.0011; Figure 5C). Next, all HNSCC patients were also classified by perineural invasion. Patients with high-risk scores exhibited poorer OS than those with low-risk scores in perineural invasion group (P <0.0001, Figure 5E) as well as in non-perineural invasion group (P <0.0001, Figure 5F). These results suggest that the novel signature derived from ADME genes is an independent prognostic predictor in HNSCC patients.

Functional and Immune Infiltration Analyses in Different Risk Groups

GO and KEGG pathway analyses were carried out to explore biological functions and pathways using the dysregulated genes in the two groups. Interestingly, GO analysis revealed that immune-related functions were enriched, including lymphocyte activation, T cell activation, leukocyte migration, leukocyte cell-cell adhesion, cytokine-mediated signaling pathway and leukocyte differentiation (Figure 6A). Likewise, KEGG analysis showed that the differentially expressed genes were enriched in immune-related pathways, such as Th1, Th2 and Th17 cell differentiation; cell adhesion molecules; viral protein interaction with cytokine and cytokine receptor; allgraft rejection; and the intestinal immune network for IgA production (Figure 6B). Moreover, GSVA analysis showed low-risk group was obviously enriched in immune-related pathways, including inflammatory response, interferon gamma response, IL2/STAT5 signaling and IL6/JAK/STAT3 signaling (Figure S3).

FIGURE 6

Figure 6 Tumor immune features and immune cell infiltration analyses in high-risk and low-risk groups. Gene ontology (GO) (A) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway (B) analyses were performed on dysregulated genes in different risk groups. (C) The infiltration levels of 22 immune cells in high-risk and low-risk groups calculated by The Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) analysis (Wilcoxon test, ns: not significant; **P < 0.01; ***P < 0.001; ***P < 0.0001). (D–G) Linear correlation between risk score and naive B cell infiltration levels (D), CD8 T cell infiltration levels (E), follicular helper T cell infiltration levels (F), M2 macrophage cell infiltration (G). (H) Boxplot of ESTIMATE score in high and low risk groups (Wilcoxon test, ***P < 0.001; ***P < 0.0001). (I) Boxplot of xCell-immune score, CYT score, MHC score in high and low risk groups (Wilcoxon test, ***P < 0.001; ***P < 0.0001). (J) Boxplot of the score of CD8 T effector and immune checkpoint in high and low risk groups (Wilcoxon test, ***P < 0.0001). (K) Correlation of risk score with other immune features by Pearson correlation analysis.

To further elucidate the relationship between the risk score and immune cell infiltration in the tumor microenvironment, bioinformatic analysis was applied using the CIBERSORT algorithm. The results indicated that immune cell infiltration levels were significantly different between the high-risk and low-risk groups (Figure 6C). The following analysis revealed that the infiltration levels of naive B cells, CD8 T cells and follicular helper T cells were negatively correlated with the risk score, while the infiltration levels of M2macrophages were positively correlated with the risk score (Figures 6D–G, all P < 0.05).

Additionally, we performed other immune cell infiltration analyses and the results showed that the immune scores (stromal score, immune score, ESTIMATE score, xcell score, CYT score, MHC score, CD8 T effector score and immune checkpoint score) in low-risk group were higher compared with those in high-risk groups (Figures 6H–J). ADME-related risk score was negatively correlated with these immune infiltration scores (Figure 6K). The higher infiltration levels of immune cells, especially B cells and T cells, in the low-risk group contributed to a better immune response against cancer cells, which may explain why the low-risk group had better prognosis than the high-risk group in HNSCC.

Immunotherapy and Chemotherapy Response Analyses

Our previous results revealed the significance of ADME-related risk score in evaluating tumor immune microenvironment. To further explore the relationship between the ADME score and benefit of immunotherapy, we collected two external cohorts (IMvigor210 and TCGA-SKCM) which received immunotherapy and performed following analyses. The results showed that patients with low-risk scores had better prognosis (Figures 7A, D) and response to immunotherapy (Figures 7B, E), indicating that the low-risk group was more likely to benefit from immunotherapy. ROC curves also proved the efficacy of the ADME-related risk score in predicting responsiveness to immunotherapy (Figures 7C, F).

FIGURE 7

Figure 7 The role of risk scores in the prediction of therapeutic outcomes. (A) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the IMvigor210 cohort. Log-rank test, P <0.0001. (B) Rate of CR (Complete Response)/PR (Partial Response) and PD (Progressive Disease)/SD (Stable Disease) in high and low risk groups in the IMvigor210 cohort. (C) ROC curve of risk score in IMvigor210 cohort on OS. (D) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the TCGA-SKCM cohort with Immunotherapy. Log-rank test, P =0.042. (E) Rate of CR/PR and PD/SD in high and low risk groups in the TCGA-SKCM cohort with Immunotherapy. (F) ROC curve of risk score in TCGA-SKCM cohort with Immunotherapy on OS. (G) Kaplan–Meier curves for the OS of high-risk and low-risk groups in the TCGA-HNSC cohort with Chemotherapy. Log-rank test, P =0.00025. (H) Rate of CR/PR and PD/SD in high and low risk groups in the TCGA-HNSC cohort with Chemotherapy. (I) ROC curve of risk score in TCGA-HNSC cohort with Chemotherapy on OS. (J) The estimated IC50s of clinical chemotherapeutic and targeted drugs of HNSCC in high-risk and low-risk groups. From left to right, cisplatin, docetaxel, paclitaxel, 5-fluorouracil, gemcitabine, cetuximab.

In addition, the predictive value of the ADME-related signature for chemotherapy was evaluated using the TCGA-HNSC cohort. Consistent with the results in the cohorts received immunotherapy, patients with low-risk scores showed longer survival and better response (Figures 7G, H). ROC analysis confirmed the predictive efficacy of the ADME signature for chemotherapy response (Figure 7I). Subsequently, we performed the “pRRophetic” algorithm predict the sensitivity to the 6 chemotherapeutic and targeted drugs. The results indicated that IC50s of cisplatin, docetaxel, paclitaxel, 5-fluorouracil and cetuximab were significantly different between the high-risk and low-risk groups (Figure 7J). In conclusion, we identified a novel ADME-related gene signature to predict immunotherapy and chemotherapy response of patients in HNSCC.

Discussion

HNSCC, as the sixth most prevalent malignant tumor, exhibits poor clinical outcomes because of local recurrence and metastasis (28). Identification of novel predictive biomarkers in HNSCC is needed for developing individualized therapy (29, 30). In the present study, a novel ADME-related prognostic model was constructed based on the TCGA training cohort and its robustness was validated by the internal TCGA validation cohort and an external GEO cohort (GSE102349). The risk scores were calculated by the LASSO algorithm to predict the patients’ prognosis in HNSCC. The analysis indicated that patients in high-risk group exhibited shorter OS and DFS, both in the training and validation cohorts. ROC analysis confirmed the predictive efficacy of the 14-gene signature for survival of the two cohorts. The above findings suggest that our ADME-related risk score was a valuable predictor for the prognosis of HNSCC patients.

More than half of ADME genes (160/298, 53.6%) are dysregulated in HNSCC according to our analysis, indicating that ADME genes may play important roles in cancer development and progression. Here, we identified 14 ADME genes to establish the prognostic model and these genes were classified into 4 groups: phase I drug-metabolizing enzymes (CYP2D6, CYP8B1, CYP7B1, ALDH2, CYP4F12), phase II drug-metabolizing enzymes (UGT1A10), transporters (ABCB1, ABCB4, ABCB8, SLCO1B1, SLCO1B3, SULT1B1) and modifiers (PPARG, SOD1).

Recent research has reported that these ADME genes participate in tumorigenesis. Depletion of CYP2D6 influences the genes involved in EMT, oncogenesis and immune-related pathways and acts as a biomarker for drug resistance in non-small cell lung cancer (31). CYP8B1 is identified as a prognostic gene for survival analysis in hepatocellular carcinoma (32). ABCB1 is regulated by N6-methyladenosine-induced ERRγ and triggers chemoresistance in cancer cells (33). ABCB8 mediates doxorubicin resistance by protecting the mitochondrial genome in melanoma (34). SLCO1B1 polymorphisms influence the estrogenic response to aromatase inhibitor treatment in breast cancer (35). PPARG activates lipid signaling pathways, and high levels of PPARG/FASN confer a poor prognosis in prostate cancer (36). SOD1 carries a housekeeping function that maintains ROS levels below a threshold, supporting oncogene-dependent proliferation in mammary tumor formation (37). Although some ADME genes were studied in HNSCC (38–41), their biological functions and molecular mechanisms remain unclear.

To explore the biological roles of ADME genes in HNSCC, functional analyses were conducted based on the differentially expressed genes between high-risk and low-risk groups. Interestingly, GO, KEGG and GSVA analyses revealed that many important immune-related functions and pathways were enriched, including lymphocyte activation; T cell activation; leukocyte cell–cell adhesion; Th1, Th2 and Th17 cell differentiation; cell adhesion molecules; inflammatory response; and interferon gamma response. The immune system plays a vital role in the evolution and progression of HNSCC by regulating the tumor microenvironment (42, 43). According to our results, ADME genes may affect metastasis, angiogenesis or growth of HNSCC through immune-related mechanisms.

Subsequently, CIBERSORT, ESTIMATE and other analyses were performed to investigate the immune cell infiltration in the tumor microenvironment between the low- and high-risk groups. The results indicated that the two groups had significantly distinct immune cell infiltration characteristics. Interestingly, we found that infiltration levels of immune cells, especially B cells and T cells, were higher in the low-risk group than that the high-risk group. The above findings suggest that the ADME-related model could estimate immune status of HNSCC patients. The higher infiltration levels of immune-activating cells contribute to a better immune response against cancer cells, which may explain the improved prognosis observed in the low-risk group.

To further evaluate the relationship between the ADME signature and benefit of immunotherapy and chemotherapy, we collected two external cohorts (IMvigor210 and TCGA-SKCM) which received immunotherapy and obtained the drug information of TCGA-HNSC. We found that patients with low-risk scores showed better prognosis and response to immunotherapy and chemotherapy compared to those with high scores. High- risk group were more sensitive to the most used chemotherapy and targeted therapy regimens (cisplatin, docetaxel, paclitaxel, 5-fluorouracil and cetuximab) in HNSCC.

There are several potential limitations in our study. First, only six patients received immunotherapy in the TCGA cohort. Thus, we explored the predictive efficacy of the ADME signature for immunotherapy response in other cancer types. Additionally, studies in different and larger populations are necessary for validating our findings in the future.

In this study, a novel ADME-related prognostic signature of HNSCC was constructed using a TCGA cohort and validated with two validation cohorts. This risk model performed well in predicting survival of patients. Furthermore, we found the relationship of ADME-related score and immune cell infiltration in the TME. The following analysis confirmed the response prediction of the ADME score for immunotherapy and chemotherapy. Our study provides a promising prognostic signature to guide individualized therapy for HNSCC patients. In addition, targeting ADME genes may reverse TME cell infiltration and transform ‘cold tumors’ into ‘hot tumors’, providing a novel insight into potential immunotherapeutic strategies and drug combination strategies for HNSCC.

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

Conceptualization, XT, DW and XW. Data curation, XT, DW and XW. Formal analysis, XT, RLi and YW. Funding acquisition, XT and XW. Software, Supervision, RLi. Validation, FZ and RLv. Writing – original draft, XT and RLi. Writing – review & editing, XW. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Guangdong Basic and Applied Basic Research Foundation (2019A1515110076, 2019A1515011427), the National Natural Science Foundation of China (82003214, 81903134, 82172673), and China Postdoctoral Science Foundation (2019TQ0368).

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/fimmu.2022.905635/full#supplementary-material

References

1. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2020. CA Cancer J Clin (2020) 70:7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bray F, Ferlay J, Soerjomataram I, RL S, LA T, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68:394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Ferlay J, Colombet M, Soerjomataram I, Mathers C, Parkin DM, Pineros M, et al. Estimating the Global Cancer Incidence and Mortality in 2018: GLOBOCAN Sources and Methods. Int J Cancer (2019) 144:1941–53. doi: 10.1002/ijc.31937

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Shield KD, Ferlay J, Jemal A, Sankaranarayanan R, Chaturvedi AK, Bray F, et al. The Global Incidence of Lip, Oral Cavity, and Pharyngeal Cancers by Subsite in 2012. CA Cancer J Clin (2017) 67:51–64. doi: 10.3322/caac.21384

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Tao L, Zhang P, Qin C, SY C, Zhang C, Chen Z, et al. Recent Progresses in the Exploration of Machine Learning Methods as in-Silico ADME Prediction Tools. Adv Drug Deliv Rev (2015) 86:83–100. doi: 10.1016/j.addr.2015.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Gleeson MP, Hersey A, Montanari D, Overington J. Probing the Links Between In Vitro Potency, ADMET and Physicochemical Parameters. Nat Rev Drug Discovery (2011) 10:197–208. doi: 10.1038/nrd3367

CrossRef Full Text | Google Scholar

7. Klein K, Tremmel R, Winter S, Fehr S, Battke F, Scheurenbrand T, et al. A New Panel-Based Next-Generation Sequencing Method for ADME Genes Reveals Novel Associations of Common and Rare Variants With Expression in a Human Liver Cohort. Front Genet (2019) 10:7. doi: 10.3389/fgene.2019.00007

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hu DG, Marri S, RA M, PI M, Meech R. Deregulation of the Genes That Are Involved in Drug Absorption, Distribution, Metabolism, and Excretion in Hepatocellular Carcinoma. J Pharmacol Exp Ther (2019) 368:363–81. doi: 10.1124/jpet.118.255018

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hovelson DH, Xue Z, Zawistowski M, MG E, EC H, SL S, et al. Characterization of ADME Gene Variation in 21 Populations by Exome Sequencing. Pharmacogenet Genomics (2017) 27:89–100. doi: 10.1097/FPC.0000000000000260

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Jittikoon J, Mahasirimongkol S, Charoenyingwattana A, Chaikledkaew U, Tragulpiankit P, Mangmool S, et al. Comparison of Genetic Variation in Drug ADME-Related Genes in Thais With Caucasian, African and Asian HapMap Populations. J Hum Genet (2016) 61:119–27. doi: 10.1038/jhg.2015.115

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Arbitrio M, Di Martino MT, Scionti F, Barbieri V, Pensabene L, Tagliaferri P. Pharmacogenomic Profiling of ADME Gene Variants: Current Challenges and Validation Perspectives. High Throughput (2018) 7:40. doi: 10.3390/ht7040040

CrossRef Full Text | Google Scholar

12. Hu DG, Mackenzie PI, McKinnon RA, Meech R. Genetic Polymorphisms of Human UDP-Glucuronosyltransferase (UGT) Genes and Cancer Risk. Drug Metab Rev (2016) 48:47–69. doi: 10.3109/03602532.2015.1131292

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Drozdzik M, Oswald S. Expression and Regulation of Drug Transporters and Metabolizing Enzymes in the Human Gastrointestinal Tract. Curr Med Chem (2016) 23:4468–89. doi: 10.2174/0929867323666161024154457

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Fisel P, Schaeffeler E, Schwab M. DNA Methylation of ADME Genes. Clin Pharmacol Ther (2016) 99:512–27. doi: 10.1002/cpt.343

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wu KC, Lin CJ. The Regulation of Drug-Metabolizing Enzymes and Membrane Transporters by Inflammation: Evidences in Inflammatory Diseases and Age-Related Disorders. J Food Drug Anal (2019) 27:48–59. doi: 10.1016/j.jfda.2018.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

16. De Iuliis F, Salerno G, Taglieri L, Scarpa S. Are Pharmacogenomic Biomarkers an Effective Tool to Predict Taxane Toxicity and Outcome in Breast Cancer Patients? Literature Review. Cancer Chemother Pharmacol (2015) 76:679–90. doi: 10.1007/s00280-015-2818-4

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Suthandiram S, Gan GG, Zain SM, Bee PC, Lian LH, Chang KM, et al. Effect of Polymorphisms Within Methotrexate Pathway Genes on Methotrexate Toxicity and Plasma Levels in Adults With Hematological Malignancies. Pharmacogenomics (2014) 15:1479–94. doi: 10.2217/pgs.14.97

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhang C, Ma Q, Shi Y, Li X, Wang M, Wang J, et al. A Novel 5-Fluorouracil-Resistant Human Esophageal Squamous Cell Carcinoma Cell Line Eca-109/5-FU With Significant Drug Resistance-Related Characteristics. Oncol Rep (2017) 37:2942–54. doi: 10.3892/or.2017.5539

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hu DG, Mackenzie PI, Nair PC, McKinnon RA, Meech R. The Expression Profiles of ADME Genes in Human Cancers and Their Associations With Clinical Outcomes. Cancers (2020) 12:3369. doi: 10.3390/cancers12113369

CrossRef Full Text | Google Scholar

20. Szklarczyk D, Franceschini A, Kuhn M, Simonovic M, Roth A, Minguez P, et al. The STRING Database in 2011: Functional Interaction Networks of Proteins, Globally Integrated and Scored. Nucleic Acids Res (2011) 39:D561–8. doi: 10.1093/nar/gkq973

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Aran D, Hu Z, Butte AJ. Xcell: Digitally Portraying the Tissue Cellular Heterogeneity Landscape. Genome Biol (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and Genetic Properties of Tumors Associated With Local Immune Cytolytic Activity. Cell (2015) 160:48–61. doi: 10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Lauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, et al. Mutational and Putative Neoantigen Load Predict Clinical Benefit of Adoptive T Cell Therapy in Melanoma. Nat Commun (2017) 8:1738. doi: 10.1038/s41467-017-01460-0

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFbeta Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature (2018) 554:544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Geeleher P, Cox N, Huang RS. Prrophetic: An R Package for Prediction of Clinical Chemotherapeutic Response From Tumor Gene Expression Levels. PloS One (2014) 9:e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, et al. Cancer Incidence and Mortality Worldwide: Sources, Methods and Major Patterns in GLOBOCAN 2012. Int J Cancer (2015) 136:E359–86. doi: 10.1002/ijc.29210

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Leemans CR, Snijders P, Brakenhoff RH. The Molecular Landscape of Head and Neck Cancer. Nat Rev Cancer (2018) 18:269–82. doi: 10.1038/nrc.2018.11

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Genome Atlas NetworkComprehensive Genomic Characterization of Head and Neck Squamous Cell Carcinomas. Nature (2015) 517:576–82. doi: 10.1038/nature14129

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Luo R, Ge C, Xiao X, Song J, Miao S, Tang Y, et al. Identification of Genetic Variations Associated With Drug Resistance in non-Small Cell Lung Cancer Patients Undergoing Systemic Treatment. Brief Bioinform (2021) 22:bbab187. doi: 10.1093/bib/bbab187

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Wang X, Liao X, Yang C, Huang K, Yu T, Yu L, et al. Identification of Prognostic Biomarkers for Patients With Hepatocellular Carcinoma After Hepatectomy. Oncol Rep (2019) 41:1586–602. doi: 10.3892/or.2019.6953

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Chen Z, Wu L, Zhou J, Lin X, Peng Y, Ge L, et al. N6-Methyladenosine-Induced ERRgamma Triggers Chemoresistance of Cancer Cells Through Upregulation of ABCB1 and Metabolic Reprogramming. Theranostics (2020) 10:3382–96. doi: 10.7150/thno.40144

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Elliott AM, Al-Hajj MA. ABCB8 Mediates Doxorubicin Resistance in Melanoma Cells by Protecting the Mitochondrial Genome. Mol Cancer Res (2009) 7:79–87. doi: 10.1158/1541-7786.MCR-08-0235

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Dempsey JM, Kidwell KM, Gersch CL, Pesch AM, Desta Z, Storniolo AM, et al. Effects of SLCO1B1 Polymorphisms on Plasma Estrogen Concentrations in Women With Breast Cancer Receiving Aromatase Inhibitors Exemestane and Letrozole. Pharmacogenomics (2019) 20:571–80. doi: 10.2217/pgs-2019-0020

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ahmad I, Mui E, Galbraith L, Patel R, Tan EH, Salji M, et al. Sleeping Beauty Screen Reveals Pparg Activation in Metastatic Prostate Cancer. Proc Natl Acad Sci U.S.A. (2016) 113:8290–5. doi: 10.1073/pnas.1601571113

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gomez ML, Shah N, Kenny TC, Jenkins EJ, Germain D. SOD1 is Essential for Oncogene-Driven Mammary Tumor Formation But Dispensable for Normal Development and Proliferation. Oncogene (2019) 38:5751–65. doi: 10.1038/s41388-019-0839-x

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Han J, Wang L, Yang Y, Zhang J. Meta-Analyses of the Effect of CYP1A1 and CYP2D6 Polymorphisms on the Risk of Head and Neck Squamous Cell Carcinoma. Oncol Res Treat (2014) 37:406–11. doi: 10.1159/000363428

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wang Y, Xu Y, Hua Q, Jiang Y, Liu P, Zhang W, et al. Novel Prognostic Model Based on Immune Signature for Head and Neck Squamous Cell Carcinoma. BioMed Res Int (2020) 2020:4725314. doi: 10.1155/2020/4725314

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Lian M, Chen J, Shen X, Hou L, Fang J. Pparg may Promote Chemosensitivity of Hypopharyngeal Squamous Cell Carcinoma. Ppar Res (2020) 2020:6452182. doi: 10.1155/2020/6452182

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Liu Y, Zha L, Li B, Zhang L, Yu T, Li L. Correlation Between Superoxide Dismutase 1 and 2 Polymorphisms and Susceptibility to Oral Squamous Cell Carcinoma. Exp Ther Med (2014) 7:171–8. doi: 10.3892/etm.2013.1375

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ferris RL. Immunology and Immunotherapy of Head and Neck Cancer. J Clin Oncol (2015) 33:3293–304. doi: 10.1200/JCO.2015.61.1509

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Miyauchi S, Kim SS, Pang J, Gold KA, Gutkind JS, Califano JA, et al. Immune Modulation of Head and Neck Squamous Cell Carcinoma and the Tumor Microenvironment by Conventional Therapeutics. Clin Cancer Res (2019) 25:4211–23. doi: 10.1158/1078-0432.CCR-18-0871

PubMed Abstract | CrossRef Full Text | Google Scholar