A Follow-Up of the Multicenter Collaborative Study on HIV-1 Drug Resistance and Tropism Testing Using 454 Ultra Deep Pyrosequencing.

BACKGROUND: Ultra deep sequencing is of increasing use not only in research but also in diagnostics. For implementation of ultra deep sequencing assays in clinical laboratories for routine diagnostics, intra- and inter-laboratory testing are of the utmost importance.
METHODS: A multicenter study was conducted to validate an updated assay design for 454 Life Sciences' GS FLX Titanium system targeting protease/reverse transcriptase (RTP) and env (V3) regions to identify HIV-1 drug-resistance mutations and determine co-receptor use with high sensitivity. The study included 30 HIV-1 subtype B and 6 subtype non-B samples with viral titers (VT) of 3,940-447,400 copies/mL, two dilution series (52,129-1,340 and 25,130-734 copies/mL), and triplicate samples. Amplicons spanning PR codons 10-99, RT codons 1-251 and the entire V3 region were generated using barcoded primers. Analysis was performed using the GS Amplicon Variant Analyzer and geno2pheno for tropism. For comparison, population sequencing was performed using the ViroSeq HIV-1 genotyping system.
RESULTS: The median sequencing depth across the 11 sites was 1,829 reads per position for RTP (IQR 592-3,488) and 2,410 for V3 (IQR 786-3,695). 10 preselected drug resistant variants were measured across sites and showed high inter-laboratory correlation across all sites with data (P<0.001). The triplicate samples of a plasmid mixture confirmed the high inter-laboratory consistency (mean% ± stdev: 4.6 ±0.5, 4.8 ±0.4, 4.9 ±0.3) and revealed good intra-laboratory consistency (mean% range ± stdev range: 4.2-5.2 ± 0.04-0.65). In the two dilutions series, no variants >20% were missed, variants 2-10% were detected at most sites (even at low VT), and variants 1-2% were detected by some sites. All mutations detected by population sequencing were also detected by UDS.
CONCLUSIONS: This assay design results in an accurate and reproducible approach to analyze HIV-1 mutant spectra, even at variant frequencies well below those routinely detectable by population sequencing.

Introduction

HIV-1 drug resistance testing is currently recommended and widely implemented as the standard of care prior to starting antiretroviral therapy (ART) [1]. The standard method used worldwide for drug resistance testing is based on population sequencing, which does not allow for the detection or reporting of minority HIV-1 drug-resistance mutations [2]. Several methods are available for detection of minority mutations [3]. Tests based on allele specific PCR (AS-PCR) can be used for high sensitivity detection of minority variants; however, this technology can only detect one mutation per assay and becomes very labor intensive for analysis of multiple resistant variants [4, 5]. Next Generation Sequencing assays, having the capability to analyze multiple drug resistance mutations in a single assay, are being used with increasing frequency [6]. Ultra deep sequencing (UDS) has been retrospectively assessed in clinical trials such as the MOTIVATE and A4001029 studies where viral tropism was determined phenotypically prior to maraviroc administration [7]. The presence of ≥2% non-R5 HIV by UDS was found to be a good predictor of virologic response to maraviroc [7].

UDS has been indicated in previous studies to be an important tool for drug resistance detection due to its ability to detect and quantify minority drug resistance mutations. A systematic review and pooled analysis showed the risk of virological failure to be significantly increased in ART-naïve patients receiving an NNRTI-based first line regimen who harbor pre-existing minority NNRTI-resistant HIV-1 variants [8] recently confirmed by a multi-cohort European case-control study [9]. The latter and all of the studies analyzed in this pooled analysis; however, were retrospective in nature. Prospective studies have not yet been completed to unequivocally show the clinical relevance of minority HIV-1 drug-resistance mutations. Furthermore, clinical cut-offs for the significance of such minority HIV-1 drug-resistance mutations have yet to be defined. More sensitive assays and studies are required to answer these questions. Here, we describe an assay with the potential to help answer these questions.

We previously described a multicenter study demonstrating high inter-laboratory correlation of UDS-detected minority HIV-1 drug-resistance mutations exclusively in HIV-1 subtype B samples. This initial study was completed using the GS FLX sequencing platform [10] and focused on the detection of minority HIV-1 drug-resistance mutations in the RTP region of HIV-1 subtype B viruses. The follow-up study described here employed an updated primer design consisting of 8 amplicons between 400 and 600 bp (including sequencing adaptors and multiplex identifier (MID); also commonly referred to as a barcode) that target multiple clades of HIV-1 samples. The new assay design takes advantage of longer read lengths and also includes the V3 loop of the env gene to provide viral co-receptor use information (Fig A in S1 File).

