Universal human papillomavirus typing by whole genome sequencing following target enrichment: evaluation of assay reproducibility and limit of detection

Human Papillomaviruses (HPV) are a group of double-stranded DNA viruses that can cause genital warts and are in the causal pathway for anogenital and oropharyngeal cancers [1]. There are more than 200 HPV types in 5 of the 45 genera within Papillomaviridae (Alphapapillomavirus, Betapapillomavirus, Gammapapillomavirus, Mupapillomavirus and Nupapillomavirus). Most assays focus on detection and typing of a subset of Alphapapillomavirus HPV types recognized to be clinically important. Conventional typing assays do not capture information on additional viral characteristics, such as integration and variant status within HPV type, therefore, systematic evaluation of the clinical importance of this information is lacking. While next generation sequencing (NGS) could address this limitation, current HPV NGS methods rely on sequencing PCR amplicons targeting a limited region of the genome for genotyping [2–7] or for identifying the integration/variant status [8–16] of a restricted number of types. To overcome these limitations, we recently developed a target enriched whole genome sequencing method (eWGS) designed to identify all known and potentially novel HPV types in a given sample [17]. The resulting whole genome sequence data is useful to address variant and integration status.

Our original eWGS report provided details on the method [target enrichment with RNA baits (based on Agilent SureSelect technology), library preparation, and sequencing using Illumina HiSeq 2500 platform] as well as initial performance metrics for HPV type determination such as genome coverage and uniformity, but reproducibility and limit of detection (LOD) were not addressed. Given the complexity of the workflow for NGS methods, determination of these important characteristics is difficult. A few amplicon-based HPV NGS methods have reported reproducibility in terms of concordance at the level of final type determination; however, these studies did not give detailed measures of variability at the level of base quality, coverage, and mapped reads [2, 5, 6]. Moreover, methods evaluating reproducibility and sensitivity of amplicon-based NGS methods targeting a limited region of the genome are not directly applicable to our conceptually different whole-genome enriched NGS method. This study addresses reproducibility and LOD for type determination based on results in defined HPV samples using our eWGS method. We used overall quality of sequencing reads, number of reads mapped to reference genomes, average depth of coverage, and fraction of genome covered by mapped reads to measure reproducibility and results on samples with decreasing copy number to determine LOD. We find that our eWGS method is highly reproducible for HPV type determination with an LOD of 25 copies/reaction even under the scenario of infection with multiple types.

This study provides an in-depth evaluation of the reproducibility and LOD for HPV genotyping with our recently described eWGS method. The approach, using defined samples varying in complexity and copy number, analyzed in replicate and duplicate assays, is applicable to most NGS methods. Importantly, the study design considered reproducibility of the enrichment and library preparation steps as well as sequencing. We chose to use reproducibility of HPV genotyping using reads mapped under the most stringent L1S1 mapping conditions, while varying cut-off criteria for number of reads, depth of coverage, and fraction of reference genome covered. Varying input amounts of HPV genomes (1–625 copies/reaction) and the number of types in the samples provided additional insights into assay robustness and LOD.

eWGS showed high reproducibility in cluster generation between flow cell lanes (CV: 5.71 and 0.56% for experiment 1 and 2, respectively) and between experiments (CV: 13.4%). Consistency of cluster generation reflects consistency during the complex workflow of library preparation and DNA determination of recovered products added to the flow cell. The eWGS results from the 4 replicates demonstrates consistency in the number and quality of reads (mean CV: 14.7% for total number of reads; 0.5% for base quality score; 0.61% for bases with Q score ≥ 30). The number of reads mapped to HPV genomes were highly correlated between the two flow-cell lanes (R² = 1) in both experiments and between the two experiments (R² = 0.99). Importantly, eWGS results in terms of number of mapped reads (mean CV: 10.1%; range: 2.5–23.8%) and fraction of reference covered (mean CV: 0.5%; range: 0–1.8%) at copies/reaction ranging from 125 to 625 was consistent for all 9 HPV plasmids. Using all three cut-off criteria, HPV genotype calling was fully reproducible at 625 copies for all 9 HPV plasmids, while at 125 copies only 6 out of 9 types were called reproducibly among the 4 replicates. HPV 6, 16 and 31 did not meet the depth of coverage criteria ≥20 at 125 copies in all replicates. These results, taken together, suggest that despite the complex nature of its workflow eWGS is robust and highly reproducible for HPV genotyping on a whole genome level at 625 copies/reaction using all three cut-off criteria, with no significant type-specific differences.

