metagenomics is defined as the study of genomic content of microbial communities in their natural environments , bypassing the need for isolation and laboratory cultivation of individual species .
its importance arises from the fact that over 99% of the species yet to be discovered are resistant to cultivation .
this limitation imposed by cultivation of isolated clones has severely skewed our view of microbial diversity .
metagenomics promises to enable scientists to study the full diversity of the microbial world , their functions and evolution , in their natural environments .
next generation sequencing ( ngs ) technologies generate data more efficiently , economically , and with a greater depth than ever before .
ngs has opened up an array of possibilities for many applications including whole - genome sequencing , epigenetics , and metagenomics . of these , the characterization of diversity of heterogeneous microbial environments , metagenomes , has recently gained significant interest .
although a host of methods for whole - genome assembly have been developed , reconstruction of individual clones from metagenomes still remains a challenge . as compared to existing technologies , reads produced by ngs are typically shorter and more error - prone .
the growth in the size of the datasets is fast outpacing the computational power needed to analyze it .
the computational method we present here aims to quantify the microbial diversity within a metagenome based on a set of deep sequencing reads . in single genome sequencing
, we can be certain that all extracted dna fragments belong to the same genome .
however , in majority of metagenomic samples , it is not possible to isolate and culture individual clones .
it is further complicated by the fact that the data comes from heterogeneous microbial communities , where the number of species as well as their relative abundance is unknown .
sequence data is usually partial and fragmentary , as environmental sequence sampling rarely produces all the sequences required for assembly .
moreover , metagenomic datasets are beset with increased amounts of polymorphism and horizontal gene transfer .
sequences from closely related species will most likely have homologous sequences shared between them , hindering their separation . as a result , the reconstruction of a whole genome is generally not possible . in the light of new data ,
an additional step in metagenomics that is not required in single genome assembly is binning the reads belonging to different species that is the need to associate the reads with its source organism . clustering methods aim to identify the species present in the sample , classify the reads by their species of origin , and quantify the abundance of each of these species .
most of the existing classification methods are supervised and depend on the availability of reference data for training [ 49 ] .
a metagenomic dataset may , however , contain reads from unexplored phyla which can not be labeled into one of the existing classes .
poor performance on the short fragments is mostly due to the high dimensionality and sparsity associated with the data . to overcome the limitation imposed by the length
however , there is the danger of assembling reads from different species together , thereby creating interspecies chimeras .
moreover , the abundances of different species can be potentially skewed such that the within - species variance overwhelms the between - species variance . in this paper
, we develop a method for clustering the short metagenome reads that addresses the challenges posed by the nature of metagenomic data .
we formulate an unsupervised two - way multispecies , multidimensional mixture model to represent reads from a metagenome .
we model the distribution of word counts along a genome as a gaussian for shorter , frequent words and as a poisson for longer words that are rare .
we employ either a mixture of gaussians or a mixture of poissons to model reads within each bin .
the proposed model is used to cluster metagenomic reads by their species of origin and to characterize the abundance of each species .
it is a composition - based method that seeks to distinguish between genomes based on their characteristic dna compositional pattern .
our method handles the high - dimensionality and sparsity associated with the data by grouping the set of words comprising the reads , to regularize the parameters in the mixture model .
this implies that , for every group , only one statistic of the words in this group is needed to cluster reads .
we show that a high clustering accuracy can be obtained at a much lower dimension .
we provide a framework that complements existing similarity - based methods . later in the paper
, we evaluate the applicability of the multidimensional mixture of distributions and its ability to estimate the parameters of genome abundance accurately , for simulated and real metagenomes .
we compare the performance of our method with likelybin and scimm , two other unsupervised composition - based method .
also , we demonstrate the robustness of our method to changes in the relative abundance of different species .
  the last decade has seen an explosion in the number of computational methods developed to analyze metagenomic data .
