- Methodology article
- Open Access

# Microbial comparative pan-genomics using binomial mixture models

- Lars Snipen
^{1}Email author, - Trygve Almøy
^{1}and - David W Ussery
^{2}

**10**:385

https://doi.org/10.1186/1471-2164-10-385

© Snipen et al; licensee BioMed Central Ltd. 2009

**Received: **24 April 2009

**Accepted: **19 August 2009

**Published: **19 August 2009

## Abstract

### Background

The size of the core- and pan-genome of bacterial species is a topic of increasing interest due to the growing number of sequenced prokaryote genomes, many from the same species. Attempts to estimate these quantities have been made, using regression methods or mixture models. We extend the latter approach by using statistical ideas developed for capture-recapture problems in ecology and epidemiology.

### Results

We estimate core- and pan-genome sizes for 16 different bacterial species. The results reveal a complex dependency structure for most species, manifested as heterogeneous detection probabilities. Estimated pan-genome sizes range from small (around 2600 gene families) in *Buchnera aphidicola* to large (around 43000 gene families) in *Escherichia coli*. Results for *Echerichia coli* show that as more data become available, a larger diversity is estimated, indicating an extensive pool of rarely occurring genes in the population.

### Conclusion

Analyzing pan-genomics data with binomial mixture models is a way to handle dependencies between genomes, which we find is always present. A bottleneck in the estimation procedure is the annotation of rarely occurring genes.

## Keywords

## Background

One of the consequences of the explosion in numbers of fully sequenced and annotated microbial genomes is that we are now facing the challenges of comparative pan-genomics [1]. The microbial pan-genome, as defined by [2], is the number of essentially different genes found within a population at a specified taxonomic level, usually within a species, though this can be extended to higher levels, such as genus. As multiple genomes of the same species are sequenced, one can construct the pan-genome, and begin to compare pan-genomes from different species.

Having a set of fully sequenced and annotated genomes from several strains within a species, one is interested in two sets of genes. The first is the set of core genes, *i.e*. the genes found in every strain within a species. The size and content of the core genome is interesting for characterizing the genomic essence of the species. The other set is the pan-genome, which is the total number of different genes found in all strains within the species. The size of this pan-genome, relative to the number of genes found in a typical strain, is an indicator of the plasticity of the species, and could be reflective of its potential for adaptation in a diverse environment.

The true core- and pan-genome sizes, here denoted *γ* and *η* respectively, will most likely remain unknown for any species, since it is impossible to sequence and annotate all existing strains. Thus, we have to rely on estimates based on existing data. The problem of estimating the size of the core- and pan-genome was first approached by [2]. They used an exponential function to explain the number of new genes introduced by each new sequenced genome, and by extrapolating this they came up with some estimates of the pan-genome size. The core-genome size was also estimated in a similar way. Modified versions of this approach have later been used by others. For example the number of new *Escherichia coli* genes contributed by each additional genome sequenced was first estimated to be rather large – 440 genes by [3]. More recent estimates, based on 17 different isolates from a wide variety of strains, brought the number of expected novel genes per new genome to be around 300, with approximately 13,000 genes estimated to be in the total *E. coli* pan-genome [4]. Based on comparison of 32 *E. coli* genome sequences, we have previously estimated the number to be around 80 novel genes per genome, with a pan-genome size of just under 10,000 genes [5].

One of the implications of early pan-genome estimates is that some bacterial species might have an "infinite" pan-genome [2, 6]. This is a dramatic statement, especially since it can be largely due to a bias from their use of an exponential model, which inherently assumes the pan-genome can be divided into two groups: The core-genes always present in all genomes, and the dispensable genes, equally likely to occur in any genome. The latter part of this assumption is often far from reality, which we will show in this paper. This was also recognized by [7], who was the first to introduce a mixture model to estimate the core- and pan-genome size. Unfortunately, they also imposed some rather heavy restrictions in their model, making their pan-genome estimates biased towards larger values.

We will, however, extend the good idea of [7] in this paper, and by avoiding their heavy restrictions hopefully come up with more realistic estimates of core- and pan-genome sizes.

