A fresh tissue specimen was collected from a patient with IIIb squamous cell carcinoma who had undergone surgery at Anyang Cancer Hospital (Henan province, China) in 2009. The patient was 56 years old and had HPV16-positive cancer. This study received ethical approval from the Institutional Review Board of the hospital. Individual informed consent had been collected from this study participant.

Genomic DNA, provided by Peking University Cancer Hospital and Institute (Beijing, China), was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. Double-stranded (ds) DNA was quantified using the Nanodrop 2000 and Qubit dsDNA HS Assay Kit (both from Thermo Fisher Scientific, Inc., Waltham, MA, United States). The average fragment size of DNA [>5 kilobase pair (kbp)] was measured (identified by comparison to DL2000 PLUS DNA Ladder, Life Technologies, Carlsbad, CA, United States) on a 1.0% agarose gel in 1X TAE buffer on the Bio-Rad CHEF DRII system.

Sheared DNA was used to make genomic DNA libraries using the NEBNext^® UltraII^TM DNA Library Prep Kit for Illumina^® (E7370L) and NEBNext^® Multiplex Oligos for Illumina^® [E6609L; New England Biolabs (NEB), Ipswich, MA, United States] as recommended by the manufacturer. In short, 200 ng DNA was sheared in a 50 μL volume using the Covaris M220 Focused-ultrasonicator (Covaris, Cambridge, MA, United States) for 340 s with a duty cycle of 10%, cycles per burst of 200, and peak power of 75. Then, fragmented DNA was end-repaired and A-tailed using the NEBNext^® UltraII^TM Mix and buffer, and ligation of DNA adapters was done using the NEBNext^® UltraII^TM Ligation Master Mix and Enhancer. Then, 0.9 × AMPure XP beads (Beckman Coulter Inc., Skyesville, MD, United States) were used to remove all traces of adapter dimers. The adapter-ligated DNA was amplified with six cycles of PCR using the NEBNext^® UltraII^TM Q5 Master Mix. The generated library was quantified with the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific) using the Qubit 3.0 fluorometer (Invitrogen, Carlsbad, CA, United States), and run on the Qsep1^TM biosystem (BiOptic Inc., NTC, Taiwan) for quality analysis. The final size of the electrophoresis fragment was about 320 bp.

Human papillomavirus probes were designed by Integrated DNA Technologies (Integrated DNA Technologies, Coralville, LA, United States), according to the full-length genome of HPV16. Overall NGS target enrichment with xGen Hybridization and a Wash Kit was conducted by Integrated DNA Technologies, which is briefly described below and performed as detailed by the plate standard protocol for xGen^® hybridization capture of DNA libraries from Integrated DNA Technologies. In brief, whole-genomic libraries were hybridized with HPV probes (Integrated DNA Technologies), absorbed onto the beads via Dynabeads MyOne Streptavidin T1 (Thermo Fisher Scientific), after which the uncaptured DNA fragments were removed by washing. Then, the eluted fragments containing the targeted gene were enriched by 15 cycles of PCR to generate libraries for sequencing. The captured library was quantified with the Qubit dsDNA HS Assay Kit using a Qubit 3.0 fluorometer (both from Thermo Fisher Scientific) and run on the Agilent 2100 TapeStation (Agilent Technologies, Santa Clara, CA, United States) for quality analysis prior to sequencing. DNA libraries were sequenced using the NextSeq platform (Illumina) with 150 bp paired-end reads.

