article: Estimating the time since infection (TI) in newly diagnosed HIV-1 patients is challenging, but important to understand the epidemiology of the infection. Here we explore the utility of virus diversity estimated by next-generation sequencing (NGS) as novel biomarker by using a recent genome-wide longitudinal dataset obtained from 11 untreated HIV-1-infected patients with known dates of infection. The results were validated on a second dataset from 31 patients. Virus diversity increased linearly with time, particularly at 3rd codon positions, with little inter-patient variation. The precision of the TI estimate improved with increasing sequencing depth, showing that diversity in NGS data yields superior estimates to the number of ambiguous sites in Sanger sequences, which is one of the alternative biomarkers. The full advantage of deep NGS was utilized with continuous diversity measures such as average pairwise distance or site entropy, rather than the fraction of polymorphic sites. The precision depended on the genomic region and codon position and was highest when 3rd codon positions in the entire pol gene were used. For these data, TI estimates had a mean absolute error of around 1 year. The error increased only slightly from around 0. 6 years at a TI of 6 months to around 1. 1 years at 6 years. Our results show that virus diversity determined by NGS can be used to estimate time since HIV-1 infection many years after the infection, in contrast to most alternative biomarkers. We provide the regression coefficients as well as web tool for TI estimation.
  At diagnosis, most HIV-1 infected patients have an established HIV-1 infection of unknown duration. This uncertainty complicates inference about the epidemiology of HIV-1. Consequently, there is limited information about the true incidence of HIV-1, the number of hidden, undiagnosed infected persons, the magnitude of the problem referred to as “late presentation” and other important aspects of HIV-1 spread. Several biomarkers that classify patients as recently or long-term infected have been used to estimate HIV-1 incidence in populations [1–7]. These biomarkers can be divided into three main categories: (i) serological incidence tests, (ii) CD4+ T-lymphocyte (CD4) -based estimates and (iii) sequence-based estimates. Importantly, these biomarkers usually do not determine the time since infection (TI), which limits their utility. Serological incidence assays are based on knowledge about the development and maturation of HIV-1 antibody responses (reviewed in [1,4–6,8]). Among the serological assays, the BED assay and the LAg avidity assay are the most widely used [4,8–10]. CD4 counts are determined as part of routine clinical care, a CD4 count below 350 cells/μL (or an AIDS-defining illness) at diagnosis is defined as late presentation [11,12]. However, CD4 count is an imprecise measure of TI because its rate of decline is quite variable [13–16]. Sequence-based methods focus on the increase in intra-patient HIV-1 sequence diversity following infection [17]. Kouyos et al [18] showed that time since infection correlated with the fraction of polymorphic nucleotides in partial HIV-1 pol gene sequences determined by Sanger sequencing. Others have later reported similar findings [19–21]. This idea was expanded to other measures of sequence diversity, such as mean Hamming distance [6,7, 22] and high-resolution melting (HRM) [23]. These studies (except HRM) used sequences generated by traditional Sanger population sequencing often performed as part of routine HIV-1 resistance testing. Here, we have investigated the utility of estimating time since HIV-1 infection using genetic diversity in whole genome deep sequencing data generated by next-generation sequencing (NGS) on the Illumina platform [24]. We show that sequence diversity is a useful biomarker that grows approximately linearly with time during the first 8 years of infection. We found that the pol gene was best suited to calculate TI because diversity, mostly at third positions, accumulated more steadily in pol than in other genomic regions. Inclusion of intra-patient single nucleotide variants (iSNVs, also referred to as “polymorphisms”) down to the detection limit of NGS (i. e. 0. 3%) improved the accuracy of TI estimations as did exclusion of 1st and 2nd codon position sites. NGS provided more accurate estimates of TI than Sanger sequencing, which at best detects iSNVs down to 25% [25–27].
  The study included two sets of HIV-1 whole-genome sequence data; a training dataset and a validation dataset. The training dataset consisted of sequence data from 11 patients who were diagnosed in Sweden between 1990 and 2003 (Table 1 and S6 Table). Details about the patients and the samples have been published [24] and are also available online at hiv. biozentrum. unibas. ch. The patients were selected based on the following criteria: 1) A relatively well-defined time of infection based on a laboratory-confirmed primary HIV-1 infection (PHI) or a negative HIV antibody test less than two years before the first positive test; 2) No antiretroviral therapy (ART) during a minimum of approximately five years following diagnosis; and 3) Availability of biobank plasma samples covering this time period. As previously described 6–12 plasma samples per patient were retrieved from biobanks and used for full-genome HIV-1 RNA deep sequencing [24]. The validation dataset consisted of data from 31 patients who were diagnosed in Sweden between 2003 and 2010 (S7 Table). The patients were selected using similar criteria as the 11 patients in the training dataset, but the follow-up time without ART was shorter (median 2. 9 years, range 1. 4–6. 2 years). For each patient, one plasma sample collected early during follow-up and one plasma sample collected late during follow-up was deep sequenced using the same method as for the 11 patients. The study was conducted according to the Declaration of Helsinki. Ethical approval was granted by the Regional Ethical Review Board in Stockholm, Sweden (registration no. 2012-505 and 2014-646 for the validation dataset and 2002-367,2004-797,2007-1533 and 2011-1854 for the training dataset). Patients participating in the study gave written and oral informed consent to participate. For each patient in the training and validation datasets the “true” time of infection (TI) was estimated using a hierarchical scheme based on clinical and laboratory findings as previously described [24,28]: If possible information from several methods to determine TI were combined. The true TIs were considered to be without measurement error. Comprehensive information is provided in the S6 and S7 Tables. Plasma HIV-1 RNA levels were measured using the Cobas AmpliPrep sample preparation system followed by analysis using the Cobas Amplicor HIV-1 monitor version 1. 5 or the Cobas TaqMan HIV-1 v1. 0 or v2. 0 (Roche Molecular Systems, Basel, Switzerland). CD4+ T-lymphocyte (CD4) cells were enumerated using flow cytometry. As part of determination of the true TIs, BED testing was performed on the first plasma sample from each study subject using the Aware BED EIA HIV-1 Incidence Test (Calypte Biomedical Corporation, Portland, OR, USA) according to the manufacturer’s instructions. Whole-genome deep-sequencing of virus RNA populations in plasma samples obtained before start of therapy was performed as previously described [24,32]. In short, total RNA in plasma was extracted using RNeasy Lipid Tissue Mini Kit (Qiagen Cat No. 74804) and amplified using a one-step RT-PCR with outer primers for six overlapping regions and Superscript III One-Step RT-PCR with Platinum Taq High Fidelity High Enzyme Mix (Invitrogen, Carlsbad, California, US). An optimized Illumina Nextera XT library preparation protocol was used together with a kit from the same supplier to build DNA libraries, which were sequenced on the Illumina MiSeq instrument with 2x250bp or 2x300bp sequencing kits (MS-102-2003/MS-10-3003). Sequencing reads are available in the European Nucleotide Archive under accession number PRJEB9618 and processed data are available at hiv. biozentrum. unibas. ch. The samples of the validation dataset were sequenced and processed using the same protocol and analysis pipeline [24,32]. Patient-specific consensus sequences were constructed using an iterative mapping procedure. All reads were then remapped against this reference to calculate iSNV frequencies (i. e. pile-ups or tables how often each nucleotide was observed at every position of the genome). Sequencing and mapping/assembly was successful for 56 of the 62 samples. Filtered short-reads were submitted to ENA and are available under study accession number PRJEB21629. iSNV frequency counts at each position of pol and gag are available as part of the analysis code repository at github. com/neherlab/HIV_time_of_infection. All analyses were done in Python using the libraries numpy, biopython, and matplotlib [33–35]. These iSNV frequency tables were then used to calculate average pairwise distances, average alignment entropies, or the number of sites with variation above a cutoff xc. We have used three different diversity measures: fraction of polymorphic sites, average pairwise distance per length, and site entropy. All of these measures can be straight-forwardly calculated from the frequencies of different nucleotides xiα at site i = 1…L and α ∈ {A, C, G, T} along the genome. Prior to calculation the nucleotide frequencies for each site were normalized to sum to unity (i. e. ignoring gaps or positions not called by the sequencer.) For all methods, we introduce a cutoff xc. Sites at which the sum of all minor variants is smaller than xc contribute zero to the diversity measure. This cutoff serves to filter sequencing errors or rare variation that cannot be reproducibly measured across samples. When using the fraction of polymorphic sites as diversity measure, xc serves as the value above which sites are considered “polymorphic”. Specifically, the fraction of polymorphic sites is defined as D A = 1 L ∑ i = 1 L Θ (1 - x i m - x c) (1) where x i m is the frequency of the dominant nucleotide at position i, and Θ (x) is 1 when x > 0 and 0 otherwise (i. e. Θ (1 - x i m - x c) = 1 when 1 - x i m > x c and 0 otherwise). DA is thus the fraction of sites at which the dominant nucleotide is less frequent than 1 − xc. The average pairwise distance is the probability that two randomly drawn sequences have different nucleotides at a specified position, averaged over all positions. It can be calculated from the xiα as D H = 1 L ∑ i = 1 L Θ (1 - x i m - x c) [ ∑ α x i α (1 - x i α) ]. (2) The quantity defined by Eq (2) is the conventional Nei-Li nucleotide diversity [36] (∑α xiα (1 − xiα) ) averaged over the sites. We refer to it as “average pairwise distance” whenever it is necessary to distinguish it from the other diversity measures introduced here, but otherwise call it simply “diversity”. This diversity measure is similar to the fraction of polymorphic sites with the important difference that the contribution of each site is weighted by a frequency dependent factor. The average site entropy is defined by D E = - 1 L ∑ i = 1 L Θ (1 - x i m - x c) [ ∑ α x i α log (x i α) ] (3) and differs from the average pairwise distance by the weighing function used. The entropy puts more weight on sites with rare variation. This can increase the information of the measure about TI, but can also be detrimental if too much weight is put on frequencies dominated by sequencing error. We evaluate and discuss the merits of the different measures below. We use the average pairwise distance as default diversity measure. Given a diversity measure D, we model TI by TI ^ = s × D + t 0, (4) where s is the conversion factor between diversity and time, and t0 is the intercept value intended to accommodate possible non-linearity of diversity at small times. To estimate values of s, t0 we minimized the average prediction error for the available data points in respect to these two parameters (see S1 Appendix for more details) The error in estimating s, t0 was calculated by randomly sampling the patients (bootstrapping over the patients). To test the accuracy of the TI inference we used ten of the eleven patients as training data (to determine the slope and the intercept) treating the eleventh patient as test data. This procedure was repeated for every patient. Leaving out one patient at a time, rather than one sample at a time, gives more accurate confidence intervals as different samples from the same patients are correlated. We included in statistical analysis only samples where more than 50% of the sites in the averaging window were successfully sequenced to minimal coverage of 100. Additional validation was provided by applying the slopes and intercepts obtained by analyzing the data from the training set of 11 patients to the validation dataset of 31 other patients with known times of infection.
  The training dataset consisted of recently published longitudinal full-genome deep sequencing data from 11 HIV-1 infected patients who were diagnosed in Sweden between 1990 and 2003 (6–12 samples per patient) and had a relatively well-defined time of infection (TI) [24]. The patient characteristics are summarized in Table 1 and fully described in the S6 Table. Nine of the eleven patients were MSM infected with HIV-1 of genetic subtype B. TI was hierarchically estimated using clinical and laboratory data (see Methods). Here, we take this estimate as the true time of infection and investigate how accurately TI can be estimated from sequence diversity in one sample. We will refer to this estimate as the estimated time since infection (ETI). To validate the findings we used a second dataset consisting of similar sequence data from 31 additional patients (two samples per patient). The patient characteristics in validation dataset was more diverse than for the 11 patients training dataset. Thus, 16 of the 31 patients were infected with non-B-subtypes of HIV-1 and 11 patients belonged to other transmission groups than MSM. See methods and (S7 Table). All three diversity measures described in Materials and Methods grew linearly with time in the eleven patients. Fig 1 shows average pairwise distance in the pol separately for each codon position. Most diversity in pol is synonymous and accumulates at 3rd codon positions, while diversity at 1st and 2nd codon positions remained low throughout. This pattern was less pronounced in other genes [24,37] (see also S1 and S2 Figs). In env, in particular, frequent selective sweeps result in a saturation of diversity later in infection [17,24]. We quantified the fraction of variation in diversity measures that could be explained by a linear regression of sample date vs. diversity using the coefficient of determination (r2), see the S3 and S4 Figs. In all patients for whom early samples were available, a linear regression explained between 70% and 90% of variation if rare iSNVs below 20% population frequency were included, that is the cutoff xc was below 0. 2. The coefficient of determination was much smaller when only iSNVs between 20% and 80% were included, that is the xc cutoff is larger than 0. 2. This decrease is due to increased noise as fewer and fewer sites contribute to the diversity measures. Furthermore, as seen from S4 Fig, the diversity at 3rd codon positions (at which mutations are mostly synonymous) exhibited higher r2, whereas the trajectories at the 1st and 2nd codon positions saturated quickly after the infection, Fig 1. Thus, in the following we limit the analysis mainly to sites in 3rd codon positions (whenever we are dealing with a whole gene, i. e. when the reading frame is known.) However, the results reported below show that inclusion of 1st and 2nd codon positions only has a small deleterious effect on the accuracy of the TI estimates. Accurate estimation of TI requires averaging diversity across many sites. To investigate which regions of the genome yields the most accurate estimates and how many sites should be averaged, we calculated the average prediction error for averaging windows of different length and in different regions of the genome. In Fig 2 the mean absolute error (MAE) for the estimated TIs is shown as a function of the genome window position for different window sizes. We found that the most precise TI estimates were obtained for windows with a length of 2000–3000bp. The precision was highest if the window covered the pol gene and significantly lower if the window included env. Smaller windows contain fewer sites and therefore gave less precise estimates. Larger windows also gave less precise TI estimates because they necessarily include regions in which diversity saturates (such as env, S2 Fig) as well as regions that often were sequenced less deeply in our dataset (again env, see [24]). The precision of the TI estimates obtained for windows corresponding to particular genes are shown in Fig 2 by lines indicating the position of the gene and the corresponding average absolute error of the TI estimate. The dashed lines correspond to estimates using only 3rd codon positions of the corresponding gene. Genetic diversity in NGS data can be quantified by different measures and we investigated the performance of three related measures—average pairwise distance, site entropy, and fraction of polymorphic sites. As discussed above, these measures put different emphasis on iSNVs at different frequency. Even though the average pairwise distance and site entropy formally do not require a cutoff on minority iSNVs, in practice a cutoff is necessary to remove low-level experimental errors. We therefore introduced an iSNV cutoff in calculating the diversity based on average pairwise distance and site entropy, xc, taking into account only the data with iSNVs above xc. For diversity measures based on polymorphic sites, the cutoff value corresponds to the iSNV level from which the site is considered polymorphic. This emulates ambiguous base calls by Sanger sequencers, which (at best) can identify minority iSNVs above a threshold of around 25% [25–27]. Given the high ability of NGS to detect low level iSNVs one can consider the dependence of the TI estimation on the cutoff value, as shown for all three diversity measures in Fig 3. All three diversity measures performed equally well and increasingly better at cutoff values xc down to around 10%. This indicates that diversity in NGS data allows more accurate estimates of TI than ambiguous base calls in Sanger sequences and that this primarily is due to better sensitivity and accuracy in detection of low level iSNVs. For cutoffs below 10% the error of TI estimates based on counting polymorphic sites was greater than those based on the other two measures since rare iSNVs which are sensitive to sequencing errors and amplification biases contribute as much as common iSNVs. Indeed, at xc → 0 this measure includes all sites and becomes completely insensitive to differences in iSNV levels. The other two measures do not suffer from this problem, because they put different weights on sites with different iSNV levels. Therefore the prediction errors for estimates based 3rd codon positions and average pairwise distance or site entropy continued to decrease all the way down to xc = 0. 003, which represents the cutoff for sequencing errors for our NGS method [24]. The noticeably non-monotonic behavior of the predictions based on average pairwise distance and site entropy when all codon positions are taken into account is due to the saturation of 1st and 2nd codon positions (i. e. non-synonymous) diversity. The 1st and 2nd codon positions tend to be conserved since they result in mostly non-synonymous mutations. Depending on the fitness cost associated with the mutation, diversity at these sites saturates at different levels [37]. As the threshold xc is lowered, sites with higher and higher fitness costs contribute to the diversity measures and the effect of the saturation behavior becomes more pronounced. Note that the non-monotonic dependence is not consistently reproduced in other genes (see supplementary S5 Fig). Similarly, the fact that diversity measures based on all codon positions perform somewhat better at xc > 10% was not consistently reproduced in other genes. In gag, diversity at 3rd positions was better at estimating TI for most xc values (S5 Fig). In the following analyses we opted for using the average pairwise distance measure, taking into account only the diversity at the 3rd codon positions, with low (xc = 0. 003) iSNV cutoff. Average pairwise distance was chosen because the results were virtually indistinguishable from those produced using site entropy, but easier to calculate and interpret. Since 1st and 2nd codon positions contribute very little time-dependent diversity and are affected by purifying selection and selective sweeps, we recommend to restrict the diversity measures to 3rd codon positions. The results above indicated that more than 50% of the estimated TIs fell within a window of about one year centered at the actual TI. Fig 4 (Left) shows a more direct analysis of the distribution of TI prediction errors. The distribution is tightly peaked around zero, but has a left tail corresponding to samples estimated to have been obtained shorter after infection than they actually were drawn, i. e. diversity being lower than expected. Most of these samples were from p6, who throughout infection had lower diversity than other patients. In addition to biological reasons for low diversity, amplification problems and low RNA template numbers (i. e. low virus levels) can explain underestimation of diversity. Some samples were estimated to have been drawn later after infection than the true duration of infection. In particular early samples from p10 and p3 had markedly higher diversity than expected. For both patients, we have evidence that their infections were established by more than one virion resulting in carry-over of diversity from the donor. This excess diversity gradually decreased in p10 and, somewhat slower, in p3. Next, we investigated how the prediction error depended on the time since infection. Fig 5 shows the average absolute error of the estimated TI versus the true TI, averaged over n = 25 consecutive data points. This average error (see for details S2 Appendix) was surprisingly stable over TIs, and only increased slightly from around 0. 6 years to around 1. 1 years as the age of infection increased from 6 months to 6 years. This increase can be attributed to bigger statistical fluctuations of diversity later after infection due to factors such as variations in the rate of diversification or differences in number and strength of selective sweeps that reduce diversity. In Table 2 we list the values of slope and intercept that can be used to estimate the infection date from the known diversity calculated as average pairwise distances for 3rd codon positions in pol gene. The supplementary materials contain similar tables for the two other diversity measures (S1 and S2 Tables), as well as for the case when all codon positions are taken into account (S3, S4 and S5 Tables). As the iSNV resolution may vary between different sequencing methods and facilities, we list the values of the parameters for different cutoffs, implying that all the frequencies below the cutoff value are set to zero and the corresponding sites do not contribute to the diversity measure. Note that the slope (and intercept) increases with increasing iSNV cutoffs because fewer and fewer sites contribute to diversity. In addition to the two parameter model, we also investigated the performance of a model with the slope as the single parameter, i. e. no intercept (t0 = 0). This model has a slightly higher absolute prediction error. However, for low values of the cutoff xc ≈ 5% these models agree and for cutoffs below 20% the two models perform equally well. In order to make our data and the method of TI inference more accessible for practical use, we added a web application to the web site containing the processed patient data, accessible at hiv. biozentrum. unibas. ch/ETI/. Given a diversity value (determined according to Eq (2) ), the web application allows to determine the time since infection for an user-selected iSNV cutoff value (xc) and genetic region, along with the appropriate error estimates. The results are presented in accessible graphical form, but also as a slope and intercept pair, equivalent to Table 2. The user can specify whether the diversity was calculated using all codon positions or only the 1st, 2nd, or 3rd. Specifying a codon position is only supported when the region used to evaluate diversity is fully contained in one gene. We validated the regression coefficients on a dataset from 31 patients with known infection dates and NGS data available at two time points. The distribution of the diversity values for these patients closely resembles that for our training dataset of 11 patients (see S8 and S9 Figs). We inferred the time since infection for the 31 validation patients using the regression coefficients obtained for the eleven patients of the training set; the results are summarized in Fig 6, which also shows (in gray) the data points of the training set (same as in Fig 4). The new data exhibit the same behavior as the training set: the TI estimates are centered around the true TI and the accuracy of the estimate is about one year, as can be seen from the histogram in the left panel of Fig 6. In order to make the histograms for the training and validation data comparable, we included for the former only the points of the infection time less than 5 years. Some outliers are present also in validation data, particularly data points with overestimated TIs, i. e. having higher diversity than expected. As mentioned above, these data points probably often represent infections established by more than one founder virus. One patient had substantial overestimation of time since infection, which might be due to superinfection from different donors rather than multiple founders from a single donor. However, it should be stressed that almost all outliers are still within +/- 2 years from the true TI. Neither overestimation nor underestimation of TI was clearly related to genetic subtype of the virus, transmission route, virus levels or template numbers (S11 and S12 Figs).
  In conclusion, we show that sequence diversity determined by NGS can be used to estimate time since HIV-1 infection with a precision that is better than most alternative biomarkers. Importantly, TI can be estimated many years after infection, whereas most alternative methods only categorize patients as being recently or long-term infected or are less precise. We found that TI was most accurately estimated using 3rd codon positions in pol sequences with a xc = 0. 003 cutoff for iSNVs and that average pairwise distances was the preferred distance measure. Samples with low virus levels and infections established by multiple virions can give rise to misleading levels of virus diversity. Algorithms based on NGS diversity in combinations with other biomarkers may prove very useful.
