Gene Expression and DNA Methylation in Human Papillomavirus Positive and Negative Head and Neck Squamous Cell Carcinomas

Abstract

High-risk human papillomaviruses (HPV) are important agents, responsible for a large percentage of the 745,000 cases of head and neck squamous cell carcinomas (HNSCC), which were identified worldwide in 2020. In addition to being virally induced, tobacco and heavy alcohol consumption are believed to cause DNA damage contributing to the high number of HNSCC cases. Gene expression and DNA methylation differ between HNSCC based on HPV status. We used publicly available gene expression and DNA methylation profiles from the Cancer Genome Atlas and compared HPV positive and HPV negative HNSCC groups. We used differential gene expression analysis, differential methylation analysis, and a combination of these two analyses to identify the differences. Differential expression analysis identified 1854 differentially expressed genes, including PCNA, TNFRSF14, TRAF1, TRAF2, BCL2, and BIRC3. SYCP2 was identified as one of the top deregulated genes in the differential methylation analysis and in the combined differential expression and methylation analyses. Additionally, pathway and ontology analyses identified the extracellular matrix and receptor interaction pathway as the most altered between HPV negative and HPV positive HNSCC groups. Combining gene expression and DNA methylation can help in elucidating the genes involved in HPV positive HNSCC tumorigenesis, such as SYCP2 and TAF7L.

Keywords: HPV; cancer; gene expression; head and neck cancer; methylation.

Figure 1

Clinical data of the TCGA HNSCC HPV- positive and negative patients. Patients were filtered according to the HPV status using the information regarding p16 expression and in situ hybridization information (only patients with information present were included in the study; n = 114, HPVP = 41, HPVN = 73). (A) Distribution of age at cancer diagnosis between two groups of patients, HPVP and HPVP; (B) distribution of patients in different age groups, and sidewise comparison of age groups and HPV status; (C) representation of male and female patients, HPVP and HPVN; (D) representation of different races, independent of sex or HPV status; (E) distribution of different anatomical sites where cancer originated; (F,G) a closer look at the specific location of HPVN patients (F) and HPVP (G).

Figure 1

Clinical data of the TCGA HNSCC HPV- positive and negative patients. Patients were filtered according to the HPV status using the information regarding p16 expression and in situ hybridization information (only patients with information present were included in the study; n = 114, HPVP = 41, HPVN = 73). (A) Distribution of age at cancer diagnosis between two groups of patients, HPVP and HPVP; (B) distribution of patients in different age groups, and sidewise comparison of age groups and HPV status; (C) representation of male and female patients, HPVP and HPVN; (D) representation of different races, independent of sex or HPV status; (E) distribution of different anatomical sites where cancer originated; (F,G) a closer look at the specific location of HPVN patients (F) and HPVP (G).

Figure 2

Clustering of the TCGA HNSCC HPVP and HPVN samples and genes. PCA and heatmap clustering shows distinct patient groups as we classify them. (A) PCA shows that patients classify in two separate groups for the most part, confirming that separation in HPVP and HPVN groups by p16 expression and in situ hybridization was a valid parameter; (B–E) are heatmaps of the most variable genes (B,C) and most abundant transcripts (D,E) among n = 114 HNSCC samples; (B) shows the top 500 most variable genes, while a closer look at the top 30 most variable genes is shown in (C); top 500 transcripts with the highest mean values are depicted in (D) with a zoomed-in perspective to the top 30 in (E) HPVN samples (labeled in black) and HPVP (in red colored numbers).

Figure 2

Clustering of the TCGA HNSCC HPVP and HPVN samples and genes. PCA and heatmap clustering shows distinct patient groups as we classify them. (A) PCA shows that patients classify in two separate groups for the most part, confirming that separation in HPVP and HPVN groups by p16 expression and in situ hybridization was a valid parameter; (B–E) are heatmaps of the most variable genes (B,C) and most abundant transcripts (D,E) among n = 114 HNSCC samples; (B) shows the top 500 most variable genes, while a closer look at the top 30 most variable genes is shown in (C); top 500 transcripts with the highest mean values are depicted in (D) with a zoomed-in perspective to the top 30 in (E) HPVN samples (labeled in black) and HPVP (in red colored numbers).

Figure 2

Clustering of the TCGA HNSCC HPVP and HPVN samples and genes. PCA and heatmap clustering shows distinct patient groups as we classify them. (A) PCA shows that patients classify in two separate groups for the most part, confirming that separation in HPVP and HPVN groups by p16 expression and in situ hybridization was a valid parameter; (B–E) are heatmaps of the most variable genes (B,C) and most abundant transcripts (D,E) among n = 114 HNSCC samples; (B) shows the top 500 most variable genes, while a closer look at the top 30 most variable genes is shown in (C); top 500 transcripts with the highest mean values are depicted in (D) with a zoomed-in perspective to the top 30 in (E) HPVN samples (labeled in black) and HPVP (in red colored numbers).

Figure 3

Volcano plot depicting differentially expressed genes. Differential expression analysis identified 1854 DEGs, 941 were downregulated and 913 were upregulated between HPVP and HPVN HNSCC patient groups (using HPVN as a baseline for comparison). Purple represents DEGs, blue is statistically significant according to the p-value, green is statistically significant according to the logFC, while grey is not statistically significant.

Figure 4

Differential methylation in HNSCC patients. Represented above is the methylation profile in HNSCC. (A) Mean methylation between HPVN and HPVP HNSCC patient samples; (B) volcano plot showing the hypomethylated genes in green and hypermethylated genes in red. HPVN samples are used as a baseline. We used $\underset{¯}{\beta }$ ≥ 0.25 and p ≤ 10⁻⁵; (C,D) show a Starburst plot that combined differential gene expression data with differential methylation data. HPVN is used as a baseline. We used $\underset{¯}{\beta }$ ≥ 0.25, FDR_expression ≤ 10⁻⁵, FDR_{DNAmethylation} ≤ 10⁻⁵ │logFC│ ≥ 1 in (C) and more stringent parameters $\underset{¯}{\beta }$ ≥ 0.25, FDR_expression ≤ 10⁻⁵, FDR_{DNAmethylation} ≤ 10⁻⁵ │logFC│ ≥ 3 in (D).

Figure 5

Workflow of the TCGA HPV-related HNSCC data. A schematic representation of the stepwise workflow of TCGA data analysis. * Clinical data analysis was performed only on RNA-seq patients’ data, and not on DNA methylation data.