The sheared DNA was used to make genomic DNA libraries using the NEBNext^® UltraII^TM DNA Library Prep kit for Illumina^® (E7370L) and NEBNext^® Multiplex Oligos for Illumina^® (E6609L), as recommended by the manufacturer, followed by DNA shearing with the M220 Focused-ultrasonicator (Covaris). Then, 200 ng genomic DNA was sheared into longer fragment sizes of about 800 bp in a 50 μL volume on the Covaris M220 (Covaris) for 50 s with a duty cycle of 20%, cycles per burst of 200, and peak power of 50 according to the manufacturer’s instructions. Then, the fragmented DNA was end-repaired and A-tailed using the NEBNext^® UltraII^TM Mix and buffer, and ligation of adapters was done using the NEBNext^® UltraII^TM Ligation Master Mix and Enhancer, after which 0.9X AMPure XP beads (Beckman Coulter) were used to remove all traces of adapter dimers. The adapter-ligated DNA was amplified with six cycles of PCR using the NEBNext^® UltraII^TM Q5 Master Mix. The generated library was quantified with the Qubit dsDNA HS Assay Kit using the Qubit 3.0 fluorometer (both from Thermo Fisher Scientific) and run on the Qsep1^TM biosystem (BiOptic) for quality analysis. The final size of the electrophoresis fragment ranged from 250 to 1500 bp. The overall capture experiment with xGen Hybridization and Wash Kit was conducted by Integrated DNA Technologies as detailed by the standard protocol for xGen^® hybridization capture of DNA libraries from Integrated DNA Technologies with slight modifications. Briefly, the HYB program extends for 4 to 6 h that incubating the tubes at 65°C with a heated lid set at 105°C. The captured library was quantified with the Qubit dsDNA HS Assay Kit using the Qubit 3.0 fluorometer (both from Thermo Fisher Scientific) and run on the Agilent 2100 Tape Station (Agilent Technologies) for quality analysis prior to nanopore sequencing. DNA libraries were sequenced using MinION^TM MK1B device from Oxford Nanopore Technologies (ONT; Oxford, United Kingdom).

The capture of HPV 16 probe-template duplexes was done using the Ligation Sequencing Kit (SQK-LSK108) with Native Barcoding Expansion (EXP-NBD103) from ONT, and performed as detailed by the 1D Native barcoding genomic DNA standard protocol from ONT. Then, 310 ng purified amplicon DNA was dA-tailed using the NEBNext Ultra II End Repair/dA Tailing module (E7546S; NEB) at a temperature of 20°C for 30 min and at 65°C for 30 min using the thermal cycle. Then the DNA was purified with AMPure XP beads (A63881; Beckman Coulter), and 60 μL (1×) of re-suspended AMPure XP beads was added by pipetting and thoroughly mixed on a rotator mixer for 5 min. The beads were separated on a magnet to remove the supernatant. The beads were kept on the magnet and washed twice using 200 μL of 70% ethanol without disturbing the pellet, followed by re-suspension of the pellet in 25 μL nuclease-free water. Then, 275 ng end-prepared DNA was ligated with native barcode NB04 from ONT using the NEB Blunt/TA Ligase Master Mix (M0367S; NEB), and then incubated at 25°C for 15 min. Following the barcode ligation reaction, the DNA was cleaned again with 1X Agencourt AMPure XP beads (Beckman Coulter) and eluted in 26 μL. The resulting DNA was pooled with another unrelated barcoded library by equivalent weight. Then, 543 ng pooled barcode DNA was used to perform the adapter ligation step, and 20 μL Barcode Adapter Mix, 10 μL Quick T4 DNA Ligase, and 20 μL of NEB Next Quick Ligation Reaction Buffer (5X; E6056S; NEB) were added in that order to the pooled barcoded 50 μL DNA. The reaction was incubated at room temperature for 15 min, and then cleaned using 0.4X Agencourt AMPure XP beads (Beckman Coulter). Then, the beads were washed twice with 140 μL Adapter Bead Binding, re-suspended in 15 μL Elution Buffer, and incubated for 10 min at room temperature before pelleting in a magnetic rack. The prepared library was combined with 35 μL running buffer with Fuel Mix, 25.5 μL Library Loading beads and 12 μL DNA library (∼158 ng) and loaded into the SpotON sample port of the R9.4 flowcell. All nanopore sequencing runs were conducted using the MinION^TM MK1B device following the recommended sequencing protocols by ONT. The captured library was sequenced individually using the FLO-MIN-106 flow cell. MinKNOW software (v18.07.18) was used to control the MinION device following a 24-h run script, and sequencing data were collected in real-time and processed into basecalls using the Metrichor^TM agent.