As expected, we observed a positive correlation between the number of mapped reads and copy number/reaction (β = 13.9, p < 0.0001), as well as the fraction of reference genome covered. Variability among replicates for both measures increased as HPV copy number/reaction decreased. The positive correlation of mapped reads and fraction of genome covered with target concentration as we observed with bait-based eWGS for HPV agrees with a recently reported study on whole genome sequencing of hepatitis C viral genomes following similar target enrichment [23]. However, PCR-based methods either targeting whole genome or amplicon sequencing did not show a copy number dependence with the number of mapped reads [2, 23]. We found that eWGS method for HPV genotyping outperformed an NGS method based on amplicon sequencing [2] with less variability among the 9 vaccine types in the pooled sample for the number of mapped reads at 25 copies/reaction (CV: eWGS, 29.6%; amplicon sequencing, 125.3%). Even at higher input (50–500 copies/reaction) of the 9 vaccine types in the pool, variability in the number of HPV reads with amplicon sequencing remained high (CV: 99.1 to 131.0%) [5]. Amplicon sequencing of replicates of clinical samples for HPV detection also showed high variability in the number of reads between libraries prepared to assess reproducibility (CV: mean 81%; range 2–141%) [2].

Because it is common for multiple HPV types to be present in clinical samples, competition among types for PCR amplification in assays relying on consensus primers could result in a bias in detection [24–26]. The eWGS method should minimize detection bias. We evaluated this by comparing the number of reads for each HPV type at 625 copies in mixed type plasmid pools (simulating multiple infections) and single plasmid pools. The numbers of reads did not vary significantly (p = 0.16–0.99) on 8 of the 9 plasmids, suggesting eWGS has minimal bias for HPV typing. The data on reproducibility among the replicates (Table 1) with HPV copies down to 25 copies/reaction (CV for number of reads =29.6%) also indicates low bias for HPV typing.

Reported LODs for NGS assays were determined with cut-offs based on target specific reads and did not consider how LOB could be used in determining the LOD. In this study, we used the number of reads mapped to HPV reference genomes as a parameter in relation to the LOB to calculate LOD [21, 22, 27]. Because we used defined samples in this study, reads mapping to HPV types not included in the sample served as a good measure of false positive reads or “noise”. Using this calculation, 25 copies/reaction was determined to be the LOD for the 9 types with eWGS method. At an LOD of 25 copies/reaction, eWGS averaged 322 ± 95.4 mapped reads for HPV genotype calling, higher than the 1–5 mapped reads used as threshold by some amplicon sequencing methods [2, 5]. This LOD compares favorably with PCR methods for HPV detection and typing that have reported sensitivity ranging from 50 to 500 copies/reaction [5, 20] and sensitivity at 100 copies/reaction for another HPV genotyping method using hybrid capture combined with isothermal whole genome amplification and Luminex detection [19]. There are a few studies [23, 28–32] addressing the determination of LOD using NGS methods for viral detection but the units of LOD (genome copies/ml, pfu/ml, IU/ml) reported in these studies are not comparable to genome copy number/reaction used in this study for HPV typing.

Our study is not without limitations, and additional work will be required to evaluate eWGS performance using clinical samples, and for other variables such as inter-operators, instrument, and laboratory variation. Determination of analytical sensitivity of an HPV detection assay is challenging because biologic samples vary widely in the proportion of normal and infected/neoplastic cells. Therefore, the best option for HPV is the determination of sensitivity denoted as copies or genomes/reaction as reported previously for NGS or PCR or hybrid capture based assays [2, 5, 18, 19], and in WHO’s proficiency study of HPV genotyping [20]. Assay performance is likely to vary when applied to standard testing conditions; the challenge of determining true results in undefined samples can only be addressed indirectly using inter-assay comparisons. Cost reductions could be achieved through automating the steps of library production and by increasing the number of samples that are pooled for the sequencing reaction. In our continued work (data not shown), we tested a total of 18 HPV types from mixing 4–7 types/reaction, and found that the number of samples/sequencing lane can be increased from 16 to 32 while maintaining the level of sensitivity (625 copies/reaction) for type determination (based on all three cut-off criteria combined) for all 18 HPV types. Additional work is needed to determine the maximum number of samples that can be pooled per sequencing reaction without compromising sensitivity. The focus of our current study was to verify sensitivity and reproducibility of HPV detection at the level of type. We are developing a bioinformatics pipeline for HPV genotyping that can be modified to identify variant and integration status from the same capture DNA sequencing data. HPV integration may be inferred from identifying single end or paired end reads that map to both the human and viral reference genomes as well as by identifying missing segments of HPV genome, as noted in our first publication [17]. In summary, we report an in-depth evaluation of the reproducibility of eWGS method for HPV genotyping employing a variety of metrics that include overall quality of reads, number of reads mapped, depth of coverage, and fraction of reference genome covered. Results indicate eWGS is highly reproducible for HPV genotyping at 625 copies /reaction using multiple cut-off criteria, with the possibility of reducing LOD to 25 copies/reaction, if HPV detection is based on the widely used number of mapped read criteria alone.

Our data suggest that eWGS method reduces type-competition, and has sensitivity competitive with widely used consensus PCR methods for HPV genotyping. In addition to genotyping, eWGS has the potential to provide highly reproducible and less biased sequence data for variant determination and identifying integration status, but may require additional studies to determine LODs specific to these applications of eWGS. The protocol used in this study, involving defined samples varying in complexity and copy number, analyzed in replicate and duplicate assays, is applicable to most WGS methods.