literature abounds in methods for classifying ( as opposed to clustering ) metagenome reads into taxon - specific bins [ 5 , 7 , 8 ] .
current approaches to metagenomics clustering can be classified into two main categories : similarity - based and composition - based methods .
the similarity - based approaches align the reads to close phylogenetic neighbors and hence depend on the availability of closely related genomes in existing databases [ 6 , 11 , 12 ] .
megan , a metagenome analysis software system , is a representative example of this kind .
it uses sequence homology to assign reads to common ancestors based on best match as given by blast ( basic local alignment search tool ) . as
most of the extant databases are highly biased in their representation of true diversity , such methods fail to find homologs for reads derived from novel species .
these methods rely on the intrinsic features of the reads such as oligonucleotide distributions [ 7 , 8 , 14 ] , codon usage preference , and gc composition   to differentiate between reads belonging to different species .
a significant limitation of most composition - based methods developed so far is that they do not perform well on reads shorter than 500  bp .
phylopythia   is a supervised composition - based classification method that trains a support vector machine to classify sequences of length greater than 1  kbp .
phymm uses interpolated markov models to characterize variable length dna sequences by their phylogenetic group .
its accuracy of assignment drops drastically ( to just 7.1% at genus level ) for short reads and reads from unknown species .
nasser et al .   demonstrated that a k - means based fuzzy classifier , trained using a maximal order markov chain , can separate fragments that are about 1  kbp long at the phylum level with a high accuracy .
compostbin is a semisupervised algorithm for grouping fragments that uses a novel weighted pca ( principal component analysis ) and a normalized cut clustering algorithm to classify the sequences .
they have demonstrated an error rate bounded by 10% , when guided by information from phylogenetic markers , on datasets of low complexity .
however , the accuracy of this method on reads less than 1  kbp has not been shown .
recently , chan et al . developed a semisupervised seeded growing self - organizing map ( s - gsom )   to cluster metagenomic sequences .
it extracts 813  kbp of flanking sequences of highly conserved 16s rrna from the metagenome and uses them as seeds to assign the remaining reads using composition - based clustering .
the caveat with soms is that it was shown to work well only on dna fragments that are longer than 8  kbp and lose much accuracy for reads with length below 1  kbp . all the above supervised methods
a metagenomic dataset may , however , contain reads from unexplored phyla which can not be labeled into one of the existing classes .
the accuracy of these methods on dataset containing reads from unknown species is yet to be demonstrated .
likelybin is an unsupervised method that clusters metagenomic sequences via a monte carlo markov chain approach .
scimm is a recently developed model - based approach to sequence clustering where interpolated markov models represent clusters and optimization is performed using a variant of the k - means algorithm . in the results section , we compare the accuracy of our proposed method with likelybin and scimm on datasets of different divergences .
however , if it is known a priori that the reads differ widely in their abundances , then we recommend using abundance bin over other binning methods .
  one of the most common genome signatures is the frequency of occurrence of words ( or oligomers ) in a dna sequence . in our method , we model each cluster , containing reads from a species , as a function of probability distributions of words comprising them .
