Background

With the considerable growth of available nucleotide sequence data over the last decade, integrated and flexible analytical tools have become a necessity. In particular, in the field of population genetics, there is a strong need for automated and reliable procedures to conduct repeatable and rapid polymorphism analyses, coalescent simulations, data manipulation and estimation of demographic parameters under a variety of scenarios.

Results

In this context, we present EggLib (Evolutionary Genetics and Genomics Library), a flexible and powerful C++/Python software package providing efficient and easy to use computational tools for sequence data management and extensive population genetic analyses on nucleotide sequence data. EggLib is a multifaceted project involving several integrated modules: an underlying computationally efficient C++ library (which can be used independently in pure C++ applications); two C++ programs; a Python package providing, among other features, a high level Python interface to the C++ library; and theegglibscript which provides direct access to pre-programmed Python applications.

Conclusions

EggLib has been designed aiming to be both efficient and easy to use. A wide array of methods are implemented, including file format conversion, sequence alignment edition, coalescent simulations, neutrality tests and estimation of demographic parameters by Approximate Bayesian Computation (ABC). Classes implementing different demographic scenarios for ABC analyses can easily be developed by the user and included to the package. EggLib source code is distributed freely under the GNU General Public License (GPL) from its websitehttp://egglib.sourceforge.net/ where a full documentation and a manual can also be found and downloaded.

C++ library

egglib-cppis a fully object-oriented C++ library, meaning that all code is organized in classes, allowing to organize programs in a modular form.egglib-cppimplements tasks related to sequence data storage, simulation and polymorphism analysis. The main classes of the library pertain respectively to aligned and non-aligned set of sequences (AlignandContainer, respectively), polymorphism analysis (NucleotideDiversity, MicrosatelliteDiversity, HaplotypeDiversity, Fstatistics, HFStatistics) and a coalescent simulator allowing recombination. These classes constitute the backbone of the whole package.

egglib-cppis available for use in native C++ applications as an independent C++ library package. The functionalities of the C++ library however are available for use in Python applications through the high-level Python interface that is described in the next section.

Python package

egglib-pyis the Python package of EggLib and fulfills several goals, that are reflected by the seven modules it contains (Figure1). The first module,binding, provides a Python interface ofegglib-cppthrough a C++-to-Python binding. In most applications, it will not be necessary to handlebindingdirectly, since the other modules are built on top ofbinding.

The module data contains data storage classes that are likely to be central to most usage of EggLib. The classes dedicated to the management of sequence data (ContainerandAlign) inherit the C++ implementation of their counterparts inegglib-cppbut also incorporates a wide range of interface and extension methods. As a result, the Python versions ofContainerandAlignprovide a wide range of functionality transparent with respect to the underlying implementation, such as FASTA import/export, introspection, data access and modification, filtering or extracting. In addition,Alignprovides several methods for polymorphism analysis. Additional pure Python classes allow to handle microsatellite data (with import/export functions), annotated sequences (incorporating a GenBank Flat File Format parser/formatter) and phylogenetic trees (incorporating a Newick parser/formatter). Similarly to sequence sets, these classes support a wide array of data access, manipulation and edition operations as methods.

The modulesimulimplements coalescent simulations. Since the underlying coalescent simulator is highly flexible, model specifications are passed through two classes (one holding options relative to the demographic model and the other holding options relative to the mutation model) rather than a long, tedious and error-prone argument list. This object-oriented design allows to readily specify complex scenarios.

The moduletoolsincludes pure Python components for sequence data manipulation (such as coding sequence translation under various models, open reading frame prediction or alignment concatenation) and extra utilities. The module wrappers provides interfaces to popular applications frequently used by population geneticists such as BLAST + [6], ClustalW [7], MUSCLE [8], PhyML [9], and codeml [10].

Thefitmodelmodule comprises all the classes pertaining to the adjustment of demographic models using Approximate Bayesian Computation (ABC), an increasingly used methodology for demographic inference [11]. Briefly, the principle of ABC is: 1) assume a demographic model that is determined by a set of parameters to be estimated, 2) draw random parameter values from a prior distribution, 3) for each set of parameters, perform a simulation under the assumed demographic model, 4) compare a set of summary statistics computed from the simulated data set to an observed data set, and 5) determine the posterior (estimated) distribution of parameters based on the fit of simulated summary statistics to the observed summary statistics [12,13]. An in-depth description of ABC foundations and methodologies is available in [11]. Compared to existing ABC software, the aim of EggLib is to provide the user with maximal freedom for designing demographic models, statistical priors, and sets of summary statistics.fitmodelhas pre-defined models, priors and statistics sets that can be replaced by user-defined classes leveraging all potentialities of EggLib (and beyond). In contrast, low-level analytical steps are implemented in C++ using the GNU Scientific Library in order to maximize performance. Since modern ABC analyses potentially generates very large data sets, files are not fully imported in memory, allowing to accomplish this step using standard workstation computers.

