ve-SEQ: Robust, unbiased enrichment for streamlined detection and whole-genome sequencing of HCV and other highly diverse pathogens

Sample collection and preparation

Samples for optimization of sequencing methods were acquired from HCV Research UK (http://www.hcvresearchuk.org/), whose clinical samples were used with informed consent, conforming to the ethical guidelines of the 1975 Declaration of Helsinki. Study protocols were approved by the NRES Committee East Midlands, Derby (Ethics reference 11/EM/0323). Samples for resistance testing were obtained from patients enrolled and consented as part of the OxBRC Prospective Cohort Study in Hepatitis C (Ethics reference 09/H0604/20) at the Oxford University Hospitals NHS Trust.

Patient plasma was collected from EDTA blood tubes by centrifugation for 10 minutes at 600g in a Heraeus Megafuge, and stored at -80°C. RNA was isolated from 500µl plasma volumes using the NucliSENS magnetic extraction system (bioMerieux) and collected in 30µl of kit elution buffer for storage in aliquots at -80°C.

Sequencing library construction, enrichment and sequencing

Libraries were prepared for Illumina sequencing using the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (New England Biolabs) with 5µl sample (maximum 10ng total RNA) and previously published modifications of the manufacturer’s guidelines (v2.0)¹³, briefly: fragmentation for 5 or 12 minutes at 94°C, omission of Actinomycin D at first-strand reverse transcription, library amplification for 15–18 PCR cycles using custom indexed primers¹⁹ and post-PCR clean-up with 0.85× volume Ampure XP (Beckman Coulter).

Libraries were quantified using Quant-iT™ PicoGreen^® dsDNA Assay Kit (Invitrogen) and analysed using Agilent TapeStation with D1K High Sensitivity kit (Agilent) for equimolar pooling, then re-normalized by qPCR using the KAPA SYBR^® FAST qPCR Kit (Kapa Biosystems) for sequencing. Metagenomic virus RNA-Seq libraries were sequenced with 100b paired-end reads on the Illumina HiSeq 2500 with v3 Rapid chemistry.

A 500ng aliquot of the pooled library was enriched using the xGen^® Lockdown^® protocol from IDT (Rapid Protocol for DNA Probe Hybridization and Target Capture Using an Illumina TruSeq^® or Ion Torrent^® Library (v1.0), Integrated DNA Technologies) with equimolar-pooled 120nt DNA oligonucleotide probes (IDT) followed by a 12-cycle, modified, on-bead, post-enrichment PCR re-amplification. The cleaned post-enrichment ve-Seq library was normalized with the aid of qPCR and sequenced with 100b paired-end reads on a single run of the Illumina MiSeq using v2 chemistry.

Sequence data analysis

De-multiplexed sequence read-pairs were trimmed of low-quality bases using QUASR v7.01²⁰ and adapter sequences with CutAdapt version 1.7.1²¹ and subsequently discarded if either read had less than 50b remaining sequence or if both reads matched the human reference sequence using Bowtie version 2.2.4²². The remaining read pool was screened against a BLASTn database containing all 165 ICTV (International Committee on the Taxonomy of Viruses) HCV genomes (http://talk.ictvonline.org/ictv_wikis/m/files_flavi/default.aspx) both to choose an appropriate reference and to select those reads which formed a majority population for de novo assembly with Vicuna v1.3²³ and finishing with V-FAT v1.0 (http://www.broadinstitute.org/scientific-community/science/projects/viral-genomics/v-fat). Reads were mapped back to the assembly using Mosaik v2.2.28²⁴, variants were called by V-Phaser v2.0²⁵ and intra-host diversity was explored with V-Profiler v1.0²⁶.

Virus RNA-seq detection of RNA viruses in plasma

We first evaluated the performance of a conventional, “metagenomic” approach to virus whole-genome sequencing¹³. Indexed sequencing libraries were constructed in duplicate from plasma RNA of 29 subjects infected with diverse HCV subtypes (1a, 1b, 2a, 2b, 3a, 4a and 4d) and a 3.5-log range of VLs (2,200–4.9 million IU/mL; 1 IU = 2.7 copies on the instrument we use) and sequenced on a single Illumina HiSeq 2500 Rapid run, producing a median of 8.0 million reads per sample (range 6.0–24.9 million), of which 0.37% originated from HCV (range 0.03%–2.8%) (Supplementary Table S1). There was a linear relationship between HCV VL and the yield of HCV reads with high mapping quality (Figure 1). Mapping the HCV reads for each sample to the closest available reference (either a database reference or a de novo assembly of the same reads) produced patterns of peaks and troughs in sequence coverage along the genome that showed some similarity between samples of different subtypes and were highly reproducible between library and sequencing technical replicates; we therefore infer patterns of coverage are caused mainly by genomic features such as secondary structure and melting temperature²⁷.

Figure 1. HCV metagenomic sequence yield is proportional to viral load.

The yield of reads that map to any HCV genome and the probability of successful de novo assembly of a complete genome sequence both depend on viral load (VL). Samples were prepared as replicate libraries that were sequenced simultaneously with consistent yield. Blue circles: successful de novo assembly (>90% complete genome length recovered); red circles: incomplete genome assembly. a. With standard mapping criteria, up to 2.8% of reads match HCV and a background 0.02–0.1% of low-complexity human-derived sequences overwhelms the HCV signal in low-VL samples. Linear trend is plotted for samples with VL > 10⁵ IU/ml. b. Under stringent mapping criteria (mapping Q > 40), lower complexity human and HCV reads are excluded and yield is proportional to VL (slope of linear trend in log-log space not significantly different from 1) across the VL range.

In its standard form, metagenomic sequencing of a batch of up to 96 samples costs <£100 per sample. In this experiment, a VL of approximately 2×10⁵ IU/mL was sufficient to attain a mean read depth across the genome of ~30 and a high probability of successful de novo assembly, but higher read depths are necessary for precise characterization of minor variants. Results are better with high-VL samples, and measures to increase library complexity and improve release of virus during RNA isolation may improve variant-calling sensitivity, but the low efficiency of metagenomic sequencing poses a fundamental problem.

ve-Seq: Probe-based enrichment increases HCV sequence yield

When the sequence of interest comprises only a small fraction of the starting material, probe-based sequence capture, as used in exome sequencing, can dramatically increase sequencing efficiency^17,28. Anticipating the challenge posed by the extreme diversity of HCV, we drew on a representative genome sequence from each of four common genotypes (1a, 2b, 3a and 4a) to construct a combined panel of biotinylated DNA oligonucleotides (xGen^® Lockdown^® probes, IDT) comprising four sets of 155–157 probes, each a 120 nt sequence fragment overlapping the next by 60 nt, and excluding the 3’ poly-(U) tract to avoid enrichment of low-complexity non-HCV sequences.

We enriched the previously-sequenced pool of libraries for HCV sequences by solution hybridization with the 4-genotype probe panel and sequenced it on the Illumina MiSeq platform. This yielded a greater-than 16 × increase in the total number of HCV reads produced, even with an output of ~14 × fewer reads than the previous metagenomic sequencing on the higher-output HiSeq (Supplementary Table S1). HCV sequence content reached 86% in the enriched pool (range 1–98% among samples), equivalent to a median 1,660 (range 10–75,700) genomic average read depth or >10³-fold enrichment for samples with mid-range VL (Supplementary Figure S1); and hit saturation point (near-100% HCV reads) for samples with higher starting HCV content. Although probe panels can be expensive to synthesize, they can be used for many (hundreds of) pooled captures, so the lower sequencing costs in ve-SEQ more than account for the extra costs of the enrichment step.

Probe-target dissimilarity reduces enrichment efficiency

We used a single-genome, subtype 1a subset of the 4-genotype probe panel to investigate the effect of varying probe-target sequence identity on ve-SEQ enrichment success (Figure 2). When a sample is enriched with probes derived from that sample’s consensus sequence, there is no detectable bias in read depth with genomic position (i.e. coverage across the genome for enriched data follows a pattern almost identical to unenriched data, albeit at much higher read depth). When a non-identical sample of the same subtype is enriched, coverage patterns coincide, but are not identical. When a sample from the same genotype but a different subtype to the probe panel is enriched, large sections of the genome are adequately sequenced, but the most divergent regions are covered poorly and whole-genome assembly fails for samples with low viral load (Supplementary Table S1). When the sample and the enrichment probe set come from different genotypes, only the most conserved parts of the genome are adequately represented with ve-SEQ data and read depth is essentially zero for divergent regions.

Figure 2. Enrichment efficiency decreases with phylogenetic distance.

Read depth across the genome before (blue, left axis) and after (red, right axis) enrichment with a single-sequence subtype 1a probe set. a. The HCV genome comprises 5’ and 3’ untranslated regions (UTRs) and a large central segment encoding a single polyprotein that is cleaved into ten proteins. b. A subtype 1a sample enriched with probes derived from its own consensus sequence yields coverage patterns across the genome essentially identical to metagenomic sequencing. c. A distinct subtype 1a sample produces highly similar but non-identical patterns of pre- and post-enrichment genomic coverage. d. A subtype 1b sample yields low read depths at loci that are relatively divergent from the 1a probe sequence (E1, E2, NS2 and NS5a). e. Sequence capture of a sample from a different genotype, 3a, is poor across large segments of the genome. f. Heat map representing average diversity (calculated as Shannon entropy) among 165 HCV reference genomes. Nucleotide diversity varies dramatically across the genome and tracks drops in enrichment efficiency between phylogenetically distinct probe-target combinations.

In order to rationalize our approach to probe choice and enable the design of an efficient, comprehensive HCV enrichment probe set, we analysed the relationship between probe-target similarity and the relative efficiency of enrichment (Figure 3). Noting a strong inflection point, we deduced that a minimum 80% identity between a 120 nt segment of sample sequence and its closest matching probe was sufficient to ensure near-maximal enrichment, assuming that each sequencing library molecule interacted with a single probe molecule and ignoring the potential effects of bridging capture (i.e. successful sequencing of a poorly matching fragment effected by hybridization of an adjacent target sequence on the same library molecule to a better-matching probe). The 20% divergence cutoff for successful enrichment falls between the mean inter-subtype (<15%) and inter-genotype (30–35%) divergence levels, explaining why enrichment with a subtype-mismatched probe set leads to only localized bias, while genotype-mismatch results in failure across most of the genome. It also follows from this analysis that when enrichment is performing well, there should be no detectable bias in the representation of single nucleotide variant alleles such as RAVs.

Figure 3. Enrichment efficiency is directly related to probe-target identity.

A set of 10 HCV samples with highest VL was sequenced before and after enrichment with a single-genome, subtype 1a probe set, and for each sample the relative read depth for each probe window was plotted against the maximum identity between target and any probe. Read depth ratio was normalized by giving the most efficiently enriched probe position (in the highly conserved 5’ UTR) a value of 1. Maximal enrichment is observed where probe-target identity exceeds approximately 80% and enrichment decreases dramatically as identity falls below 80%.

Design of a comprehensive probe set for HCV

As is evident from the previous section, a probe panel based on just four subtype-representative sequences cannot perfectly capture HCV global diversity. Exploiting the observation that some regions of the HCV genome (e.g. the 5’UTR) are well-enough conserved to not require multiple probe sets, together with the 20% divergence cutoff for efficient capture, we implemented an algorithm for efficient probe set design that would facilitate a comprehensive HCV enrichment panel as well as, in principle, efficient probe sets for other organisms.

We started with the 4-genotype probe panel and added extra probes to improve coverage for already-included subtypes 1a, 2b, 3a and 4a as well as the extra subtypes 1b, 2a, 2c, 5a and 6a, using a database of 482 reference whole-genome sequences. First we calculated a consensus sequence for each subtype. Then, starting with the existing probe set and the first genome in the most common subtype (1b), we identified genomic regions with less than 80% identity to any of the probes already in the panel. For each such region the subtype consensus sequence was considered as a potential probe but only used if it was ≥80% identical to the genomic sequence it replaced; otherwise the genomic sequence fragment was added as a new probe. The process was repeated for each 1b reference sequence and then similarly for each subtype.

In contrast to the naïve design of probe sets with the standard IDT approach that requires 155–157 probes per HCV target genome, we were able to augment our 4-genome probe panel to represent the known diversity of nine subtypes spanning six of the seven recognized genotypes with only another 491 probes (1,116 total). Our algorithm substantially and automatically reduces redundancy: a completely naïve approach that simply encoded every genome in the reference set, without accounting for similarity between genomes, would have dictated a prohibitively expensive set of ~75,000 probes. In contrast, if we had instead started from scratch, we estimate that our simple algorithm could have produced an equally effective combined panel for nine subtypes with as few as 955 probes. In informal testing, a typical sample from the newly added subtype 1b achieved near-zero bias even though its exact sequence was not encoded in the probe set but was instead covered by reference to recorded sequence diversity (Supplementary Figure S2), and a sample from subtype 4d, not included in the revised probe set, achieved adequate although imperfect enrichment (Supplementary Figure S3), consistent with previous subtype-mismatched captures. Although feasible and relatively inexpensive, we have deferred the addition of probes for remaining rare subtypes.

Detection of resistance-associated variants in clinical samples

To explore the potential utility of high-depth RAV data in predicting the clinical effectiveness of HCV treatment, we used ve-SEQ to analyse retrospectively plasma samples collected from 33 genotype 1-infected patients before NS3-targeting DAA therapy with Boceprevir (14 patients) or Telaprevir (19 patients) (Supplementary Table S2). We obtained whole-genome sequences for all samples, with a mean read depth of 4600 across the NS3 gene. We first confirmed that our sequence data (28 subtype 1a and 5 subtype 1b) matched clinical subtyping data where the latter was available.

Mutations in the NS3 gene, denoted T54S and V55I, were detected in patient P23, in whom Boceprevir treatment failed to suppress HCV. Only one other patient had relevant baseline resistance: P6 possessed a single T54S mutation, yet cleared infection with 48 weeks of BCP. Additionally, Simeprevir RAVs Q80K/R were detected in five patients with genotype 1a virus, consistent with the reported prevalence of these mutations in PI-naïve patients²⁹. Variants associated with NS5A inhibitor resistance were detected in 11 patients, including nine with combinations of two or more RAVs, previously associated with higher relapse rates than Lidipesvir/Sofosbuvir³⁰.

In samples taken after treatment cessation, five patients carried both V36M and R155K NS3 variants, associated with drug resistance but also reduced virus fitness in the absence of treatment^31,32, including three patients illustrated in Figure 4. RAVs V36M and R155K were each detected independently of the other (in P30 and P33, respectively) and virus sampled in P27 during treatment revealed approximately 2-fold more V36M variants than R155K, confirming that V36M alone was sufficient to confer resistance on individual genomes. Telaprevir had failed to suppress virus in subject P24 by week 4 when V36M and R155K variants circulated in approximately half of virus. It is therefore not surprising that a subsequent treatment attempt also failed, providing a real-world clinical example of where sequencing might have prevented futile retreatment. Six weeks after the second treatment attempt had failed, the R155K mutation had reverted to the wild-type arginine residue in all sequence reads. Partial reversion was also observed in P18, although in this instance, reversion of V36M occurred some 20 or more weeks after the cessation of treatment and R155K was still present in 100% of variants 1 year later.

Figure 4. Detection of resistance-associated variants after DAA treatment failure.

VL and RAV status for three patients who failed to achieve sustained virological response after Telaprevir-based therapy. Grey shading: duration of therapy (weeks starting at time 0); squares: VL measurements; inverted triangles: samples sequenced using the comprehensive probe panel (open: no Telaprevir RAVs detected, black: RAVs and supporting read proportions, where <100%).