Integrative Genomics and Transcriptomics Analysis Reveals Potential Mechanisms for Favorable Prognosis of Patients with HPV-Positive Head and Neck Carcinomas

Abstract

Patients with HPV-positive head neck squamous cell carcinomas (HNSCC) usually have a better prognosis than the HPV-negative cases while the underlying mechanism remains far from being well understood. We investigated this issue by an integrative analysis of clinically-annotated multi-omics HNSCC data released by the Cancer Genome Atlas. As confirmatory results, we found: (1) Co-occurrence of mutant TP53 and HPV infection was rare; (2) Regardless of HPV status, HNSCCs of wild-type TP53 implied a good survival chance for patients and had fewer genome-wide somatic mutations than those with a mutation burden on the gene. Our analysis further led to some novel observations. They included: (1) The genes involved in "DNA mismatch repair" pathway were up-regulated in HPV-positive tumors compared to normal tissue samples and HPV-negative cases, and thus constituted a strong predictive signature for the identification of HPV infection; (2) HPV infection could disrupt some regulatory miRNA-mRNA correlations operational in the HPV-negative tumors. In light of these results, we proposed a hypothesis for the favorable clinical outcomes of HPV-positive HNSCC patients. That is, the replication of HPV genome and/or its invasion into the genomes of cancer cells may enhance DNA repair mechanisms, which in turn limit the accumulation of lethal somatic mutations.

Figure 1. Kaplan–Meier (K-M) survival curves for patient groups defined by the proposed tumor stratification.

Green (subtype A1, n = 65): tp53&cdkn2a-mut_HPV−. Blue (subtype A2, n = 147): tp53-mut_HPV−. Red (subtype C, n = 42): tp53-wild_HPV−. Yellow (subtype D, n = 40): tp53-wild_HPV+.

Figure 2. Cumulative distributions of somatic mutations present in HNSCC samples.

For each gene catalogue (i.e. HUGO, MutSig or COSMIC), the sample proportion (y) corresponding to a specific mutation burden (x) is calculated by dividing the number of samples with mutation burden ≥x by the total number of samples. A specific mutation burden (x) is quantified by the number of mutated genes.

Figure 3. KEGG pathways over-represented by the significant genes identified in five comparisons.

In this figure, an over-representation relationship (adj.p < 0.01) is highlighted in grey. The pathway clusters are determined by a hierarchical cluster analysis (Manhattan distance and Ward method) with a 0/1 matrix (i.e. M) as the input. In the matrix, rows and columns represent pathways and comparisons, respectively. When the i^th pathway is over represented (adj.p < 0.01) by the significant genes identified from the j^th comparison, the element m_ij of M is 1. Otherwise, it is 0.

Figure 4. Cancer subtype-specific expression profiles of the genes in MMR pathway.

Bisque: Normal tissue. Green: tp53&cdkn2a-mut_HPV−. Blue: tp53-mut_HPV−. Red: tp53-wild_HPV−. Purple: tp53-wild_HPV+. The expression profiles of CDKN2A and TP53 genes are depicted in the last two plots as a reference. The 13 significant genes in the comparison of “CTR-CD” (purple vs red) are marked with stars. Of them, the genes in subsets of (MSH2, MSH6, PCNA, RFC1-5), (EXOL1), (RPA2, POLD1, POLD3) and (LIG1) are involved in mismatch recognition, the excision of mismatched DNA, and DNA re-synthesis and ligation, respectively (http://www.genome.jp/kegg-bin/show_pathway?ko03430).

Figure 5. Evaluation of the expression profile of MMR genes as a prognostic signature for HPV-positive HNSCCs.

Left: Illustration of the proposed SVD-based classification algorithm. SVD-u1 and SVD-u2 represent the first and second left vectors of the Singular Value Decomposition of the transpose of the row-centered expression matrix of 22 MMR genes. The score w_i = f (u_1i, u_2i) for the i^th tumor represents the distance from the corresponding data point to the center of the quarter circle. The coordinates of the center are determined by the minimums of SVD-u1 and SVD-u2. Right: Demonstration of the predictive strength of the score w_i as an independent predictive variable, compared with the other individual predictors, including expression level of CDKN2A, SVD-u1 or SVD-u2.

Figure 6. Identifying (tumor) subtype-specific miRNA-mRNA correlation-network modules by the hierarchical clustering algorithm.

Red: the top positive correlations. Yellow: the top negative correlations. Orange: pseudo or unconsidered correlations. For the sake of clarification, the heatmap for each subtype was further refined by removing the columns corresponding to the miRNAs not involved in the identified modules.

Figure 7. miRNA target site enrichment analysis for six miRNA-mRNA modules.

The plots in the top and middle rows depict the four semi-canonical regulatory modules identified in the tumors of subtype A2, A1 or C, respectively. The plots in the bottom row depict the two major negative-connection modules identified in HPV-positive tumors (i.e. subtype D).

Figure 8. A heuristic model for the potential carcinogenesis of HNSCCs and the genetically defined progressive relationships between the different subtypes.

The figure was created by the first author, Dr. Wensheng Zhang, of this paper.