the inherent basis of this method is that the set of reads sequenced from a species have a characteristic genome signature that distinguishes it from reads belonging to other species .
the distribution of word counts along a genome can be approximated as a gaussian for shorter , frequent words and as a poisson for longer words that are rare .
we propose an unsupervised multidimensional naive bayes poisson mixture model and derive an expectation maximization algorithm for the same . the corresponding algorithm for gaussian mixture model
however , with the increase in the length of words , the dimensionality of the data increases exponentially , while the word counts become sparse . to tackle high dimensionality and sparsity of word counts
in essence , the proposed method provides a general statistical framework for associating each read with its species of origin , based on its genome signatures .  
a genome signature is a compositional parameter reflecting the relative abundance of different words along the genome . in general
some words that are deemed to be biologically significant are very common in a genome , while others may never be encountered .
composition - based methods use genome signatures to ascertain the origin of the dna reads .
the underlying basis is that the distribution of words in a dna is specific to each species and undergoes only slight variations along the genome . by establishing the dictionary of words used by a species and their frequency of occurrence
literature abounds in methods that study the statistical distribution of the word locations along a sequence and word frequencies [ 21 , 25 ] .
the exact distribution of count of words is known under the hypothesis that the letters are independent ( bernoulli ) or under the markov model . however , in practice , it is extremely time consuming to compute the exact distribution for long sequences or for frequent words .
distribution of word counts along a genome can be approximately modeled as a gaussian distribution for short words ( that are more frequent ) , or a poisson distribution for longer words ( that are rare ) .
as different words occur with different frequencies along the genome , each word follows its own distribution .
thus , reads belonging to a species can be modeled as a multidimensional distribution of words ( one dimension for each word ) comprising them .  
figure 1 illustrates the distribution of dimer and pentamer counts across reads sampled from the genome of haemophilus influenzae ( for the purpose of clarity , only a few distributions are shown ) .
we see that count of each dimer ( a short word ) across the reads tends to a gaussian distribution with a different mean and standard deviation and that of a pentamer tends to a poisson distribution .
hence , the problem of clustering metagenomic reads can be cast as a multidimensional mixture of gaussians ( or poissons for longer words ) where distribution of each word is modeled as a gaussian ( or poisson ) . in other words
, this corresponds to the multidimensional naive bayes model , where each dimension is modeled as a unimodal gaussian ( or poisson ) distribution .
such a general statistical model takes into account the compositional heterogeneity of words along the genome .   in this paper
, we formulate an unsupervised multidimensional poisson mixture model for clustering reads within a metagenome by their species of origin .
we propose to model the reads from a species as a multidimensional distribution of the words comprising them .
therefore , each cluster is represented by the distribution of word counts within the species .
mixture models cover the data well , that is , dominant patterns in the data are captured by the component distributions .
they allow better approximations of the true distributions , and their parameters are relatively easy to estimate .
an additional advantage of using generative models is that they are flexible and can handle a large number of classes .
for instance , a mixture of poissons can be multimodal , while a poisson distribution is always unimodal .
we begin with a metagenome , x = { x1 , x2 ,  , xn } , containing n reads from m species .
let m be the proportion of species m in the dataset , with m=1m = 1 .
we assume that x is observed and is governed by some density function p(x |  ) with parameter . our goal is to cluster the reads by their species of origin , based on the frequency of words that appear in the reads .
for every species m , we want to determine m , its proportion in the dataset , and  , the parameter governing the distribution of words within the reads .
we assume that yi = m for m  1,  m if the ith read belongs to the mth species .
we call ( x , y ) the complete dataset . for a word of length l , we obtain p = 4 different words ( combinations of a , c , t , g ) , denoted by w = { w1 , w2 ,  , wp}. each read xi is represented by a p - dimensional feature vector , xi = { xi1 , xi2 ,  , xip } , where xij is the count of word wj in read xi .
we model the distribution of words within every species m by a multidimensional poisson distribution , say m = { m1 , m2 ,  , mp}. that is , given that read xi belongs to species m , the distribution of each word wj is poisson with parameter mj , where m = 1,2 ,  ,
m and j = 1,2 ,  , p , 
 ( 1)p(wj    mj)=(wj    mj)=emjmjxijxij!. 
 we assume independence between features of read vector .
the probability of a read xi , given it belongs to species m , is , 
 ( 2)p(xi    yi = m,)=p(xi    m)=j=1p(xij    mj )  . 
 at first glance , it might seem imprudent to represent a read as a collection of words comprising it , because it leads to the loss in information about the sequencing read . strictly speaking ,