summary: HIV-1 establishes a chronic infection, which may last for many years before the infected person is diagnosed. The resulting uncertainty in the date of infection leads to difficulties in estimating the number of infected but undiagnosed persons as well as the number of new infections, which is necessary for developing appropriate public health policies and interventions. Such estimates would be much easier if the time since HIV-1 infection for newly diagnosed cases could be accurately estimated. Three types of biomarkers have been shown to contain information about the time since HIV-1 infection, but unfortunately, they only distinguish between recent and long-term infections (concentration of HIV-1-specific antibodies) or are imprecise (immune status as measured by levels of CD4+ T-lymphocytes and viral sequence diversity measured by polymorphisms in Sanger sequences). In this paper, we show that recent advances in sequencing technologies, i. e. the development of next generation sequencing, enable significantly more precise determination of the time since HIV-1 infection, even many years after the infection event. This is a significant advance which could translate into more effective HIV-1 prevention.
keywords: distance measurement
  sequencing techniques
  medicine and health sciences
  pathology and laboratory medicine
  evolutionary biology
  engineering and technology
  pathogens
  population genetics
  microbiology
  biomarkers
  retroviruses
  viruses
  immunodeficiency viruses
  next-generation sequencing
  rna viruses
  genome analysis
  human genetics
  molecular biology techniques
  population biology
  thermodynamics
  research and analysis methods
  genomics
  entropy
  genetic polymorphism
  medical microbiology
  hiv
  microbial pathogens
  hiv-1
  molecular biology
  nucleotide sequencing
  physics
  biochemistry
  viral pathogens
  transcriptome analysis
  genetics
  biology and life sciences
  physical sciences
  computational biology
  lentivirus
  dna sequencing
  organisms
  measurement
title: Estimating time of HIV-1 infection from next-generation sequence diversity