## Results

### Algorithm

#### Gene families

For a given species *G* different genomes have been sequenced and annotated. The first step in any pan-genome analysis is to come up with a list of gene families in the current sample. A deeper analysis of this problem is not the focus of this paper, and we have at this stage taken the approach used by [7] and [5]. First an all-against-all BLASTing (blastp) is performed, and only alignments with at least 50% identity along at least 50% of both sequences are considered. Two sequences belong to the same gene family if both their reciprocal alignments fulfill the 50-50-cutoff rule. The results of this procedure is typically stored in a *pan-matrix* M= {*m*_{
ij
}} where each row corresponds to a gene family and each column to a genome. If gene family *i* has at least one member in genome *j* then *m*_{
ij
}= 1, else *m*_{
ij
}= 0.

#### Mixture model

The pan-genome size, *η*, is the number of gene families found in all strains, also including the gene families not yet observed in the *G* genomes sequenced so far. Summing row *i* in M we get the number of genomes in which gene family *i* has been observed. Tabulating all these row-sums gives us the number of gene-families observed in 1,..., *G* genomes, which we denote *y*_{1},...*y*_{
G
}. The sample pan-genome size is
, while *y*_{
G
}is usually listed as the sample core-genome size. The true pan-genome size also includes *y*_{0}, the number of gene families observed in zero genomes so far. Hence *η* = *n* + *y*_{0} and estimating *η* is equivalent to predicting *y*_{0}.

*y*

_{0}we need a model that relates

*y*

_{0}to

*y*

_{1},...,

*y*

_{ G }. Consider y= (

*y*

_{0},

*y*

_{1},...,

*y*

_{ G }). Since the total sum of gene families,

*η*, is constant y is a multiniomial vector if we assume independence between gene families,

*i.e*. y~ Mult(θ,

*η*). The multinomial probabilities θ= (

*θ*

_{0},...,

*θ*

_{ G }) are the probabilities of a gene family to be detected in 0,...,

*G*genomes, respectively. The expected value of

*y*

_{0}is

*E*(

*y*

_{0}) =

*ηθ*

_{0}due to the multinomial model. Also, a similar argument leads to

*E*(

*n*) =

*η*(1 -

*θ*

_{0}). Combined they lead to

*n*as an estimate of

*E*(

*n*) we can predict

*y*

_{0}if we can estimate

*θ*

_{0}. This estimate can be found by assuming some degree of smoothness across the multinomial probabilities. One way of obtaining this is by using a binomial mixture model. This means we assume

*detection probability ρ*

_{ k }. Thus, the multinomial probabilities are expressed as a combination of

*K*binomial probability mass functions (PMF). The shape and location of these binomial PMFs will determine how

*θ*

_{ g }are related to each other, and more specifically how

*θ*

_{0}relates to

*θ*

_{ g },

*g*= 1,...,

*G*. Figure 1 illustrates this idea for a three component model, i.e. we use the combination of three binomial PMFs to describe the 11 multinomial probabilities. Component

*k*in this mixture model may be interpreted as a class of gene families with probability

*ρ*

_{ k }of being detected (probability of "success") in a genome. If

*ρ*

_{ k }is low, these genes are typically rarely observed in the sequenced genomes, and vice versa. A binomial mixture like this was also used by [7].

It is natural to reserve one of the mixture components for the class of core genes. Core genes are special, since these genes should always be present in all genomes, and it is natural to assign them detection probability 1.0, as was also done by [7]. We define the first component as the core component, hence *ρ*_{1} = 1.0.

#### Estimation

The parameters of the binomial mixture model cannot be estimated directly from y, again because *y*_{0} is missing. This led [7] to impose some heavy restrictions on their model, which is not necessary. A commonly used approach for such models is to estimate parameters maximizing the zero-truncated log-likelihood [8].

*n*the vector y

_{+}= (

*y*

_{1},...,

*y*

_{ G }) is also a multinomial, with probability

*θ*

_{ g }