The nanopore sequencing results were basecalled using the EPI2ME interface (v. 2.59.1896509). For passed 1D reads, quality scores (Q-score ≥ 7) and length distributions were evaluated using FastQC. Obtained fastq files were converted to fasta files using the FASTX-toolkit¹. These fasta files were aligned to the Illumina adapter database by LAST (version: lastal 956; Frith et al., 2010; Kiełbasa et al., 2011), which alignment parameters were mismatch cost: 1, gap extension cost: 1, gap existence cost: 1, and minimum score for gapped alignments: 45. Then these reads were divided into two group reads, which we defined that the only library (OL) structure reads and multiple library (ML) structure reads (Figure 2B). The ML reads were split to the OL reads by the Illumina adapter sequences. All OL reads were aligned to the human genome (GRCh37/hg19) and HPV genome (NC_001526.2) databases using Blat (version: v36; Kent, 2002), with parameters were stepSize = 5, repMatch = 2,253, minScore = 20, and minIdentity = 0. We also defined the highest score alignment result as the best OL reads alignment result to find the breakpoint of HPV integration. We used Samtools (version: v1.9) software to analyze the HPV16 genome sequencing depth. In each one OL reads, we defined the breakpoint of HPV integration which the gap and overlap between the best human alignment result and the best HPV alignment result less than 10 bp. The breakpoint was annotated to gene function using ANNOVAR² and an in-house program with human reference and HPV16 reference, respectively. The sequencing reads generated by Illumina sequencer were analyzed using the SEGF pipeline (Xu et al., 2018). A diagram of the integration site distribution was made using Illustrator software. Gene function analysis was applied using the DAVID online analysis tool (Huang et al., 2009a, b).

To validate the selected HPV integration into the human genome detected by nanopore sequencing, primers were designed, one of which was derived from the human genome at the potential site of integration, and the other of which was against HPV sequences suspected of being near the site of integration within the HPV genome. Both primers were designed 100∼400 bp away from the detected integration sites based on the nanopore targeted-capture reads, and 13 HPV integration sites (n ≥ 5) were selected. PCR was performed using each primer set for DNA that was previously used for HPV capture. As a template, 1 μL (30 ng) genomic DNA solution was used in the subsequent PCR. The PCR reaction mix was prepared in a total volume of 20 μL containing 4 μL 5X Phoenix Hot Start Taq Reaction Buffer, 2 μL dNTP (2.5 mM), 0.5 μL of each forward and reverse primer (10 mM), 0.2 μL Phoenix Hot Start Taq DNA Polymerase (500 U), and 11.8 μL nuclease-free water (not DEPC-Treated). The PCR conditions were as follows: 5 min at 95°C; 35 cycles of 30 s at 94°C; 60 s at 50–60°C for annealing; and 60 s at 72°C; followed by 72°C for 1 min. PCR products were analyzed by electrophoresis on 2.0% agarose gels, purified, and Sanger sequenced.

The captured library was constructed and sequenced by ONT MinION. The sequencing run generated 381,475 (149.1 Mb) sequencing reads from 977 active pores. To obtain high-quality reads, the raw reads were filtered using the Metrichore 1D base calling program and kept for further analysis if a Qscore ≥ 7 was obtained. In total, 333,028 1D reads (130.2 Mb) were retained with the read lengths ranging from 67 to 7,404 bp, with a mean read length of 502.8 bp, and the quality score ranged from 7 to 17, with a mean value of 12.8. The distribution diagram of the read lengths and quality scores are shown in Figure 1 and other sequence details are summarized in Table 1. The captured Illumina library was sequenced by Illumina sequencing machine and 11,093,630 raw reads (1.67 Gb raw bases) were generated. After data cleaning, 11,093,535 clean reads (1.64 Gb) were obtained with an average length of 151 bp. The summary of sequencing results from the two platforms are listed in Table 1 and the sequence coverage of HPV16 genome by Illumina and nanopore sequencer was shown in Supplementary Figure S1. It showed that HPV16 full genome was covered by both sequencing platforms and Illumina result provided obviously higher sequence depth than nanopore. The sequence depth of the most HPV16 genome regions are even both on Illumina and nanopore results, excepting a small region around 4,000 bp location in Illumina. Although it showed a significant decrease of the sequence depth in this region, it provided more than 100× sequencing depth of it.

