Abstract
Human papillomavirus (HPV) plays a major role in oncogenesis and circular extrachromosomal DNA (ecDNA) is found in many cancers. However, the relationship between HPV and circular ecDNA in human cancer is not understood.
Forty-four primary tumor tissue samples were obtained from a cohort of patients with HPV-positive oropharynx squamous cell carcinoma (OPSCC). Twenty-eight additional HPV oropharyngeal cancer (HPVOPC) tumors from The Cancer Genome Atlas (TCGA) project were analyzed as a separate validation cohort. Genomic, transcriptomic, proteomic, computational, and functional analyses of HPVOPC were applied to these datasets.
Our analysis revealed circular, oncogenic DNA in nearly all HPVOPC, with circular human and human–viral hybrid ecDNA present in over a third of HPVOPC and viral circular DNA in remaining tumors. Hybrid ecDNA highly express fusion transcripts from HPV promoters and HPV oncogenes linked to downstream human transcripts that drive oncogenic transformation and immune evasion, and splice multiple, diverse human acceptors to a canonical SA880 viral donor site. HPVOPC have high E6*I expression with specific viral oncogene expression pattern related to viral or hybrid ecDNA composition.
Nonchromosomal circular oncogenic DNA is a dominant feature of HPVOPC, revealing an unanticipated link between HPV and ecDNA that leverages the power of extrachromosomal inheritance to drive HPV and somatic oncogene expression.
The incidence of human papillomavirus (HPV)-driven oropharynx cancer has increased in recent years, and in the United States, its frequency is the second-fastest growing of solid tumors. We reveal that human–viral extrachromosomal circular DNA is a strong driver of oncogenic HPV transcription, and that circular HPV genomes may exist in rearranged, noncanonical states in the cancer genome. These circular oncogenic DNA structures enable expression of oncogenic elements in nearly all HPVOPC. The formation and maintenance of oncogenic circular DNA are tasks that are both unique to cancer, and thus represent potential therapeutic targets. Furthermore, such elements are not constricted by the rules of Mendelian inheritance, and they enable more rapid tumor evolution, even in the face of targeted therapies. As a result, detecting and profiling circular DNA in cancers presents an important potential prognostic indicator for patient outcome.
Introduction
Oropharynx cancer has become the second-fastest growing cause of cancer death and the third-fastest growing in frequency among solid organ cancers in the United States (1, 2). The main histology is oropharynx squamous cell carcinoma (OPSCC), which is driven by high-risk human papillomavirus (HPV) type 16 (3, 4). The annual number of HPV-related oropharynx carcinoma (HPVOPC) cases has already surpassed the number of cervical cancer cases in the United States in 2009, and by 2030, approximately half of all head and neck cancers in the United States are predicted to be HPV related (3). Although HPVOPC exhibits an improved clinical prognosis compared with HPV-negative OPSCC, 20% to 35% of tumors exhibit an aggressive course despite multimodality therapy (5). A major hurdle to understanding HPV-mediated oncogenesis is an incomplete understanding of the role of viral and viral–human hybrid transcripts and viral–human DNA integration.
Extrachromosomal DNA (ecDNA) has recently been shown to play a critical role in human cancer (6–9). Because of its nonchromosomal mechanism of inheritance, ecDNA can drive high copy number while promoting intratumoral heterogeneity, promoting accelerated tumor evolution and drug resistance (6, 10). Moreover, chromatin rewiring on ecDNA allows for higher accessibility and increased expression of oncogenes (6, 7, 11). More recent reports have conjectured that hybrid human–virus ecDNA formation could be a possible mechanism for increased copy number of the HPV oncogenes E6 and E7 (12–16). Our previous studies have demonstrated that the HPVOPC cell line UPCI:SCC090 features hybrid human–viral circular ecDNA containing FOXE1 and HPV16 through conventional and long-read whole-genome sequencing (WGS), and which we verified in vitro using fluorescence in situ hybridization (FISH; ref. 11).
Given these data, we hypothesized that the genetic structure and viral gene expression in primary HPVOPC as well as the expression of human viral hybrid transcripts may be related to ecDNA. We combined WGS, conventional RNA sequencing (RNA-seq), and long-read RNA-seq to analyze HPV and human viral hybrid genomic and transcriptomic structure in the context of HPVOPC (17). Analysis of ecDNA and associated transcript structure clarified HPV transcript structure and the role of viral, human, and hybrid ecDNA in enhancing expression of diverse and oncogenic viral, human, and hybrid transcripts with functional validation.
Materials and Methods
Patient samples
Forty-four primary tumor tissue samples were obtained from a cohort of patients with HPV-positive OPSCC from the Johns Hopkins Tissue Core [Institutional Review Board (IRB) protocol #NA_00-36235] and Moores Cancer Center Biorepository and Tissue technology shared at University of California, San Diego Human Research Protections Program (IRB-approved protocol HRPP# 181755). Pathology of the primary tumors confirmed by two independent pathologists and tumor tissue was microdissected to yield at least 80% tumor purity. HPV tumor status was determined by in situ hybridization or p16 immunohistochemistry. In equivocal cases, HPV16 E6 and E7 viral oncoproteins were detected via PCR for confirmation. WGS using paired-end Illumina sequencing along with conventional RNA-seq was acquired for 37 samples. Full clinical characteristics of the cohort are presented in Supplementary Table S1.
WGS
DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen) for high-quality extraction per the manufacturer's instructions. DNA samples from tumor were quantified using a Qubit (Thermo Fisher Scientific). Greater than 1 μg of each sample was prepared using a sonication-based library construction and enrichment method per the Beijing Genomics Institute (BGI) as described previously (18).
DNA was isolated from 0.35-mm thick frozen tissue cuts digested in 1% SDS (Sigma-Aldrich) and 50 μg/mL proteinase K (Invitrogen) solution at 48°C for 48 hours. The DNA was purified by phenol–chloroform extraction and ethanol precipitation. DNA was resuspended in LoTE buffer, and the DNA concentration was quantified using the NanoDrop spectrophotometer. Sequencing was performed with the Illumina Hiseq Xten 151PE strategy with 350-bp insert library. The pipeline steps included preparation of HPV reference genome file, performance of quality control on BAM files, extraction of unmapped read pairs, conversion of unmapped read pairs to FASTQ format, and alignment of unmapped read pairs to the HPV reference genomes (accession number: AY686584.1).
RNA preparation
Frozen tissue specimens were cut into 0.35-mm thick sections and RNA was extracted according to the Qiagen RNeasy Plus Mini Kit (Qiagen). RNA concentration was verified using a NanoDrop spectrophotometer (Thermo Fisher Scientific). The absorbance ratio of 260 to 280 nm was used to verify adequate quality, defined as >1.8. An RNA integrity number (RIN) of 7.0 or greater was required for quality assessment. RNA from all eleven tumors passed quality assessment.
cDNA library preparation, long-read RNA-seq, and alignment to HPV16 genome
Briefly, the RNA was extracted from 0.35-mm thick frozen tissue sections and a stranded RNA library was prepared using the Illumina TruSeq stranded total RNA seq poly A+ Gold Kit following the manufacturer's recommendations. Long-read RNA-seq of full-length transcripts was performed on 2 nonintegrated and 3 integrated tumors according to the PacBio Iso-Seq pipeline. Briefly, 500 nanograms of purified RNA was used to prepare cDNA using the Clontech SMARTer PCR cDNA Synthesis Kit and cDNA was then repaired. Large-scale PCR was performed using the Blue-Pippin size selection system for three sized cDNA libraries (<1.5 kb, 1.5–2.5 kb, >2.5 kb). SMRTbell templates were then purified and sequenced on the PacBio SMRT Sequencing platform. The general SMARTer IIA oligonucleotide was used to anneal to the polyA tail of transcripts during cDNA sample preparation. Junction-spanning reads covered by fewer than 5 reads were dropped from analysis.
The Spliced Transcripts Alignment to a Reference (STAR) software was used to align long-read RNA-seq reads to the HPV16 reference genome (GenBank: AY686584.1; ref. 19). Full-length transcripts were visualized with IGV for confirmation, and erroneously mapping transcripts were removed from analysis.
Short-read RNA-seq alignment and analysis
Standard (short-read) RNA-seq was performed as described previously (20). A ribosomal RNA reduction was performed and the purified RNA was fragmented, then converted to double-stranded cDNA, and the cDNA was 3′ adenylated and ligated with barcode adapters. The library was then enriched using PCR and AMPure XP bead purification. Sequencing was then performed using the HiSeq 2500 platform sequencer (Illumina), and the TruSeq Cluster Kit for 2 × 100 bp sequencing. The reads were trimmed to remove adapter sequences and low-quality reads using Trim Galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/).
The RNA sequences were aligned to the HPV16 genome (GenBank: AY686584.1) and hg19 assembly using MapSplice2 version 2.0.1.9. Integration and expression of HPV genes were identified by taking reads in RNA-seq data aligned to a combined database of human reference genome and high-risk HPV16, HPV33, and HPV35 reference genomes using MapSplice (https://github.com/favorov/viruses-in-sequencing; ref. 21). MapSplice was run with the default command line arguments. RNA-seq reads were aligned to the HPV16 genome and reads spanning canonical HPV16 splice junctions were extracted. For splicing analysis, RNA-seq reads were normalized by dividing by the total number of junction-spanning reads in the sample. Junction-spanning reads were discarded if the junction constituted <1% of all junctions in the sample.
Quantification of RNA-seq expression was performed following the HISAT, StringTie, and Ballgown pipeline (22). Briefly, HISAT2 was used to align the RNA-seq reads to the hg19+HPV reference genomes. StringTie was run on each individual alignment to identify the assembled transcripts within the sample. All identified transcripts across all samples were merged using StringTie to create a consistent set of reference transcripts across the entire dataset. Abundances for each transcript for each sample was re-estimated using StringTie, and Ballgown was run on the resulting output to obtain read counts, coverage, and expression data across all samples.
Detection of focal amplifications with AmpliconArchitect
WGS reads were aligned to hg19 and 337 viral genomes using BWA-MEM (23). AmpliconArchitect (AA) seed detection was performed with CNVKit (24). Copy number amplification regions matching low complexity, repetitive, or poorly mappable genomic regions were filtered using the AA database. Some unfiltered regions corresponding to repetitive genomic regions still existed after this step and were shared across multiple samples. We removed any such regions existing in 10% or more of samples. Remaining copy number seed intervals larger than 10 kbp and with estimated copy number (CN) > 4.3 were used as input to AA. Resulting amplicons generated by AA were examined for the presence of breakpoint graph cycles containing solely amplified human DNA (human ecDNA) and viral breakpoint edges linking HPV to an amplified cyclic human DNA structure (human–viral, ‘hybrid’, ecDNA).
Integrative Genome Viewer confirmation
Putative full-length transcripts and splice junctions were visualized using Integrative Genome Viewer (IGV/Broad Institute, version 2.4.5; ref. 25). From long-read and short-read RNA-seq data, BAM files were loaded into IGV and visualized at the start and end of each junction. Long-read isoforms from IGV were individually verified and isoforms with mapping error were removed from analysis.
Integration, fusion, and splicing analysis
WGS and RNA-seq reads were aligned to hg19 and viral genomes using BWA-MEM. ViFi was run on each aligned BAM file to detect viral integration and transcription fusion location. Several transcription fusion events did not have a proximal viral integration event. Closer inspection of these fusion transcripts revealed that they had much lower support (mean 74 RNA-seq reads supporting the fusion event) compared to fusion transcripts that were proximal to a viral integration event (mean 955 RNA-seq reads supporting the fusion event). As such, we removed all fusion transcripts without supporting genomic integration.
Samples were classified according to the presence of viral integration and transcript fusion events. Samples containing a viral integration were classified as Hybrid-DNA. Samples that also contained a fusion transcript were classified as Hybrid-RNA.
In-depth splicing analysis was performed by taking the reads aligned to the HPV16 reference genome and counting the number of splice events detected by the alignment denoted by HPV16 donor and acceptor site pair (i.e., SDx-SAy). For the cases in which the splice junction started with an HPV16 donor site and ended in the human genome, it was denoted by the HPV16 donor site to human splice event (i.e., SDx-Human). From these data, we generated a splicing matrix where each row was a sample and each column was a splicing event, and the entries in the matrix were the total number of times that splicing event was observed within the sample. We performed principal component analysis on the splicing matrix in order to examine how the samples clustered.
Unsupervised hierarchical clustering and splicing analysis
Conventional RNA-seq reads spanning canonical HPV16 splice donor and acceptor sites (SD226, SA409, SA526, SA742, SD880, SD1302, SA2582, SA2709, SA3358, SD3632, and SA5639) from the institutional and TCGA cohorts were extracted and normalized to the number of reads specific to that tumor. Tumors with fewer than 500 mapped RNA-seq reads were excluded from inclusion in the heatmap (T14 – 2 41, T12 – 4, T30 – 2). RNA-seq reads spanning HPV16 splice donor to human acceptor sites were also extracted and normalized to the number of reads specific to that tumor. To calculate the proportion of E6 gene reads with truncation to E6*I in each tumor, the following formula was applied to each tumor:
We also quantified the relative frequency of HPV16 SD880 to human splicing events with the following:
Insertion of HPV16 into Human Genome Analysis
Integration breakpoints and intragenomic viral breakpoints were identified with ViFi and AA. For DNA-based breakpoints both ViFi and AA were used to identify integration sites. To detect intragenomic viral breakpoints, AA alone was used. For RNA-seq data ViFi and MapSplice2 were used to identify splicing and human–viral chimeric sequences.
Analysis of noncanonical HPV16 structures
Canonical and noncanonical circular viral genome structure status was determined by AA analysis. Tumor samples that did not have hybrid ecDNA were classified as noncanonical if they contained a cyclic AA graph decomposition and >100 bp of rearranged genomic content (including indels), while canonical circular status was assigned if no such large rearrangements were present and cyclic AA graph decomposition of virus was present.
Quantification of splice acceptor cluster ranges for hybrid splicing events
For measurement of splice cluster ranges, splice clusters were defined via k-means clustering. Any group of donor, sample, and chromosome with fewer than 10 samples was excluded from consideration. Clustering was repeated 30 times at a given threshold to account for random seeding, with the optimal clustering at a given k threshold determined via silhouette score. The number of clusters was increased until a loss in performance was observed, and the number of clusters was confirmed by visual inspection. Clusters containing fewer than five observations were then filtered out, and the range of splicing sites was then computed on each remaining cluster.
Hybrid RNA-seq reads from both Institutional and TCGA cohorts were extracted and mapped to specific viral splice donor and human acceptor sites for each tumor. The location of the human acceptor (chromosome and nucleotide) for each read was then determined using split reads, which were identified as having a primary alignment to HPV and a secondary alignment to the human genome. Histograms were created to map the distribution of splice acceptor sites for a given HPV16 donor (i.e. SD226, SD880, etc.) across the human genome for each tumor. Samples with viral read counts <10 were removed from analysis.
Functional studies
Proliferation of HCT116 and NOKSI cells was investigated in the presence of empty vector (negative control), as well as E6E7 (positive controls), and parent/daughter constructs. Cells were seeded in 96-well plates at a density of 3,000 cells/well for NOKSI and 4,000 cells/well for HCT116. Individual vectors composed of daughter constructs were then transfected by X-tremeGENE9 (Roche). Proliferation was measured as a ratio of relative absorbance 2 days after transfection versus day of transfection.
The effect of FOXE1 siRNA was investigated on cell line SCC090. Cells were seeded at a density of 2,000 cells/well. SiFOXE1 (Santa Cruz Biotechnology) was added at a concentration of 10 nmol/L. Percent viability was measured using % viability = (absorbance of siRNA)/(absorbance of vehicle) × 100. For proliferation experiments, each datapoint is the average of five replicates with standard error represented by error bars, and all experiments were repeated at least three times demonstrating consistent results.
HCT116 and SCC090 cells were obtained from ATCC; NOKSI was provided as a gift by the Silvio Gutkind Lab (Department of Pharmacology, University of California, San Diego). Cell lines were used for between 4 to 20 passages after thawing from frozen stock. Mycoplasma testing was conducted monthly using the MycoAlert-Plus Mycoplasma Detection Kit (Lonza).
Data and materials availability
AA amplicon visualizations, breakpoint graphs, and amplicon decompositions have been uploaded to FigShare (doi: 10.6084/m9.figshare.13520087; URL: dx.doi.org/10.6084/m9.figshare.13520087).
Plasmids generated in this study are available upon request.
The codebases utilized in this study are available at:
Results
Forty-four primary tumors were acquired from a cohort of patients with HPV-positive OPSCC. HPV tumor status was determined by in situ hybridization or p16 immunohistochemistry (see Supplementary Table S1 for clinical information). In equivocal cases, HPV16 E6 and E7 viral DNA were detected via PCR for confirmation. Forty of 44 samples were HPV16-positive, three were HPV33-positive, and one was HPV35-positive. WGS using paired-end Illumina reads at mean coverage of 30× along with RNA-seq was acquired for 38 HPV16 samples (20), and long-read RNA-seq of full-length transcripts was generated for 5 samples using PacBio Iso-Seq technology (26). Twenty-eight additional HPVOPC tumors from TCGA project were analyzed as a separate validation cohort. WGS and RNA-seq data were mapped to the hg19 reference and analyzed using ViFi (27).
ecDNA that carry oncogenes are common in HPVOPC
As we had previously demonstrated the presence of ecDNA in an HPVOPC cell line (11), we hypothesized that ecDNA may be present in primary HPVOPC. A recently developed method, AA analyzes whole genome sequences to predict ecDNA with 85% precision and 83% sensitivity, as well as reconstruct the fine structure of the amplicons (8). We applied AA to the 28 samples from the HPVOPC cohort (Supplementary Fig. S1; ref. 11). Remarkably, we found six hybrid viral–human ecDNA, and another six with human-only ecDNA, (one tumor exhibited both hybrid and human ecDNA – T14; Fig. 1A) in our institutional cohort. In addition, 18 tumors contained only HPV viral circular DNA (vcDNA). One tumor was not classified due to low viral copy number (T26). HPV vcDNA was present in an intact form including a complete, nonrearranged (canonical) form in 16 tumors. Interestingly, another 15 tumors contained a noncanonical truncated vcDNA with deletions mostly in the L1 and L2 region (Supplementary Fig. S2), suggesting that a significant fraction of HPVOPC tumors contain HPV genomes, which have undergone substantial genomic rearrangement prior to enrichment in copy number.
To test that the prediction of ecDNA in HNSCC samples was not specific to the institutional cohort, we analyzed WGS samples from the HNSC data in TCGA. Ten samples contained viral–human hybrid ecDNA and eight contained human-only ecDNA (four tumors exhibited both; Fig. 1B). vcDNA was also prevalent, with two tumors containing canonical vcDNA while 13 other tumors contained noncanonical truncated HPV vcDNA. One tumor was not classified due to low viral copy number (CV-7406).
ecDNA do not carry centromeres and therefore segregate independently, allowing tumors to rapidly modulate copy numbers of genes on ecDNA, specifically when the genes provide a growth or proliferative advantage (28). Consistent with this hypothesis, seven of the ecDNA+ tumors (both hybrid and human-only) in the institutional cohort carried oncogenic protein-coding genes or ncRNA (Supplementary Table S1), including many known oncogenes: EGFR, SEC61G, VOPP1, and VSTM2A on chr7 in T1; DUSP4 and KIF138 on chr8 as also CD93 in T29; and CST1 and THBD on chr20 in T14. The TCGA cohort revealed similar findings, with protein-coding genes found in five hybrid ecDNA-carrying samples, and three human-only ecDNA samples (Supplementary Table S2). The structures were largely sample specific. However, we did observe 2 samples with ecDNA segments on chromosome 11 carrying six oncogenes (ANO1, CTTN, FADD, MIR548K, PPFIA1, SHANK2). Tumor TCGA-CQ-5323 contained an ecDNA with 13 cancer-associated genes, including ANO1, CCND1, CPT1A, CTTN, TRPC4AP, FADD, IGHMP2, MIR548K, MRPL21, ORAOV1, PPFIA1, and SHANK2 (Supplementary Table S2). In TCGA-CV-5443, a hybrid ecDNA amplicon containing and amplifying the immune regulating ligand PD-L1 was identified (Fig. 1D). PD-L1 amplification occurs in a subset of lung, kidney, bladder, and head and neck cancers. The PD-1 checkpoint is the most common immunotherapeutic target in solid tumors currently and PD-1/PD-L1–directed immunotherapy has gained FDA approval for first-line treatment of unresectable or metastatic head and neck cancer (29).
Overlaying RNA-seq data on to the hybrid ecDNA structure in the institutional cohort showed hybrid RNA combining HPV16 E6, E7, E1, E4, L1, and L2 with a multitude of human sequences containing genes EGFL7, TBCD (chr17; T1), SOX2-OT (chr3; T19), TTC33 (chr5; T41; see Fig. 1C), PVT1 (chr8; T14), LINC01363 (chr1; T47), and TBC1D16 (chr17; T49), As one interesting example, we identified an ecDNA in T41 amplicon that connected viral promoter to multiple exons in TTC33 (Fig. 1C). Hybrid splicing was additionally confirmed using long-read Iso-Seq data (Fig. 1C) consistent with the circular hybrid ecDNA structure. The TTC33 gene (tetraticopeptide repeat domain 33) has been implicated as an mRNA chimera in breast, ovarian, stomach, colon, kidney, and uterine cancer (30). These interesting patterns suggested a possible rewiring of the regulatory circuitry in hybrid ecDNA.
Hybrid ecDNA is associated with increased human gene expression
Prior data have shown that ecDNA mediate increased expression of oncogenic human transcripts contained within ecDNA structure. To examine the effect of viral DNA genomic integration and viral integration in ecDNA, we examined expression of both viral and human transcripts in these contexts. For each gene on an ecDNA, we computed the ratio of its expression (in FPKM units) in the target sample to its mean expression value in all samples where the gene was not on an ecDNA (Materials and Methods) and called it the “FPKM-ratio”. Genes associated with the 38 ecDNA amplicons in the institutional cohort and the TCGA cohort were upregulated nearly 150×, with a mean FPKM-ratio of 149.8 (SD 1,015) and median 4.26 (IQR 1.67–8.92; Fig. 2A and B; Supplementary Table S2). Thirty-three of 38 (86.8%) of these genes were oncogenes or associated with oncogenic phenotypes. For example, oncogene TNFSF4 on chromosome 1 in TCGA-CR-6473 was upregulated to FPKM-ratio of 116.58. TNFSF4 has been reported to be upregulated in brain metastases (31). EGFR on chromosome 7 in T48 was also upregulated with FPKM-ratio 30.7. PVT1 transcripts are found at a high level in lung cancer and contribute to VEGFC expression (32). CD274/PDL1, associated with immune checkpoint activation, was one of the most upregulated genes by tumor TCGA-CV-5443 in the TCGA cohort with FPKM-ratio 24.4.
Importantly, human genes associated with hybrid transcriptomes showed increased expression for both the institutional and TCGA cohorts (Fig. 2C). The increased expression was most pronounced in tumors with hybrid transcripts, but also increased expression was noted for all genes located on hybrid ecDNA. To further explore this relationship, we assembled and annotated human reads overlapping with human viral splice junctions, and spatially defined expression along genomic fusions. We noted that strand-specific expression of human genes downstream of viral sequences in fusion DNA structures could exceed 30-fold compared to surrounding genes. For example, T41 shows dramatic increase in TTC33 expression downstream of HPV-human hybrid sequence (Fig. 2D). A similar phenomenon was noted in the context of other hybrid ecDNA structures (Supplementary Fig. S3).
Viral transcripts show diverse isoforms in hybrid ecDNA and hybrid transcript expression
Unsupervised hierarchical clustering based on frequency of HPV16 splicing junctions in mapped RNA reads showed a number of distinct patterns. First, we observed that SD226-SA409 represent a significant portion of junction-spanning reads in the HPV16 cohort, occurring in every sample with high frequency 36.3%, range 18.3% to 63.5%, and SD 9.7% of junction-spanning reads. The use of this junction creates a shortened form of E6, called E6*I, which results in a premature stop codon (33). The mean fraction of reads demonstrating truncation of E6 to E6*I was 81.1% (SD 12.0%; min 44%.5, max 94.7%) across the cohort of HPV16-positive tumors (see Materials and Methods), demonstrating that E6*I is more commonly expressed than full-length E6 across all HPVOPC. Our results are consistent with previous results identifying E6*I in cervical dysplasia, cervical cancer, and HPVOPC (33–37). We calculated the proportion of E6 transcripts that were E6*I in the institutional cohort and found that tumors with either form of ecDNA had reduced E6*I production compared with the non-ecDNA tumors [0.72 (0.15) vs. 0.82 (0.09); P = 0.0197 by t test; mean (SD)]. Human ecDNA tumors had reduced E6*I production compared with the nonhuman ecDNA tumors [0.69 (0.13) vs. 0.81 (0.11); P = 0.0288 by t test; mean (SD)]. Hybrid ecDNA tumors did not have reduced E6*I production compared with the nonhuman ecDNA tumors [0.75 (0.18) vs. 0.80 (0.11); P = 0.42 by t test; mean (SD)]. To validate these findings, we examined the TCGA cohort of 28 HPVOPC tumors. Similar to the institutional cohort, the majority of E6 was truncated to E6*I (mean 0.89, SD 0.05, median 0.89, IQR 0.87–0.92), although we did not detect a difference in proportion of E6*I based on ecDNA status. This does confirm that, contrary to the classic model of HPV carcinogenesis, E6*I, rather than E6, is the most common viral transcript in HPVOPC.
Second, although splicing of the 5′ SD880 splice donor site to the 3′ SA3358 site had been described as the most frequent splicing event in HPV16 cervical cancers and in cell lines (38–40), we found that 10 of 37 tumors (27%) preferentially spliced from SD880 to a human splice acceptor site instead of the canonical SA3358 (Fig. 1A). In these 10 tumors, 47% (mean; SD = 11%) of reads spanning the SD880 splice site spliced to a human locus, and only 3% (mean; SD = 3%) of SD880 reads spliced to the canonical SA3358 receptor (P < 0.001 by t test). This questions previous findings that efficient usage of SA3358 is necessary for production of E6, E7, E4, E5, L1, and possibly L2 proteins (38). However, we did not detect preferential splicing of SD880 to a human acceptor to be associated with ecDNA status (40% of any ecDNA tumors preferentially spliced SD880 to human versus 22% non-ecDNA; P = 0.28; 40% of hybrid ecDNA tumors preferentially spliced SD880 to human versus 25% nonhybrid ecDNA; P = 0.482; 40% of human ecDNA tumors preferentially spliced SD880 to human versus 25% non-ecDNA; P = 0.482). We also analyzed splicing patterns of SD880 in TCGA and found that 10 of 28 (35%) preferentially spliced to a human splice acceptor rather than canonical SA3358. Tumors with either form of ecDNA were significantly more likely to splice to a human acceptor from SD880 (8/14; 57% vs. 2/14; 14%, P = 0.018).
We also noted a strong association of splicing patterns depending on hybrid DNA or RNA status in the institutional cohort, irrespective of ecDNA (Fig. 1A and B). Hybrid-DNA tumors (n = 18) exhibited a greater fraction of RNA reads covering SD880-human junctions than nonhybrid-DNA tumors [n = 19; 0.26 (0.25) vs. 0.02 (0.01); P < 0.001; and fewer RNA reads covering the SD880-SA3358 junction [0.21 (0.21) vs. 0.42 (0.12); P < 0.001 by t test]. Similarly, hybrid-RNA tumors (n = 12) exhibited a greater fraction of RNA reads covering SD880–human junctions than nonhybrid-RNA tumors [n = 25; 0.38 (0.21) vs. 0.01 (0.05); P < 0.001 by t test] and fewer RNA reads covering the SD880–SA3358 junction [0.09 (0.12) vs. 0.43 (0.12); P < 0.001 by t test]. Principal component analysis of splicing patterns in hybrid-RNA tumors also demonstrated a distinct subset based on splicing signature, in which HPV splicing pattern most closely relates to the presence of viral–human hybrid transcripts (Fig. 3A; Supplementary Fig. S4). We observed selective enrichment of E6/E7 regions in both WGS and RNA data, which is pronounced in hybrid-RNA tumors compared to nonhybrid-RNA tumors (Fig. 3B), as well as depletion of L2 in hybrid samples.
We defined aggregate splice donors in HPV and hybrid transcripts, and noted that SD226, SD880, and other known splice sites in the early region of HPV16 genome were strongly preserved and limited to a single donor canonical nucleotide. However, splice acceptors in hybrid transcripts showed broad variation. In the institutional cohort, the mean variation of the SD880 splice acceptor was 11,060 nucleotides (SD 37,217), but tighter for SD226 (mean 885, SD 2,290) and SD1302 (mean 48, SD 58). In TCGA, the degree of variation was similar (Fig. 3C). We also noted that splicing patterns varied depending on hybrid DNA or RNA status in TCGA, irrespective of ecDNA. Sixty-four percent (18/28) exhibited a hybrid genome and 53% (15/28) exhibited hybrid transcriptomes, and hybrid DNA was a prerequisite for hybrid RNA (P < 0.001). Similarly, SD226-SA409 represented a significant portion 40.0% (range 19.1%–73.2%, SD 17.1%) of junction-spanning reads; TCGA tumors also preferentially expressed E6*I compared with full-length E6 [89.2% E6*I reads (SD 5.3%)], and hybrid genome tumors (n = 18) exhibited significantly higher percentage of splicing events from SD880 to a human locus [11.0 (10.0%) vs. 2.9 (9.0%); P = 0.0508] and fewer reads covering the SD880-SA3358 junction [22.2 (23.9%) vs. 50.0 (19.2%); P = 0.005]. Hybrid transcriptome tumors (n = 15) also exhibited significantly higher fraction of splicing events from SD880 to a human locus [13.2 (10.0%) vs. 2.2 (8.0%); P = 0.0041] and fewer reads covering the SD880-SA3358 junction [16.5 (21%) vs. 50.0 (17.4%); P = 0.001; Supplementary Tables S3 and S4].
Novel HPV transcript structures related to ecDNA are found in HPVOPC
Long-read polyA RNA-seq provides direct sequencing of full-length transcripts and can avoid artifacts introduced by short-read transcript assembly. To provide a more precise understanding of HPV transcripts in HPVOPC, we performed long-read RNA whole-genome polyA transcript sequencing on a subset of two tumors without hybrid transcripts (T2 and T38), one human ecDNA tumor with hybrid transcript expression (T19), and one hybrid ecDNA tumor with hybrid transcripts (T45; Fig. 4A; Supplementary Table S5).
Long-read sequencing confirmed a divergent transcript structure of hybrid ecDNA and vcDNA HPV transcripts. For example, T2 and T38 vcDNA (nonhybrid DNA/transcriptome tumors) exhibited 0% of conventional RNA-seq reads covering SD880 spliced to a human junction. In contrast, T19 (human ecDNA/hybrid transcripts) exhibited 19.5% and T45 (hybrid ecDNA/hybrid-RNA tumors) 100% of conventional RNA-seq reads of SD880 to be spliced to a human splice acceptor. Of interest, even though long-read RNA-seq is not quantitative in nature, we observed that both T19, a nonhybrid ecDNA tumor that expressed hybrid transcripts, and T45, a hybrid ecDNA/hybrid transcript tumor, essentially displayed no full-length transcripts that mapped to HPV alone; the transcripts were all hybrid. In T19, 36/40,474 (0.1%) long reads mapped to HPV16 and in hybrid ecDNA tumor T45 17/55,814 (0.03%) long reads mapped to HPV16. Conversely, the two nonhybrid vcDNA tumors T2 and T38 carried more full-length HPV-only transcripts. In T2, 515/19,220 (2.6%) long reads mapped to HPV and in T38 405/38,026 (1.1%) long reads mapped to HPV (P < 0.001 between nonhybrid and hybrid tumors).
Long-read data confirmed the presence of canonical splicing events seen in conventional RNA sequencing. The most common full-length transcript in nonhybrid tumors was 1,476 nt long, beginning at the p97 promoter with splicing at SD226-SA409 and SD880-SA3358 extending to the early polyA tail, with coding potential for the E6 oncoprotein variant E6*I defined by SD226-SA409, full-length E7, full-length E4, and full-length E5 (Fig. 4B and C). This transcript was observed in 55% of full-length reads mapped to HPV16 in T2 and 65% of full-length reads mapped to HPV16 in T38, and included splice junctions commonly observed in RNA-seq data, including SD226-SA409.
The long-read results and the RNA-seq data from the institutional and TCGA cohorts suggested that the predominant form of E6 in full-length transcripts found in HPVOPC was not the full-length E6 isoform, but truncated version E6*I defined by SD226-SA409. Additional isoforms E6*II, and E6*III (41, 42), were also noted in long-read transcripts. Full-length E7 was also common and present in the majority of HPV coding transcripts (91%). Finally, E1⁁E4, which is the result of SD880-SA3358 splicing was observed in 11/23 (48%) of distinct full-length isoforms.
ecDNA hybrid HPV transcripts are functionally active
We selected a tumor (T41) for functional characterization of transcripts due to presence of fusion RNA reads, in addition to the fusion WGS reads and observation of up to 512× increased expression of segments associated on this ecDNA structure. After confirming the presence of AA-predicted hybrid junctions using RT-PCR and sequencing (Supplementary Fig. S5), we cloned the entire transcript (E6-E7-E1-TTC33*-E5*) as well as daughter constructs into a pcDNA 3.1(+)-myc-His A vector (Genscript, Inc.; Fig. 5A and B), as well as the most common full-length HPV16 in the form of component HPV gene transcripts, into the same backbone. To explore the functional effects of these transcripts, we transfected these constructs into HPV null diploid HCT116 (p53 and Rb wt MMR-deficient colorectal carcinoma) cells as well as an HPV-null normal oral keratinocyte NOKSI (spontaneously immortalized oral keratinocyte) cell lines that respond to HPV E6/E7 gene expression with enhanced proliferation, to provide an assessment of the effects of transcripts from this primary tumor (43, 44). In HCT116 cells, E6*I, E6E7, E6*I/E7, E6*I/E7/E4/E5, E7, TTC33, E6E7E1, and the entire hybrid transcript from T41 induced significant growth compared with the empty vector (P = 0.02, 0.02, 6 × 10−3, 0.03, 3 × 10−4, 2 × 10−4, 7 × 10−3, 0.01, respectively; Student t test; Fig. 5C). Similarly, in NOKSI cells, E6*I, E6E7, E6E7E1, and the entire intact hybrid transcript from T41 increased proliferation (P = 0.02, 0.01, 0.03, and 0.04, respectively; Student t test; Fig. 5D).
We have previously reported on the presence of ecDNA in an HPVOPC cell line UPCI;SCC090, demonstrating reconstruction of a complex hybrid structure (>100 kbp) containing the oncogene FOXE1 as well as highly expressing HPV16 sequence, and shown the presence of FOXE1 in both chromosomal HSRs as well as in ecDNA using FISH probes in metaphase imaging (11). To examine the functional contribution of FOXE1 in an ecDNA context, we treated SCC090 cells with siRNA for FOXE1 demonstrating significant inhibition of growth, indicating that overexpression of FOXE1 via ecDNA-mediated mechanisms is a major driver of growth in SCC090 cells (Fig. 5E).
To define the potential for hybrid ecDNA in HPVOPC to drive protein expression related to immune evasion, data derived via reverse-phase protein arrays (RPPA) corresponding to head and neck squamous cell carcinoma samples included in TCGA were extracted from The Cancer Proteome Atlas (TCPA; ref. 45). Nine tumor samples identified as being HPV+ in the TCGA cohort had RPPA-derived expression data available in TCPA. Mean PD-L1 protein expression in the tumor with PD-L1 present in a hybrid ecDNA structure (TCGA-CV-5443) was increased 7.6× relative to the mean of the other eight TCGA tumors (one sample t test, P < 0.001; Fig. 5F).
We analyzed the presence of hybrid DNA, hybrid RNA, hybrid ecDNA, and human circular ecDNA as compared with clinical features in the institutional cohort and in TCGA. In the institutional cohort there was no significant correlation with aggressive pathologic features (perineural invasion, lymphovascular invasion, or extranodal extension), poorly differentiated tumors, and smoking status. Hybrid-RNA status was more likely in patients with a drinking history (11/11 vs. 18/26; P = 0.038), and in patients who had a drinking history and had more than 10 pack-years of smoking (6/11 vs. 3/23, P = 0.010). There was no association with the above groups and AJCC v7 staging. We then compared overall and recurrence-free survival based on presence of hybrid DNA, hybrid RNA, hybrid circular ecDNA, and human circular ecDNA. There were 10 deaths, with the median time to death being 121 months. There were 11 recurrences and/or deaths, with the median time to event being 104 months. We found no significant difference in overall or recurrence-free survival. In a multivariable proportional hazards model adjusting for age, stage, and treatment modality (surgery, radiation, and/or chemotherapy), we found no significant relationship between recurrence-free survival and hybrid DNA (P = 0.105), hybrid RNA (P = 0.540), hybrid ecDNA (P = 0.486), or human circular ecDNA (P = 0.282). Nor did we find a significant relationship between overall survival and hybrid DNA (P = 0.150), hybrid RNA (P = 0.525), hybrid ecDNA (P = 0.667), or human circular ecDNA (P = 0.195).
In TCGA, there were 7 deaths in the 28 patients for whom survival data was available. Multivariable survival analysis was not possible due to the size of the study population, however univariable analysis showed that neither patients with hybrid DNA (P = 0.172) nor human circular ecDNA (P = 0.649) were more likely to die. Patients with hybrid ecDNA had decreased likelihood of death (0/7 deaths vs. 7/21 deaths in nonhybrid ecDNA; P = 0.078), as did patients with hybrid RNA (1/15 deaths vs. 6/13 in nonhybrid RNA; P = 0.016).
It is critical to note that neither our cohort nor TCGA was powered to detect an expected difference in survival. We were unable to perform multivariable analyses to adjust for additional clinical factors due to the small cohort and infrequency of failure events.
Discussion
Intrachromosomal oncogene transcription and amplification have been the dominant paradigm for oncogene-mediated transformation for decades. Similarly, HPV viral integration into human chromosomes and expression of oncoproteins E6 and E7 have been the classic mechanism for HPV-mediated oncogenesis (4). The recent definition of ecDNA as a driver of oncogene amplification and overexpression that is found in nearly half of all human cancers has altered the paradigm of the role of extrachromosomal circular DNA in carcinogenesis (6, 11). Recently nonintegrated HPV has been reported in HPVOPC, indicating that HPV may exert oncogenic effects while maintaining status as vcDNA (46). The data we present describe an unexpected interaction of the HPV genome with cancer-associated ecDNA. Specifically, our results suggest that HPVOPC frequently employ ecDNA in human as well as human–viral hybrid forms that leverage ecDNA mediated viral gene transcription to drive high expression of hybrid viral–human oncogene transcripts. Furthermore, HPVOPC ecDNA structures express a broad variety of human and hybrid cancer-related transcripts including functional oncogenic noncoding RNA and immune evasion checkpoint proteins, in addition to traditional oncogenes. HPVOPC hybrid transcripts drive high expression of E6*I as well as other noncanonical HPV transcripts. This confirms and expands the spectrum of oncogenic HPV gene expression beyond the traditional expression of E6 and E7. In addition, HPV vcDNA is found in HPVOPC in deletion structures that lack coding for viral capsid proteins without evidence of integration into host DNA, as well as in the traditional full-length episomal state. Most strikingly, nearly all HPVOPC contain transcriptionally active, circular, and oncogenic DNA outside of host chromosomes, including vcDNA, human ecDNA, or viral–human hybrid ecDNA.
As noted, HPVOPC tumors express hybrid human viral sequences with increased human and viral gene expression driven by HPV promoters from integrated and ecDNA structures. These hybrid transcripts often include human genes implicated in carcinogenesis, and ecDNA hybrid transcripts are often expressed at a higher level than hybrid transcripts expressed from genomic HPV integration. Our data show that E6 and E7 expression is slightly enhanced in tumors in the context of HPV integration into chromosomal loci, and E1, E2, E4, and E5 expression is only slightly decreased. However, hybrid genomes incorporated into circular ecDNA show with significant upregulation of E6 and E7 as well as dramatic upregulation of human genes within 10kb up and down-stream of the recombined region, with loss of E2, E4, and E5 expression. Taken together, these data indicate that the mechanism for the classic model of HPV carcinogenesis characterized by E6 and E7 overexpression is often driven by the formation of viral–human hybrid ecDNA structures, facilitating amplification of E6 and E7 and associated human genes.
The near universal presence of oncogenic, circular DNA that expresses oncogenic transcripts, whether viral, human, or hybrid, indicates that HPVOPC creates a genomic context that is generally permissive for the stable maintenance of nonchromosomal, transcriptionally active, circular DNA. Traditionally, the E2 protein has been shown to be the regulator of episomal maintenance for high-risk HPV (47, 48). However, we describe hybrid ecDNA structures that exclude E2 in tumors that have dominant or exclusive expression of E6 and E7. It is possible that high-risk HPV gene products other than E2 may have effects on genomic maintenance and chromosomal structure, allowing for perturbations in DNA homeostasis that facilitate stability, replication, and heritability of circular, nonchromosomal DNA structures. The evolution and progression of benign HPV infection to HPVOPC, therefore, requires ongoing maintenance of stable ecDNA or vcDNA as part of carcinogenesis. In this context, these data show three main viral genomic/transcriptomic HPVOPC pathways: (i) nonintegrated HPV vcDNA in a canonical or noncanonical form with broad expression of HPV gene products; (ii) chromosomal integration of HPV with retained expression of HPV gene products and overexpression of human hybrid transcripts; and (iii) formation of viral–human hybrid ecDNA as well as integrated viral DNA with dramatic amplification and overexpression of E6 and E7 as part of human viral fusion transcripts, with dramatic reduction in early E2, E4, and E5 HPV gene product expression. These pathways are independently supported by recent data showing HPV integration does not necessarily result in high levels of E6 and E7 transcripts, and that tumors with nonintegrated HPV16 actually overexpress E6 compared with tumors with integrated or mixed genomes (49, 50).
We identified the most common viral HPVOPC mRNAs as polycistronic transcripts typically >1,000 nt in length, and that the most common full-length transcript in viral ecDNA HPVOPC is 1,476 nt long, beginning at the p97 promoter with splicing at SD226-SA409 and SD880-SA3358 extending to the early polyA tail, with coding potential for E6*I, full-length E7, E1⁁E4, and full-length E5. The function of the E6*I protein (a truncated protein with 43 residues) remains incomplete, although E6*I is found predominantly in high-risk HPV strains and not in low-risk strains (51). Studies suggest it may have opposing functions to full-length E6 with respect to p53 and, through procaspase 8, the cellular response to the TNF-family of cytokines (52, 53). Meanwhile, the E6*I RNA (which contains the E7 ORF) functions as a source of E7 mRNA and increases the efficiency of E7 protein translation (36, 54). E6*I has been observed in E6-expressing keratinocytes to upregulate IL6, a key cytokine in tumorigenesis and inflammation which functions via the JAK–STAT pathway (55). In addition, E6*I expressing cells demonstrated increased p53 levels as well as increased reactive oxygen species levels, which has the effect of increasing DNA damage in cells (56). E6*I RNA levels were significantly higher in cervical samples from patients who had higher grades of dysplasia (57). In our cell line systems, we were able to demonstrate that the most common transcript expressing E6*I, full-length E7, E1⁁E4, and full-length E5 was able to induce proliferation, indicating that the net biologic effect of coordinated expression of these genes can support a malignant phenotype. Conventional and long-read RNA-seq of HPVOPC also suggests that mechanisms of oncogenesis independent of E6 exist, as we found the truncated form E6*I to be more common than full-length E6 in mRNAs. This finding is of interest because the crystal structure of the ternary complex that inactivates p53 is comprised of full-length E6, not E6*I (58). We found E5 protein was commonly expressed, and E5 is increasingly being recognized as a driver of HPV-mediated tumorigenesis (59) and mediates resistance to checkpoint inhibitor blockade in head and neck SCC and can be targeted (60). Finally, we have recently published data demonstrating a similar proliferative cellular phenotype in systems where E2/E4/E5 are expressed in comparison with E6/E7 expression, indicating that E6/E7 gene products may not be as critical for HPV-mediated carcinogenesis (61).
We have noted a striking propensity for canonical splicing junctions to be preserved in HPV viral transcripts, and that these canonical junctions serve as donors for acceptors to diverse human hybrid transcripts, resulting in consistent expression of specific HPV transcripts. Previous studies have identified hybrid viral–human transcripts in HPV anogenital lesions with the HPV splice donor being in the E1 ORF, but these findings have not yet been described in HPVOPC (62). We found that a unifying feature of tumors with hybrid genomes was preferential use of a human splice acceptor site for SD880, rather than the SA3358 HPV site. This is of interest because SA3358 has been documented as the primary splice acceptor involved in transcripts coding for E4, E5, E6, and E7 (38). Together, these data show that alternatively spliced forms of HPV genes, including E6*I, are in fact more highly expressed than conventional transcripts. These transcripts may be key to HPV carcinogenesis with potential value as therapeutic targets. In addition, the proliferation data we provide as well as prior reports of E2, E4, and E5 cooperative ability to support carcinogenesis indicate that coordinated expression of multiple transcripts may be necessary to reproduce phenotypic characteristics of malignancy (61). Similarly, we have found that the overexpressed human transcripts that are associated with human and hybrid ecDNA are quite diverse, comprising traditional oncogenes, for example, CCND1, EGFR, oncogenic long noncoding RNAs like PVT1, as well as immune modulatory genes, including PD-L1. In addition, circular HPV DNA molecules have often been assumed to represent traditional intact, full-length forms found in intact virus. However, we found that noncanonical HPV genome, including deletion in the L1/L2 region, was noted in the majority of tumors that contained HPV vcDNA, and that full-length intact HPV genome was not found in these tumors with noncanonical HPV vcDNA.
These data do have limitations and present opportunities for further investigation. The effect of viral promoters and genes that facilitate high expression of hybrid transcripts and human oncogene expression in hybrid ecDNA structures is presumably facilitated by chromatin modification that is found in human ecDNA (7). However, the additional effect of viral promoters in an ecDNA context may facilitate chromatin structural effects. Although we have previously noted that human ecDNA is associated with poorer prognosis for human cancers in general, we have not seen worse prognosis associated with human ecDNA in HPVOPC, perhaps due to the small sample sizes available for ecDNA characterization, as well as the general favorable prognosis of HPVOPC.
The understanding of the role of ecDNA in HPVOPC has direct implications for development of HPVOPC therapy as ecDNA-mediated amplification has been recognized as a potential mechanism of therapeutic resistance and specific, overexpressed oncogenes on ecDNA may be targeted with precision medicine approaches (6). ecDNA may also have interactions with chromosomes, affecting stability and integration that may also provide therapeutic opportunities for HPVOPC (63), including the presence of ecDNA as a marker for susceptibility to DNA repair inhibition. Indeed, HPVOPC has been shown to be sensitive to Wee-1 and CHK1/2 inhibition, and SCC090 cells that we have identified as ecDNA driven in this manuscript have shown dramatic apoptotic response to ionizing radiation in combination with the Chk1/2 inhibitor prexasertib (64–66). Nonchromosomal circular oncogenic DNA is present and transcriptionally active in nearly all HPVOPC in the form of ecDNA and vcDNA, and the mechanisms of circular DNA formation and maintenance themselves may be potential therapeutic targets. It is possible that the HPV gene products that facilitate general mechanisms of nonchromosomal circular DNA formation may be targeted in all HPVOPC, to target vcDNA as well as human and hybrid ecDNA that transcribe a variety of oncogenic products. Finally, human genes overexpressed via ecDNA in individual tumors have key driver oncogenic roles as well as roles in immune evasion that may serve as targets for personalized therapy.
Authors' Disclosures
N. Nguyen reports other support from Boundless Bio, Inc., outside the submitted work. A. Finegersh reports grants from American Head and Neck Society during the conduct of the study. P.S. Mischel reports being a co-founder of Boundless Bio, Inc., a company focused on developing new treatments for patients with ecDNA-driven cancers; P.S. Mischel also reports equity interest in and is a member of the scientific advisory board in Boundless Bio, Inc. V. Bafna reports grants, personal fees, non-financial support, and other support from Boundless Bio, Inc., and non-financial support and other support from Abterra, Inc. outside the submitted work; V. Bafna is also a co-founder, consultant, SAB member, and has equity interest in Boundless Bio, inc. and Abterra, Inc. No disclosures were reported by the other authors.
Authors' Contributions
J. Pang: Conceptualization, data curation, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. N. Nguyen: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. J. Luebeck: Conceptualization, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. L. Ball: Formal analysis, writing–original draft, writing–review and editing. A. Finegersh: Formal analysis, writing–original draft. S. Ren: Data curation, formal analysis, investigation, visualization. T. Nakagawa: Data curation, investigation. M. Flagg: Formal analysis, writing–original draft. S. Sadat: Data curation. P.S. Mischel: Writing–original draft, writing–review and editing. G. Xu: Data curation, visualization. K. Fisch: Data curation, writing–original draft, writing–review and editing. T. Guo: Writing–original draft, writing–review and editing. G. Cahill: Writing–original draft, writing–review and editing. B. Panuganti: Data curation. V. Bafna: Conceptualization, resources, supervision, funding acquisition, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. J. Califano: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.
Acknowledgments
The following funding sources are noted: NIH/NIDCD T32 DC000028 (to J. Pang), NIH/NIDCR R01DE023347-04 and NIH 1R01 CA204264-01 (to J. Califano), NIH/CTSA UL1TR001442 (to K. Fisch), and R01GM114362 (to V. Bafna).
In addition, we acknowledge Amanda Birmingham for assistance with data processing within the UCSD Center for Computational Biology.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.