Functional variants of human papillomavirus type 16 demonstrate host genome integration and transcriptional alterations corresponding to their unique cancer epidemiology
© The Author(s). 2016
Received: 19 May 2016
Accepted: 25 October 2016
Published: 2 November 2016
The Erratum to this article has been published in BMC Genomics 2016 17:1018
Human papillomaviruses (HPVs) are a worldwide burden as they are a widespread group of tumour viruses in humans. Having a tropism for mucosal tissues, high-risk HPVs are detected in nearly all cervical cancers. HPV16 is the most common high-risk type but not all women infected with high-risk HPV develop a malignant tumour. Likely relevant, HPV genomes are polymorphic and some HPV16 single nucleotide polymorphisms (SNPs) are under evolutionary constraint instigating variable oncogenicity and immunogenicity in the infected host.
To investigate the tumourigenicity of two common HPV16 variants, we used our recently developed, three-dimensional organotypic model reminiscent of the natural HPV infectious cycle and conducted various “omics” and bioinformatics approaches. Based on epidemiological studies we chose to examine the HPV16 Asian-American (AA) and HPV16 European Prototype (EP) variants. They differ by three non-synonymous SNPs in the transforming and virus-encoded E6 oncogene where AAE6 is classified as a high- and EPE6 as a low-risk variant. Remarkably, the high-risk AAE6 variant genome integrated into the host DNA, while the low-risk EPE6 variant genome remained episomal as evidenced by highly sensitive Capt-HPV sequencing. RNA-seq experiments showed that the truncated form of AAE6, integrated in chromosome 5q32, produced a local gene over-expression and a large variety of viral-human fusion transcripts, including long distance spliced transcripts. In addition, differential enrichment of host cell pathways was observed between both HPV16 E6 variant-containing epithelia. Finally, in the high-risk variant, we detected a molecular signature of host chromosomal instability, a common property of cancer cells.
We show how naturally occurring SNPs in the HPV16 E6 oncogene cause significant changes in the outcome of HPV infections and subsequent viral and host transcriptome alterations prone to drive carcinogenesis. Host genome instability is closely linked to viral integration into the host genome of HPV-infected cells, which is a key phenomenon for malignant cellular transformation and the reason for uncontrolled E6 oncogene expression. In particular, the finding of variant-specific integration potential represents a new paradigm in HPV variant biology.
KeywordsHuman papillomavirus HPV16 E6 oncogene variants Organotypic rafts Viral integration Transcriptomics Pathogen-host relationship
Approximately 20 % of human cancers are caused by infectious agents , including >500,000 patients diagnosed annually with human papillomavirus (HPV) associated cancers. Oncogenic HPV, denoted as “high-risk”, is the primary risk factor for cervical cancer due to its exclusive tropism for mucosal tissues [2, 3]. Upon persistent infections of the cervical mucosa, oncogenic HPVs can cause progression from low- to high-grade cervical intraepithelial neoplasias that, without ablative treatment, may develop into invasive carcinomas. At the molecular level HPV is a double-stranded DNA virus and, to date, the sequences of over 200 types have been described . The ~8 kbp genome of HPV contains 8 functional open reading frames (ORFs) that encode 5 early gene products (E1, E2, E5, E6 and E7) and 3 late gene products (E4, L1 and L2). While E1 and E2 are involved in DNA replication and transcriptional regulation of the viral genome , HPV’s potent tumourigenicity is primarily due to E6 , E7 , and E5 . L1 and L2 are structural proteins that self-assemble to form icosahedral capsids , while the fused product of ORFs E1 and E4 (E1^E4) is most abundant in the productive viral life cycle, coinciding with the onset of viral DNA amplification .
Among the HPV types, HPV16 (a member of species Alphapapillomavirus 9) is the most prevalent in cervical cancers. Intriguingly, and perhaps related to its prevalence, the HPV16 genome is polymorphic. Evolutionary analyses have revealed that the worldwide diversity of HPV16 genomes evolved for over 200,000 years , leading to five phylogenetic branches representing isolates from Africa, Europe, Asia and the Americas . Furthermore, each branch can be further dissected into intratypic single nucleotide polymorphisms (SNPs) or variants differing in their host persistence and frequency of detection in human pre-cancers and cancers (reviewed in ). The tumourigenic differences of these SNPs have been ascribed largely to those within the E6 oncogene [14–17]. The Asian-American (AAE6) and European Prototype (EPE6) are common HPV16 genome variants which differ by six SNPs in their E6 genes, three of which are non-synonymous, leading to the 151-residue AAE6 protein differing by three amino-acids: Q14H, H78Y, and L83V  (with residue 14 and 83 being under Darwinian constraint ).
Epidemiological studies showed that the AAE6 genome variant is a higher risk factor for dysplasia as well as an earlier onset of invasive tumours than EPE6 [20–26]. As well, AAE6 has a greater transforming, migratory, and invasive potential than EPE6 when retrovirally transduced into primary human keratinocytes during recent long-term in vitro immortalization studies [27–30]. These results suggested that coding changes in E6 have strong mechanistic and functional consequences for infection and thus contribute to marked differences in cancer risk of HPV16 variants.
Here, to further advance our mechanistic understanding of the impact of these common but epidemiologically and clinically important E6 SNPs, we conducted an “-omics” analysis on the NIKS-based organotypic epithelia containing the HPV16 variants AAE6 and EPE6 (Fig. 1). Modern deep sequencing techniques have been used to study HPV [35–39], but only recently in the context of intratypic variants , and not using an organotypic epithelial model with full viral variant genomes. Instead, our complete approach allowed a comparison of these variants with regards to their integration capacity and subsequent transcriptional consequences in close to in vivo conditions, resulting in viral integration and a molecular signature of host chromosomal instability for AAE6 only.
Results and discussion
Viral integration in the HPV16 AAE6 but not EPE6 epithelium
To permit the viral life cycle in a raft culture system, we transfected the keratinocytes, prior to rafting, with complete viral genomes containing either the HPV16 EPE6 or AAE6 variant. A similar technique was used in a recent study to successfully study varicella zoster virus , providing a keratinocyte model and a “global” perspective of all changes in host transcription in response to a pathogen. As illustrated in Fig. 1, over a 14 day differentiation process, we observed that the NIKS were normal epithelia whereas NIKS with HPV16 EPE6 exhibited a mild dysplasia and NIKS with HPV16 AAE6 exhibited a moderate dysplasia.
Based on the Dr.VIS (Viral Integration Site) v2.0 database of HPV16 integration sites , this exact region (5q32) of integration is not frequent, but potentially recurrent as it was found in 2 out of 878 previously documented sites. The nearest fragile site was 13 Mb upstream of this integration site: FRA5C, 5q31.1. Since repeated regions might be prone to genome rearrangements and therefore prone to HPV integration, we scanned the adjacent regions using the UCSC hg19 genome browser RepeatMasker track for human repeat elements and found a nearby 158 bp long interspersed nuclear element (LINE): L1MB5 located from Chr5 nt position 149,347,143 to 149,347,300. Indeed, L1MB5-derived sequences have been documented as breakpoints, such as in the human genes HPRT , CYP2C , and in proximity of genes containing the ubiquitin ligase Mib-herc2 domain, which mediates Notch signalling . Strikingly, this domain contains the Hect region, homologous to the E6-associated protein carboxyl terminus, raising the question of whether or not the underlying homology could play a role in this target site selection. Another, non-exclusive hypothesis is that the frequent hypo-methylation of LINE elements plays a role to facilitate access to the chromosomal DNA and associated genomic instability [47, 48]. Altogether, our three-dimensional organotypic cultures demonstrated that the HPV16 AAE6 variant had integrated into the host genome while the EPE6 variant remained episomal, suggesting an increased propensity towards integration due to AAE6. A previous study of HPV16 integration propensity with respect to the variants did not demonstrate a statistically significant difference (P-value = 0.28, two-tailed Fisher’s exact test) between EPE6 (3 episomal and 20 integrated cases) and the E-T350G variant (6 episomal and 16 integrated cases, responsible for one of the residue changes also found in AAE6: L83V) . Only one tumour sample in their set contained the AA variant, therefore precluding a formal analysis of its propensity to integrate, but notably it was in integrated form.
The HPV16 AAE6 epithelium has a unique transcriptional profile
To quantitatively account for sample variance, we also performed differential expression analysis of the viral gene counts using DESeq . DESeq software tests for differential expression in library size-corrected count data using a negative binomial distribution model. In agreement with our previous RT-qPCR results , we found significantly more E6 (24.05 fold higher, P < 10−10) and E7 (17.30 fold higher, P < 10−10) counts in triplicate AAE6 rafts in comparison to triplicate EPE6 rafts (Fig. 3c). Taken together, analyses of viral transcriptome data revealed that the AAE6 viral transcriptome significantly differs from that of EPE6 in a manner that is indicative of integration, with increased E6 and E7 levels [53–55]. Evidently, AAE6 transcriptome profiles are lacking E2 and have increased E6/E7 oncogene expression, perhaps due to loss of transcriptional repression by E2. We therefore reasoned that the increased levels of E6/E7 expression between the variants were ultimately due to their viral integration status, as we hypothesized in our phenotypic study, and confirmed by Capt-HPV, leading to a significant effect on the host transcriptome .
Nature of viral-human fusion transcripts detected in HPV16 AAE6 epithelium
Integration loci detected by ViralFusionSeq
Mapped human transcript†
HPV transcript breakpoint(s)‡
Solute carrier family 26, member 2
Cyclic GMP- Phosphodiesterase 6A alpha subunit
In accordance with the structure of the HPV integration, the transcript breakpoints mapped to either the E1 or L1 HPV16 ORF. Alternative splicing was detected with the viral nucleotide position at the fusion site of one class of the viral-human fusion transcripts (Fig. 3a): nt 880 (splice donor, SD) in the E1 gene . This is the same SD site for the E1^E4 splice transcript typically expressed in the late stage of the viral life cycle , and previously shown to be expressed in our EPE6 epithelia . HPV16 viral-human fusion transcripts are often detected with a breakpoint at this natural splice donor site [56, 67, 68], and the coverage plot for AAE6 shows decreased coverage for transcripts downstream of this E1 SD site, supporting the hypothesis of alternative splicing. With respect to the L1 breakpoints, the typical L1 splice acceptor (SA) site is at nt 5639 , but notably in our study, the viral-human fusion transcripts here had a putative downstream SA site at nt 5778. Interestingly, the coverage plot of the viral transcriptome shows nt 5778 as the site where L1 coverage begins to be detected in AAE6 rafts (Fig. 3a), so we reasoned that this discrepancy in SA site could be due to either a cryptic SA site in the HPV16 W12E genome (although not found previously in the literature) or simply due to integration truncating the upstream region of L1.
Next, we mapped the human portion of the fusion transcripts using VFS’s clipped-seq (CS) and read-pair (RP) methods. Confirmed by both these methods, two fusions mapped to the human chromosome location 5q32, occurring within the solute carrier family 26 (anion exchanger), member 2 (SLC26A2) and phosphodiesterase 6A, cGMP-specific, rod, alpha (PDE6A) human ORFs (Table 1). Strikingly, along with detection of fusion transcripts with these genes, we detected a significant increase in the expression of human genes from this region in AAE6 epithelia compared to normal epithelia, namely SLC26A2 (114.19 fold increase, P = 2.14 × 10−173) and colony-stimulating factor 1 receptor (CSF1R, 407.82 fold increase, P = 4.70 × 10−112, which was only detected as RP fusion reads by VFS, and not confirmed by CS). This observation is in agreement with others who have found that, in numerous cervical carcinomas across multiple high-risk HPV types, HPV integration leads to an increase in the expression of genes adjacent to integration loci . To explain the molecular basis of this cis-effect, it has been proposed to be the result of viral promotor-driven expression or somatic genome amplification at the integration site [70, 71]. In the present case, this last hypothesis is unlikely because the AAE6 integration produced a clean 11 bp deletion of the target region that led to two co-linear viral-human junctions (2J-COL), which is not associated with gene amplification .
Functional human fusion proteins can be formed due to chromosomal translocations in cancer cells . The elucidation of novel protein-coding viral-human fusion transcripts is particularly intriguing due to their potentially functional roles within host cells. Using immunofluorescence for the expressed portion of the SLC26A2 protein in formalin-fixed and paraffin embedded (FFPE) rafts, we determined that SLC26A2 protein expression was aberrantly high in AAE6 compared to EPE6, supposedly as a result of its viral-human fusion and increased transcription (Fig. 2b). This translated fusion protein contains exon 3 of the transmembrane protein SLC26A2, previously known as diastrophic dysplasia sulfate transporter (DTDST) , which encodes the carboxy-terminal cytoplasmic sulfate transporter and anti-sigma factor (STAS) domain . We cannot find any evidence in the literature of this unique viral-human fusion protein in other HPV-integrated samples. Overall, these chimeric molecules are unique for each sample and to the specific integration site, with presently unknown effect on host cell functions, an aspect to be further researched due to its importance for understanding mechanisms of tumourigenesis as well as in the emerging field of personalized medicine.
The HPV16 AAE6 epithelium reveals a signature of chromosomal instability conducive to host genome integration
Integration of HPV DNA into the host genome is considered to be a key factor for cervical cancer development [67, 75, 76], but the cellular events that initiate the integration process (and selection of insertion sites) remain to be better understood. A reasonable hypothesis is that the integration is triggered by a rare and stochastic target site event, such as a replicative fork stalling or an accidental chromosome double-strand break, leading to an ultimate use of the viral DNA for repair via recombination, template switching (FoSTeS) and/or microhomology-mediated break-induced replication (MMBIR) ([42, 71, 77], and references within each). Indeed, infections with pathogens can cause chromosomal instability by inactivating the host DNA damage response . For HPV, this has been linked to the expression of both HPV16 E6 and E7 oncoproteins, affecting the infected cell’s genome integrity [79–82]. A model of early carcinogenesis due to HPV16 E6 and E7 suggests that this chromosomal instability is caused by uncontrolled proliferation, leading to an insufficient nucleotide pool that cannot support normal replication . Alternatively, E6 alone, through the inactivation of p53, can promote chromosomal instability, at least during early onset of carcinogenesis . Presently, HPV16 AAE6 demonstrated enhanced integration propensity over EPE6 and exhibited increased E6 and E7 oncogene expression, which is in accordance with elevated E6 and E7 levels reported in other studies [53–55]. This enhanced integration ability is based on AAE6’s greater proliferation ability, leading to chromosomal instability. The underlying mechanism of its increased cell growth is the result of a deregulated sugar metabolism (Warburg effect), as we reported previously  and currently under study (Cuninghame et al., in preparation: unpublished observations).
HPV16 AAE6 epithelium exhibits a proliferating phenotype as a consequence of viral integration into the host genome
We have systematically characterized the viral integration process of a common high-risk HPV16 variant and its consequences for the affected host cell. This and earlier work lend themselves to propose a model of increased tumourigenicity in human keratinocyte epithelia where AAE6’s enhanced ability to proliferate leads to chromosomal instability. In such an environment, the host genome may be susceptible to viral integration subsequently increasing E6/E7 oncogene expression and ultimately driving additional tumourigenic changes. Previously, we performed phenotypic studies of the EPE6 and AAE6 variants in a 3D raft model of early carcinogenesis  and determined the functional differences of these variants in longitudinal monolayer cell cultures [27–30]. While necessary for studying the viral life cycle, limitations of the current organotypic model are the lack of immune components, vasculature, and the complexity of tissue heterogeneity that arises. Our current study builds on the foundation of these investigations. We have applied a wide range of molecular analyses, creating a framework which can benefit future virus-host interaction studies with various organotypic cell culture models. A variant-specific integration is worth reporting and should be further investigated, with additional samples from independent donors, as it represents a new paradigm in HPV variant biology. Here we report a viable integration mechanism in a robust viral life cycle model for AAE6. The findings of the current and other studies reported by us [27–31], and others [89–91], are consistent with cancer epidemiology studies demonstrating that the HPV16 AA variant is a higher risk factor for high-grade intraepithelial neoplasia and progression to invasive cervical cancer [22–24, 89]. In the future, HPV variant genotyping could be used as a clinical prognostic factor for patient-centered health services, while the role of individual host genomics on integration, including characterization of integration sites, will be important to consider for personalized medicine approaches.
As described by us previously , we used the Normal/Near-Diploid Immortalized Keratinocytes (NIKS) cell line  to establish 3D organotypic epithelia cultures. These spontaneously immortalized cells were originally derived from neonatal human foreskin and are non-tumourigenic, though contain an additional long arm piece of chromosome 8 (8q). In monolayer they are grown on mitomycin-C-treated Swiss mouse J2/3T3 fibroblast feeder layers , while primary human foreskin fibroblasts (ATCC CRL-2097) are incorporated into the dermal equivalent of organotypic NIKS cultures .
Detection of integrated papillomavirus sequences by next-generation DNA-Seq: Capt-HPV
DNA-Seq was used to confirm the presence and location of the viral integration sites in the human genome using DNA extracted from formalin-fixed paraffin embedded (FFPE) samples which had been prepared previously . DNA was extracted using the DNeasy Blood and Tissue Kit (QIAGEN, Cat# 69504) with the recommended pre-treatment for FFPE samples and the optional RNase treatment. To overcome the limitations of traditional techniques, such as DIPS-PCR (Detection of Integrated Papillomavirus Sequences by ligation mediated PCR), we used an unbiased and state-of-the-art next-generation DNA sequencing technique for detecting HPV viral integration sequences in our samples . Library preparation, sequence capture, and high-throughput sequencing was carried out at the Institut Curie on an Illumina MiSeq platform with a V2 Nano chip (~1 × 106 total reads) with 2 × 151 base pair read length. Analysis of sequencing data was performed using the Galaxy platform [93–95], with the primary goal of detecting the viral-human junction site locations. Packages used were FASTQ Groomer , Bowtie2 , Picard MarkDuplicates , SAMtools BAM-to-SAM and Filter SAM .
RNA-Seq library preparation and sequencing
Isolation of high-quality total RNA from the epithelium of organotypic keratinocyte cultures containing full-length HPV16 E6 variant genomes, European Prototype (EPE6) and Asian-American (AAE6), was described previously . Our keratinocyte model was grown for 14 days to allow simultaneous epithelial differentiation and occurrence of an active viral life cycle. Total RNA for EPE6, AAE6, and HPV16 negative cultures (NIKS), three organotypic raft cultures (n = 3) each, were sent for library preparation and sequencing at The Centre for Applied Genomics, Hospital for Sick Children, Toronto, Canada. RNA-Seq libraries were prepared by Illumina TruSeq® RNA Sample Preparation kit followed by sequencing using an Illumina HiSeq® 2500 platform with Illumina v3 chemistry. One lane of multiplexed, paired-end, 2 × 101 base pair sequencing was performed with nine samples: yielding an average of 40.4 million total reads (~20 to 25 million fragments) per sample (Additional file 1: Table S2).
Viral variant read alignment, mapping, and coverage plotting
The human papillomavirus type 16 W12E isolate genome [GenBank: AF125673] [54, 99] was used as a viral reference sequence since it was the parental sequence modified by site-directed mutagenesis to generate the EPE6 and AAE6 viral genomes used in this study . Only the three non-synonymous nucleotide changes differentiated EPE6 and AAE6 genomes: EPE6 was made by mutating the parental W12E genome at G350T while AAE6 was mutated at G145T and C335T. Prior to alignment and mapping, Bowtie2  was used to build a reference index for HPV16 using the AF125673 W12E isolate RefSeq. TopHat2  was used for alignment to our viral RefSeq. Variant-specific non-synonymous SNPs were confirmed by variant calling with SAMtools . The Broad Institute’s Integrative Genomics Viewer (IGV)  was used to visualize alignment coverage for each sample. Gene-level counts of the HPV16 W12E ORF’s were generated using SAMtools , and normalized with library-size correction factors using the Bioconductor project DESeq  in the statistical environment R . DESeq was also used for differential viral gene expression analysis. DESeq uses a default false discovery rate (FDR) of 10 % for its binomial statistical inference tests to determine differentially expressed genes. Clustered heatmaps of normalized viral gene counts were generated using the gplots package .
Identification of viral-human fusion transcripts
ViralFusionSeq (VFS)  was used, with default parameters, to identify any viral-human fusion transcripts in each of our sample RNA-Seq datasets. As with viral alignment by TopHat2 (described above), the W12E genome was used as a reference sequence for VFS. Briefly, VFS is a Perl script that searches in high-throughput sequencing data (RNA or DNA-Seq) for viral-human fusion transcripts, which are present as a result of viral integration events into host DNA. This software uses read pair (RP) and clipped sequences (CS) to accurately discover and identify viral-fusion sequences . Additionally, VFS is able to reconstruct fusion transcripts by a targeted de novo assembly process. These methods allow us to identify, with single-base resolution, viral-human fusion transcripts present within our epithelial cultures. Viral-human fusion transcripts were compared to known HPV16 integration sites and fusion transcripts with assistance from the database of disease related viral integration sites (Dr. VIS v2.0, ).
We sought to perform protein-level confirmation of highly expressed viral-human fusion transcripts containing exons from human targets SLC26A2 and CSF1R. SLC26A2 protein expression was detected in raft cultures by immunofluorescence, as described previously . Based on the viral-human fusion RNA-Seq data, the primary antibody (rabbit polyclonal, 1:500 dilution, Bethyl Laboratories Inc., Cat. No. A304-467A) was chosen to have specificity for translated exon 3 (epitope between amino acid residue 689 and 739). Although also highly up-regulated, no suitable commercial antibody was found for CSF1R exons 20 to 22.
Human read alignment, mapping, and count generation
Read alignment, mapping, and count generation for the human reference genome (hg19, UCSC nomenclature for GRCh37) was performed by The Centre for Applied Genomics, Hospital for Sick Children, Toronto, Canada. TopHat2  was used for RefSeq while gene- and exon-level counts were generated using HTSeq . Number of reads and percentage of human RefSeq reads defined as aligned, exon, and exon-exon are reported in Additional file 1: Table S2 for each sample analyzed.
Differential expression analysis of human transcriptome
Differential analysis of pair-wise human gene-level counts between NIKS and EPE6, NIKS and AAE6, and EPE6 and AAE6 were performed using the Bioconductor project DESeq  package implemented in the statistical environment R . Raw gene counts from HTSeq were first normalized by estimating the sample library sizes (Additional file 1: Table S3) and applying the size-factor correction to all counts within a given sample. A dispersion plot was made to visualize the variance estimation step prior to differential expression inference (Additional file 2: Figure S4). A clustered heatmap with hierarchical dendrograms was used to show overall sample and biological replicate clustering: the gene expression profile of AAE6 samples was distinct from EPE6 and NIKS (control) samples (Additional file 2: Figure S5). Although EPE6 replicate 3 and NIKS replicate 1 cluster outside of their specific sample group, viral RNA-Seq analysis has confirmed these sample ID’s are correct, and that their grouping is likely a result of the minor host transcriptomic difference between NIKS and EPE6 cultures. DESeq uses a default false discovery rate (FDR) of 10 % for its binomial statistical inference tests to determine differentially expressed genes. However, for downstream analyses of down- and up-regulated genes we used a more stringent adjusted P-value cut-off of 10−5.
CIN70 scoring and micronuclei detection
Host chromosomal instability was assessed, using normalized human gene count data from our RNA-Seq experiments, by calculating a CIN70 gene expression signature score  for EPE6 and AAE6 relative to NIKS epithelia. For each of the 70 genes, a normalized human gene count ratio was calculated for all EPE6 and AAE6 samples relative to the average of the NIKS samples. Relative ratio values were then averaged for all 70 genes in each sample and a Welch’s T-test, for unequal variance, was used to determine whether there was a statistically significant difference in host chromosomal instability signature between EPE6 and AAE6 epithelia. We used a significance level of P < 0.05. As a morphological assessment of chromosomal instability we screened haematoxylin and eosin-stained sections from formalin-fixed and paraffin-embedded NIKS, EPE6, and AAE6 epithelia for micronuclei (MN). These aberrant nuclei structures  were detected using light microscopy with high-magnification (at least 400×).
Gene set enrichment analysis and networks
Enrichment of host biological processes of differentially expressed human genes was determined using the Gene Ontology (GO) Term Enrichment Service hosted on the AmiGO 2 website . Only biological processes were included. Terms were considered significantly enriched if the Bonferroni-corrected P-value was less than 0.05. To aid in the visual interpretation of down- and up-regulated gene sets, co-expression networks were constructed with Cytoscape software . Pearson correlation coefficients were calculated for each gene-gene pairwise comparison in highly significant down- and up-regulated genes between AAE6 and EPE6 (Additional file 6 for down- and up-regulated gene-gene pairwise comparisons, respectively). Pearson correlation coefficient cut-offs used for networking were selected strategically to produce small distinct clusters of genes, since setting the threshold too low results in all nodes connected, and setting the threshold too high results in a lack of clusters.
Thank you to Dr. Allyson Holmes at the Institut Curie for her valuable feedback and collaboration on the DNA-Seq experiments. Special thanks go to Melissa Togtema for her insightful comments while preparing the manuscript as well as Darryl Willick for his help in setting up and maintaining the Galaxy platform hosted at the Lakehead University High Performance Computing Centre (LUHPCC).
This work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) grants to IZ (#355858-2008, #435891-2013, #RGPIN-2015-03855), NSERC Alexander Graham Bell Canada Graduate Scholarship-Doctoral (CGS-D) to RJ (#454402-2014), NSERC Alexander Graham Bell Canada Graduate Scholarship-Masters (CGS-M) to SC (#442618-2013), and an NSERC Undergraduate Student Research Award (USRA) to JB (#483630-2015). The funding bodies had no role in study design, data collection, data analysis and interpretation, or preparation of the manuscript.
Availability of data and materials
Raw sequence data used in this article can be accessed via the Sequence Read Archive (SRA), study accession number SRP055094 (http://www.ncbi.nlm.nih.gov/sra/SRP055094) and National Center for Biotechnology Information (NCBI) BioProject, accession number PRJNA275642 (http://www.ncbi.nlm.nih.gov/bioproject/PRJNA275642). Remaining supporting data can be accessed as Additional files, while software and tools used have been cited throughout the Methods section.
This interdisciplinary study was initially conceived by IZ and RJ refined the bioinformatics portion in collaboration with BR, WF, SL, and AN. RJ, PL, and IZ designed and carried out the 3D organotypic skin culturing experiments. IZ and PL contributed reagents, materials, and methods for culturing experiments. RJ, BR, SC and JB performed RNA-Seq and follow-up data analyses. AN contributed reagents, materials, and methods for DNA sequencing. SL and RJ performed DNA-Seq and follow-up analyses. RJ, BR, SL, SC, JB, WF, PL, AN, and IZ contributed to data interpretation. All authors contributed to writing the paper with RJ being the lead author and IZ having considerable input into the writing. All authors have read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Bouvard V, Baan R, Straif K, Grosse Y, Secretan B, Ghissassi FE, et al. A review of human carcinogens—Part B: biological agents. Lancet Oncol. 2009;10:321–2.View ArticlePubMedGoogle Scholar
- zur Hausen H. Papillomavirus infections—a major cause of human cancers. BBA-Rev Cancer. 1996;1288:F55–78.Google Scholar
- zur Hausen H. Papillomaviruses and cancer: from basic studies to clinical investigations. Nat Rev Cancer. 2002;2:342–50.View ArticlePubMedGoogle Scholar
- Kocjan BJ, Bzhalava D, Forslund O, Dillner J, Poljak M. Molecular methods for identification and characterization of novel papillomaviruses. Clin Microbiol Infect. 2015;21:808–16.View ArticlePubMedGoogle Scholar
- Doorbar J, Quint W, Banks L, Bravo IG, Stoler M, Broker TR, et al. The biology and life-cycle of human papillomaviruses. Vaccine. 2012;30:F55–70.View ArticlePubMedGoogle Scholar
- Vande Pol SB, Klingelhutz AJ. Papillomavirus E6 oncoproteins. Virology. 2013;445:115–37.View ArticlePubMedPubMed CentralGoogle Scholar
- Roman A, Münger K. The papillomavirus E7 proteins. Virology. 2013;445:138–68.View ArticlePubMedPubMed CentralGoogle Scholar
- Maufort JP, Williams SM, Pitot HC, Lambert PF. Human papillomavirus 16 E5 oncogene contributes to two stages of skin carcinogenesis. Cancer Res. 2007;67:6106–12.View ArticlePubMedPubMed CentralGoogle Scholar
- Conway MJ, Meyers C. Replication and assembly of human papillomaviruses. J Dent Res. 2009;88:307–17.View ArticlePubMedPubMed CentralGoogle Scholar
- Middleton K, Peh W, Southern S, Griffin H, Sotlar K, Nakahara T, et al. Organization of human papillomavirus productive cycle during neoplastic progression provides a basis for selection of diagnostic markers. J Virol. 2003;77:10186–201.View ArticlePubMedPubMed CentralGoogle Scholar
- Bernard HU. The clinical importance of the nomenclature, evolution and taxonomy of human papillomaviruses. J Clin Virol. 2005;32:1–6.View ArticleGoogle Scholar
- Yamada T, Manos MM, Peto J, Greer CE, Munoz N, Bosch FX, et al. Human papillomavirus type 16 sequence variation in cervical cancers: a worldwide perspective. J Virol. 1997;71:2463–72.PubMedPubMed CentralGoogle Scholar
- Burk RD, Harari A, Chen Z. Human papillomavirus genome variants. Virology. 2013;445:232–43.View ArticlePubMedPubMed CentralGoogle Scholar
- Grodzki M, Besson G, Clavel C, Arslan A, Franceschi S, Birembaut P, et al. Increased risk for cervical disease progression of French women infected with the human papillomavirus type 16 E6-350G variant. Cancer Epidemiol Biomarkers Prev. 2006;15:820–2.View ArticlePubMedGoogle Scholar
- Zehbe I, Wilander E, Delius H, Tommasino M. Human papillomavirus 16 E6 variants are more prevalent in invasive cervical carcinoma than the prototype. Cancer Res. 1998;58:829–33.PubMedGoogle Scholar
- Zehbe I, Voglino G, Delius H, Wilander E, Tommasino M. Risk of cervical cancer and geographical variations of human papillomavirus 16 E6 polymorphisms. Lancet. 1998;352:1441–2.View ArticlePubMedGoogle Scholar
- Zehbe I, Voglino G, Wilander E, Delius H, Marongiu A, Edler L, et al. p53 codon 72 polymorphism and various human papillomavirus 16 E6 genotypes are risk factors for cervical cancer development. Cancer Res. 2001;61:608–11.PubMedGoogle Scholar
- Cornet I, Gheit T, Franceschi S, Vignat J, Burk RD, Sylla BS, et al. Human papillomavirus type 16 genetic variants: phylogeny and classification based on E6 and LCR. J Virol. 2012;86:6855–61.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen Z, Terai M, Fu L, Herrero R, DeSalle R, Burk RD. Diversifying selection in human papillomavirus type 16 lineages based on complete genome analyses. J Virol. 2005;79:7014–23.View ArticlePubMedPubMed CentralGoogle Scholar
- Xi LF, Koutsky LA, Galloway DA, Kiviat NB, Kuypers J, Hughes JP, et al. Genomic variation of human papillomavirus type 16 and risk for high grade cervical intraepithelial neoplasia. J Natl Cancer Inst. 1997;89:796–802.View ArticlePubMedGoogle Scholar
- Villa LL, Sichero L, Rahal P, Caballero O, Ferenczy A, Rohan T, et al. Molecular variants of human papillomavirus types 16 and 18 preferentially associated with cervical neoplasia. J Gen Virol. 2000;81:2959–68.View ArticlePubMedGoogle Scholar
- Berumen J, Ordonez RM, Lazcano E, Salmeron J, Galvan SC, Estrada RA, et al. Asian American variant of human papillomavirus 16 and risk for cervical cancer: a case–control study. J Natl Cancer Inst. 2001;93:1325–30.View ArticlePubMedGoogle Scholar
- Xi LF, Koutsky LA, Hildesheim A, Galloway DA, Wheeler CM, Winer RL, et al. Risk for high-grade cervical intraepithelial neoplasia associated with variants of human papillomavirus types 16 and 18. Cancer Epidemiol Biomarkers Prev. 2007;16:4–10.View ArticlePubMedGoogle Scholar
- Zuna RE, Moore WE, Shanesmith RP, Dunn ST, Wang SS, Schiffman M, et al. Association of HPV16 E6 variants with diagnostic severity in cervical cytology samples of 354 women in a US population. Int J Cancer. 2009;125:2609–13.View ArticlePubMedPubMed CentralGoogle Scholar
- Schiffman M, Rodriguez AC, Chen Z, Wacholder S, Herrero R, Hildesheim A, et al. A population-based prospective study of carcinogenic human papillomavirus variant lineages, viral persistence, and cervical neoplasia. Cancer Res. 2010;70:3159–69.View ArticlePubMedPubMed CentralGoogle Scholar
- Freitas LB, Chen Z, Muqui EF, Boldrini NAT, Miranda AE, Spano LC, et al. Human Papillomavirus 16 Non-European Variants Are Preferentially Associated with High-Grade Cervical Lesions. PLoS One. 2014;9:e100746.View ArticlePubMedPubMed CentralGoogle Scholar
- Zehbe I, Richard C, DeCarlo CA, Shai A, Lambert PF, Lichtig H, et al. Human papillomavirus 16 E6 variants differ in their dysregulation of human keratinocyte differentiation and apoptosis. Virology. 2009;383:69–77.View ArticlePubMedGoogle Scholar
- Richard C, Lanner C, Naryzhny S, Sherman L, Lee H, Lambert PF, et al. The immortalizing and transforming ability of two common human papillomavirus 16 E6 variants with different prevalences in cervical cancer. Oncogene. 2010;29:3435–45.View ArticlePubMedGoogle Scholar
- Niccoli S, Abraham S, Richard C, Zehbe I. The Asian-American E6 variant protein of human papillomavirus 16 alone is sufficient to promote immortalization, transformation, and migration of primary human foreskin keratinocytes. J Virol. 2012;86:12384–96.View ArticlePubMedPubMed CentralGoogle Scholar
- Togtema M, Jackson R, Richard C, Niccoli S, Zehbe I. The human papillomavirus 16 European-T350G E6 variant can immortalize but not transform keratinocytes in the absence of E7. Virology. 2015;485:274–82.View ArticlePubMedGoogle Scholar
- Jackson R, Togtema M, Lambert PF, Zehbe I. Tumourigenesis Driven by the Human Papillomavirus Type 16 Asian-American E6 Variant in a Three-Dimensional Keratinocyte Model. PLoS One. 2014;9:e101540.View ArticlePubMedPubMed CentralGoogle Scholar
- Allen-Hoffmann BL, Schlosser SJ, Ivarie CA, Sattler CA, Meisner LF, O’Connor SL. Normal growth and differentiation in a spontaneously immortalized near-diploid human keratinocyte cell line. NIKS J Invest Dermatol. 2000;114:444–55.View ArticlePubMedGoogle Scholar
- Schütze DM, Snijders PJ, Bosch L, Kramer D, Meijer CJ, Steenbergen RD. Differential In Vitro Immortalization Capacity of Eleven, Probable High-Risk Human Papillomavirus Types. J Virol. 2014;88:1714–24.View ArticlePubMedPubMed CentralGoogle Scholar
- Poreba E, Broniarczyk JK, Gozdzicka-Jozefiak A. Epigenetic mechanisms in virus-induced tumorigenesis. Clin Epigenetics. 2011;2:233–47.View ArticlePubMedPubMed CentralGoogle Scholar
- Mine KL, Shulzhenko N, Yambartsev A, Rochman M, Sanson GF, Lando M, et al. Gene network reconstruction reveals cell cycle and antiviral genes as major drivers of cervical cancer. Nat Commun. 2013;4:1806.View ArticlePubMedPubMed CentralGoogle Scholar
- Khoury JD, Tannir NM, Williams MD, Chen Y, Yao H, Zhang J, et al. The Landscape of DNA Virus Associations Across Human Malignant Cancers Using RNA-Seq: An Analysis of 3775 Cases. J Virol. 2013;87:8916–26.View ArticlePubMedPubMed CentralGoogle Scholar
- Bryant D, Onions T, Raybould R, Flynn Á, Tristram A, Meyrick S, et al. mRNA sequencing of novel cell lines from human papillomavirus type‐16 related vulval intraepithelial neoplasia: Consequences of expression of HPV16 E4 and E5. J Med Virol. 2014;86:1534–41.View ArticlePubMedGoogle Scholar
- Chandrani P, Kulkarni V, Iyer P, Upadhyay P, Chaubal R, Das P, et al. NGS-based approach to determine the presence of HPV and their sites of integration in human cancer genome. Br J Cancer. 2015;112:1958–65.View ArticlePubMedPubMed CentralGoogle Scholar
- Cullen M, Boland J, Schiffman M, Zhang X, Wentzensen N, Yang Q, et al. Deep sequencing of HPV16 genomes: A new high-throughput tool for exploring the carcinogenicity and natural history of HPV16 infection. Papillomavirus Research. 2015;1:3–11.View ArticlePubMedPubMed CentralGoogle Scholar
- Lavezzo E, Masi G, Toppo S, Franchin E, Gazzola V, Sinigaglia A, et al. Characterization of Intra-Type Variants of Oncogenic Human Papillomaviruses by Next-Generation Deep Sequencing of the E6/E7 Region. Viruses. 2016;8:79.View ArticlePubMedPubMed CentralGoogle Scholar
- Jones M, Dry IR, Frampton D, Singh M, Kanda RK, Yee MB, et al. RNA-seq analysis of host and viral gene expression highlights interaction between varicella zoster virus and keratinocyte differentiation. PLoS Pathog. 2014;10:e1003896.View ArticlePubMedPubMed CentralGoogle Scholar
- Holmes A, Lameiras S, Jeannot E, Marie Y, Castera L, Sastre-Garau X, et al. Mechanistic signatures of HPV insertions in cervical carcinomas. Genome Med. 2016;1:16004.Google Scholar
- Yang X, Li M, Liu Q, Zhang Y, Qian J, Wan X, et al. Dr.VIS v2.0: an updated database of human disease-related viral integration sites in the era of high-throughput deep sequencing. Nucl Acids Res. 2015;43:D887–92.View ArticlePubMedGoogle Scholar
- Williams M, Rainville IR, Nicklas JA. Use of inverse PCR to amplify and sequence breakpoints of HPRT deletion and translocation mutations. Environ Mol Mutagen. 2002;39:22–32.View ArticlePubMedGoogle Scholar
- Zhou S. Cytochrome P450 2D6: structure, function, regulation and polymorphism. CRC Press; 2016 Feb 24Google Scholar
- del Rosario RC, Rayan NA, Prabhakar S. Noncoding origins of anthropoid traits and a new null model of transposon functionalization. Genome Res. 2014;24:1469–84.View ArticlePubMedPubMed CentralGoogle Scholar
- Richards KL, Zhang B, Baggerly KA, Colella S, Lang JC, Schuller DE, et al. Genome-wide hypomethylation in head and neck cancer is more pronounced in HPV-negative tumors and is associated with genomic instability. PLoS One. 2009;4:e4941.View ArticlePubMedPubMed CentralGoogle Scholar
- Baba Y, Watanabe M, Murata A, Shigaki H, Miyake K, Ishimoto T, et al. LINE-1 hypomethylation, DNA copy number alterations, and CDK6 amplification in esophageal squamous cell carcinoma. Clin Cancer Res. 2014;20:1114–24.View ArticlePubMedGoogle Scholar
- Xu B, Chotewutmontri S, Wolf S, Klos U, Schmitz M, Dürst M, et al. Multiplex identification of human papillomavirus 16 DNA integration sites in cervical carcinomas. PLoS One. 2013;8:e66693.View ArticlePubMedPubMed CentralGoogle Scholar
- Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.View ArticlePubMedPubMed CentralGoogle Scholar
- Durst M, Bosch FX, Glitz D, Schneider A, zur Hausen H. Inverse relationship between human papillomavirus (HPV) type 16 early gene expression and cell differentiation in nude mouse epithelial cysts and tumors induced by HPV-positive human cell lines. J Virol. 1991;65:796–804.PubMedPubMed CentralGoogle Scholar
- Jeon S, Allen-Hoffmann BL, Lambert PF. Integration of human papillomavirus type 16 into the human genome correlates with a selective growth advantage of cells. J Virol. 1995;69:2989–97.PubMedPubMed CentralGoogle Scholar
- Daniel B, Rangarajan A, Geetasree M, Elizabeth V, Krishna S. The link between integration and expression of human papillomavirus type 16 genomes and cellular changes in the evolution of cervical intraepithelial neoplastic lesions. J Gen Virol. 1997;78:1095–101.View ArticlePubMedGoogle Scholar
- Lace MJ, Anson JR, Klussmann JP, Wang DH, Smith EM, Haugen TH, et al. Human papillomavirus type 16 (HPV-16) genomes integrated in head and neck cancers and in HPV-16-immortalized human keratinocyte clones express chimeric virus-cell mRNAs similar to those found in cervical cancers. J Virol. 2011;85:1645–54.View ArticlePubMedGoogle Scholar
- Hawkins TB, Dantzer J, Peters B, Dinauer M, Mockaitis K, Mooney S, et al. Identifying viral integration sites using SeqMap 2.0. Bioinformatics. 2011;27:720–2.View ArticlePubMedPubMed CentralGoogle Scholar
- Westermann AJ, Gorski SA, Vogel J. Dual RNA-seq of pathogen and host. Nat Rev Microbiol. 2012;10:618–30.View ArticlePubMedGoogle Scholar
- Bonfert T, Csaba G, Zimmer R, Friedel CC. Mining RNA–Seq Data for Infections and Contaminations. PLoS One. 2013;8:e73071.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen Y, Yao H, Thompson EJ, Tannir NM, Weinstein JN, Su X. VirusSeq: software to identify viruses and their integration sites using next-generation sequencing of human cancer tissue. Bioinformatics. 2013;29:266–7.View ArticlePubMedGoogle Scholar
- Li JW, Wan R, Yu CS, Wong N, Chan TF. ViralFusionSeq: accurately discover viral integration events and reconstruct fusion transcripts at single-base resolution. Bioinformatics. 2013;29:649–51.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang Q, Jia P, Zhao Z. VirusFinder: Software for Efficient and Accurate Detection of Viruses and Their Integration Sites in Host Genomes through Next Generation Sequencing Data. PLoS One. 2013;8:e64465.View ArticlePubMedPubMed CentralGoogle Scholar
- Katz JP, Pipas JM. SummonChimera infers integrated viral genomes with nucleotide precision from NGS data. BMC Bioinformatics. 2014;15:348.View ArticlePubMedPubMed CentralGoogle Scholar
- Lau CC, Sun T, Ching AKK, He M, Li JW, Wong AM, et al. Viral-Human Chimeric Transcript Predisposes Risk to Liver Cancer Development and Progression. Cancer Cell. 2014;25:1–15.View ArticleGoogle Scholar
- Johansson C, Schwartz S. Regulation of human papillomavirus gene expression by splicing and polyadenylation. Nat Rev Microbiol. 2013;11:239–51.View ArticlePubMedGoogle Scholar
- Doorbar J. The papillomavirus life cycle. J Clin Virol. 2005;32:7–15.View ArticleGoogle Scholar
- Wentzensen N, Ridder R, Klaes R, Vinokurova S, Schaefer U, Doeberitz M. Characterization of viral-cellular fusion transcripts in a large series of HPV16 and 18 positive anogenital lesions. Oncogene. 2002;21:419–26.View ArticlePubMedGoogle Scholar
- Kraus I, Driesch C, Vinokurova S, Hovig E, Schneider A, von Knebel Doeberitz M, et al. The majority of viral-cellular fusion transcripts in cervical carcinomas cotranscribe cellular sequences of known or predicted genes. Cancer Res. 2008;68:2514–22.View ArticlePubMedGoogle Scholar
- Ojesina AI, Lichtenstein L, Freeman SS, Pedamallu CS, Imaz-Rosshandler I, Pugh TJ, et al. Landscape of genomic alterations in cervical carcinomas. Nature. 2014;506:371–5.View ArticlePubMedGoogle Scholar
- Peter M, Stransky N, Couturier J, Hupé P, Barillot E, de Cremoux P, et al. Frequent genomic structural alterations at HPV insertion sites in cervical carcinoma. J Pathol. 2010;221:320–30.View ArticlePubMedGoogle Scholar
- Akagi K, Li J, Broutian TR, Padilla-Nash H, Xiao W, Jiang B, et al. Genome-wide analysis of HPV integration in human cancers reveals recurrent, focal genomic instability. Genome Res. 2014;24:185–99.View ArticlePubMedPubMed CentralGoogle Scholar
- Rabbitts TH. Chromosomal translocations in human cancer. Nature. 1994;372:143–9.View ArticlePubMedGoogle Scholar
- Hästbacka J, Kerrebrock A, Mokkala K, Clines G, Lovett M, Kaitila I, et al. Identification of the Finnish founder mutation for diastrophic dysplasia (DTD). Eur J Human Genet. 1999;7:664–7.View ArticleGoogle Scholar
- Sharma AK, Rigby AC, Alper SL. STAS domain structure and function. Cell Physiol Biochem. 2011;28:407–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Pett M, Coleman N. Integration of high‐risk human papillomavirus: a key event in cervical carcinogenesis? J Pathol. 2007;212:356–67.View ArticlePubMedGoogle Scholar
- Bodelon C, Vinokurova S, Sampson JN, den Boon JA, Walker JL, Horswill MA, et al. Chromosomal copy number alterations and HPV integration in cervical precancer and invasive cancer. Carcinogenesis. 2016;37:188–96.View ArticlePubMedGoogle Scholar
- Hu Z, Zhu D, Wang W, Li W, Jia W, Zeng X, et al. Genome-wide profiling of HPV integration in cervical cancer identifies clustered genomic hot spots and a potential microhomology-mediated integration mechanism. Nat Genet. 2015;47:158–63.View ArticlePubMedGoogle Scholar
- Weitzman MD, Weitzman JB. What’s the damage? The impact of pathogens on pathways that maintain host genome integrity. Cell Host Microbe. 2014;15:283–94.View ArticlePubMedPubMed CentralGoogle Scholar
- White AE, Livanos EM, Tlsty TD. Differential disruption of genomic integrity and cell cycle regulation in normal human fibroblasts by the HPV oncoproteins. Genes & Dev. 1994;8:666–77.View ArticleGoogle Scholar
- Kessis TD, Connolly DC, Hedrick L, Cho KR. Expression of HPV16 E6 or E7 increases integration of foreign DNA. Oncogene. 1996;13:427–31.PubMedGoogle Scholar
- Duensing S, Lee LY, Duensing A, Basile J, Piboonniyom SO, Gonzalez S, et al. The human papillomavirus type 16 E6 and E7 oncoproteins cooperate to induce mitotic defects and genomic instability by uncoupling centrosome duplication from the cell division cycle. Proc Natl Acad Sci U S A. 2000;97:10002–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Duensing S, Münger K. The human papillomavirus type 16 E6 and E7 oncoproteins independently induce numerical and structural chromosome instability. Cancer Res. 2002;62:7075–82.PubMedGoogle Scholar
- Bester AC, Roniger M, Oren YS, Im MM, Sarni D, Chaoat M, et al. Nucleotide deficiency promotes genomic instability in early stages of cancer development. Cell. 2011;145:435–46.View ArticlePubMedPubMed CentralGoogle Scholar
- Havre PA, Yuan J, Hedrick L, Cho KR, Glazer PM. p53 inactivation by HPV16 E6 results in increased mutagenesis in human cells. Cancer Res. 1995;55:4420–4.PubMedGoogle Scholar
- Carter SL, Eklund AC, Kohane IS, Harris LN, Szallasi Z. A signature of chromosomal instability inferred from gene expression profiles predicts clinical outcome in multiple human cancers. Nat Genet. 2006;38:1043–8.View ArticlePubMedGoogle Scholar
- How C, Bruce J, So J, Pintilie M, Haibe-Kains B, Hui A, et al. Chromosomal instability as a prognostic marker in cervical cancer. BMC Cancer. 2015;15:1.View ArticleGoogle Scholar
- Samanta S, Dey P, Nijhawan R. Micronucleus in Cervical Intraepithelial Lesions and Carcinoma. Acta Cytol. 2011;55:42–7.View ArticlePubMedGoogle Scholar
- Zhang CZ, Spektor A, Cornils H, Francis JM, Jackson EK, Liu S, et al. Chromothripsis from DNA damage in micronuclei. Nature. 2015;522:179–84.View ArticlePubMedPubMed CentralGoogle Scholar
- Sichero L, Sobrinho JS, Villa LL. Oncogenic potential diverge among human papillomavirus type 16 natural variants. Virology. 2012;432:127–32.View ArticlePubMedGoogle Scholar
- Hochmann J, Sobrinho JS, Villa LL, Sichero L. The Asian-American variant of human papillomavirus type 16 exhibits higher activation of MAPK and PI3K/AKT signaling pathways, transformation, migration and invasion of primary human keratinocytes. Virology. 2016;492:145.View ArticlePubMedGoogle Scholar
- Zacapala-Gómez AE, Del Moral-Hernández O, Villegas-Sepúlveda N, Hidalgo-Miranda A, Romero-Córdoba SL, Beltrán-Anaya FO, et al. Changes in global gene expression profiles induced by HPV 16 E6 oncoprotein variants in cervical carcinoma C33-A cells. Virology. 2016;488:187–95.View ArticlePubMedGoogle Scholar
- Muller E, Brault B, Holmes A, Legros A, Jeannot E, Campitelli M, et al. Genetic profiles of cervical tumors by high-throughput sequencing for personalized medical care. Cancer Med. 2015;4:1484–93.View ArticlePubMedPubMed CentralGoogle Scholar
- Giardine B, Riemer C, Hardison RC, Burhans R, Elnitski L, Shah P, et al. Galaxy: a platform for interactive large-scale genome analysis. Genome Res. 2005;15:1451–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Goecks J, Nekrutenko A, Taylor J. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol. 2010;11:R86.View ArticlePubMedPubMed CentralGoogle Scholar
- Blankenberg D, Kuster GV, Coraor N, Ananda G, Lazarus R, Mangan M, et al. Galaxy: a web‐based genome analysis tool for experimentalists. Curr Protoc Mol Biol. 2010;Chapter 19:Unit 19.10.1-21Google Scholar
- Blankenberg D, Gordon A, Von Kuster G, Coraor N, Taylor J, Nekrutenko A. Manipulation of FASTQ data with Galaxy. Bioinformatics. 2010;26:1783–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Picard Tools. http://broadinstitute.github.io/picard/. Accessed 13 May 2016
- Flores ER, Allen-Hoffmann BL, Lee D, Sattler CA, Lambert PF. Establishment of the human papillomavirus type 16 (HPV-16) life cycle in an immortalized human foreskin keratinocyte cell line. Virology. 1999;262:344–54.View ArticlePubMedGoogle Scholar
- Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562–78.View ArticlePubMedPubMed CentralGoogle Scholar
- R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing. 2013. http://www.R-project.org
- gplots Package for R. http://cran.r-project.org/web/packages/gplots/gplots.pdf
- Anders S, Pyl PT, Huber W. HTSeq — A Python framework to work with high-throughput sequencing data. bioRxiv. 2014. doi:https://doi.org/10.1101/002824.Google Scholar
- Carbon S, Ireland A, Mungall CJ, Shu S, Marshall B, Lewis S. AmiGO: online access to ontology and annotation data. Bioinformatics. 2009;25:288–9.View ArticlePubMedGoogle Scholar
- Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011;27:431–2.View ArticlePubMedGoogle Scholar