even if the sequence of bases in a dna is independently and identically distributed , distribution of word occurrences are not independent , due to overlaps .
bayesian networks or belief networks can be used to represent the conditional dependencies between the words comprising the reads .
although , in practice , methods for exact inference in bayesian networks are often computationally expensive .
an attractive alternative to bayesian networks is the naive bayes algorithm that assumes independence between the different features of the read .
naive bayes is known to perform well on complex models and takes time that is linear in the number of components .
in addition , lost information can be restored at later stages . in this paper , we have presented the formulation of mixture models with the assumption that the different features ( word counts ) of the read are independent of each other .
to initialize the estimation algorithm , we randomly assign each read to a cluster m. the posterior probability qi , m is set to 1 , if read i is assigned to cluster m and 0 otherwise . with the initial posterior probabilities , a maximization step ( m - step )
expectation stepwe estimate the posterior probability qi , m of read xi belonging to species m. by bayes theorem , we have 
 ( 3)p(yi = mxi,)=mp(xim)k=1mkp(xik)=qi , m , qi , mmj=1p(xij    mj )  subject  tom=1mqi , m=1 . 
 
we estimate the posterior probability qi , m of read xi belonging to species m. by bayes theorem , we have 
 ( 3)p(yi = mxi,)=mp(xim)k=1mkp(xik)=qi , m , qi , mmj=1p(xij    mj )  subject  tom=1mqi , m=1 . 

maximization stepthe m - step uses qi , m to compute the expectation of complete data log likelihood , 
 ( 4)q((t+1),(t))=ep(y    x,)[log  p(x , y    )]=m=1 m  i=1np(yi = m    xi,(t ) )  log(p(xi , yi = m    (t+1)))=m=1 m  i=1n(qi , mlog(mp(xi    m ) ) ) . 

we also take into account the constraint , which requires that m 's sum to 1 by adding a lagrange multiplier : 
 ( 5)q((t+1),(t))=m=1 m  i=1n(qi , mlog(mp(xi    m ) ) )  + (m=1mm1 ) . 

we maximize the above expression with respect to the parameters ,  = arg  maxq( ,  ) , and update the parameters , 
 ( 6)m(t+1)=i=1nqi , mn ,  mj(t+1)=i=1nqi , mxiji=1nqi , m . 
 finally , these two steps are repeated as necessary .
each iteration is guaranteed to increase the log - likelihood , and the algorithm is guaranteed to converge to a local maximum of the likelihood function . 
 
the m - step uses qi , m to compute the expectation of complete data log likelihood , 
 ( 4)q((t+1),(t))=ep(y    x,)[log  p(x , y    )]=m=1 m  i=1np(yi = m    xi,(t ) )  log(p(xi , yi = m    (t+1)))=m=1 m  i=1n(qi , mlog(mp(xi    m ) ) ) . 

we also take into account the constraint , which requires that m 's sum to 1 by adding a lagrange multiplier : 
 ( 5)q((t+1),(t))=m=1 m  i=1n(qi , mlog(mp(xi    m ) ) )  + (m=1mm1 ) . 

we maximize the above expression with respect to the parameters ,  = arg  maxq( ,  ) , and update the parameters , 
 ( 6)m(t+1)=i=1nqi , mn ,  mj(t+1)=i=1nqi , mxiji=1nqi , m . 
 finally , these two steps are repeated as necessary .
each iteration is guaranteed to increase the log - likelihood , and the algorithm is guaranteed to converge to a local maximum of the likelihood function .  
higher - order words are known to be more discriminative than shorter ones . with the increase in length of the word
firstly , the distribution of words tends to poisson and not gaussian ( by law of rare numbers ) , see figure 1 .
secondly , the length of the read vector grows exponentially ( e.g. , for l = 10 , 4  10 ) . with increase in dimensions , many words will tend to have similar distributions and hence can be clustered together into a  word group .  at the same time , the number of distinct words in any read is usually substantially smaller than the number of dimensions .
hence , the model may fail to predict the true feature distribution of different components .
however , reduction of the number of words using feature selection can not be too aggressive , otherwise the clustering accuracy will suffer .
 a supervised two - way poisson mixture model with word grouping was originally proposed by li and zha for simultaneous document classification . such a two - way clustering involves simultaneous clustering of reads as well as of words .
