BACKGROUND
next-generation sequencing  techniques  <cit>  allow for a high-throughput dna sequencing, producing from thousands to billions of sequence fragments  composed of tens to hundreds of nucleotide bases.
ngs has the potential to replace sanger sequencing for many applications, including de-novo sequencing, re-sequencing, meta-genomics and intra-host characterisation of infectious pathogens  <cit> .
de-novo sequencing implies a genome assembly problem, which is the reconstruction of a unique genome from a set of sequence fragments.
several methods and software for genome assembly have been developed  <cit> .
these methods were designed initially for sanger sequencing, and have been revised for ngs technology  <cit> , given different error rates among ngs machineries  <cit> .
re-sequencing conjugates with the problem of single nucleotide polymorphisms  discovery.
recent studies characterised snps or drug-induced mutations with ngs, considering the human immunodeficiency virus  and the hepatitis b virus   <cit> .
more specifically, re-sequencing can be useful for the characterisation of variants within a quasispecies harbouring an infected host.
many rna viruses are present in a carrier  as a swarm of highly genetically related variants, i.e.
a quasispecies, due to the error prone characteristics of the viral polymerases and high viral replication rates.
this intra-host variability represents a substrate for the selective pressure exerted by the immune system of the host or by drug exposure, which leads to the continuous evolution of viruses.
quasispecies reconstruction would allow detailed description of the composition of individual viral genomes, genetic linkage and evolutionary history.
for example, in hiv or hbv infection, the development of drug resistance is a major problem and the early diagnosis of drug-resistant variant selection might help in designing targeted therapeutic interventions.
here, we addressed the problem of reconstructing a viral quasispecies from a ngs data set, which is a relatively new topic that is not widely investigated in literature.
we aimed to reconstruct all coexistent individual variants within a population, along with their prevalence, rather than a reconstruction of a single or predominant genome.
current assembly software is not designed to accomplish this task, nor to deal easily with the reconstruction of highly variable genomes.
the huge coverage and base pair output provided by ngs enables the design of experiments to investigate and validate theoretical methods for quasispecies reconstruction.
at present, only a few methodological papers have been published presenting new algorithms for quasispecies reconstruction that are able to infer both genomes of population variants and their prevalence  <cit> .
in  <cit> , the authors proposed an algorithm based on a generative model of the sequencing process and a tailored probabilistic inference and learning procedure for model fitting.
in  <cit> , a set of methodologies was proposed both for error correction and inference about the structure of a population from a set of short sequence reads as obtained from ngs.
the authors assumed a known mapping of reads to a reference genome, defined a sliding window over the reference genome and associated each aligned read to one or more windows by trimming the reads accordingly to the windows' bounds.
sequencing errors were corrected by locally clustering reads in the windows.
a set of single variants of the quasispecies  was obtained by constructing an overlap graph of non-redundant, error-free, aligned reads, and by calculating a minimal coverage set of paths over the graph.
the frequency estimation was done with an expectation maximisation algorithm and was proven to be more efficient than a naïve procedure based on uniform read sampling.
one drawback of this methodology is that the variant reconstruction phase did not account for the relations among frequencies of distinct variants  that were overlapping consistently across the sliding windows: this may lead potentially to selection of in-silico recombinants and the procedure of haplotype frequency may be biased from the exclusion of real  paths.
after the paper, free software was released, named shorah  <cit> .
in  <cit> , a scalable assembling method for quasispecies based on a novel network flow formulation was presented, applied efficiently for the assembly of hepatitis c virus.
in  <cit> , a refinement of the original procedures presented in  <cit>  was given, substituting k-means clustering with a dirichlet process mixture for locally inferring haplotypes and correcting reads.
in this work, a set of formulae for combinatorial analysis of quasispecies genome fragments sampled by ngs was derived, and a new greedy algorithm for quasispecies reconstruction was introduced.
the formulae derivation provided some theoretical bounds explaining the difficulty in reconstructing a set of individual variants of the quasispecies, by conditioning on several parameters, such as the genome length, the fragment  size, or the overlap length between two sampled fragments.
the reconstruction algorithm was based on combinations of multinomial distributions and was designed to minimise the reconstruction of in-silico recombinants.
unlike previous approaches, the algorithm selects and reconstructs variants not only by coupling reads that have consistent overlaps, but also considering reads that have similar frequencies across the various amplicons.
for our combinatorial analysis, we assumed that the problem of re-sequencing, including reference alignment and error correction, is solved.
in other words, a set of error-free reads is available, aligned univocally to a reference or consensus sequence.
such a reference may either have been directly reconstructed using assembly software or selected from literature.
the assumption for the unique mapping of a read against the reference may not be always fulfilled when in presence of short reads and genomes with long repeats.
however, this problem can be negligible when considering coding regions of highly-variable viral pathogens targeted by inhibitors, with a few regulatory regions where the repeats usually are located.
the sequence quality may be a major concern for reconstruction algorithms and this is often an experiment- and machine-dependent problem: procedures for alignment and error corrections have been investigated elsewhere  <cit> , with different methodologies, along with protocols for sample preparation.
another critical point with ngs is the presence of contaminants that must be detected and excluded.
the problem can be solved easily when the contamination is from different organisms, with a test statistic on read/reference alignment scores during re-sequencing  <cit> .
it is harder when the ngs experiment comprises a mixture of closely related organisms, for instance when samples of patients infected with the same virus are put together in one ngs experiment  <cit> .
as a second assumption, our algorithm required a non empty set of overlapping regions, called amplicons, which cover the reference genome.
each read has to be assigned to one of these amplicons.
roche  <dig> gsflx technology has a double working modality that allows both for shotgun sequencing and for amplicon sequencing with specific primer design, although the latter option is generally more expensive.
with this technology, amplicons can be defined a priori.
in contrast, if shotgun sequencing is performed, additional data elaboration has to be made in order to define a set of amplicons: one solution is to define amplicons via sliding windows over the genome and cluster reads accordingly to their mapping region  <cit> .
the proposed reconstruction algorithm was applied to 1) simulated and error-free data; and 2) then empirical sequence data derived from blood samples from hbv-infected patients processed via the roche  <dig> gsflx titanium machine.
this second dataset was designed to assess the performance of quasispecies reconstruction in presence of sequencing errors.
RESULTS
algorithm: ngs data processing and amplicon definition
this work analysed a re-sequencing experiment of a viral quasispecies carried out using ngs machinery.
since currently the maximum read length of a ngs experiment does not exceed a few hundred of bases, we were interested in genome regions of quasispecies whose length is much larger than the average read length, i.e.
when it is not possible that a read spans entirely the genome of interest.
we required then that a reference genome is available and that reads are significantly aligned  against this reference genome.
this can be achieved by aligning each read in forward- or reverse-strand against the reference genome, using the smith-waterman-gotoh local alignment  <cit> , which is an exact algorithm, and keeping the highest alignment score.
reads then can be filtered by excluding those that do not show a significant  alignment score, as compared to a score distribution obtained from quasi-random sequences  aligned to the reference genome, as described in  <cit> .
we assumed also that sequencing errors were corrected.
the condition of error-free reads was required only for the combinatorial analysis, whilst the quasispecies reconstruction algorithm can be applied also to noisy data .
given a reference genome g and a read alignment over g, we define then a sliding window partition of g composed of w+ <dig> windows, that we call amplicons.
these amplicons cover the entire genome and two adjacent amplicons have a partial overlap, for a total of w overlaps.
amplicons do not need to be necessarily of the same length.
clearly, each amplicon size has to be smaller than the average read size, so that a read can span an amplicon entirely.
as stated in the introduction, amplicons can be designed a priori if roche  <dig> gsflx technology is used, or determined with a fixed sliding window approach from any shotgun sequencing.
after defining the amplicons, reads that spanned entirely an amplicon were trimmed so that their start/end positions corresponded exactly to the amplicon start/end.
consequently, all the reads in one amplicon had the same length.
reads that span more than one amplicon entirely are considered multiple times, whilst those that do not cover exactly at least one amplicon were discarded.
figure  <dig> illustrates with an example a three-amplicon design over a reference genome, with corresponding read assignment and trimming.
algorithm: combinatorial analysis of ngs over a quasispecies
consider a number of x variants  determined by their genomic sequences of length n, either nucleotidic or amino-acid .
the number x of variants is unknown, along with the prevalence of each variant.
we assume that the quasispecies is stable over a fixed set of variants, in a mutation-selection balance  <cit> .
in other words, the number of distinct viral variants in the quasispecies that a carrier  harbours is x, although each variant can be present with different prevalence, presumably due to different viral fitness or immune response or drug-induced selection.
after a multiple sequence alignment, a consensus sequence can be generated or one variant can be used as a reference.
we define a point difference between two aligned sequences  as the presence of two different nucleotides in one position of the alignment.
the pairwise difference between two variants d is the number of point differences divided by the genome length , i.e.
d=dij=∑k=1nn
the average pairwise difference of a variant si with respect to all the others is defined as
d=di=/, for j= <dig> ..x.
finally, the average overall pairwise difference among a set of aligned variants is then
davg=∑i=1x−1∑j=i+1xdijx <dig>
let us define now a reference sequence sref as the variant with the lowest average pairwise difference as compared to all other variants, i.e.
dref = d = min{d}, for i =  <dig> ..x.
if we define the diversity - that we regard as a probability of mutation - of the quasispecies as m = davg, then we can also approximate m by doubling the value of dref, i.e.
m/ <dig> is the average pairwise difference of our reference variant with respect to any other variant .
indeed, this approximation is correct when no identical base changes  happen in the same alignment position of two variants as compared to the reference, and this is dependent on the genome length and the mutation probability.
the approximation gets better either if the genome length increases or the mutation rate decreases .
more details on the efficacy of this approximation are given in additional file  <dig>
we previously defined a joined set of w+ <dig> amplicons, which induces an overlapping ordered coverage over the quasispecies genome space.
assume that each amplicon has a fixed length of k bases and overlaps with its neighbour over q bases for w times.
we assume also that, given three adjacent amplicons and two corresponding overlaps, these latter do not share any position in common, i.e.
there are no overlapping overlaps.
since there are no nested amplicons, we can define an amplicon identification number as its ordinal position with respect to the reference genome and the other amplicons.
thus, amplicon <dig> starts at position  <dig> and ends at position k, amplicon <dig> starts at k-q+ <dig> and ends at 2k-q, et cetera.
the overlaps are clearly at the end of each amplicon and at the beginning of the adjacent one.
each amplicon is associated with a set of reads, or sequence fragments, sampled uniformly from the quasispecies.
these reads, by definition, are significantly aligned to the reference genome, error-free and span exactly an amplicon region, being trimmed at its ends.
thus, after sampling, each count of distinct reads across each amplicon cannot exceed the number of variants x.
given two reads associated to two adjacent amplicons, we say that their overlapping region is consistent if the two reads share in that region the same characters .
we aim to calculate the probability that  the overlapping region of two adjacent reads is consistent and  at least one overlapping region across the amplicons is consistent.
for point , we first define the probability that i mutations are present in a sequence fragment of q length over a genome of length n  as
p=
note that the diversity m/ <dig> here is multiplied by n , obtaining the expected number of changes .
the probability that two random regions of q length over a genome of length n have both i mutations is
p=) <dig>
where the terms of eq.
<dig> are the square root of the terms of eq.
<dig>
the probability that these random regions share exactly i mutations at the same positions  is
p=pos)|q,n,m)=)2⋅i
where the term i accounts for the 4-letter alphabet since we are considering nucleotides.
in the binary case, the term has to be dropped.
in the general case, for an alphabet of size σ, it would correspond to )i.
thus, the probability that the two sequence fragments are the same is
p=∑i=0i=nm/2)2⋅i
in our context, fragment <dig> and fragment <dig> refer to the overlapping region of two distinct reads in two adjacent amplicons.
for point , let's define the set a = { | ai ∈ n, a1+a2+...+aw+aw+ <dig> = nm/2}, as the space of frequency distributions where nm/ <dig> mutations can distribute either in w overlaps or in the remaining  parts, grouped in the additional variable w+ <dig>  given a generic element a =  ∈ a, each ai contains a certain number of mutations and the sum is the total number of mutations.
of note, the formula that gives the number of elements of the space a, as a function of n, m and w is
|a|==
and corresponds to the number of combinations with repetitions of w+ <dig> elements of nm/ <dig> class.
the probability that nm/ <dig> mutations distribute over the overlaps and the non-overlapping parts in a mode  is
p|q,n,m,w)=...
thus, for two vectors a =  ∈ a and b =  ∈ a, at least one overlapping region  will be consistent if, excluding the non-overlapping part, either  both a and b have the same element set to zero  or  both have one or more identical elements in the same overlap and within this overlap the mutations are in the same sites .
for case , let p be the probability  for a generic distribution ai = ) ∈ a, where at least one element aij is equal to zero.
define pij as the joint probability between two distributions, i.e.
pij = pp.
the sum of all joint probabilities ∑pij, where ∀ i ∃ j, k | aik = ajk =  <dig>  k ≠ w+ <dig>  yields the probability of a consistent overlap.
for case  we show how to calculate the probability associated to two distributions a and b, when they share at least one identical element, different from  <dig>  otherwise the case reduces to .
consider the two products
π1=...π2=...
and choose one of them, say π <dig>
if the two distributions have j identical elements  in the same sites , naming them  <dig>   <dig>  ..., j, we can write the following
α1=π1−π1⋅a1α2=α1−α1⋅a2αj=αj−1αj−1⋅aj
any αj can be interpreted as the number of combinations  that do not present the same elements in the same positions.
finally
p=π <dig>
is the probability for the two generic distributions a and b to have at least one identical overlap.
note that eq.
<dig> is valid under the constraint q > = nm/ <dig> and  > = nm/ <dig>  the sum of all joint probabilities ∑ p ij, where ∀ i ∃ j, k | aik = bjk, k ≠ w+ <dig>  aik ≠  <dig>  yields the probability of a consistent overlap.
eq.
<dig> is computationally intensive: for a small value of w and n it is possible to calculate it exactly, but for larger values , it is preferable to rely on numerical simulations.
algorithm: reconstruction of the quasispecies
from the definition in the above paragraph, we have a set of x variants  over a quasispecies, with a genome length of n and a mutation probability m. each variant has an associated prevalence p, p, ..., p, such that p+p+...+p =  <dig>  by using ngs machinery, we are able to sample  uniformly a large number of variant sequence fragments from the quasispecies population.
upon the definition of amplicons, we obtain w+ <dig> population samples, each one of length k, where k > n and an amplicon overlap of q sites.
previous studies investigated the probability of covering all bases of a single genome by shotgun sampling  <cit>  and the probability of covering all bases of different variants in a quasispecies  <cit> .
nowadays, ngs machineries are able to cover with high support a quasispecies of genomes of a few kilobases length.
we define the multinomial distribution ci = , i =  <dig> ... , where the generic element cij contains the number of identical reads  found in the amplicon i; thus, we have w+ <dig> available distributions.
since x is unknown, we assume initially that x is the maximum number of distinct reads that can be found in one amplicon, and we order the cik decreasingly, assuming that, given two distributions ci and cj, each cik and cjk correspond to a sample from the same variant.
based on the samples and the corresponding distributions, we aim to reconstruct the genomes associated to the unknown variants and their number, which eventually may be different from the initial value of x.
the objectives may be easier to reach if all the amplicons were designed such that the variants were different in all the overlaps, if the number of reads sequenced and covering the amplicons was sufficiently large and if the reads were completely error-free.
the problem becomes more challenging in presence of ambiguous overlaps , non-uniform or biased sampling, and uncorrected read errors.
for the latter  scenario, we design a set of algorithms in order to reconstruct a consistent set of variants that explains the ci distributions.
in the trivial case of a unique amplicon over the whole genome length, for a sufficiently large sample size, we may estimate variant probabilities from the distribution c <dig> as p = c1i/∑jc1j.
with multiple amplicons, depending on n, m, q, k and sample size, the distributions ci vary: if the hypothesis of error-free reads was fulfilled, the equations of the previous paragraph permit to calculate some confidence bounds.
in the real case, we expect that the multinomial distributions calculated for the amplicons are related, but we have to account for the uncertainty coming from the sampling process, cases of ambiguous overlaps and uncorrected read errors.
having a set of c <dig>  ..., cw+ <dig> distributions, we may be interested to find which is the most probable distribution under a given set of parameters or a model, i.e.
which distribution explains better the entire data.
this would be useful when applying a reconstruction algorithm, as explained in the next paragraphs.
if the probability of an event x dependent on parameter set θ  is written p, then the likelihood of the parameters given the data is l. in our case, θ corresponds to one of the cis and x is the set of remaining distributions x = {cj | j =  <dig> .. w+ <dig>  j ≠ i}.
we aim to find i such that l is the maximum.
however, since the derivation of l may be difficult, we use a minimum chi-square criterion  <cit> .
for each ci, i =  <dig> .. w+ <dig>  calculate and sum the chi-square statistic associated with all other cjs, and pick up the index i for which the sum of chi-square statistics is the minimum.
we may exclude as candidate model any ci for which |ci| <max{ |cj| j =  <dig> ..w+ <dig> }.
we define now a procedure that reconstructs a set of candidate variants of the quasispecies: the procedure takes into account both read distributions over the amplicons and calculation of consistent overlaps.
the algorithm is as follows:
<dig>  construct a matrix m = , i =  <dig> .. x, j =  <dig> .. w+ <dig>  where the columns represent the absolute frequency  distributions of distinct reads in the amplicons and each row contains distinct read representatives with their associated frequencies.
thus, the generic element mij is the number of distinct reads in amplicon j that correspond to a hypothetical variant i. each column of the matrix is ordered decreasingly.
since x is estimated as the maximum number of distinct reads found considering each amplicon, in amplicons where the number of distinct reads is less than x, missing values are all set to  <dig>
<dig>  choose a guide distribution among the amplicon distributions , say the one corresponding to amplicon g∈ { <dig>   <dig>  ..., w+1}.
<dig>  for each mgj ∈ m, j =  <dig> .. x, check iteratively if mgj is consistent with any other mik, i ≠ g, k =  <dig> .. x. if there is more than one consistent overlap, choose the index k whose absolute difference with the actual j is the lowest .
<dig> .
when a consistent set of distinct reads is obtained, i.e.
one variant is reconstructed with corresponding read-amplicon indices {ĵ <dig>  ..., ĵ}, subtract the number of distinct reads corresponding to the mgĵ value from the other mjĵ elements and update them in m. if some of the subtractions lead to negative values, set them to zero.
<dig>  if there is not a column of m with all zero elements  or if one variant has been constructed or the scan through amplicon distributions has not ended, go to point  <dig>
<dig>  output the variants reconstructed.
in the beginning, the algorithm counts all distinct reads for each amplicon.
distinct read representatives are ordered decreasingly by their frequency, creating w+ <dig> multinomial distributions of size x, which is the maximum number of distinct reads seen considering each amplicon.
if less than x distinct reads are found in one amplicon, the remaining elements of the corresponding multinomial distribution are set to zero-frequency.
a guide multinomial distribution is chosen either at random or by the minimum chi-square criterion.
the first read representative of the guide distribution corresponding to the count mg, <dig> is compared with the read representatives of an adjacent amplicon , starting from the first read  checking if there is a consistent overlap between the two read representatives.
if a consistent overlap is found, then a partial variant is reconstructed  and the step is repeated on another adjacent amplicon , until the whole set of amplicons is analysed.
if a consistent overlap between two reads is not found, for instance between read representatives corresponding to mg, <dig> and mg+ <dig> , then the procedure checks for a consistent overlap between the current read at position mg, <dig> and the next read in the same adjacent amplicon g+ <dig>  which is mg+ <dig>  in this case.
every time that a variant is reconstructed spanning all the amplicons, the algorithm subtracts the frequency count of the current read in the guide distribution  from the counts of all reads in the other amplicons that concurred to the variant reconstruction.
read frequencies that go to zero or below zero are considered exhausted and are not further evaluated.
negative values might appear due to variations generated by the ngs in the total read counts across the amplicons.
consider the trivial example of a quasispecies composed by a unique variant and two amplicons , where in the first amplicon the total read count is  <dig> and in the second is  <dig>  if the guide distribution corresponds to the second amplicon, the subtraction leads to - <dig>  the algorithm stops when all reads have been examined or if one amplicon distribution has all zero-frequencies.
a step-by-step example of the reconstruction algorithm is given in the additional file  <dig>  the computational complexity of the algorithm is o, which grows exponentially with the number of amplicons.
clearly, it is desirable to have a limited number of overlaps, e.g.
of amplicons, in order to decrease the computational burden.
note that initially the algorithm assumes that the number of variants in the quasispecies  is given by the maximum number of distinct reads observed across all amplicons , but the final number or reconstructed variants can be different, depending on the frequency distribution and overlap consistence.
indeed, in  <cit>  it was shown that in some cases the number of variants is higher than the number of distinct reads.
using our algorithm, the number of variants would be exactly x if the multinomial distributions mi were allowing for just one consistent overlap between elements in the same row ) and if the frequency subtraction was always exhausting all the mij.
in order to evaluate the effectiveness of the reconstruction algorithm and of the guide distribution choice policy, we designed and executed multiple simulation experiments over fixed parameters , varying mutation and sample size.
functions for the goodness of fit were  the prevalence of variants reconstructed correctly,  the number of false in-silico recombinants, and  number of reconstructed variants over the set of full consistent paths.
we obtained distributions of these loss functions executing multiple simulation runs and compared them via parametric test statistics.
testing: combinatorial analysis
in the methods section we derived two main formulae that provide theoretical bounds for the probabilities of consistent overlaps.
in particular, eq.
<dig> describes the probability that one overlap is consistent, given the genome length, the number of amplicons and the overlapping region size.
more generally, eq.
<dig> calculates the probability that at least one overlap is consistent.
by fixing n, q and w as above, and by simulating eq.
<dig> with  <dig> million of iterations, we calculated that the probabilities of at least one consistent overlap for m/ <dig> = { <dig> %, 1%,  <dig> %, 2%,  <dig> %,  <dig> %} are p = { <dig> ,  <dig> ,  <dig> ,  <dig> ,  <dig> ,  <dig> }, respectively.
for instance, at an m/ <dig> rate of  <dig> % there is 40% of chance that at least one overlap is consistent.
this gives a description of how much it could be difficult to reconstruct exact variants when the diversity is low.
testing: reconstruction algorithm on simulated data
the reconstruction algorithm, along with the investigation of guide distribution choice, was evaluated using simulated data.
a quasispecies composed by x =  <dig> variants was designed, considering a genome of length n =  <dig>  over a 4-letter alphabet.
variant prevalence was the following: p = 2%; p = 4%; p = 5%; p = 7%; p = 9%; p = 11%; p = 13%; p = 14%; p = 17%; p = 18%.
the amplicons consisted of w+ <dig> =  <dig> regions, each one of length k =  <dig> and overlap q =  <dig>  different uniform mutation probabilities were considered, specifically: m/ <dig> = { <dig> %, 1%,  <dig> %, 2%,  <dig> %,  <dig> %}.
we tested either a random guide distribution or a guide distribution chosen by maximum likelihood.
we executed ngs simulations for a sample of  <dig>  reads.
the reads were error-free and uniformly distributed along the genome.
figure  <dig> reports simulation results over a set of  <dig> independent runs, shuffling the mutational sites.
with  <dig>  read samples, the method reconstructed on average exactly the  <dig> variants at values of m/ <dig> around 2%.
by decreasing m/ <dig> to 1%, on average more than a half of the original variants were reconstructed, but there was higher prevalence of in-silico recombinants.
as it concerns the sole reconstruction of correct variants, comparison of the usage of a random guide distribution vs. one based on maximum likelihood did not yield significant differences.
however, the maximum likelihood policy reconstructed, on average, a lower number of in-silico recombinants.
note that, since the multinomial distributions are ordered decreasingly, we expect to reconstruct variants from the most prevalent to the less prevalent.
another way to evaluate the robustness of the algorithm is by looking at the number of potential variants  as a function of the per-site mutation probability, as depicted in figure  <dig>
in our simulation study, for an m/ <dig> = 1% on average there would be ≈ <dig>  paths, i.e.
candidate variants.
our algorithm on average chose 5- <dig> out of  <dig> correctly and did not reconstruct more than  <dig> variants.
by increasing m to  <dig> %, the number of paths would be still fairly high, i.e.
≈ <dig> : in this case the algorithm on average reconstructed > 80% variants correctly and the total number did not exceed  <dig>
using the same sets of simulated data , we compared our algorithm with the shorah  program; however it should be noted that shorah has not been designed to work on amplicons, but rather on shotgun modality.
although the current release provides the possibility to vary sliding window and the step size parameters, we could not reproduce exactly our amplicon settings, since the sliding window procedure is designed to cover multiple times each base over a uniform  fragment sequencing.
however, the average number of total reconstructions yielded by shorah was comparable to our method, across different runs and m values.
on average, at m/ <dig> =  <dig> %, the percentage of correct reconstruction was > 70% over different runs.
figure  <dig> depicts a phylogenetic tree constructed by pooling the original quasispecies together with the reconstructed variants from shorah and our method, over a single simulation run at m/ <dig> =  <dig> %.
seven shorah variants clustered significantly  with the original variants, over a total number of  <dig> reconstructions.
interestingly, our method reconstructed  <dig> variants .
a figure indicating recombination patterns is available in the additional file  <dig>
testing: reconstruction algorithm on real data
the algorithm was also applied to real ngs data.
we designed an experiment amplifying hbv sequences from  <dig> infected patient using a roche  <dig> gsflx titanium machine based on the amplicon sequencing modality.
patients' samples were processed in the same plate using barcodes  <cit> .
three amplicons were defined with specific primers, each one with a length of { <dig>   <dig>  394} bases and with two overlaps of length { <dig>  109}.
see the additional file  <dig> for experiment details.
one patient was infected with a genotype a virus  and four with a genotype d .
overall, average  read length was  <dig>   bases.
the same hbv reference sequence  was used for read alignment and individual genome re-sequencing of each patient.
we selected only reads that were significantly aligned with the reference .
three-percent of reads was discarded.
the average diversity m/ <dig> was  <dig> %.
according to the amplicon coverage, we reduced the amplicon lengths to { <dig>   <dig>  290} and overlaps to { <dig>  90} bases.
finally, we selected those reads that covered entirely one amplicon region with a gap percentage below 5%.
for each amplicon, exactly  <dig>  reads for patients were retained, selecting them at random, without replacement, from the previous set of filtered sequences.
all reads from the different patients were pooled together in a unique file, thus obtaining  <dig>  reads per patient and  <dig>  reads in total, with a fixed read/amplicon/patient ratio.
we were able to reconstruct virus consensus genomes from each individual using the read alignment, but we did not know a-priori the composition of the viral quasispecies of the patients.
however, for each read we knew the corresponding patient.
the purpose of this experiment was to see if the reconstruction algorithms were able to reconstruct a swarm of variants closely related to each patient's virus consensus genome, without mixing the population and without creating incorrect, populations.
both shorah  and the reconstruction algorithm were run on this joined data set, considering - as a simple error correction procedure - only reads with a frequency > =  <dig>  requiring that at least one read was seen in reverse-strand and another in forward-strand.
shorah identified  <dig> distinct variants, with a median  prevalence of  <dig>  .
the number of shorah variants with prevalence above the 95th percentile of the overall distribution was  <dig>  our reconstruction algorithm reconstructed  <dig> unique variants.
we executed a phylogenetic analysis pooling together the set of reconstructed genomes, the  <dig> shorah variants, the  <dig> unique variants obtained with our algorithms, and two additional outgroups .
the phylogenetic tree was estimated via a neighbour-joining method and the logdet distance, assessing node support with  <dig>  bootstrap runs.
all the variants reconstructed with our algorithm clustered with the corresponding patients, and in four cases out of five the phylogenetic clusters had a support > 75%.
the same held when looking at the shorah variants, although a considerable number of variants clustered apart from the patients.
figure  <dig> depicts the phylogenetic tree.
of note, in patient # <dig>  two variants reconstructed with our algorithm were indeed recombinants between patient # <dig> and patient # <dig>
discussion
in this paper we addressed the problem of quasispecies determination and variant reconstruction by using ngs machinery.
original assumptions were:  to have a uniform random sampling of the population,  a reference genome,  a unique, error-free, alignment of each read against the reference, and  a sliding window partition of the reference into a set of amplicons.
we derived first a set of formulae in order to analyse the probability of consistent overlaps given two sequence fragments over a set of amplicons.
we showed that many factors, including diversity and overlap length, can affect the chance to detect spurious consistent overlaps.
we introduced then the concept of multinomial distribution as a model for the classification of distinct reads and relative prevalence within amplicons.
upon this, we designed a greedy algorithm that reconstructs a set of paths through the whole set of amplicons , coupling elements of different multinomial distributions, and trying to minimise the chance to reconstruct in-silico recombinants.
the algorithm is based on a "guide distribution" policy that can be either random or based on maximum-likelihood.
with a practical example , we highlighted the reasons for which any quasispecies reconstruction procedure should consider read frequencies in order to avoid the estimation of false variants.
in fact, our reconstruction algorithm tends to select variants not only looking at the consistent overlaps , but also considering reads that have similar frequencies across the various amplicons.
simulation results proved that, exploring a fixed set of parameters, our method was able to select a compact and correct set of variants even at low diversities.
at m/ <dig> =  <dig> %, the algorithm was able to reconstruct on average > 80% of correct variants, with an estimated number of variants close to the real value .
we also executed a test on real ngs data, prone to contain sequencing errors, considering a mixed population of hbv-infected patients with a low average diversity.
in this case our algorithm was able to distinguish variants corresponding to different patients, with a minimal evidence of in-silico recombination.
in addition, the algorithm did not generate variants that could be different from the sequenced population.
in its current definition, though, our model possesses several limitations.
first, a reference genome and a sliding window amplicon partition are needed: thus, the variant reconstruction method is suitable only for quasispecies for which there is at least one available genome.
however, de-novo quasispecies determination can be easily achieved by pre-processing ngs data with existing whole genome assembly software and obtaining a usable reference genome.
another important critical point is that we assume a uniform distribution of diversity along the genome, which is an ideal hypothesis.
one solution may be the design of amplicons and overlaps of different lengths.
the reconstruction algorithm works even with size-variable amplicons and overlaps, but the formulae introduced in the preliminary combinatorial analysis should be derived again, taking into account these modifications.
another issue is the assumption of a unique mapping of each read with respect to the reference genome, which may be not always fulfilled when in presence of long repeats .
however, this problem does not affect the reconstruction algorithm once the read mapping is given along with the sliding window amplicon setup.
several approaches have been proposed in literature  <cit>  and may be applied to ngs.
as future refinements of the reconstruction algorithm we foresee the estimation of exact variant prevalence, since currently we report variants just in decreasing prevalence order: one idea is to calculate average and standard errors of distinct read frequencies from the various multinomial distributions joined during the reconstruction phase; another approach could be to estimate the prevalence a-posteriori, using expectation-maximisation as it was done in  <cit> .
a broader perspective would be to relax the need for a reference genome and to estimate the quasispecies independently from the read mapping and the amplicon definition, but under these general settings the theoretical results here obtained would be hardly reusable.
CONCLUSIONS
the presented combinatorial analysis and the reconstruction algorithm may be a fundamental step towards the characterisation of quasispecies using ngs.
immediate applications can be found in analysing genomes of infectious pathogens, like viruses, currently targeted by inhibitors and developing resistance.
the investigation of in-depth, intra-host, viral resistance evolutionary mechanisms and interactions among mutations is crucial in order to design effective treatment strategies, even at early disease stages, and to maximise further treatment options.
authors' contributions
all authors read and approved the final manuscript.
lp carried out combinatorial analysis and algorithmic specifications.
mp performed simulations, evaluation on real data, phylogenetic analysis, statistical comparisons with other methods, and manuscript writing.
gu revised the mathematical methods.
ab provided software installations and runs.
ia, gr, dv, and mcs provided expertise in next-generation machinery and amplicon sequencing, and performed laboratory experiments.
mrc leaded the research and supervised authors' contributions.
supplementary material
additional file 1
this file includes additional details on:  average population diversity estimation;  step-by-step example of the quasispecies reconstruction algorithm;  information on the sample preparation and laboratory protocols for the experiment on roche glsflx platform;  figure of recombination patterns for a reconstructed quasispecies given a simulation experiment.
click here for file
acknowledgements
this work has been partly supported by the dynanets eu fet open project , and by grants from italian ministry of health .
we would like to thank dr. rebecca gray  for revising the english language.

 next-generation sequencing techniques allow for a high-throughput dna sequencing, producing from thousands to billions of sequence fragments composed of tens to hundreds of nucleotide bases. ngs has the potential to replace sanger sequencing for many applications, including de-novo sequencing, re-sequencing, meta-genomics and intra-host characterisation of infectious pathogens. de-novo sequencing implies a genome assembly problem, which is the reconstruction of a unique genome from a set of sequence fragments. several methods and software for genome assembly have been developed. these methods were designed initially for sanger sequencing, and have been revised for ngs technology, given different error rates among ngs machineries. re-sequencing conjugates with the problem of single nucleotide polymorphisms discovery. recent studies characterised snps or drug-induced mutations with ngs, considering the human immunodeficiency virus and the hepatitis b virus. more specifically, re-sequencing can be useful for the characterisation of variants within a quasispecies harbouring an infected host. many rna viruses are present in a carrier as a swarm of highly genetically related variants, i.e. a quasispecies, due to the error prone characteristics of the viral polymerases and high viral replication rates. this intra-host variability represents a substrate for the selective pressure exerted by the immune system of the host or by drug exposure, which leads to the continuous evolution of viruses. quasispecies reconstruction would allow detailed description of the composition of individual viral genomes, genetic linkage and evolutionary history. for example, in hiv or hbv infection, the development of drug resistance is a major problem and the early diagnosis of drug-resistant variant selection might help in designing targeted therapeutic interventions. here, we addressed the problem of reconstructing a viral quasispecies from a ngs data set, which is a relatively new topic that is not widely investigated in literature. we aimed to reconstruct all coexistent individual variants within a population, along with their prevalence, rather than a reconstruction of a single or predominant genome. current assembly software is not designed to accomplish this task, nor to deal easily with the reconstruction of highly variable genomes. the huge coverage and base pair output provided by ngs enables the design of experiments to investigate and validate theoretical methods for quasispecies reconstruction. at present, only a few methodological papers have been published presenting new algorithms for quasispecies reconstruction that are able to infer both genomes of population variants and their prevalence. in, the authors proposed an algorithm based on a generative model of the sequencing process and a tailored probabilistic inference and learning procedure for model fitting. in, a set of methodologies was proposed both for error correction and inference about the structure of a population from a set of short sequence reads as obtained from ngs. the authors assumed a known mapping of reads to a reference genome, defined a sliding window over the reference genome and associated each aligned read to one or more windows by trimming the reads accordingly to the windows' bounds. sequencing errors were corrected by locally clustering reads in the windows. a set of single variants of the quasispecies was obtained by constructing an overlap graph of non-redundant, error-free, aligned reads, and by calculating a minimal coverage set of paths over the graph. the frequency estimation was done with an expectation maximisation algorithm and was proven to be more efficient than a naïve procedure based on uniform read sampling. one drawback of this methodology is that the variant reconstruction phase did not account for the relations among frequencies of distinct variants that were overlapping consistently across the sliding windows: this may lead potentially to selection of in-silico recombinants and the procedure of haplotype frequency may be biased from the exclusion of real paths. after the paper, free software was released, named shorah. in, a scalable assembling method for quasispecies based on a novel network flow formulation was presented, applied efficiently for the assembly of hepatitis c virus. in, a refinement of the original procedures presented in was given, substituting k-means clustering with a dirichlet process mixture for locally inferring haplotypes and correcting reads. in this work, a set of formulae for combinatorial analysis of quasispecies genome fragments sampled by ngs was derived, and a new greedy algorithm for quasispecies reconstruction was introduced. the formulae derivation provided some theoretical bounds explaining the difficulty in reconstructing a set of individual variants of the quasispecies, by conditioning on several parameters, such as the genome length, the fragment size, or the overlap length between two sampled fragments. the reconstruction algorithm was based on combinations of multinomial distributions and was designed to minimise the reconstruction of in-silico recombinants. unlike previous approaches, the algorithm selects and reconstructs variants not only by coupling reads that have consistent overlaps, but also considering reads that have similar frequencies across the various amplicons. for our combinatorial analysis, we assumed that the problem of re-sequencing, including reference alignment and error correction, is solved. in other words, a set of error-free reads is available, aligned univocally to a reference or consensus sequence. such a reference may either have been directly reconstructed using assembly software or selected from literature. the assumption for the unique mapping of a read against the reference may not be always fulfilled when in presence of short reads and genomes with long repeats. however, this problem can be negligible when considering coding regions of highly-variable viral pathogens targeted by inhibitors, with a few regulatory regions where the repeats usually are located. the sequence quality may be a major concern for reconstruction algorithms and this is often an experiment- and machine-dependent problem: procedures for alignment and error corrections have been investigated elsewhere, with different methodologies, along with protocols for sample preparation. another critical point with ngs is the presence of contaminants that must be detected and excluded. the problem can be solved easily when the contamination is from different organisms, with a test statistic on read/reference alignment scores during re-sequencing. it is harder when the ngs experiment comprises a mixture of closely related organisms, for instance when samples of patients infected with the same virus are put together in one ngs experiment. as a second assumption, our algorithm required a non empty set of overlapping regions, called amplicons, which cover the reference genome. each read has to be assigned to one of these amplicons. roche gsflx technology has a double working modality that allows both for shotgun sequencing and for amplicon sequencing with specific primer design, although the latter option is generally more expensive. with this technology, amplicons can be defined a priori. in contrast, if shotgun sequencing is performed, additional data elaboration has to be made in order to define a set of amplicons: one solution is to define amplicons via sliding windows over the genome and cluster reads accordingly to their mapping region. the proposed reconstruction algorithm was applied to 1) simulated and error-free data; and 2) then empirical sequence data derived from blood samples from hbv-infected patients processed via the roche gsflx titanium machine. this second dataset was designed to assess the performance of quasispecies reconstruction in presence of sequencing errors.  algorithm: ngs data processing and amplicon definition this work analysed a re-sequencing experiment of a viral quasispecies carried out using ngs machinery. since currently the maximum read length of a ngs experiment does not exceed a few hundred of bases, we were interested in genome regions of quasispecies whose length is much larger than the average read length, i.e. when it is not possible that a read spans entirely the genome of interest. we required then that a reference genome is available and that reads are significantly aligned against this reference genome. this can be achieved by aligning each read in forward- or reverse-strand against the reference genome, using the smith-waterman-gotoh local alignment, which is an exact algorithm, and keeping the highest alignment score. reads then can be filtered by excluding those that do not show a significant alignment score, as compared to a score distribution obtained from quasi-random sequences aligned to the reference genome, as described in. we assumed also that sequencing errors were corrected. the condition of error-free reads was required only for the combinatorial analysis, whilst the quasispecies reconstruction algorithm can be applied also to noisy data. given a reference genome g and a read alignment over g, we define then a sliding window partition of g composed of w+ windows, that we call amplicons. these amplicons cover the entire genome and two adjacent amplicons have a partial overlap, for a total of w overlaps. amplicons do not need to be necessarily of the same length. clearly, each amplicon size has to be smaller than the average read size, so that a read can span an amplicon entirely. as stated in the introduction, amplicons can be designed a priori if roche gsflx technology is used, or determined with a fixed sliding window approach from any shotgun sequencing. after defining the amplicons, reads that spanned entirely an amplicon were trimmed so that their start/end positions corresponded exactly to the amplicon start/end. consequently, all the reads in one amplicon had the same length. reads that span more than one amplicon entirely are considered multiple times, whilst those that do not cover exactly at least one amplicon were discarded. figure illustrates with an example a three-amplicon design over a reference genome, with corresponding read assignment and trimming. algorithm: combinatorial analysis of ngs over a quasispecies consider a number of x variants determined by their genomic sequences of length n, either nucleotidic or amino-acid. the number x of variants is unknown, along with the prevalence of each variant. we assume that the quasispecies is stable over a fixed set of variants, in a mutation-selection balance. in other words, the number of distinct viral variants in the quasispecies that a carrier harbours is x, although each variant can be present with different prevalence, presumably due to different viral fitness or immune response or drug-induced selection. after a multiple sequence alignment, a consensus sequence can be generated or one variant can be used as a reference. we define a point difference between two aligned sequences as the presence of two different nucleotides in one position of the alignment. the pairwise difference between two variants d is the number of point differences divided by the genome length, i.e. d=dij=∑k=1nn the average pairwise difference of a variant si with respect to all the others is defined as d=di=/, for j=..x. finally, the average overall pairwise difference among a set of aligned variants is then davg=∑i=1x−1∑j=i+1xdijx let us define now a reference sequence sref as the variant with the lowest average pairwise difference as compared to all other variants, i.e. dref = d = min{d}, for i =..x. if we define the diversity - that we regard as a probability of mutation - of the quasispecies as m = davg, then we can also approximate m by doubling the value of dref, i.e. m/ is the average pairwise difference of our reference variant with respect to any other variant. indeed, this approximation is correct when no identical base changes happen in the same alignment position of two variants as compared to the reference, and this is dependent on the genome length and the mutation probability. the approximation gets better either if the genome length increases or the mutation rate decreases. more details on the efficacy of this approximation are given in additional file we previously defined a joined set of w+ amplicons, which induces an overlapping ordered coverage over the quasispecies genome space. assume that each amplicon has a fixed length of k bases and overlaps with its neighbour over q bases for w times. we assume also that, given three adjacent amplicons and two corresponding overlaps, these latter do not share any position in common, i.e. there are no overlapping overlaps. since there are no nested amplicons, we can define an amplicon identification number as its ordinal position with respect to the reference genome and the other amplicons. thus, amplicon starts at position and ends at position k, amplicon starts at k-q+ and ends at 2k-q, et cetera. the overlaps are clearly at the end of each amplicon and at the beginning of the adjacent one. each amplicon is associated with a set of reads, or sequence fragments, sampled uniformly from the quasispecies. these reads, by definition, are significantly aligned to the reference genome, error-free and span exactly an amplicon region, being trimmed at its ends. thus, after sampling, each count of distinct reads across each amplicon cannot exceed the number of variants x. given two reads associated to two adjacent amplicons, we say that their overlapping region is consistent if the two reads share in that region the same characters. we aim to calculate the probability that the overlapping region of two adjacent reads is consistent and at least one overlapping region across the amplicons is consistent. for point, we first define the probability that i mutations are present in a sequence fragment of q length over a genome of length n as p= note that the diversity m/ here is multiplied by n, obtaining the expected number of changes. the probability that two random regions of q length over a genome of length n have both i mutations is p=) where the terms of eq. are the square root of the terms of eq. the probability that these random regions share exactly i mutations at the same positions is p=pos)|q,n,m)=)2⋅i where the term i accounts for the 4-letter alphabet since we are considering nucleotides. in the binary case, the term has to be dropped. in the general case, for an alphabet of size σ, it would correspond to )i. thus, the probability that the two sequence fragments are the same is p=∑i=0i=nm/2)2⋅i in our context, fragment and fragment refer to the overlapping region of two distinct reads in two adjacent amplicons. for point, let's define the set a = { | ai ∈ n, a1+a2+...+aw+aw+ = nm/2}, as the space of frequency distributions where nm/ mutations can distribute either in w overlaps or in the remaining parts, grouped in the additional variable w+ given a generic element a = ∈ a, each ai contains a certain number of mutations and the sum is the total number of mutations. of note, the formula that gives the number of elements of the space a, as a function of n, m and w is |a|== and corresponds to the number of combinations with repetitions of w+ elements of nm/ class. the probability that nm/ mutations distribute over the overlaps and the non-overlapping parts in a mode is p|q,n,m,w)=... thus, for two vectors a = ∈ a and b = ∈ a, at least one overlapping region will be consistent if, excluding the non-overlapping part, either both a and b have the same element set to zero or both have one or more identical elements in the same overlap and within this overlap the mutations are in the same sites. for case, let p be the probability for a generic distribution ai = ) ∈ a, where at least one element aij is equal to zero. define pij as the joint probability between two distributions, i.e. pij = pp. the sum of all joint probabilities ∑pij, where ∀ i ∃ j, k | aik = ajk = k ≠ w+ yields the probability of a consistent overlap. for case we show how to calculate the probability associated to two distributions a and b, when they share at least one identical element, different from otherwise the case reduces to. consider the two products π1=...π2=... and choose one of them, say π if the two distributions have j identical elements in the same sites, naming them..., j, we can write the following α1=π1−π1⋅a1α2=α1−α1⋅a2αj=αj−1αj−1⋅aj any αj can be interpreted as the number of combinations that do not present the same elements in the same positions. finally p=π is the probability for the two generic distributions a and b to have at least one identical overlap. note that eq. is valid under the constraint q > = nm/ and > = nm/ the sum of all joint probabilities ∑ p ij, where ∀ i ∃ j, k | aik = bjk, k ≠ w+ aik ≠ yields the probability of a consistent overlap. eq. is computationally intensive: for a small value of w and n it is possible to calculate it exactly, but for larger values, it is preferable to rely on numerical simulations. algorithm: reconstruction of the quasispecies from the definition in the above paragraph, we have a set of x variants over a quasispecies, with a genome length of n and a mutation probability m. each variant has an associated prevalence p, p,..., p, such that p+p+...+p = by using ngs machinery, we are able to sample uniformly a large number of variant sequence fragments from the quasispecies population. upon the definition of amplicons, we obtain w+ population samples, each one of length k, where k > n and an amplicon overlap of q sites. previous studies investigated the probability of covering all bases of a single genome by shotgun sampling and the probability of covering all bases of different variants in a quasispecies. nowadays, ngs machineries are able to cover with high support a quasispecies of genomes of a few kilobases length. we define the multinomial distribution ci =, i =..., where the generic element cij contains the number of identical reads found in the amplicon i; thus, we have w+ available distributions. since x is unknown, we assume initially that x is the maximum number of distinct reads that can be found in one amplicon, and we order the cik decreasingly, assuming that, given two distributions ci and cj, each cik and cjk correspond to a sample from the same variant. based on the samples and the corresponding distributions, we aim to reconstruct the genomes associated to the unknown variants and their number, which eventually may be different from the initial value of x. the objectives may be easier to reach if all the amplicons were designed such that the variants were different in all the overlaps, if the number of reads sequenced and covering the amplicons was sufficiently large and if the reads were completely error-free. the problem becomes more challenging in presence of ambiguous overlaps, non-uniform or biased sampling, and uncorrected read errors. for the latter scenario, we design a set of algorithms in order to reconstruct a consistent set of variants that explains the ci distributions. in the trivial case of a unique amplicon over the whole genome length, for a sufficiently large sample size, we may estimate variant probabilities from the distribution c as p = c1i/∑jc1j. with multiple amplicons, depending on n, m, q, k and sample size, the distributions ci vary: if the hypothesis of error-free reads was fulfilled, the equations of the previous paragraph permit to calculate some confidence bounds. in the real case, we expect that the multinomial distributions calculated for the amplicons are related, but we have to account for the uncertainty coming from the sampling process, cases of ambiguous overlaps and uncorrected read errors. having a set of c..., cw+ distributions, we may be interested to find which is the most probable distribution under a given set of parameters or a model, i.e. which distribution explains better the entire data. this would be useful when applying a reconstruction algorithm, as explained in the next paragraphs. if the probability of an event x dependent on parameter set θ is written p, then the likelihood of the parameters given the data is l. in our case, θ corresponds to one of the cis and x is the set of remaining distributions x = {cj | j =.. w+ j ≠ i}. we aim to find i such that l is the maximum. however, since the derivation of l may be difficult, we use a minimum chi-square criterion. for each ci, i =.. w+ calculate and sum the chi-square statistic associated with all other cjs, and pick up the index i for which the sum of chi-square statistics is the minimum. we may exclude as candidate model any ci for which |ci| <max{ |cj| j =..w+ }. we define now a procedure that reconstructs a set of candidate variants of the quasispecies: the procedure takes into account both read distributions over the amplicons and calculation of consistent overlaps. the algorithm is as follows: construct a matrix m =, i =.. x, j =.. w+ where the columns represent the absolute frequency distributions of distinct reads in the amplicons and each row contains distinct read representatives with their associated frequencies. thus, the generic element mij is the number of distinct reads in amplicon j that correspond to a hypothetical variant i. each column of the matrix is ordered decreasingly. since x is estimated as the maximum number of distinct reads found considering each amplicon, in amplicons where the number of distinct reads is less than x, missing values are all set to choose a guide distribution among the amplicon distributions, say the one corresponding to amplicon g∈ {..., w+1}. for each mgj ∈ m, j =.. x, check iteratively if mgj is consistent with any other mik, i ≠ g, k =.. x. if there is more than one consistent overlap, choose the index k whose absolute difference with the actual j is the lowest.. when a consistent set of distinct reads is obtained, i.e. one variant is reconstructed with corresponding read-amplicon indices {ĵ..., ĵ}, subtract the number of distinct reads corresponding to the mgĵ value from the other mjĵ elements and update them in m. if some of the subtractions lead to negative values, set them to zero. if there is not a column of m with all zero elements or if one variant has been constructed or the scan through amplicon distributions has not ended, go to point output the variants reconstructed. in the beginning, the algorithm counts all distinct reads for each amplicon. distinct read representatives are ordered decreasingly by their frequency, creating w+ multinomial distributions of size x, which is the maximum number of distinct reads seen considering each amplicon. if less than x distinct reads are found in one amplicon, the remaining elements of the corresponding multinomial distribution are set to zero-frequency. a guide multinomial distribution is chosen either at random or by the minimum chi-square criterion. the first read representative of the guide distribution corresponding to the count mg, is compared with the read representatives of an adjacent amplicon, starting from the first read checking if there is a consistent overlap between the two read representatives. if a consistent overlap is found, then a partial variant is reconstructed and the step is repeated on another adjacent amplicon, until the whole set of amplicons is analysed. if a consistent overlap between two reads is not found, for instance between read representatives corresponding to mg, and mg+, then the procedure checks for a consistent overlap between the current read at position mg, and the next read in the same adjacent amplicon g+ which is mg+ in this case. every time that a variant is reconstructed spanning all the amplicons, the algorithm subtracts the frequency count of the current read in the guide distribution from the counts of all reads in the other amplicons that concurred to the variant reconstruction. read frequencies that go to zero or below zero are considered exhausted and are not further evaluated. negative values might appear due to variations generated by the ngs in the total read counts across the amplicons. consider the trivial example of a quasispecies composed by a unique variant and two amplicons, where in the first amplicon the total read count is and in the second is if the guide distribution corresponds to the second amplicon, the subtraction leads to - the algorithm stops when all reads have been examined or if one amplicon distribution has all zero-frequencies. a step-by-step example of the reconstruction algorithm is given in the additional file the computational complexity of the algorithm is o, which grows exponentially with the number of amplicons. clearly, it is desirable to have a limited number of overlaps, e.g. of amplicons, in order to decrease the computational burden. note that initially the algorithm assumes that the number of variants in the quasispecies is given by the maximum number of distinct reads observed across all amplicons, but the final number or reconstructed variants can be different, depending on the frequency distribution and overlap consistence. indeed, in it was shown that in some cases the number of variants is higher than the number of distinct reads. using our algorithm, the number of variants would be exactly x if the multinomial distributions mi were allowing for just one consistent overlap between elements in the same row ) and if the frequency subtraction was always exhausting all the mij. in order to evaluate the effectiveness of the reconstruction algorithm and of the guide distribution choice policy, we designed and executed multiple simulation experiments over fixed parameters, varying mutation and sample size. functions for the goodness of fit were the prevalence of variants reconstructed correctly, the number of false in-silico recombinants, and number of reconstructed variants over the set of full consistent paths. we obtained distributions of these loss functions executing multiple simulation runs and compared them via parametric test statistics. testing: combinatorial analysis in the methods section we derived two main formulae that provide theoretical bounds for the probabilities of consistent overlaps. in particular, eq. describes the probability that one overlap is consistent, given the genome length, the number of amplicons and the overlapping region size. more generally, eq. calculates the probability that at least one overlap is consistent. by fixing n, q and w as above, and by simulating eq. with million of iterations, we calculated that the probabilities of at least one consistent overlap for m/ = { %, 1%, %, 2%, %, %} are p = {,,,,, }, respectively. for instance, at an m/ rate of % there is 40% of chance that at least one overlap is consistent. this gives a description of how much it could be difficult to reconstruct exact variants when the diversity is low. testing: reconstruction algorithm on simulated data the reconstruction algorithm, along with the investigation of guide distribution choice, was evaluated using simulated data. a quasispecies composed by x = variants was designed, considering a genome of length n = over a 4-letter alphabet. variant prevalence was the following: p = 2%; p = 4%; p = 5%; p = 7%; p = 9%; p = 11%; p = 13%; p = 14%; p = 17%; p = 18%. the amplicons consisted of w+ = regions, each one of length k = and overlap q = different uniform mutation probabilities were considered, specifically: m/ = { %, 1%, %, 2%, %, %}. we tested either a random guide distribution or a guide distribution chosen by maximum likelihood. we executed ngs simulations for a sample of reads. the reads were error-free and uniformly distributed along the genome. figure reports simulation results over a set of independent runs, shuffling the mutational sites. with read samples, the method reconstructed on average exactly the variants at values of m/ around 2%. by decreasing m/ to 1%, on average more than a half of the original variants were reconstructed, but there was higher prevalence of in-silico recombinants. as it concerns the sole reconstruction of correct variants, comparison of the usage of a random guide distribution vs. one based on maximum likelihood did not yield significant differences. however, the maximum likelihood policy reconstructed, on average, a lower number of in-silico recombinants. note that, since the multinomial distributions are ordered decreasingly, we expect to reconstruct variants from the most prevalent to the less prevalent. another way to evaluate the robustness of the algorithm is by looking at the number of potential variants as a function of the per-site mutation probability, as depicted in figure in our simulation study, for an m/ = 1% on average there would be ≈ paths, i.e. candidate variants. our algorithm on average chose 5- out of correctly and did not reconstruct more than variants. by increasing m to %, the number of paths would be still fairly high, i.e. ≈: in this case the algorithm on average reconstructed > 80% variants correctly and the total number did not exceed using the same sets of simulated data, we compared our algorithm with the shorah program; however it should be noted that shorah has not been designed to work on amplicons, but rather on shotgun modality. although the current release provides the possibility to vary sliding window and the step size parameters, we could not reproduce exactly our amplicon settings, since the sliding window procedure is designed to cover multiple times each base over a uniform fragment sequencing. however, the average number of total reconstructions yielded by shorah was comparable to our method, across different runs and m values. on average, at m/ = %, the percentage of correct reconstruction was > 70% over different runs. figure depicts a phylogenetic tree constructed by pooling the original quasispecies together with the reconstructed variants from shorah and our method, over a single simulation run at m/ = %. seven shorah variants clustered significantly with the original variants, over a total number of reconstructions. interestingly, our method reconstructed variants. a figure indicating recombination patterns is available in the additional file testing: reconstruction algorithm on real data the algorithm was also applied to real ngs data. we designed an experiment amplifying hbv sequences from infected patient using a roche gsflx titanium machine based on the amplicon sequencing modality. patients' samples were processed in the same plate using barcodes. three amplicons were defined with specific primers, each one with a length of { 394} bases and with two overlaps of length { 109}. see the additional file for experiment details. one patient was infected with a genotype a virus and four with a genotype d. overall, average read length was bases. the same hbv reference sequence was used for read alignment and individual genome re-sequencing of each patient. we selected only reads that were significantly aligned with the reference. three-percent of reads was discarded. the average diversity m/ was %. according to the amplicon coverage, we reduced the amplicon lengths to { 290} and overlaps to { 90} bases. finally, we selected those reads that covered entirely one amplicon region with a gap percentage below 5%. for each amplicon, exactly reads for patients were retained, selecting them at random, without replacement, from the previous set of filtered sequences. all reads from the different patients were pooled together in a unique file, thus obtaining reads per patient and reads in total, with a fixed read/amplicon/patient ratio. we were able to reconstruct virus consensus genomes from each individual using the read alignment, but we did not know a-priori the composition of the viral quasispecies of the patients. however, for each read we knew the corresponding patient. the purpose of this experiment was to see if the reconstruction algorithms were able to reconstruct a swarm of variants closely related to each patient's virus consensus genome, without mixing the population and without creating incorrect, populations. both shorah and the reconstruction algorithm were run on this joined data set, considering - as a simple error correction procedure - only reads with a frequency > = requiring that at least one read was seen in reverse-strand and another in forward-strand. shorah identified distinct variants, with a median prevalence of. the number of shorah variants with prevalence above the 95th percentile of the overall distribution was our reconstruction algorithm reconstructed unique variants. we executed a phylogenetic analysis pooling together the set of reconstructed genomes, the shorah variants, the unique variants obtained with our algorithms, and two additional outgroups. the phylogenetic tree was estimated via a neighbour-joining method and the logdet distance, assessing node support with bootstrap runs. all the variants reconstructed with our algorithm clustered with the corresponding patients, and in four cases out of five the phylogenetic clusters had a support > 75%. the same held when looking at the shorah variants, although a considerable number of variants clustered apart from the patients. figure depicts the phylogenetic tree. of note, in patient # two variants reconstructed with our algorithm were indeed recombinants between patient # and patient # discussion in this paper we addressed the problem of quasispecies determination and variant reconstruction by using ngs machinery. original assumptions were: to have a uniform random sampling of the population, a reference genome, a unique, error-free, alignment of each read against the reference, and a sliding window partition of the reference into a set of amplicons. we derived first a set of formulae in order to analyse the probability of consistent overlaps given two sequence fragments over a set of amplicons. we showed that many factors, including diversity and overlap length, can affect the chance to detect spurious consistent overlaps. we introduced then the concept of multinomial distribution as a model for the classification of distinct reads and relative prevalence within amplicons. upon this, we designed a greedy algorithm that reconstructs a set of paths through the whole set of amplicons, coupling elements of different multinomial distributions, and trying to minimise the chance to reconstruct in-silico recombinants. the algorithm is based on a "guide distribution" policy that can be either random or based on maximum-likelihood. with a practical example, we highlighted the reasons for which any quasispecies reconstruction procedure should consider read frequencies in order to avoid the estimation of false variants. in fact, our reconstruction algorithm tends to select variants not only looking at the consistent overlaps, but also considering reads that have similar frequencies across the various amplicons. simulation results proved that, exploring a fixed set of parameters, our method was able to select a compact and correct set of variants even at low diversities. at m/ = %, the algorithm was able to reconstruct on average > 80% of correct variants, with an estimated number of variants close to the real value. we also executed a test on real ngs data, prone to contain sequencing errors, considering a mixed population of hbv-infected patients with a low average diversity. in this case our algorithm was able to distinguish variants corresponding to different patients, with a minimal evidence of in-silico recombination. in addition, the algorithm did not generate variants that could be different from the sequenced population. in its current definition, though, our model possesses several limitations. first, a reference genome and a sliding window amplicon partition are needed: thus, the variant reconstruction method is suitable only for quasispecies for which there is at least one available genome. however, de-novo quasispecies determination can be easily achieved by pre-processing ngs data with existing whole genome assembly software and obtaining a usable reference genome. another important critical point is that we assume a uniform distribution of diversity along the genome, which is an ideal hypothesis. one solution may be the design of amplicons and overlaps of different lengths. the reconstruction algorithm works even with size-variable amplicons and overlaps, but the formulae introduced in the preliminary combinatorial analysis should be derived again, taking into account these modifications. another issue is the assumption of a unique mapping of each read with respect to the reference genome, which may be not always fulfilled when in presence of long repeats. however, this problem does not affect the reconstruction algorithm once the read mapping is given along with the sliding window amplicon setup. several approaches have been proposed in literature and may be applied to ngs. as future refinements of the reconstruction algorithm we foresee the estimation of exact variant prevalence, since currently we report variants just in decreasing prevalence order: one idea is to calculate average and standard errors of distinct read frequencies from the various multinomial distributions joined during the reconstruction phase; another approach could be to estimate the prevalence a-posteriori, using expectation-maximisation as it was done in. a broader perspective would be to relax the need for a reference genome and to estimate the quasispecies independently from the read mapping and the amplicon definition, but under these general settings the theoretical results here obtained would be hardly reusable.  the presented combinatorial analysis and the reconstruction algorithm may be a fundamental step towards the characterisation of quasispecies using ngs. immediate applications can be found in analysing genomes of infectious pathogens, like viruses, currently targeted by inhibitors and developing resistance. the investigation of in-depth, intra-host, viral resistance evolutionary mechanisms and interactions among mutations is crucial in order to design effective treatment strategies, even at early disease stages, and to maximise further treatment options. authors' contributions all authors read and approved the final manuscript. lp carried out combinatorial analysis and algorithmic specifications. mp performed simulations, evaluation on real data, phylogenetic analysis, statistical comparisons with other methods, and manuscript writing. gu revised the mathematical methods. ab provided software installations and runs. ia, gr, dv, and mcs provided expertise in next-generation machinery and amplicon sequencing, and performed laboratory experiments. mrc leaded the research and supervised authors' contributions. supplementary material additional file 1 this file includes additional details on: average population diversity estimation; step-by-step example of the quasispecies reconstruction algorithm; information on the sample preparation and laboratory protocols for the experiment on roche glsflx platform; figure of recombination patterns for a reconstructed quasispecies given a simulation experiment. click here for file acknowledgements this work has been partly supported by the dynanets eu fet open project, and by grants from italian ministry of health. we would like to thank dr. rebecca gray for revising the english language.