The last module,utils, contains components supporting the interactive commands described hereafter.

Interactive commands

A program provided in theegglib-pydistribution allows to run directly (from a command terminal) a set of pre-programmed commands. These commands are only a subset of what could be achieved with Python programs using EggLib, but they provide a set of immediately available applications. Commands broadly fall into five categories: 1) BLAST-based tools, 2) primer-designing tools, 3) data file conversion or edition, 4) tree manipulation, and 5) ABC estimation of demographic parameters. The latter are the most elaborate. In particular, the commandabc_sampleperforms the steps of coalescent simulation and computation of summary statistics (see the short description of the ABC above), andabc_fitperforms the step of estimation of the posterior parameter distribution. In addition, several commands allow to compute marginal or joint posterior distributions, generate graphical plots (using Matplotlib; [14]) and perform posterior simulations using the fitted model as a null model.

Documentation

The documentation of the C++ classes was generated using Doxygen [15] and that of all Python code was generated using Sphinx [16]. Both Doxygen and Sphinx generate navigable HTML documentation. In addition, a general introduction to EggLib, a manual and description pages have been generated using Sphinx. The whole documentation contains the manual and documentation of both the C++ and Python parts and is available for browsing fromhttp://egglib.sourceforge.net/ and for downloading from the project download page.

Feature overview

An overview of the different type of services provided by Egglib is given in Table1 and shows whether those features are implemented in other frequently used software. Whereas no general class of features is exclusive to EggLib, the point of EggLib is to bring together most tasks routinely performed in population genetics analyses within a single framework, whenever possible as built-in features (which are efficient and convenient to use). EggLib also brings specific features, such as missing data management and several coalescent simulation options (mutation bias, explicit position of markers, diploid model with selfing). Missing data (and alignment gaps) are a recurrent concern of empirical studies. EggLib can perform nucleotide diversity analyses allowing a given proportion of missing data (the statistics are computed on the remaining data). The power of this approach to detect polymorphic sites that would be otherwise ignored is depicted in Figure2.

Usage of egglib-py

The programming interface of the Python package EggLib was designed to be intuitive, simple to use, and to allow fast development of scripts automating population genetics analyses. This was done by providing high-level interface layers above components implemented in C++ and internalizing much of the complexity. We present a simple example to demonstrate how a data set comprising an arbitrary number of loci can be analyzed in a compact and readable fashion by combining Python and EggLib simple syntax (Figure3). The example's comments describe what each block achieves (a full documentation of EggLib's classAlignand simul module is available in the online reference manual). Here, we will point out the parts of the code exploiting EggLib's potentialities. Line 16 (align = egglib.Align(locus)) creates an alignment instance. The user is only required to specify the name of the FASTA file containing the alignment (locus). Line 19 (pol = align.polymorphism()) performs a polymorphism analysis with default settings, that correspond to the standard approach. One of the options (not shown here, see the reference manual) allows to support missing data (see above and Figure2). The returned value,pol, is a dictionary (associative array), that allows straightforward access to computed statistics. Whenever several populations and/or outgroup sequences are present in the alignment, between-population and outgroup-based statistics will be automatically computed. Finally, lines 44-46 demonstrate the usage of the coalescent simulator. In this example, the simplest possible model is used: a single constant-sized population with an infinite-site model of mutation. Three steps are performed: creation of aCoalesceParamSetinstance (specifying the number of samples; line 44), creation of aFiniteAlleleMutatorinstance (specifying the mutator type and the rate of mutation; line 45), and, finally, call to thecoalescefunction that returns a list of Align instances (line 46). The advantage of this three-step syntax for configuring coalescent simulation is that it can accommodate both simple models (as the one used here) and more complex scenarios exploiting all potentialities of the coalescent simulator.

User-defined ABC model

Several commands accessible from the command line utilities of EggLib allow one to perform ABC analysis using command-line tools. However, the set of pre-defined models cannot be exhaustive and one of our aims is to allow using all EggLib functionalities to design any possible demographic model.

The model presented in Figure4 is an arbitrary example of model that is not available in thefitmodelmodule. This model is depicted at the top of the figure. It has five different parameters:THETA(θ in the picture),DATE, SIZE, MIGR1andMIGR2. This model can be viewed as a double, simultaneous domestication from two partially isolated stocks (time runs from top to bottom).DATEis the age of the domestication event, the migration parameters specify exchange ratesMIGR1andMIGR2between pairs of populations andSIZEgives the relative size of cultivated populations.