the clusters means are regularized by dividing the words into groups and constraining the parameters for the words within the same group to be identical .
the grouping of the words is not predetermined but optimized as part of the model estimation .
this implies that , for every group , only one statistic for all the words in this group is needed to cluster reads .
for instance , in figure 1 , we observe the distribution of pentamers falls into three distinct group .
therefore , words following similar distributions can be clustered together into a  word group .
 we extend our formulation to an equivalent two - way unsupervised poisson mixture model in order to simultaneously cluster word features and classify reads and derive an expectation maximization algorithm to estimate its parameters .  
note that we make a distinction on the use of  cluster  to refer to binning of reads belonging to the same species and  group  to refer to binning of words within read in a cluster .
recall that the genome signature is similar between closely related species and dissimilar between nonrelated species .
therefore , we can assume that , within each cluster , words in different reads have equal poisson parameters , while , for reads in different clusters , words may follow different poisson distributions .
for simplicity , we assume that all clusters have the same number of word groups .
it is trivial to extend to the case where different clusters may have different number of word groups .
, l denote the word groups . we define a group assignment function c(m , j )  1,2 ,  , l ,
which denotes the group to which word wj belongs in class m. words in the same word group will have the identical parameters , that is , mk = mj = m , l , if c(m , k ) = c(m , j ) .
the likelihood of xi is now 
 ( 7)p(xi    m)=j=1pp(xij    mj)=j=1pp(xij    m , c(m , j ) ) . 

we can derive an em algorithm similar to the one outlined above to estimate the poisson parameters m , l where m  1 ,  , m , l  1 ,  ,
m , j  1 ,  , p and the prior mixture components m , for m  1 ,  , m. we initialize by setting each value of the group assignment function c(m , j ) randomly to a number in 1 ,  ,
l. we start with the same word group partition for all the clusters , that is , c(m , j)s are initially identical over m. we update the parameters as follows : 
 ( 8)m(t+1)=i=1nqi , mn ,  m , l(t+1)=i=1nqi , mj1xijmli=1nqi , m ,  where  c(m , j)=l . 
 once ml is fixed , the word cluster index c(m , j ) can be found by doing a linear search over all components : 
 ( 9)c(t+1)(m , j)=argmaxli=1qi , m(xijlogm , l(t+1)m , l(t+1 )  ) . 
 

xp ) are independent poisson variables with parameters 1 , 2 ,  , p , respectively , then the conditional distribution of ( x1 , x2 ,  , xp ) given that x1
+ x2 + +xp = n is multinomial with parameters j/ , where  = j , that is , mult(n ,  ) , where  = ( 1/ , 2/ ,  , p/ ) . 
  if ( x1 , x2 ,  , xp ) are independent poisson variables with parameters 1 , 2 ,  , p , respectively , then the conditional distribution of ( x1 , x2 ,  , xp ) given that x1 + x2 + +xp = n is multinomial with parameters j/ , where  = j , that is , mult(n ,  ) , where  = ( 1/ , 2/ ,  , p/ ) .
the above theorem implies that the unconditional distribution ( x1 , x2 ,  , xp ) can be factored into a product of two distributions : a poisson for the overall total and a multinomial distribution of x , x ~ mult(n ,  ) .
therefore , the likelihood - based inferences about  are the same whether we regard x1 , x2 ,  , xp as sampled from p independent poissons or from a single multinomial .
here , n refers to the length of the reads and our interest lies in the proportion of words in the reads .
any estimates , tests , inferences about the proportions will be the same whether we regard n as random or fixed .  
we can now derive the naive bayes mixture of multinomials as standardized mixture of poissons .
we assume that the distribution of words within the reads of a species is governed by the parameters of a multinomial distribution  = ( 1 , 2 ,  , m ) , where each m is the parameter for species m and
therefore , the likelihood of the data will be 
 ( 10)p(xi    yi = m)=p(xi    m)=ni!j=1pxij!j=1pmjxij , 
 the sum of the probabilities satisfies the constraint j=1mj = 1 . the em algorithm for naive bayes mixture of multinomials