*/*(1 -

*θ*

_{0}) for element

*g*= 1,...,

*G*. Thus, the zero-truncated log-likelihood is

where *θ*_{0},..., *θ*_{
G
}depend on π and ρ as described in (2) and (3), and *C* is a constant independent of these parameters. Thus, for some choice of *K*, we estimate π and ρ by maximizing the criterion in (4), which only involves (*y*_{1},..., *y*_{
G
}). This can be done with some iterative optimization algorithm. These estimates, denoted
and
for *k* = 1,..., *K*, are used in (3) and (2) to get the estimates of *θ*_{0} and this is in turn plugged into (1) to compute the corresponding prediction
.

*K*in the binomial mixture, i.e how many binomial PMF do we need to approximate the distribution of the observed data (

*y*

_{1},...,

*y*

_{ G }). Since our criterion in (4) is a log-likelihood function for the data, we have adopted the Bayesian Information Criterion (BIC) to select the proper model complexity [9], a choice also supported by [10]. Hence, we look for a

*K*where

is minimized, where (2*K* - 2) is the number of free parameters in the model since the sum of mixing proportions is always 1.0 and the core component has a fixed detection probability *ρ*_{1}.

*K*the estimated core- and pan-genome sizes are

where is the estimated mixing proportion for component 1, the core component.

We have observed that the pan-size estimate may be heavily influenced by the chosen number of components, a generic property discussed by [10]. In order to stabilize the estimates, [10] propose a bagging-based estimator, which we have adopted. This is a bootstrap procedure that will smooth the estimate over various choices of components, and making the final estimate more stable.

Notice that this corresponds to *y*_{0} being predicted from *y*_{1} and *y*_{2} only.

### Implementation

All computations, including the parsing of BLAST results, setting up the pan-matrix and performing all estimations have been implemented in R [12] and is freely available from the corresponding author. An R-package for microbial pan-genomics is under construction and will be made available as soon as it is operational.

### Testing

#### Estimating core- and pan-sizes

*Francisella tularensis*currently has the largest fraction covered, at 73%; this seems reasonable, in that the total pan-genome for an intracellular organism would be expected to be relatively small, compared to environmental isolates. The bacterial species with the smallest fraction of the estimated total pan-genome covered is that of

*E. coli*, with a mere 30% covered so far, based on 22 genomes completely sequenced.

*E. coli*, the size of the estimated total pan-genome is about 43,000 gene families – or nearly twice the size of the human genome. On the other hand, for

*B. aphidicola*, the total pan-genome is estimated to be about 2600 genes.

#### Distribution of gene families

*E. coli*has only a fraction of the pan-genome covered by the observed data, with one quite large component that is red (very small detection probability). On the other hand,

*F. tularensis*has most of the pan-genome already covered by the data examined; this is also the case for

*Coxiella burnetii*and

*Yersinia pestis*.

#### Effect of growing data set

*E. coli*, we have already 22 fully sequenced genomes. Still, the coverage, defined as sample pan-genome size divided by estimated pan-genome size, is as low as 30%. An interesting question is of course how many more genomes do we need to sequence in order to have a coverage of, say, 90% of the

*E. coli*pan-genome? Upon examination of this question, we discovered that this number appears to grow as more genomes are sequenced. That is, with only a few genomes sequenced, it might appear say that 100 genomes might be enough to cover the estimated pan-genome. However, even with only 22 genomes sequenced, now it looks as though perhaps around 220 additional

*E. coli*genomes would be needed. In coming up with this estimate, we find that, as more

*E. coli*genomes are sequenced, the total estimated diversity increases, resulting in a steep increase in the estimate of the pan-genome total size, as shown in Figure 5.

#### Effect of gene prediction

The use of a mixture model makes it apparent that the estimate of pan-genome size must depend on how many gene families we observe in few genomes. Especially those gene families observed in only one genome, are most likely important. These genes are often referred to as ORFans. Upon inspection of the data, we found that the annotation "hypothetical protein" is severely over-represented among the ORFans in all 16 species (Fisher exact test p-values less than 10^{-10}). Thus, false positives from the gene prediction, *i.e*. predicted gene who are not actually genes, are most likely influencing the number of ORFans most since false positives typically are "hypothetical proteins". This makes the number of ORFans uncertain, and estimation of pan-genome size even more difficult.

*E. coli*data, which is the largest data set. First, we removed 10% and 50% of the shortest hypothetical proteins in the data set, because we believe these are the most uncertain predictions. A pan-genome size was estimated for these reduced data sets. Next, we also made a completely new prediction of genes for all 22 genomes using the Easygene tool ([14, 15]), and made another estimate from these data as well. The results are displayed in Table 1. The number of ORFans drops dramatically consistent with the idea that perhaps a large fraction of the ORFans are due to artifacts of gene finding. The pan-size estimates also tend to decrease as an effect of this, but the mixture model estimates show some variability.

Effect of gene predictions

Data set | Observed | ORFans | Chao | Bin. mix. |
---|---|---|---|---|

Original NCBI | 12599 | 5438 | 26614 | 42640 |

Reduced 10% | 11273 | 4470 | 22549 | 32528 |

Reduced 50% | 9336 | 3272 | 17083 | 27456 |

Easygene | 9211 | 3121 | 17041 | 29818 |

## Discussion

The use of a binomial mixture model for estimating the pan-genome size was introduced by [7], but the use of mixture models for population size estimation is by no way new, *e.g*. [8, 10, 16]. The estimation of a population size has a long history in ecology, under the names of capture-recapture problems (*e.g*. [17]), or in epidemiology, called multiple record systems (*e.g*. [18]). Mixture models are suitable when we are faced with a larger number of recaptures/records/genomes and heterogeneous detection probabilities, which is exactly the case for pan-genomics.

From our results in Figure 2 we notice that for none of the species the optimal mixture model has 2 components. This would be expected if the gene pool could be divided into core-genes and dispensable genes, as implicitly assumed by [2, 6]. There is always at least a third group, and frequently even more. This observation corresponds to the results shown by [19], where they find that for bacteria and archaea in general, genes could be divided into three classes; core (always occurring), shell (moderately occurring) and cloud (rarely occurring).

A reason for this heterogeneity in detection probabilities may be skewed sampling. If some of the sequenced genomes are sampled in the same "corner" of the population, the genes characteristic for this "corner" will occur more frequently than they should. Another reason may be that some genes are simply frequently occurring in the population, reflecting a divergence from a fairly recent ancestor. In this perspective, it must be expected that there is a large number of true detection probabilities, which is at least partly supported by the fact that the more genomes we consider the more components we estimate (see Figure 2).

The fact that microbial genomic diversity is caused by both vertical mutations and horizontal transfer makes it also plausible to expect heterogenous detection probabilities.

From Figure 2 we also see that even for 22 genomes (*E. coli*) we only estimate 6 components. In [7] a mixture of 7 components were used for a data set of 8 genomes, which seems to be a too complex model. Using too complex mixture models will tend to over-estimate the pan-genome size, since it makes the estimate of the smallest detection probability artificially small.

In Figure 3 we see that a larger sample pan-genome tends to result in a larger estimated pan-genome.

This is due to the fact that larger data sets allow more complex models, and more complex models allow more extreme estimates. Uncertainties, as indicated by the rough confidence intervals, also tend to grow when estimates grow, which is reasonable.

In Figure 4 we have constructed a way to plot the estimated mixture models for comparative pan-genomics. In this picture the actual size of the core- and pan-genome is not visible, but we focus instead on the relative distribution of detection probabilities. Some species, typically have a large proportion of stable genes (blue area), while at the other end of the scale we find those with little overlap between genomes. A larger number of components indicates a more complex pan-genome with respect to heterogeneity in detection probabilities.

From the results in Figure 3 we can compute the coverage for each species, which is simply the size of the sample pan-genome divided by the estimated pan-genome size. Ideally, we should expect this to increase as the number of genomes increase, because the sample pan size should approach the true pan size. There is no such tendency in our results. We even observe that two of the largest data sets (*S. enterica* and *E. coli*) have two of the smallest coverages. Figure 4 also clearly demonstrates that, at least for *E. coli*, as more genomes become available the pan-genome estimates get even higher. This is typical for a population with a large fraction of ORFans. Since ORFans have a small detection probability, only a few of them will show up in every genome. Hence, it requires a substantial number of genomes before we can estimate their true abundance. In this perspective, the binomial mixture model will tend to under-estimate the true pan-size for smaller data sets.

In Table 1 we show that there are effects of possible false positive predicted genes on the estimates of pan-genome size. By removing hypothetical proteins from the data set, the number of ORFans drops. This again leads to a decreased pan-size estimates. Predicting new genes with Easygene gives the largest reduction in ORFans, but the effect on the mixture model estimated pan-size is less. This is due to the fact that the mixture model depends on the entire data distribution, not only the ORFans.

Our approach assume a closed pan-genome, i.e. *η* is a parameter. In an open pan-genome, the total number of genes is not fixed, and in a very long term perspective this is most likely the case, assuming new genes form and old genes disappear. However, in a reasonably short time window, the number of genes available to any population must be limited, and can be assumed constant. Wether genes are shared vertically or horizontally within the population should have no impact on the closedness of the gene pool.

A recent publication [20] has suggested alternative ways of estimating pan-genome size, based on power-laws and regression. Our, more probabilistic approach, is fundamentally different, and more in line with existing methods in capture-recapture modelling. However, as suggested by the results in Table 1, a major problem in pan-genome size estimation is the fact that the data themselves are estimates, and thus the uncertainty in the computation of gene families will influence the results, sometimes severely. In order to improve the estimation of bacterial genomic diversity, future efforts should probably be focused on this aspect.

## Conclusion

We have shown how to use binomial mixture models to estimate microbial core- and pan-genome size, and the vast literature on capture-recapture methods should be further exploited in microbial pangenomics, as it has been in closely related fields like metagenomics [21]. Our results indicate that pan-genomes for bacterial species are in general large compared to the size of individual genomes. Especially for *E. coli*, who has the largest number of completely sequenced and annotated genomes so far, we find that the pan-genome is significantly larger than the human genome. We also show that our pan-size estimates are most likely too moderate since the addition of new genomes tend to push them upwards. In order to improve reliability of estimates, more focus should be devoted to the computation of gene families.

## Declarations

### Acknowledgements

We wish to thank Carsten Friis, Centre for Biological Sequence Analysis, Technical University of Denmark, for his assistance in performing the Easygene predictions.

## Authors’ Affiliations

## References

- Read TD, Ussery DW: Opening the pan-genomics box. Current Opinion in Microbiology. 2006, 9: 10.1016/j.mib.2006.08.010.Google Scholar
- Tettelin H, Masignani V, Cieslewicz MJ, Donati C, Medini D, Ward NL, Angiuoli SV, Crabtree J, Jones AJ, Durkin AS, DeBoy RT, Davidsen TM, Mora M, Scarselli M, y Ros IM, Peterson JD, Hauser CR, Sundaram JP, Nelson WC, Madupu R, Brinkac LM, Dodson RJ, Rosovitz MJ, Sullivan SA, Daugherty SC, Haft DH, Selengut J, Gwinn ML, Zhou L, Zafar N, Khouri H, Radune D, Dimitrov G, Watkins K, OConnor KJB, Smith S, Utterback TR, White O, Rubens EC, Grandi G, Madoff LC, Kasper DL, Telford JL, Wessels MR, Rappuoli R, Fraser CM: Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: Implications for the microbial pan-genome. PNAS. 2005, 102: 13950-13955. 10.1073/pnas.0506758102.PubMed CentralView ArticlePubMedGoogle Scholar
- Chen S, Hung C, Xu J, Reigstad C, Magrini V, Sabo A, Blasiar D, Bieri T, Meyer R, Ozersky P, Armstrong J, Fulton R, Latreille J, Spieth J, Hooton T, Merdis E, Hultgren S, Gordon J: Identification of genes subject to positive selection in uropathogenic strains of
*Escherichia coli*: A comparative genomics approach. PNAS. 2006, 103 (15): 5977-5982. 10.1073/pnas.0600938103.PubMed CentralView ArticlePubMedGoogle Scholar - Rasko D, Rosovitz GMJ, Myers , Mongodin E, Fricke W, Gajer P, Crabtree J, Sebaihia M, Thomson N, Chaudhuri R, Henderson I, Sperandio V, Ravel J: The Pangenome Structure of
*Escherichia coli*: Comparative Genomic Analysis of*E. coli*Commensal and Pathogenic Isolates. Journal of Bacteriology. 2008, 190 (20): 6881-6893. 10.1128/JB.00619-08.PubMed CentralView ArticlePubMedGoogle Scholar - Willenbrock H, Hallin PF, Wassenaar TM, Ussery DW: Characterization of probiotic Escherichia coli isolates with a novel pan-genome microarray. Genome Biology. 2007, 8: 10.1186/gb-2007-8-12-r267.Google Scholar
- Medini D, Donati C, Tettelin H, Masignani V, R R: The microbial pan-genome. Current Opinion in Genetics & Development. 2005, 15: 589-594. 10.1016/j.gde.2005.09.006.View ArticleGoogle Scholar
- Hogg JS, Hu FZ, Janto B, Boissy R, Hayes J, Keefe R, Post JC, Erlich GD: Characterization and modelling of the Haemophilus influenzae core- and supra-genomes based on the complete genomic sequences of Rd and 12 clinical nontypeable strains. Genome Biology. 2007, 8: R103-10.1186/gb-2007-8-6-r103.PubMed CentralView ArticlePubMedGoogle Scholar
- Bunge J, Barger K: Parametric Models for Estimating the Number of Classes. Biometrical Journal. 2008, 50 (6):Google Scholar
- Schwarz G: Estimating the Dimension of a Model. The Annals of Statistics. 1978, 6 (2): 461-464. 10.1214/aos/1176344136.View ArticleGoogle Scholar
- Kuhnert R, Del Rio Villas VJ, Gallagher J, Böhning D: A Bagging -Based Correction for the Mixture Model Estimator of Population Size. Biometrical Journal. 2008, 50 (6):Google Scholar
- Chao A: Estimating the population size for capture-recapture data with unequal catchability. Biometrics. 1987, 43: 783-791. 10.2307/2531532.View ArticlePubMedGoogle Scholar
- The R project. [http://www.r-project.org/]
- NCBI. [http://www.ncbi.nlm.nih.gov/Genomes/]
- Schou-Larsen T, Krogh A: EasyGene – a prokaryotic gene finder that ranks ORFs by statistical significance. BMC Bioinformatics. 2003, 4 (21):Google Scholar
- Nielsen P, Krogh A: Large-scale prokaryotic gene prediction and comparison to genome annotation. Bioinformatics. 2005, 21: 4322-4329. 10.1093/bioinformatics/bti701.View ArticlePubMedGoogle Scholar
- Wang JZ, Lindsay BG: A penalized Nonparametric Maximum Likelihood Approach to Species Richness Estimation. Journal of the American Statistical Association. 2005, 100 (471):Google Scholar
- Seber GAF: The Estimation of Animal Abundance and Related Parameters. 1973, Hafner PressGoogle Scholar
- El-Khoratzaty MN, Imrey PB, Koch GG, Wells B: Estimating the Total Number of Events with Data from Multiple-Record Systems: A Review of Methodological Strategies. Internation Statistical Review. 1977, 45: 129-157. 10.2307/1403312.View ArticleGoogle Scholar
- Koonin E, Wolf Y: Genomics of bacteria and archaea: the emerging dynamic view of the prokaryotic world. Nucleic Acids Research. 2008, 36 (21):Google Scholar
- Tettelin H, Riley D, Cattuto C, Medini D: Comparative genomics: the bacterial pan-genome. Current Opinions in Microbiology. 2008, 12: 472-477. 10.1016/j.mib.2008.09.006.View ArticleGoogle Scholar
- Schloss PD, Handelsman J: A statistical toolbox for metagenomics: assessing functional diversity in microbial communities. BMC Bioinformatics. 2008, 9:Google Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.