Validation based on comparison to state of the art genotyping is a pre-requisite for adoption of ultra deep sequencing in routine care settings. We investigated the inter- and intra-laboratory performance of a prototype HIV-1 primer set and sequencing protocol for ultra-deep pyrosequencing. In addition, we sought to demonstrate concordance with state of the art genotyping as well as the detection of minority HIV-1 drug-resistance mutations in samples from multiple HIV-1 clades.

Materials and Methods

Ethics statement

Human plasma used for this study was derived from anonymous blood donations obtained from healthy blood donors provided by the Blood Donation Service. Written informed consent or an ethics statement for the use of plasma from anonymous blood donations are not required for research purposes.

Study design

HIV-1 containing samples were provided to 10 study sites in Europe and one study site in South Africa. Site numbering began with the number 20 to disambiguate the new data files from those generated in the previous study. The samples were comprised of 30 HIV-1 subtype B and 6 subtype non-B samples (A2, AE, A/G, B/G, C, and G) with viral titers of 3,940–447,400 copies/mL, two dilution series (52,129–1,340 and 25,130–734 copies/mL, respectively), triplicate samples and two plasmid controls (Table A in S1 File). All samples were sent blinded to the sites. Also provided to each site were detailed laboratory protocols, four 96-well thermocycler plates with ready to use cDNA primers as well a set of four 96-well thermocycler plates that contained pre-loaded barcoded PCR primers for the generation of eight overlapping amplicons. For traceability through the workflow, the barcode on the PCR primers in the first column of each plate was unique to that plate. For the concordance analysis, population sequencing of all samples was performed with the ViroSeq HIV-1 genotyping system by Dr. Däumer’s clinical laboratory.

Primer design

The overlapping amplicons span PR codons 10–99, RT codons 1–251 and the entire V3 region. The amplicons were designed to provide dual coverage for Stanford Database (version 6.0.9) Drug Resistance PR and RT mutations with a score ≥5 (any drug) and tropism determinants in V3 (Fig A and Table A in S1 File). All primers were designed using the principles previously described, with the exception of allowing for 98% conservations instead of 99% as described in the referenced paper [11]. The longer read length chemistry employed in this study reduced the required number of overlapping RTP amplicons from six to four and also allowed for placement of primers in sufficiently conserved areas outside the V3 region.

Sample characteristics

The RTP samples used in this study were cultured recombinant virus diluted in plasma from HIV-1 negative individuals. The cloned inserts included the coding regions for the complete protease and the 5’ 302 amino acids (aa) of the reverse transcriptase (nucleotide positions 2,002 to 3,457 based on the reference sequence HIV-1HXB2, GenBank accession number K03455). The V3 samples were cultured samples from patient isolates, diluted in HIV-1 negative plasma. Two control samples were generated for inter- and intra-laboratory comparison purposes (control sample A) and for amplification control (control sample B). Control sample A is a mixture of two plasmid samples obtained through the AIDS Research and Reference Reagent Program, Division of AIDS, NIAID, NIH: pM46I/L63P/V82T/I84V and pL10R/M46I/L63P/V82T/I84V from Dr. Emilio Emini [12]. These two plasmids were mixed at a 20:1 ratio to create a sample with the L10R mutation at approximately 5%. Control sample B is a preparation of plasmid pM46I/L63P/V82T/I84V.

HIV-1 RNA extraction, cDNA synthesis, amplification and pooling

1 ml of each sample was ultracentrifuged for 1.5 hrs at 4°C. 800 μl of the supernatant was removed and discarded, leaving ~200 μL of viral pellet and plasma. HIV-1 RNA was extracted using the QIAamp MinElute Virus Spin kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. MS2 carrier RNA (Roche, Mannheim, Germany) was then added to each eluate to a final concentration of 10 ng/μL. 13 μL of the RNA eluate was reverse transcribed using Transcriptor Reverse Transcriptase (Roche) with each of the “4R”, “5R” and/or “V3R” primers (provided in ready to use 96-well plates), resulting in three first strand cDNA-fragments covering the protease, reverse transcriptase and env V3 coding regions of HIV-1. The cDNA was treated with RNAseH (Roche) and 3 μL cDNA product used for each PCR reaction.