can be derived similarly , and we only give the final set of equations : 
 ( 11)m=i=1nqi , mn,mj=i=1nqi , mxiji=1nj=1pqi , mxij=i=1nqi , mxiji=1nqi , mni . 
 if we assume the length of each read to be a constant n , we get the same results as that with poisson distribution ; hence , the two distributions are equivalent in modeling the distribution of words within reads of a species . also ,
since the multinomial distribution is single distribution , we do not perform a two - way dimension reduction on it .
the algorithm has been implemented in matlab and c. the space and time complexity scale linearly with the number of reads and species .
our method converged for all the cases we tested and was robust to the choice of initial conditions .
metagenomics being a relatively new field lacks standard datasets for the purpose of testing clustering algorithms . as the  true solution  for sequence data generated from most metagenomic studies
we also apply our method to the actual acid mine drainage dataset to identify the dominant species . in order to test the accuracy of our proposed method , we used metasim to simulate synthetic metagenomes
. metasim takes as input the sequencing technology to be used ( sanger , 454 , exact ) , a set of known genomes , length of the reads , and an abundance / coverage profile which determines the relative abundance of each genome in the simulated dataset .
the genomes used for generating the synthetic metagenomes were downloaded from national center of biotechnology information ( ncbi ) .
we generated datasets with reads of lengths between 50 and 1000  bp and various abundance ratios .
in the first part of this section , we demonstrate the performance of the multidimensional gaussian mixture model on several datasets .
as the number of dimensions is relatively small , we do not perform word grouping .
next , we describe the results using the two - way poisson mixture model with a word length of 5 . the method has been implemented for word lengths from 2 to 9 . in order to calculate the clustering accuracy
, we assign each cluster to the source species that is most frequent in the cluster .
the number of species in each dataset is supplied as an input . determining the number of clusters from a statistical perspective is a difficult problem and has been addressed by .
previously , 16  s/18  s rdna have been used for phylotyping and assessing species diversity using a rarefaction curve .
tools such as metaphyler and treephyler can be used for making an educated guess of the number of species [ 34 , 35 ] .
estimating species diversity is still an active area of research , and we do not address it in this paper .
experiments in the 1960s and 1970s have shown that the dinucleotide relative abundance in a genome is a remarkably stable property [ 36 , 37 ] .
closely related organisms display more similar dinucleotide composition than do distant organisms . in , the authors proposed a measure of intergenomic difference between two sequences f and g , called the average dinucleotide relative abundance , 
 ( 12)(f , g)=116x , y|xy(f)xy(g)| , 
 where xy*(f ) = fxy*/fy*fy * and fx * denotes the frequency of x in f. a measure of intergenomic difference was obtained by comparing different genome signatures . in order to assess the robustness of our method , we test it across datasets representative of  * values ranging from 34 to 340 . in general , lower 
