Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction

Data import

Poppr allows import of data in several formats for dominant/codominant, haploid/diploid and geographic data. The R package adegenet, that defines the genind data structure that poppr utilizes, allows support for importing data natively from Structure, Genetix, Genepop, and Fstat. While these formats are very common and widely supported, these do not allow for import of geographic and/or regional data. Furthermore, adegenet will only handle diploids with this format, though manual import is possible. To aid in importing data, poppr has newly added the function read.genalex(), to read data from GenAlEx formatted text files into the genind data object of the package adegenet (Jombart, 2008; Jombart & Ahmed, 2011; Peakall & Smouse, 2006). GenAlEx is a popular add-in for Microsoft Excel that can handle data including codominant/dominant and haploid/diploid markers as well as geographic and regional data. This function further facilitates the import of haploid, geographic, and regional data.

Transferring data to new formats and manipulating data by hand, such as collapsing data into clones or subsetting data into different hierarchical levels, is tedious, creates redundancy, and can result in lost or misrepresented data. Poppr includes tools to automate such repetitive tasks. Many currently available data formats and software implementations allow analysis of only one or two levels of a population hierarchy. With poppr the user can import a single data set with an unlimited number of hierarchical levels. This is achieved by having the user combine the levels using a common delimiter (e.g., “Year_Country_City”). These combined levels are then used as the defining population factor in the input file and can easily be manipulated within R.

Data analysis

Once data is imported into R, the user can dynamically access and manipulate the population hierarchy with the function splitcombine(), subset the data set by population with popsub(), and check for cloned multilocus genotypes using mlg(). For data sets that include clones, the poppr function clonecorrect() will censor clones with respect to any level of a population hierarchy. In the case of missing data we use the commonly implemented, most parsimonious approach of treating missing states as novel alleles. This inherently makes analysis sensitive to missing data and genotyping error, but the user has tools available such as missingno() to filter out missing data at a per-individual or per-locus level. The user can also decide how uninformative loci (e.g., alleles occurring at minor frequencies; monomorphic loci; fixed heterozygous loci) are treated using the function informloci(). Thus, the user can specify a frequency for removal of uninformative loci. The user is encouraged to conduct analysis with and without missing data/uninformative loci to assess sensitivity to these issues when making inferences. A full list of functions available in poppr is provided in Table 1.

Typical analyses in poppr start with summary statistics for diversity, rarefaction, evenness, MLG counts, and calculation of distance measures such as Bruvo’s distance, providing a suitable stepwise mutation model appropriate for microsatellite markers (Bruvo et al., 2004). Poppr will define MLGs in your data set, show where they cross populations, and can produce graphs and tables of MLGs by population that can be used for further analysis with the R package vegan (Oksanen et al., 2013). Many of the diversity indices calculated by the vegan function diversity() are useful in analyzing the diversity of partially clonal populations. For this reason, poppr features a quick summary table (Table 2) that incorporates these indices along with the index of association, I_A (Brown, Feldman & Nevo, 1980; Smith et al., 1993), and its standardized form, ${\bar{r}}_d$, which accounts for the number of loci sampled (Agapow & Burt, 2001). Both measures of association can detect signatures of multilocus linkage and values significantly departing from the null model of no linkage among markers are detected via permutation analysis utilizing one of four algorithms described in Table 3 (Agapow & Burt, 2001). The user can specify the number of samples taken from the observed data set to obtain the null distribution expected for a randomly mating population. Detailed examples of these analyses can be found in the poppr manual.

Visualizations

Poppr generates bar charts of MLG counts found within each population of your data set (Fig. 1). Histograms with rug plots for I_A and ${\bar{r}}_d$ allow visual assessment of the quality of the distribution derived from resampling to see if a higher number of replications are necessary (Fig. 2). Poppr automatically produces custom minimum spanning networks for Bruvo’s or other distances using Prim’s algorithm, as implemented in the package igraph (Csardi & Nepusz, 2006), with the functions bruvo.msn() for Bruvo’s distance (Fig. 3) and poppr.msn() for any distance matrix. The combination of data structures from adegenet and igraph allow graphing that is color coded by population with vertices grouped by MLG (Jombart, 2008; Jombart & Ahmed, 2011; Csardi & Nepusz, 2006). The user can further customize the appearance of the graphs directly within R by utilizing igraph’s layout() and plot() functions. Poppr also includes visualization of dendrograms using UPGMA (Schliep, 2011) and Neighbor-Joining (Paradis, Claude & Strimmer, 2004) algorithms with bootstrap support for Bruvo’s distance using the function bruvo.boot() (Fig. 4). Neither graphing of minimum spanning networks or dendrograms with bootstrap support are currently possible for populations in any other R packages.

Performance