next-generation sequencing offers a unique opportunity for high-throughput genomics and has potential to replace sanger sequencing in many fields, including de-novo sequencing, re-sequencing, meta-genomics, and characterisation of infectious pathogens, such as viral quasispecies. here we introduce a set of formulae for the combinatorial analysis of a quasispecies, given a ngs re-sequencing experiment and an algorithm for quasispecies reconstruction. the combinatorial analysis provided a description of the difficulty to reconstruct a quasispecies, given a determined amplicon partition and a measure of population diversity. results on empirical ngs data from patients infected with hepatitis b virus, confirmed its ability to characterise different viral variants from distinct patients. the reconstruction algorithm showed good performance both considering simulated data and real data, even in presence of sequencing errors. we require that sequenced fragments are aligned against a reference genome, and that the reference genome is partitioned into a set of sliding windows. the reconstruction algorithm is based on combinations of multinomial distributions and is designed to minimise the reconstruction of false variants, called in-silico recombinants. the reconstruction algorithm was applied to error-free simulated data and reconstructed a high percentage of true variants, even at a low genetic diversity, where the chance to obtain in-silico recombinants is high. although methodologies and software for whole genome assembly and genome variation analysis have been developed and refined for ngs data, reconstructing a viral quasispecies using ngs data remains a challenge.
 next-generation sequencing offers a unique opportunity for high-throughput genomics and has potential to replace sanger sequencing in many fields, including de-novo sequencing, re-sequencing, meta-genomics, and characterisation of infectious pathogens, such as viral quasispecies. although methodologies and software for whole genome assembly and genome variation analysis have been developed and refined for ngs data, reconstructing a viral quasispecies using ngs data remains a challenge. this application would be useful for analysing intra-host evolutionary pathways in relation to immune responses and antiretroviral therapy exposures. here we introduce a set of formulae for the combinatorial analysis of a quasispecies, given a ngs re-sequencing experiment and an algorithm for quasispecies reconstruction. we require that sequenced fragments are aligned against a reference genome, and that the reference genome is partitioned into a set of sliding windows. the reconstruction algorithm is based on combinations of multinomial distributions and is designed to minimise the reconstruction of false variants, called in-silico recombinants.  the reconstruction algorithm was applied to error-free simulated data and reconstructed a high percentage of true variants, even at a low genetic diversity, where the chance to obtain in-silico recombinants is high. results on empirical ngs data from patients infected with hepatitis b virus, confirmed its ability to characterise different viral variants from distinct patients.  the combinatorial analysis provided a description of the difficulty to reconstruct a quasispecies, given a determined amplicon partition and a measure of population diversity. the reconstruction algorithm showed good performance both considering simulated data and real data, even in presence of sequencing errors. 