* values correspond to  closely related species  and higher values correspond to  distant species .  in figure 3 , we plot the performance of our proposed multidimensional poisson model over 450 datasets with  * values ranging from 34 to 340 . we observed a positive correlation between the intergenomic difference and the accuracy of our method , as also noted in .
the initial increase in the accuracy with word length is justified by the increased discriminative power of higher order words .
however , any further increase in word length has to be accompanied by dimension reduction , otherwise owing to the high dimensional and sparse nature of feature matrix , the accuracy begins to drop . in figure 4
, we compare the accuracy of our proposed multidimensional gaussian model with two other unsupervised composition - based methods likelybin   and scimm   on several datasets .
we varied the read length between 200 to 500  bp ,  * values from 60 to 300 , and the abundance ratio up to 1  :  5 .
note that the distribution of dimers tends to a gaussian . as the number of dimensions is relatively small ( 4 = 16 )
our method clearly outperforms likelybin and performs as well or better than scimm on most instances .
another point worth noting from the figure is that our method 's error rate is bounded by 10% for datasets with read length as short as 200  bp .
we analyzed the accuracy and applicability of our method on binning reads from low complexity communities , containing 35 species ( see table 1 ) . with the increase in number of species
, there was a slight degradation in performance , though the accuracy was consistently above 85% .
this is in agreement with the results from the 2 species dataset , considering that the total coverage of each species is much lower in a multispecies dataset ( reads from b. burgdorferi form only 6% of the 5th dataset ) .
next , we evaluated the robustness of our method to changes in the abundance ratio between species as well as the length of the reads .
we varied the abundance ratio from 10  :  1 to 1  :  10 in stages for the two species . from figure 5
, we note that there was only a slight drop in performance for extreme abundance ratios .
it is noteworthy to point out that estimates are good at all abundances . in order to test the usefulness of the method for analyzing data produced by the current ngs technologies ( especially solexa and solid ) that generate short reads , we tested three datasets of varying  * values for read lengths between 50 and 1000  bp . with
the decrease in read length from 1000 to 50  bp , the drop in accuracy of our method is bounded by 15% .
recall that with the increase in the length of the words and the simultaneous increase in the number of dimensions , the distribution of the words tends to a poisson and word grouping becomes necessary . in this section , we present the clustering results obtained by estimating the two - way poisson mixture model with different number of word groups l. we observed the variation in classification accuracy to be more prominent for lower values of l. therefore , in table 2 , we report the results for values of l < 50 for a 2 species dataset .
if word grouping is not performed , then clustering based on mixture model is essentially the naive bayes algorithm with each dimension modeled by a poisson distribution ( last column of table 2 ) . from the results
, we can infer that word grouping resulted in considerable increase in accuracy compared to the naive bayes algorithm .
that is , the characteristic vectors are of a much lower dimension with l  p. also , a high clustering accuracy can be achieved using no more than ml dimensions , significantly smaller than the original dimension , 1024 .
note that it is difficult to know a priori , the exact value of l that yields the best clustering .
however , among the values we tested , lower values of l provided a higher accuracy .    in table 3
, we compare the performance of our 2-way poisson mixture model with gaussian mixture model for datasets with low  * values . in real situations , it is difficult to know beforehand , the most discerning order of the word to use .
however , from our experiments , we can infer that higher - order word - based models , in general , tend to be more discriminatory than those based on lower order words .
if it is known a priori that lower - order words ( of length 2 - 3 ) are more discriminatory in the dataset , then we recommend using a gaussian mixture model . for other datasets ,
our method 's accuracy in classifying reads from the datasets composed of species across various taxonomic ranks is reported in table 4 we used the poisson mixture model without word grouping .
we can infer that the accuracy is mostly correlated to the phylogenetic distances between the species .
for example , reads from datasets containing species with taxonomic differences at the level of class were classified with a very high accuracy .   the ultimate goal of binning methods is to cluster reads in a real metagenome , by their species of origin .
clustering in real situations is error - prone and affects our final estimates of species abundance .
moreover , evaluating clustering methods on real metagenomes can be problematic as the true taxonomic composition of the data is mostly unknown .
the accuracy of unsupervised clustering methods decreases with increase in the complexity of metagenomes and for species present at very low abundances
. however , the composition of acid mine drainage metagenome has been substantially characterized , and we used this dataset to evaluate the performance of our proposed method .
the amd microbial community is reported to consist of two dominant populations ( ferroplasma sp .
group ii ) and three other less abundant ones ( ferroplasma acidarmanus type i , leptospirillum sp .
we downloaded the reads , as well as the scaffolds assembled from the reads for the 5 species of the actual amd dataset from ncbi .
only 58% of the amd reads can be mapped back to the assembled scaffolds using blast .
therefore , in order to compute the accuracy of our method , we simulated a metagenome with reads sampled from the downloaded scaffolds .
the simulated amd dataset consisted of 110,000 reads of average length 732  bp ( average read length in the actual amd dataset ) from the 5 species , in the ratio 4  :  4  :  1  :  1  :  1 .
therefore , we can simplify the problem by first separating the reads into two bins based on their abundance .
in the first stage , the reads were grouped into two bins , using abundance bin , with a resulting accuracy of 93.3% .
the bins corresponding to the abundance levels of 4 and 1 had a cluster purity of 93.2% and 98.2% , respectively . in the next stage
, we used the reads from each of the bins output by abundance bin , as an input to our proposed 2-way poisson mixture model , to further classify the reads by their species of origin .
our method clustered the reads from the bin containing dominant species into two clusters corresponding to ferroplasma sp .
the other bin consisted of very few reads from the remaining three species ferroplasma acidarmanus type i , leptospirillum sp .
our method clustered the reads from this bin into three clusters , with an accuracy of 70.34% ( with l = 10 ) .
  in this paper , we formulated an unsupervised two - way multispecies , multidimensional mixture model to represent reads from a metagenome .
we used the proposed model to cluster metagenomic reads by their species of origin and to characterize the abundance of each species .
the distribution of word counts along a genome can be approximated as a gaussian for shorter , frequent words and as a poisson for longer words that are rare .
therefore , we use a multidimensional mixture of gaussians or poissons to model the reads from each bin . an additional reason to use these distributions
this is critical for success as most metagenomic datasets contain reads from unexplored phyla which can not be labeled into one of the existing classes .
our probabilistic approach can be used to identify reads which belong to more than one species and occlude the cluster boundaries .
note that our proposed method is primarily a composition - based method that seeks to distinguish between genomes based on their characteristic dna compositional pattern .
therefore , it can not distinguish between genomes unless their dna compositions are sufficiently divergent ( see figure 4 , dataset with b. burgdorferi and c. jejuni ) .
it is unlikely that our method will be able to accurately distinguish between strains of the same species . for such datasets ,
composition - based methods must be used in conjunction with other similarity - based methods and abundance - based methods to yield better performance .
note that the two - way poisson mixture model was originally proposed for classification of documents . in this work ,
we demonstrate the relevance and applicability of such a general statistical framework for modeling metagenome reads .
we have illustrated that the proposed method can accurately classify reads from low to medium complexity datasets into taxon - specific bins , based on genome signatures .
our framework complements the existing similarity - based and abundance - based methods and hence can be combined with such methods to obtain a better performance .
we intend to develop such hybrid methods in the future that can tackle the problem of classifying sequences in complex metagenomic communities .
 a major challenge facing metagenomics is the development of tools for the characterization of functional and taxonomic content of vast amounts of short metagenome reads . 
 the efficacy of clustering methods depends on the number of reads in the dataset , the read length and relative abundances of source genomes in the microbial community . in this paper 
 , we formulate an unsupervised naive bayes multispecies , multidimensional mixture model for reads from a metagenome . 
 we use the proposed model to cluster metagenomic reads by their species of origin and to characterize the abundance of each species . 
 we model the distribution of word counts along a genome as a gaussian for shorter , frequent words and as a poisson for longer words that are rare . 
 we employ either a mixture of gaussians or mixture of poissons to model reads within each bin . 
 further , we handle the high - dimensionality and sparsity associated with the data , by grouping the set of words comprising the reads , resulting in a two - way mixture model . 
 finally , we demonstrate the accuracy and applicability of this method on simulated and real metagenomes . 
 our method can accurately cluster reads as short as 100 bps and is robust to varying abundances , divergences and read lengths . 
1. Introduction
2. Related Work
3. Methods
4. Results
5. Discussion