Most of the functions in Poppr were written and optimized for performance in R and are available for inspection and/or download at https://github.com/grunwaldlab/poppr. Algorithms of ≥ O(n²) complexity were written in the byte-compiled C language to optimize runtime performance.

For comparisons of I_A and Bruvo’s distance, we utilized the data set nancycats (237 diploid individuals genotyped at nine microsatellite loci) from the adegenet package. Calculations were run independently 10 times and then averaged. Bruvo’s distance was calculated on a machine with OSX 10.8.4 and a 2.9 GHz intel processor. The I_A and ${\bar{r}}_d$ calculations were performed on a machine with OSX 10.5.8 and a 2.4 GHz intel processor due to the inability of the software multilocus to work on any later version of OSX.

Example data sets

Along with example data sets preloaded into the adegenet package, poppr offers two data sets of clonal populations. The data set Aeut is comprised of 187 isolates of Aphanomyces euteiches genotyped over 56 AFLP loci with defined populations and subpopulations (Grünwald & Hoheisel, 2006). The partial_clone data set contains 50 simulated individuals genotyped over 10 microsatellite loci produced with the software SimuPOP v.1.0.8 (Peng & Amos, 2008). Each data set can be loaded into R using the commands data(Aeut) and data(partial_clone), respectively.

Citation of methods implemented in poppr

Several of the methods implemented in poppr are described elsewhere. Users should refer to the original publications for interpretations and citation. See Table 4 for a full list of citations. As with any R package, users should always cite the R Core Team (2013).

New functionalities

Poppr implements several new functionalities. As of this writing, aside from poppr, there exist two programs that calculate I_A: LIAN (Haubold & Hudson, 2000) and multilocus (Agapow & Burt, 2001). LIAN can calculate I_A for haploid data and is only available online or for ∗nix systems with a C compiler such as OSX and Linux (Haubold & Hudson, 2000). Multilocus implemented ${\bar{r}}_d$, a novel correction for I_A, but is no longer supported (Agapow & Burt, 2001). Multilocus will only calculate index values for one data set at a time and LIAN requires the user to structure the data set with populations in contiguous blocks to analyze multiple populations within a single file. Thus poppr provides significant improvements for calculation of linkage disequilibrium, and handles both haploid and diploid data, works on all major operating systems, and is capable of batch analysis of multiple files and multiple populations defined within a file including the possibility of clone correction and sub-setting. A comparison of the capabilities of these programs are summarized in Table 5.

To test significance for I_A and ${\bar{r}}_d$, poppr offers four permutation algorithms. Each one will randomly shuffle data at each locus, effectively unlinking the loci. The algorithm previously utilized by mutlilocus is included. The multilocus-style algorithm shuffles genotypes, maintaining the associations between alleles at each locus (Agapow & Burt, 2001). More appropriately, alleles are expected to assort independently in panmictic populations. Poppr thus provides three new algorithms for permutation that allow for independent allele assortment at each locus. The default algorithm permutes the alleles at each locus and the remaining two will randomly sample alleles from a multinomial distribution parametrically and non-parametrically (Weir, 1996). Details of these algorithms are presented in Table 3. Because the index of association is calculated using a binary measure of dissimilarity, we have also made this available as a distance measure called diss.dist(). This pairwise distance is based on the percent allelic differences.

Poppr also newly implements Bruvo’s genetic distance that utilizes a stepwise mutation model appropriate for microsatellite data (Bruvo et al., 2004). While this distance is implemented in the program GenoDive (Meirmans & Van Tienderen, 2004) and the R package polysat (Clark & Jasieniuk, 2011), there are a few caveats with these two implementations. GenoDive is closed-source, and only implemented in OSX. Both poppr and polysat are open-source and available on all platforms, but polysat, being optimized for polyploid individuals with ambiguous allelic dosage, is inappropriate for analyzing diploids. Polysat will collapse homozygous individuals into a single allele and attempt to infer the second allelic state in comparison with heterozygous individuals. Since haploid and diploid individuals show clear allelic dosage, this procedure creates a bias misrepresenting the true distance. Not only is poppr not subject to this bias, but it also newly introduces bootstrap support for this distance as shown in Fig. 4.

Performance

Poppr reduces the amount of intermediate files and repetitive tasks needed for basic population genetic analyses and implements computationally intensive functions, such as Bruvo’s distance and the index of association in C to improve performance. The polysat package calculation of Bruvo’s distance took 58.3 s on average whereas poppr’s calculation was over 190 times faster, averaging 0.3 s (Table 6). For calculation of I_A and ${\bar{r}}_d$ with 100 permutations and Nei’s genotypic diversity (Nei, 1978), multilocus required around 9.12 min on average, as compared to 13.4 s with poppr.