The sequencing data from nanopore sequencing were analyzed according to the pipeline described in the Materials and Methods, shown in Figure 2A. There were 12,951 sequence reads in the human genome and 389,839 reads in the HPV genome. Next, the potential integrated reads were analyzed using BLAST with the human and HPV genome sequences to determine the exact integration breakpoints. Without performing error correction analysis, 7,859 reads were identified in both the human and HPV sequences. To filter out the chimeric product from the library construction process, we used the number of overlapped bps as filter criteria. The HPV and human part of the same sequencing read were identified with different bp-length gaps. Only reads, whose gap was shorter than 10 bp length, were kept for further analysis. A total of 339 integration sites were finally identified (Supplementary Table S1). The location and gene information of the matched sequence in either the human or HPV genome were annotated. The read number of unique breakpoints was calculated and only those breakpoints with ≥ 2 reads identified were kept for further comparison analysis with two closed sites (HPV16:3327,chr6:7328094 and HPV16:3329,chr6:7328093)combined. Excepting 2 site also identified by Illumina with 3 reads (HPV16:4405,chr2:99438640 and HPV16:2716,chr13:74250432; Supplementary Table S2). There were 60 breakpoints in total were included. The distribution of all breakpoints in the different chromosomes is shown in Figure 3A. Chromosomes chr20, chr2, chr6 were the three chromosomes with 12, 11, and 8 breakpoints, respectively. Both chr1 and chrX had 7 and 7 integration sites each and the rest of the chromosomes had ≤5 integration sites. The identified read numbers of each integration site were highly different, with a range change from 2 to 406. There were 31.7% (19 of 60) integration sites had more than 10 reads. The top 19 abundant integration sites were distributed in chromosomes 20(6), 6(4), X (4), 11(3), and 2(2). Among the top five most abundant breakpoints, There were 3 located in chromosome 20, HPV16:2804,chr20:32516985 (406 reads), HPV16:7139,chr20:32478733 (169 reads), HPV16:4276,chr20:32502143(51 reads), and two were in chrX, HPV16:5534,chrX:20464412 (132 reads), HPV16:3163,chrX: 20462930 (50 reads). Figure 3B indicates the integration site distribution in the HPV genome. Within the 60 identified integration sites, 18 were in the L2 gene, 12 in L1, 11 in E1, 7 in E2, and <4 in the each of the other genes (E6, LCR, E5, and E7). On the human genome, there are 38 integration sites in the intergenic region, 18 in the intronic region, 2 in the ncRNA region, 1 in the exonic and 1 in the downstream region of the genes, as indicated in Figure 3C. To further look at the breakpoints in the same chromosome, in chr20, most breakpoints (10 of 12) located in the intergenic region of two genes CHMP4B, RALY-AS1, one of the other two breakpoints in intergenic region of KIF3B, ASXL1, and the other is in the intronic region of gene ATRN. With the exception of chr20, other chromosomes also have obvious cluster tendency of integration sites, like six integration sites of chrX in the intergenic region of genes RPS6KA3 and CNKSR2. Three breakpoints in chr6 were in the intronic or downstream of gene CAGE1. 10 regions of human genome were identified with more than two breakpoints clustered together, indicating there were integration sites cluster in the human genome. In addition, 14 breakpoints without any cluster tendency were observed. Illumina identified 1,718 integration sites, and only the integration sites with ≥ 3 reads were collected for further analysis. All sites are listed in Supplementary Table S2. The characteristic analyses of these integration sites are summarized in Figure 3.