Amplicons were generated from the cDNA and the two plasmid control samples (A and B) by 40 cycles of PCR with FastStart HiFi Polymerase (Roche) as depicted in Fig A in S1 File using PCR conditions previously described [10] in a total of four 96-well thermocycler plates. Primer sequences are given in Table B in S1 File. Amplicons were then purified with AMPure XP beads (Agencourt, Beckman Coulter Inc.) and quantified using the Quant-iT PicoGreen dsDNA kit (Invitrogen, Darmstadt, Germany). Amplicons with concentrations <5 ng/μL were further analyzed on a 2100 Bioanalyzer (Agilent Technologies, Böblingen, Germany). PCR products with a molar ratio of primer-dimer to amplicon above 1:3 were excluded from sample pooling. The remaining amplicons were pooled equimolarly by sample MID (6 RTP and 2 V3 amplicons for each MID), with any missing amplicon being compensated for by increasing the amount of an overlapping amplicon, where possible. The samples were then pooled by plate to make 4 final pools, 1 pool per plate corresponding to 1 pool per region of a 454 sequencing pico titer plate (PTP) fitted with a 4-lane gasket. In total, 36 samples (18 for sequencing of the RTP region and 18 for the env V3 region) were amplified in PCR plates 1 and 2, corresponding to regions 1 and 2 of the sequencing run. 3 additional samples were amplified in triplicate on plate 3 (9 samples total, sequencing both the RTP- and env V3 regions), corresponding to region 3 of the sequencing run. Amplicons for the two sample dilution series were generated on plate 4 (10 samples total, sequencing both the RTP- and env V3 regions), corresponding to region 4. Control sample A was amplified in each of the first 3 plates and included in the pools for sequencing. Control sample B was used as an amplification control on each of the four plates to ensure proper execution of the PCR step as was a negative water control to ensure no contamination of the PCR; neither of these controls was included in the pools for sequencing.

Ultra deep sequencing

Each of the 4 amplicon pools was separately subjected to clonal emulsion PCR amplification on beads using reagents that enabled sequencing in both the forward and reverse directions (Lib-A kit, 454 Life Sciences—A Roche Company, Branford). DNA-containing beads were counted on a Multisizer3 instrument (Beckman Coulter Inc.) and prepared for sequencing. Approximately 790,000 beads from each of the four pools were loaded in the four regions of a PTP fitted with a 4-lane gasket. Sequencing was performed on a GS FLX Titanium System (454 Life Sciences -A Roche Company, Branford). Resulting sequence reads were aligned to HIV-1HXB2 by the GS Amplicon Variant Analyzer (AVA) v2.3 software included with the 454 Life Sciences instrument. The AVA software demultiplexed barcoded reads and aligned them to the HIV-1HXB2 reference while simultaneously eliminating primer-derived sequences. The software then reported the prevalence of both sequence conservation and mutations relative to HIV-1HXB2. Since AVA detects variants in nucleotide space, configuration files were provided to each site that defined amino acid mutations based on triplets of nucleotides. Sff files are deposited in the European Nucleotide Archive (accession number PRJEB12165)

Data analyses

Data was analyzed for intra- and inter-laboratory concordance using the R statistical package [13] or GraphPad Prism software (GraphPad Software, San Diego, CA). The correlation of variant prevalence for sample replicates within sites and common samples across sites was also computed. In addition, the minimum, maximum and median variant levels, as reported by each site, were calculated to assess concordance with population sequencing results for all samples. For minority variant calling, 1% was chosen as cut off based on previous reports showing combined PCR and pyrosequencing error rates of 0.1–0.3% using HIV-derived clonal sequences [14-16].

Tropism analysis was performed with the geno2pheno co-receptor algorithm for each V3 sample at a false positive rate setting of 3.75% using AVA preprocessed files. Finally, standard population sequencing was performed on each sample using the ViroSeq HIV-1 genotyping system (Abbott Laboratories, Abbott Park, Illinois).

Results

Overall performance

Ten out of eleven sites delivered a complete data set. Sequencing information was obtained for 2,986 of 3,520 (85%) expected amplicons. One site was unable to deliver information for the second of their four plates (Table 1) and other sites encountered sporadic amplification failures. Table C in S1 File describes the details for each amplicon drop-out. The higher rate of drop-outs in plate 4 are due to the lower viral titers associated with the dilution series samples. Fig B I in S1 File shows the uniform distribution of the total number of high quality (HQ) reads (those reads that passed the amplicon data processing pipeline) per sequencing region for each site. Furthermore, almost all of the HQ reads could be aligned to HIV-1HXB2 with the exception of those from site 27 (region 2) and sites 26 and 28 which delivered a reduced total number of HQ reads.

Table 1

Detection across sites of 10 preselected variants.

The table lists the percentage of reads carrying the codon variant.

Drug resistance mutation	Sample	HIV-1 subtype	Study sites
			20	21	23	24	25	26	27	28	29	30	31	Median	IQR