The code at the bottom of Figure4 shows the implementation of this model within the EggLib framework. Note that ABC models should conform to a few requirements: they must be formalized as a class; they must define their name and the names of all parameters; their constructor must accept at least one argument specifying whether recombination must be implemented (and deal with it appropriately), but it can accept more arguments; and they must contain a generate method which specifies the body of the model implementation. This very piece of code can be used in conjunction tofitmodel(within a Python script), but it can also be used to add this model to the list of models available through the interactive commandabc_sample. Custom priors and sets of summary statistics can be incorporated using a similar system, although currently (as of version 2.1.2)abc_sampledoes not currently support run-time addition of priors and summary statistics (such support is planned).

Thegeneratemethod is the hook that connects the model to the rest of the ABC framework. It must take two arguments: a sample configuration and a set of parameter values drawn from the prior (thefitmodeldocumentation provides details of the exact format of these data). Thegeneratemethod must return simulated data sets (using a type defined infitmodel). Apart from these constraints, the user has full freedom with regard to what is actually done for generating the data set. Obviously, all potentialities of the coalescent simulator incorporated within EggLib are allowed. Furthermore, all forms of post-processing operations are not only possible, but easy to implement usingegglib-py. For example, one can readily include error rates or sampling or ascertainment biases and set them as model parameters.

Performance

The running time and maximum memory usage of programs performing common population genetics operations using EggLib compared with alternatives (whenever available) is shown on Tables2,3,4 and5. All tests were run on a laptop computer and (except for coalescent simulations) were repeated 10 times. Tests were performed using EggLib 2.1.2, Biopython 1.58, libsequence 1.7.4 [26], analysis 0.8.1 (containingcompute, polydNdSandrsq) [26], the version of ms updated December 11, 2009 [20], coasim-python 1.3 [21], msABC 20111219 [24] and ABCreg 2009-07-30 [25] (which were all the latest available versions at time of testing). For coalescent simulations (including in ABC), EggLib was set to use at most 4 processor cores.

EggLib is comparatively more efficient than Biopython for importing large FASTA files (Table2). TheAlignclass of EggLib is slightly more efficient thanAlignIOof BioPython for importing a large alignment. For importing data files representing the wholeOryza sativagenome, theContainerclass of EggLib is much more efficient thanSeqIOof Biopython (EggLib is able to import these two files fully in memory in a few seconds and with a limited memory overhead: the memory use is hardly larger than the file size). However, the difference between EggLib and Biopython reflects a difference in paradigm (import all file at once for EggLib, and read sequence one at a time for Biopython).

The comparison of an EggLib script for analyzing polymorphism with programs developed using the C++ libsequence library (computefor standard statistics,polydNdSfor coding sequence statistics andrsqfor linkage disequilibrium) shows that skipping unneeded statistics can significantly fasten the analysis. For analyzing a single alignment, libsequence programs are better, but for processing many alignments in a row a single loop using EggLib is more efficient. In EggLib, the linkage disequilibrium analysis is comparatively more efficient, and the coding sequence analysis (based on the wrapping of Bio++) is comparatively less efficient.

The comparison of theeggcoal, msand CoaSim simulators shows thatmsis consistently and significantly the fastest and the least memory-demanding (Table4).eggcoallies betweenmsand CoaSim. EggLib has a generic design that makes it difficult to maximize performance, explaining part of the discrepancy. However, we believe that future versions will improve performance, especially thanks to improved implementation of recombination and multithreading scheme that are currently planned.

We compared the performances of EggLib commands for ABC to the very efficient programsmsABC(for the simulation phase) andABCreg(for the analysis phase). We used two different summary statistics sets: SDZ (number of polymorphic sites, Tajima'sDand Fay and Wu'sH) and SFS (site frequency spectrum with 8 categories). The SFS was available only with EggLib. EggLib was used through the interactive commandsabc_sampleandabc_fit. We found that EggLibabc_samplewas slower thanmsABCand used more memory (chiefly because of Python-level multithreading). This is explained by the performance of the original ms program (see above), that was efficiently leveraged inmsABC. However, the overhead tends to be decreased compared with theeggcoal/mscomparison presented before, showing that the EggLib integration does not worsen performance. We therefore expect that future improvements of the coalescent simulator will bring EggLib closer to the level ofmsABC. For the analysis step, a large data file of 5,000,000 samples was imported and analyzed. We observed that this step of the ABC procedure was not limiting in running time (compared to the simulation step) but could be limiting in memory use. Therefore we followed a strategy favoring data access from file, which is relatively slower but more memory efficient. EggLib andABCregare therefore complementary regarding the speed/memory balance.