Capture-based Illumina identified 54 integration sites with ≥ 3 reads sequenced, and four sites (HPV16:7619,chr2:99404946 and HPV16:6558,chr2:99431140, HPV16:850,chr20:32515476, and HPV16:462,chr20:30942164) with less than three but overlapping with nanopore results showed in Table 2. Altogether, there are 54 integration sites were included for further analysis. A previously published paper (Liu et al., 2016a) reported 19 integration sites from the same sample we used. All integration sites were compared and the overlapping sites by either method were recorded. The numbers are shown in the Venn diagram in Figure 4A, and the detailed list of integration sites is shown in Figure 4B and Table 2. There were 13 integration sites identified by all three platforms, which indicated good repetitiveness between the different platforms. A total of 18 sites overlapped by nanopore and Illumina, three overlapped by nanopore sequencing and in the previous paper, and only one overlapped by Illumina and the published paper. Among the identified integration sites in the two platforms, some breakpoints had highly abundant read numbers, such as HPV16:2804,chr20:32516985 with 406 reads in the nanopore results and 1439 reads in the Illumina results. Other breakpoints had very few reads numbers such as HPV16:4250,chr21:97550764 with two reads in the nanopore results, whereas there were only four reads in the Illumina results. Besides the overlapping sites, each platform had their own unique integration sites. Nanopore, Illumina, and Liu et al. had 26, 22, and 2 unique integration sites, respectively. A total of 6 of the 26 nanopore integration sites were highly abundant with ≥ 6 reads sequenced, namely HPV16:7139,chr20:32478733 (169 reads), HPV16:4276, hr20:32502143 (51 reads), HPV16:4029,chr11:103891933 (21 reads), HPV16:3312, chr6:17081322 (8 reads), and HPV16:3295,chr11:31761709 (6 reads), HPV16:5311, chr1:455824 (6 reads). Within the integration sites identified by all the three methods, Illumina, nanopore and Liu paper, with relative reliable abundance, nanopore identified the most validated integration sites. All the unique integration sites and their detailed information are listed in Table 2. To test the accuracy of the new integration sites identified by our data, we performed Sanger sequencing. Together, 14 integration sites were chosen for verification, and 13 were PCR-amplified and successfully sequenced, indicating the true positive of these integration sites, which were mostly identified by both nanopore and Illumina platforms. Detailed information on the designed primers and sequencing results with chromas images are summarized in Supplementary Table S3 and Supplementary Figure S2. The verified integration sites show different coexisting situations. Eight coexisted in nanopore and Illumina integration sites datasets, which indicating the advantage of our integration detection pipeline of both platforms compared with previouse published pipeline. There were 9 integration sites with more than 10 identified reads in each of the three datasets, and four integration sites with less than 10 reads in all three datasets were verified. The high verification rate (92.86%) of integration sites also indicated the high integration sites detection accuracy. Also, four sites were verified to exist in the intronic region, and the remaining nine integration sites were found to exist in the intergenic region. The PCR gel image and the sequencing result were shown in Figure 5 and Supplementary Table S3, respectively.

All the integration sites identified by the three platforms were combined to the final unique integration sites list (Table 2) for further analysis with nanopore reads (n ≥ 2) and Illumina reads (n ≥ 3). A total of 83 unique integration sites were found, each of which were annotated, and further function classification and pathway analyses were conducted. If the integration sites located in intergenic region, the closer gene was used for analysis. The genes were clustered into eight biology processes including positive regulation of transcription from RNA polymerase II promoter (6 genes), negative regulation of transcription from RNA polymerase II promoter (6 genes), regulation of transcription from RNA polymerase II promoter (5 genes), transcription from RNA polymerase II promoter (4 genes), negative regulation of retinoic acid receptor signaling pathway (2 genes),establishment of skin barrier (2 genes), proximal/distal pattern formation (2 genes) and embryonic limb morphogenesis (2 genes) and six cytoplasm (13 genes), nucleoplasm (8 genes), extracellular exosome (8 genes), protein binding (6 genes), DNA binding (2 genes) and actin binding (2 genes). Detailed information on the function classification results is shown in Supplementary Figure S3 and Supplementary Table S4.