NNRTI-G190A	1	G	28.5	23.4	26.5	27.7	27.2	35.3	32.6	18.2	33.0	27.8	27.8	27.8	3.7
NNRTI-F227L	3	B	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
NRTI-K219R	4	B	62.6	58.3	65.6	60.6	58.6	58.2	62.1	67.0	57.2	62.8	65.0	62.1	5.5
NRTI-V75A	5	A2	2.8	1.6	1.1	1.1	0.7	1.7	0.7	2.5	0.5	1.2	1.3	1.2	0.8
PI-p-V82I	19	B	1.7	1.4	n.a.a	1.8	1.1	0.6	—b	1.4	0.9	1.4	1.3	1.4	0.3
PI-I54V	20	B	54.8	99.3	98.6	98.7	98.4	99.7	0.0	95.2	51.0	49.6	99.2	98.4	46.1
NRTI-L74I	21	B	8.7	12.2	15.7	11.4	9.1	11.0	—	10.4	9.7	11.1	11.0	11.0	1.5
NNRTI-V179D	24	AG	0.0	0.0	0.0	0.0	0.0	0.0	—	0.0	0.0	0.0	0.0	0.0	0.0
PI-G73S	26	B	2.8	3.7	4.5	4.8	4.7	4.9	—	4.1	3.4	1.5	3.6	3.9	1.2
NRTI-V118I	27	C	30.0	35.4	32.7	34.8	36.6	56.3	—	39.0	15.1	28.2	27.8	33.8	7.7

a not applicable, i.e., fewer than 15 reads in each direction

b no data returned.

Analysis of the distribution of the 8 amplicons revealed that there was a reduced number of HQ reads (across sites) for amplicons RTP4 and RTP5 for all sites when compared to the number of reads for amplicons RTP1-3, 6, V3A and V3B. Fig B II in S1 File depicts the HQ read distribution data per site for each amplicon using region 3 of the 4-region PTP as a representative example. Across all plates and sites, the median sequencing depth was 1,829 reads per position for RTP (IQR = 592–3,488) and 2,410 reads per position for V3 (IQR = 786–3,695).

Inter-laboratory concordance

To determine inter-laboratory variability of the current workflow for end users, each site was given a questionnaire to identify and report variant percentages using the AVA software. Ten preselected drug resistance mutations were chosen for this analysis, two were expected negatives (V179D, F227L), three were <5% prevalence (V75A, V82I, G73S) and the remaining five were >10% prevalence (Table 1). All readings reported in this table required a minimum of 15 reads in addition to successful amplicon generation as indicated by the amplification control on the relevant plate. A read coverage of <15 was defined as sequencing failure. Of those amplicons reported with sufficient read coverage, all but one was accurately detected and reported. The single highly inconstant variant was the I54V protease mutation on plate #2 of site 27. The site reported technical errors with this particular plate. Overall, the 10 preselected drug resistant variants, including those with non-B subtypes, showed high inter-laboratory correlation across all sites (Table 1).

Detection of a minority variant in a plasmid mix

To determine the consistency of variant detection, not only between laboratories but also across plate amplification, an engineered plasmid mixture was sequenced in triplicate at each site (Table 2). The replicate plasmid mixtures were distributed over 3 of the 4 PCR plates and contained a single low level mutation in the protease region, L10R. The mean reported frequencies for this variant across all 11 sites for plates 1, 2 and 3 were 4.7, 4.8 and 4.8% with standard deviations of 0.5, 0.6 and 0.3%, respectively. The replicate plates within each site had means ranging from 4.3 to 5.2% with standard deviations of 0.3 and 0.4%, respectively (Table 2).

Table 2

Consistency of variant detection for an engineered plasmid mix.

A plasmid mix containing a variant present at approximately 5% was sequenced in triplicate at each site. The frequency of the minority variant is given in percent. The replicates were distributed over 3 of the 4 PCR plates.

	Plates
Site	1	2	2	Median	IQR
20	4.4	4.0	4.6	4.4	0.3
21	5.2	4.7	5.2	5.2	0.3
23	5.4	5.1	4.5	5.1	0.5
24	4.4	4.3	4.9	4.4	0.3
25	4.3	5.1	5.3	5.1	0.5
26	4.2	5.0	4.6	4.6	0.4
27	4.9	5.9	4.7	4.9	0.6
28	4.5	3.7	4.3	4.3	0.4
29	4.1	5.2	4.8	4.8	0.6
30	5.6	4.7	5.4	5.4	0.4
31	5.0	5.1	5.0	5.0	0.0
Median	4.5	5.0	4.8
IQR	0.8	0.6	0.5

Concordance with population sequencing

The ViroSeq assay was performed on the sample set by a single laboratory to determine concordance with UDS. 99%-100% of drug resistance mutations detected by the population sequencing-based assay were also found by UDS by all sites (Table 3). Three categories of mutations were analyzed: i) major, defined as having a Stanford drug resistance database score >5, ii) polymorphic, defined as having a Stanford drug resistance database score <5 and iii) bulk, defined as any other difference from HIV-1HXB2 at the amino acid level. The single discrepancy between UDS and population sequencing for the polymorphic mutation may be due to low sequencing coverage.

Table 3

Concordance between amino acid mutations detected by population and ultra-deep pyrosequencing.

Included in the analysis were all non-redundant RTP samples (1–9, 19–27, 37 replicate #2, 38 replicate #2, and 39 replicate #2).

Mutation type	Population sequencing	Also detected by 454_UDS (%)
Major mutations (Stanford score >5)	72	72 (100)
Polymorphisms (Stanford score <5)	113	112 (99)
Bulk, differences from HIV-1HXB2	391	391 (100)

Of clear clinical interest, we found that several population sequencing-called mutations were incorrect due to the inability of population sequencing to phase adjacent nucleotides. Fig C in S1 File is a representative screenshot from the AVA software showing the composition of the UDS sequence reads at codon 215 in the reverse transcriptase. Population sequencing called several possible mutations that could be present at this location although the mutations could not be resolved due to the lack of phase information. This phenomenon occurs since population sequencing only provides information about the composition at each nucleotide location but not about how the nucleotides are connected in each variant. As shown in Fig C in S1 File, UDS resolved the actual variant codon compositions, making it clear that only TAT (Tyr) and ACT (Thr) were present. Population sequencing is not able to disambiguate between the 4 possible combinations TAT (Tyr), AAT (Asn), TCT (Ser), and ACT (Thr) and thus the associated software called all possible known resistance-linked variants at this position. These obvious incorrect combinations in the population sequencing data were excluded prior to the concordance analysis shown in Table 3.

Additional mutations found by UDS but not population sequencing are listed in Table 4. Most of these were minority drug resistance mutations, i.e., found in less than 20% in a given sample. However, several of these mutations were detected at levels >20% by UDS, likely reflecting the inherent variability in population sequencing sensitivity. Most of these mutations were detected and reported at similar frequencies by all sites with sufficient coverage for the respective amplicons (Table 4).

Table 4

Drug resistance mutations found only by UDS.

For each detected drug resistance mutation, the table indicates the sample, the number of sites detecting the mutation at a frequency above 0.5%, the number of sites with coverage of at least 15 reads and the median, mean, minimum and maximum variant frequencies. Details on reports of each site are given in Table D in S1 File.

Drug resistance mutation	Sample	Sites ≥0.5%a	Sites with coverageb	Median (%)	Mean (%)	Minimum (%)	Maximum (%)
NNRTI-F227L	1	11	11	1.64	1.82	0.83	5.99
NNRTI-K103N	2	11	11	10.71	9.98	5.07	15.56
NNRTI-K103R	2	9	11	1.23	1.71	0	6.19
PI-V82A	3	10	11	4.87	4.55	0	6.37
NRTI-Y115F	4	11	11	4.20	3.83	1.14	5.52
NRTI-K70E	5	11	11	6.88	6.85	5.21	9.62
NRTI-V75A	5	11	11	1.19	1.39	0.55	2.84
NNRTI-G190A	7	11	11	3.44	3.33	1.52	5.88
NNRTI-Y188C	7	11	11	19.09	19.47	16.29	21.86
NRTI-M184V	7	11	11	7.54	7.58	3.22	10.27
PI-V82A	7	10	11	2.09	2.06	0.31	3.18
NNRTI-E138G	8	5	10	0.66	1.18	0	4.77
PI-M46V	8	11	11	2.70	2.79	1.65	4.25
NNRTI-V179D	9	11	11	1.43	1.71	1.03	2.74
NRTI-T215Y	19	10	10	8.32	8.03	5.55	10.38
NRTI-T69A	19	10	10	4.65	5.01	2.97	9.83
NRTI-T69N	19	7	10	1.02	1.71	0.16	5.08
NRTI-M184I	20	10	10	6.21	6.91	5.12	10.15
NRTI-L74I	21	10	10	10.99	11.03	8.77	15.74
NNRTI-K103R	23	9	9	19.79	20.82	17.01	26.84
PI-I47V	23	9	9	4.30	4.96	2.44	9.87
NNRTI-E138G	24	9	10	6.69	7.50	0	17.54
NNRTI-K103R	25	9	10	17.97	16.74	0	26.59
NRTI-K219E	25	9	10	6.38	5.93	0	9.52
NRTI-T69N	25	9	10	22.41	22.38	0	42.57
NNRTI-F227L	26	10	10	12.18	12.14	9.31	14.17
PI-G73S	26	10	10	3.89	3.81	1.51	4.93

a 0.5% was chosen to cover the expected errors on measurements of a 1% variant

b every site was included reporting at least 15 reads at the certain position. Fewer than 15 reads were defined as unsuccessful amplification.

Reproducibility and repeatability

To determine inter- and intra-laboratory reproducibility, a set of triplicate samples was distributed to each of the sites. The triplicate samples 37, 38 and 39 had previously been sequenced and analyzed in the initial multicenter collaborative study [10]. These samples had viral titers of 230,700, 447,400 and 80,000 copies/mL, respectively. Variant percentages measured in the triplicate samples are shown for each mutation in Fig 1. Triplicate readings for mutations detected in these samples are closely clustered within each site, indicating high intra-site reproducibility (repeatability). The highest degree of variability recorded for several variants was seen in site 28. This increased variation between triplicates may have been related to the very low number of HQ aligned reads in region 3 for this site (Fig B I in S1 File) which was less than half of all other sites. In addition, this data set also indicates high inter-site reproducibility as all mutations were detected across sites with an average standard deviation of 1.0.

Fig 1

Triplicate sample analysis.

A set of three samples was measured in triplicate at each of the 11 sites. Variant percentages measured in the triplicate samples are shown for each mutation.

Variant detection in dilution series

To examine the limitations of sequencing samples with reduced viral titers, two samples with starting viral titers of 52,119 and 25,130 copies/mL were serially diluted in human plasma four times each to reach a lower viral titer of 1,340 and 734 copies/ml, respectively (Table A in S1 File). To ensure accuracy of the viral titers in the dilution series, the titer was re-determined on aliquots from each dilution. Subsequently, 1 ml of each dilution was extracted, had amplicons generated and was sequenced alongside the other HIV-1 samples. Fig 2 displays the amplicon sequencing coverage for the second dilution series sample. Drug resistance mutations in frequencies of >20% were detected by 90–100% of the sites, independent of the viral titer. Lower viral titers led to a decrease in detection of minority drug resistance mutations, especially when present at 1–2%. Of note, even these very low frequency mutations were still detected (on average) at approximately one quarter of the sites at a viral titer of 1,340 copies/mL (sample 44; Fig 2A).

Fig 2

Drug resistance mutation detection in limiting dilution series.

Two HIV-1 subtype B samples were used for the dilution series (A, samples 40–44; B, samples 45–49). The four groups shown in the histogram are based on the median frequency of each drug resistance mutation in its respective undiluted sample: >20% (white bars), 10–20% (bright grey bars), 2–10% (dark grey bars), and 1–2% (black bars). The means and standard deviations are given of the percentage of sites reporting drug-resistance mutations in these categories. The number of mutations in each category is represented by n.

Co-receptor usage prediction

Co-receptor usage calls for all V3 samples are shown in Table 5. The geno2pheno[454] algorithm (Max-Planck-Institute for Informatics, Saarbrücken, Germany) was used to analyze the UDS reads with a false positive ratio (FPR) setting of 3.75%. 14 of the samples sequenced by UDS for the V3 region were also analyzed by population sequencing. For those 14 samples, 13 (93%) co-receptor usage calls were concordant between both systems. The single discrepancy between the population sequencing and UDS based co-receptor usage predictions (sample 16) was due to a discordant result from a single site. This site reported 86.4% for V3A and 75.7% for V3B whereas all other sites reported a detection level of 0.0% X4 for both amplicons V3A and V3B and may be due to a sample swap.

Table 5

Co-receptor usage predictions.

Calculations for UDS were performed using geno2pheno with a false positive rate setting of 3.75% and reported as means with standard deviations for each sample across all sites.

V3 Sample	X4 prediction based on UDS (%)				Prediction based on population sequencing [FPR]a
	Amplicon V3-A		Amplicon V3-B
	Mean	StdDev	Mean	StdDev
10	0.1	0.2	0.6	1.8	R5 [42.3]
11	0.0	0.0	0.0	0.0	R5 [12.5]
12b	n.a.c	n.a.	n.a.	n.a.
13	0.0	0.0	0.1	0.4
14	99.8	0.4	99.7	0.3
15	0.0	0.0	0.0	0.0	R5 [93.5]
16d	6.9	22.8	8.6	27.3	R5 [35.3]
17	98.6	4.8	99.2	2.8
18	97.6	2.0	94.8	3.8	X4 [0.2]
28	82.8	12.3	84.8	12.5	X4 [0.2]
29	86.4	7.5	92.6	4.4	X4 [2.6]
30	99.7	0.3	99.2	2.1	X4 [0.1]
31	0.2	0.5	0.2	0.6	R5 [35.3]
32	0.0	0.1	0.1	0.2	R5 [76.0]
33	0.2	0.3	0.3	0.4	R5 [13.8]
34	0.1	0.3	0.1	0.2	R5 [38.8]
35	0.0	0.0	1.7	5.0	R5 [46.0]
36	0.0	0.0	0.0	0.1	R5 [38.5]
37–2	99.6	1.1	99.9	0.1
38–2	99.6	1.1	100.0	0.1
39–2	100.0	0.0	100.0	0.0

a Prediction of HIV-1 tropism by the geno2pheno [co-receptor] algorithm was done using a FPR of 10% according to the European guidelines on the clinical management of HIV-1 tropism testing [34]

b geno2pheno algorithm3 only accepted <1% of the reads for sample 12

c not applicable

d The high stdev for sample 16 resulted from a discrepancy at a single site. Without that site the %X4 would have been 0% for both V3A and V3B, agreeing with an R5 population sequencing- based prediction.

Discussion

The data presented here demonstrate that ultra deep sequencing can provide highly consistent results together with excellent cross-site correlation. For all drug resistance mutations identified by the Stanford Drug Resistance Database in the sequenced region, we demonstrate near perfect concordance with population sequencing results. In addition, 41 mutations were found by UDS in the reverse transcriptase and protease regions at frequencies below 20%. These low frequency mutations could not be detected by population sequencing due to its inherently lower sensitivity [2]. Several studies have previously reported high concordance of UDS and standard population sequencing for mutations present at levels greater than 20–25% (reviewed in [6, 17]). One such study reported similar reproducibility of the UDS results with a correlation value of r2 = 0.969 [18].

Additionally, some variants called by population sequencing were incorrect due to the inability of population sequencing to phase individual nucleotides. In the example shown in Fig C in S1 File, population sequencing analysis indicated four possible mutations at codon 215, i.e., T215 (ACT), T215Y (TAT), T215S (TCT), and T215N (AAT). In this particular case, the Serine and Asparagine codons were listed as a combinatorial possibility by the ViroSeq data analysis software, but UDS shows that both are not present among the actual quasispecies. Since T215Y is detected with both methods, it makes no clinical difference in this case, but it is easy to see how in certain instances this could cause interpretation of resistance where none is present, leading to unnecessary avoidance of a potentially useful drug regimen. Looking at the technical performance of the assay, it was clear that certain amplicons performed more robustly than others, both in PCR and sequencing. In general, shorter amplicons can be expected to amplify with higher efficiency than longer ones [19]. This phenomenon negatively affected the number of HQ reads for amplicons RTP4 and RTP5, the two longest amplicons in the set, across all sites. It was also noted that the sequencing run at site 29 yielded a disproportionally large number of HQ reads. This difference was found to be due to a custom change in their software pipeline settings allowing more HQ reads to be reported; however, this can lead to reduced accuracy and is generally not recommended. The AVA software automatically discarded the majority of the additional low accuracy reads and the variant percentages reported from site 29 were similar to those reported by other sites. Conversely, sites 26 and 28 had average read counts but produced a significantly reduced number of aligned reads. This was likely due to the presence of increased amounts of primer dimers that were introduced, along with the amplicons of interest, into the emPCR step prior to the sequencing run. The presence of primer dimers in emPCR is deleterious to the amplification of target amplicons as these shorter fragments are preferentially amplified on the beads [20]. This phenomenon results in a decreased number of amplified target amplicons available for sequencing.

Assay performance at low viral titers was investigated with two dilution series. Data gathered across all sites demonstrated that, even at low titers, mutations present at >20% of the viral population were generally detected. Minority HIV-1 drug resistance mutations in the range of 2–10% were detected at the majority of sites, and minority HIV-1 drug resistance mutations in the range of 1–2% were still detected by some sites. It should be noted that the frequency categories in Fig 2 are based on the variant percentage measured in the undiluted sample which most closely reflects the true representation of the viral species present in the sample. When the viral titer decreases, stochastic sampling and decreased PCR efficiency will cause distortion of the ratio between the different viral species after amplification, and thus the frequencies of individual variants [21].

Sequencing of the V3 region yielded reads spanning the entire V3 loop, which means that each variant can be analyzed without making any inferences about phase. The geno2pheno co-receptor usage prediction software can then simply assign each read to the appropriate predicted tropism group (ie., X4 or R5) and make an overall call based on a relevant false positive rate setting [7, 22]. As noted in the introduction, this approach has been tested in retrospective clinical trials. Theoretically, the increased sensitivity as well as the full phase information should allow for more accurate X4 detection than population based sequencing. Such detailed information is also of great interest for studies of HIV variant dynamics [23, 24]. In this study, the only available comparison data was generated by population sequencing, and we demonstrate good concordance between the two technologies. The single discrepancy noted in Table 5 between the population sequencing and the 454 UDS X4/R5 prediction is for sample 16 where the high standard deviation was caused by an outlying result from a single site. If the data from this site is excluded, the %X4 reporting is 0% for both V3A and V3B amplicons in concordance with the population sequencing based prediction.

There is an increasing awareness of errors that are artificially introduced during RT-PCR and subsequent pyrosequencing and can be falsely interpreted as real mutations [6]. Errors introduced during the cDNA synthesis can currently not be distinguished from real mutations and occur sporadically and randomly [25]. They would need to be the subject of resampling during the following PCR to account for at least 1% of the viral population; the cut-off chosen for our assay. Substitutions occurring during PCR could be addressed using, for instance, random sequence primer identifiers (ID) [26]. Again, our approach to set the cut-off for minority HIV-1 drug-resistance mutations to 1% minimizes false positive calling of most of PCR-introduced mutations except those which occur in the first PCR cycles. Another source of potential error is artificial recombination, the rate of which is highly correlated with the number of DNA templates used in the PCR reaction [14, 27, 28]. The current assay does not contain a step to normalize the transferred cDNA molecules and may therefore not be suitable to link certain drug resistance mutations at low levels, especially in very high viral load samples.

Finally, we must consider errors introduced during the process of pyrosequencing. The 1% cut-off for minority HIV-1 drug-resistance mutations is approximately a magnitude higher than the reported substitution rates introduced by pyrosequencing [29, 30]. In addition, mutations must be detected in both forward and reverse reads at the approximate same frequency to be deemed real. Insertions and deletions (indels) occur predominantly in homopolymer regions [29, 30] and artifactual variant calling can be mitigated by using appropriate platform specific analysis software such as AVA to appropriately align reads and by manually inspecting the alignments in regions with homopolymer runs. Given the very high reproducibility of variants across sites in this study, we feel confident that we have adequately addressed these sources of error for potential use in a clinical setting. In addition, we have greatly reduced the possibilities that exist for introducing human error during manual handling in the laboratory. Compared to the initial multicenter study [10], where several severe errors occurred, the new concept of providing dried down primers in plates proved very successful. In addition to “locking down” the positions of cDNA and PCR primers, we introduced a tracking MID on each plate to ensure that pools from each plate ended up in the intended region on the sequencing run. Additional improvements in the making are color coding and orientation marks on individual plates. As long as materials for individual sequencing runs are kept separate, the main source of manual error remains potential sample swaps, a problem that all clinical laboratories face and have learned to minimize.

454 ultra deep sequencing will only be supported through 2016. Nevertheless, this assay design, either directly or via primer recombination, could be transferred to other next generation sequencing platforms allowing long amplicon sequencing provided by, for instance, Illumina [31, 32], Pacific Biosciences [32], Ion Torrent [32], or Oxford Nanopore Technologies. Regardless of the technology used, our inter- and intra-laboratory testing in a multicenter collaborative study shows that this kind of validation is of the utmost importance when UDS assays are to be considered for implementation in a routine diagnostic setting.

Conclusion

The field of virology continues to study minority HIV-1 drug resistance mutations and their relevance in the drug resistance field, although the clinical impact is still controversially discussed [8, 9, 33] and there remain challenges such as the search for clinically relevant cut offs. Nevertheless, as next-generation sequencing technologies become more affordable and accessible, they will more regularly play a key role in these studies due to their accuracy, reproducibility and throughput. Thus, intra- and inter-laboratory testing is important for implementation of ultra deep sequencing assays in clinical laboratories for routine diagnostics. The updated assay design described herein for resistance testing and co-receptor usage prediction by ultra deep sequencing allows for an accurate and highly reproducible approach to analysis of the HIV-1 mutant spectra, even at frequencies well below those detected by standard population sequencing.

Fig A: Overview of amplicon locations. Reference amplicon sizes (incl. adapters and MIDs) are as follows: (A) RTP1 = 419 bp, RTP2 = 510 bp, RTP3 = 400 bp, RTP4 = 599 bp, RTP5 = 558 bp, RTP6 = 434 bp; (B) V3A = 465 bp, V3B = 423 bp. Fig B: HQ reads and aligned reads distribution. I) The number of HQ reads, aligned and not aligned, per sequencing region of the PTP for each site II) The percent composition of total sequencing reads (in brackets below the site name) obtained from each amplicon at each site using region 3 as a representative example. Fig C: AVA software screen shot of the global alignment of sample 3 from site 31 of low frequency mutation RT215. Population sequencing indicated the possibility of amino acids Threonine (ACT), Tyrosine (TAT), Serine (TCT), and Asparagine (AAT) at codon 215 in the reverse transcriptase. These could not be further resolved, as the codon nucleotides cannot be phased. The 454 pyrosequencing reads resolved the actual variant codon compositions: 52.0% TAT (Tyr), and 46.9% ACT (Thr) (only mutations present at >1% included). The variants TCT (Ser) and AAT (Asn) could not be detected. Table A: Overview of sample viral titers and clades used for this study. 30 HIV-1 subtype B and 6 subtype non-B samples were used for RTP and V3 sequencing, three HIV-1 subtype B samples were run in triplicate (samples 37, 38 and 39) and two HIV-1 subtype B samples were used for the dilution series (samples 40–44 and 45–49). The RTP, dilution samples and triplicate samples are cultured recombinant virus, whereas V3 samples are cultured samples from patient isolates. The viral titer for each is shown. Table B: Primers. Table C: Details of amplicon drop out across sites per amplicon. Lack of amplicon production was reported by each site as any amplicon measuring less than 1 ng/μL or that was comprised of primer dimer alone. All amplicon concentration measurements were taken following PCR and purification. Table D: Drug resistance mutations found only by UDS. Details on reports of each site. The summary of the data are given in table 4. Hits, number of reads reporting the variant; denom, number of total reads at this position.

(DOC)

Click here for additional data file.