Prospects

EggLib is under active development and we expect new features to be added in the future. Our current routes for improving the package include: improving the performance of the coalescent simulator thanks to a new design of the recombination process and an improved parallelization scheme; easing the definition by users of custom ABC models, sets of summary statistics and priors using automated helpers; improving the performance of the ABC framework by internalizing replications within the C++ layer and removing unnecessary steps (such asAlignconversion) without interfering with the general flexibility of the framework; putting a special effort in documentation, especially by providing tutorials besides the complete reference manual.

1. Schuster S. Next-generation sequencing transforms today's biology. Nat Methods. 2008;5:16–18. doi: 10.1038/nmeth1156. [DOI] [PubMed] [Google Scholar]
2. Stajich JE, Block D, Boulez K, Brenner SE, Chervitz SA, Dagdigian C, Fuellen G, Gilbert JGR, Korf I, Lapp H, Lehvaslaiho H, Matsalla C, Mungall CJ, Osborne BI, Pocock MR, Schattner P, Senger M, Stein LD, Stupka E, Wilkinson MD, Birney E. The bioperl toolkit: Perl modules for the life sciences. Genome Res. 2002;12:1611–1618. doi: 10.1101/gr.361602. [DOI] [PMC free article] [PubMed] [Google Scholar]
3. Holland RCG, Down TA, Pocock M, Prlic A, Huen D, James K, Foisy S, Drager A, Yates A, Heuer M, Schreiber MJ. BioJava: an open-source framework for bioinformatics. Bioinformatics. 2008;24:2096–2097. doi: 10.1093/bioinformatics/btn397. [DOI] [PMC free article] [PubMed] [Google Scholar]
4. Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, Friedberg I, Hamelryck T, Kauff F, Wilczynski B, de Hoon MJL. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009;25:1422–1423. doi: 10.1093/bioinformatics/btp163. [DOI] [PMC free article] [PubMed] [Google Scholar]
5. Bassi S. A primer on Python for life science researchers. PLoS Comput Biol. 2007;3:2052–2057. doi: 10.1371/journal.pcbi.0030199. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST plus: architecture and applications. BMC Bioinformatics. 2009;10:421. doi: 10.1186/1471-2105-10-421. [DOI] [PMC free article] [PubMed] [Google Scholar]
7. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, Thompson JD, Gibson TJ, Higgins DG. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23:2947–2948. doi: 10.1093/bioinformatics/btm404. [DOI] [PubMed] [Google Scholar]
8. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–1797. doi: 10.1093/nar/gkh340. [DOI] [PMC free article] [PubMed] [Google Scholar]
9. Guindon S, Dufayard JF, Hordijk W, Lefort V, Gascuel O. PhyML: fast and accurate phylogeny reconstruction by maximum likelihood. Infect Genet Evol. 2009;9:384–385. [Google Scholar]
10. Yang ZH. PAML 4: Phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24:1586–1591. doi: 10.1093/molbev/msm088. [DOI] [PubMed] [Google Scholar]
11. Beaumont MA. Approximate Bayesian computation in evolution and ecology. Annu Rev Ecol Evol Syst. 2010;41:379–406. doi: 10.1146/annurev-ecolsys-102209-144621. [DOI] [Google Scholar]
12. Beaumont MA, Zhang W, Balding DJ. Approximate Bayesian computation in population genetics. Genetics. 2002;162:2025–2035. doi: 10.1093/genetics/162.4.2025. [DOI] [PMC free article] [PubMed] [Google Scholar]
13. Marjoram P, Tavaré S. Modern computational approaches for analysis molecular genetics variation data. Nat Rev Genet. 2006;7:759–770. doi: 10.1038/nrg1961. [DOI] [PubMed] [Google Scholar]
14. Hunter JD. Matplotlib: a 2D graphics environment. Computing in Science & Engineering. 2007;9:90–95. [Google Scholar]
15. van Heesch D. Doxygen: generate documentation from source code.http://www.stack.nl/~dimitri/doxygen/index.html
16. Sphinx: Python documentation generation.http://sphinx.pocoo.org/
17. Knight R, Maxwell P, Birmingham A, Carnes J, Caporaso JG, Easton BC, Eaton M, Hamady M, Lindsay H, Liu ZZ, Lozupone C, McDonald D, Robeson M, Sammut R, Smit S, Wakefield MJ, Widmann J, Wikman S, Wilson S, Ying H, Huttley GA. PyCogent: a toolkit for making sense from sequence. Genome Biol. 2007;8:R171. doi: 10.1186/gb-2007-8-8-r171. [DOI] [PMC free article] [PubMed] [Google Scholar]
18. Dutheil J, Gaillard S, Bazin E, Glémin S, Ranwez V, Galtier N, Belkhir K. Bio++: a set of C++ libraries for sequence analysis, phylogenetics, molecular evolution and population genetics. BMC Bioinformatics. 2006;7:188. doi: 10.1186/1471-2105-7-188. [DOI] [PMC free article] [PubMed] [Google Scholar]
19. Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003;19:2496–2497. doi: 10.1093/bioinformatics/btg359. [DOI] [PubMed] [Google Scholar]
20. Hudson RR. Generating samples under a Wright-Fisher neutral model of genetic variation. Bioinformatics. 2002;18:337–338. doi: 10.1093/bioinformatics/18.2.337. [DOI] [PubMed] [Google Scholar]
21. Mailund T, Schierup MH, Pedersen CNS, Mechlenborg PJM, Madsen JN, Schauser L. CoaSim: a flexible environment for simulating genetic data under coalescent models. BMC Bioinformatics. 2005;6:252. doi: 10.1186/1471-2105-6-252. [DOI] [PMC free article] [PubMed] [Google Scholar]
22. Cornuet JM, Santos F, Beaumont MA, Robert CP, Marin JM, Balding DJ, Guillemaud T, Estoup A. Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation. Bioinformatics. 2008;24:2713–2719. doi: 10.1093/bioinformatics/btn514. [DOI] [PMC free article] [PubMed] [Google Scholar]
23. Wegmann D, Leuenberger C, Neuenschwander S, Excoffier L. ABCtoolbox: a versatile toolkit for approximate Bayesian computations. BMC Bioinformatics. 2010;11:116. doi: 10.1186/1471-2105-11-116. [DOI] [PMC free article] [PubMed] [Google Scholar]
24. Pavlidis P, Laurent S, Stephan W. msABC: a modification of Hudson's ms to facilitate multi-locus ABC analysis. Mol Ecol Resour. 2010;10:723–727. doi: 10.1111/j.1755-0998.2010.02832.x. [DOI] [PubMed] [Google Scholar]
25. Thornton KR. Automating approximate Bayesian computation by local linear regression. BMC Genet. 2009;10:35. doi: 10.1186/1471-2156-10-35. [DOI] [PMC free article] [PubMed] [Google Scholar]
26. Thornton K. Libsequence: a C++ class library for evolution genetic analysis. Bioinformatics. 2003;22:2325–2327. doi: 10.1093/bioinformatics/btg316. [DOI] [PubMed] [Google Scholar]
27. St. Onge KR, Källman T, Slotte T, Lascoux M, Palmé AE. Contrasting demographic history and population structure in Capsella rubella and Capsella grandiflora, two closely related species with different mating systems. Mol Ecol. 2011;20:3306–3320. doi: 10.1111/j.1365-294X.2011.05189.x. [DOI] [PubMed] [Google Scholar]
28. Li Y, Stocks M, Hemmila S, Källman T, Zhu H, Zhou Y, Chen J, Liu J, Lascoux M. Demographic histories of four spruce (Picea) species of the Qinghai-Tibetan Plateau and neighboring areas inferred from multiple nuclear loci. Mol Biol Evol. 2010;27:1001–1014. doi: 10.1093/molbev/msp301. [DOI] [PubMed] [Google Scholar]
29. Li ZH, Zhang QA, Liu JQ, Källman T, Lascoux M. The Pleistocene demography of an alpine juniper of the Qinghai-Tibetan Plateau: tabula rasa, cryptic refugia or something else? J Biogeogr. 2011;38:31–43. doi: 10.1111/j.1365-2699.2010.02400.x. [DOI] [Google Scholar]
30. De Mita S, Chantret N, Loridon K, Ronfort J, Bataillon T. Molecular adaptation in flowering and symbiotic recognition pathways: insights from patterns of polymorphism in the legume Medicago truncatula. BMC Evol Biol. 2011;11:229. doi: 10.1186/1471-2148-11-229. [DOI] [PMC free article] [PubMed] [Google Scholar]
31. Dereeper A, Nicolas S, Le Cunff L, Bacilieri R, Doligez A, Peros JP, Ruiz M, This P. SNiPlay: a web-based tool for detection, management and analysis of SNPs, Application to grapevine diversity projects. BMC Bioinformatics. 2011;12:134. doi: 10.1186/1471-2105-12-134. [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials