Package 'poppr'

Description

Poppr provides tools for population genetic analysis that include genotypic diversity measures, genetic distances with bootstrap support, native organization and handling of population hierarchies, and clone correction.

To cite poppr, please use citation("poppr"). When referring to poppr in your manuscript, please use lower case unless it occurs at the beginning of a sentence.

Details

This package relies on the adegenet package. It is built around the genind and genlight object. Genind objects store genetic information in a table of allele frequencies while genlight objects store SNP data efficiently by packing binary allele calls into single bits. Poppr has extended these object into new objects called genclone and snpclone, respectively. These objects are designed for analysis of clonal organisms as they add the @mlg slot for keeping track of multilocus genotypes and multilocus lineages.

Documentation

Documentation is available for any function by typing ?function_name in the R console. Detailed topic explanations live in the package vignettes:

Essential functions for importing and manipulating data are detailed within the Data import and manipulation vignette, details on algorithms used in poppr are within the Algorithms and equations vignette, and details for working with multilocus genotypes are in Multilocus Genotype Analysis.

Examples of analyses are available in a primer written by Niklaus J. Grünwald, Zhian N. Kamvar, and Sydney E. Everhart at https://grunwaldlab.github.io/Population_Genetics_in_R/.

Getting help

If you have a specific question or issue with poppr, feel free to contribute to the google group at https://groups.google.com/d/forum/poppr. If you find a bug and are a github user, you can submit bug reports at https://github.com/grunwaldlab/poppr/issues. Otherwise, leave a message on the groups. Personal emails are highly discouraged as they do not allow others to learn.

Functions in poppr

Below are descriptions and links to functions found in poppr. Be aware that all functions in adegenet are also available. The functions are documented as:

• function_name() (data type) - Description

Where ‘data type’ refers to the type of data that can be used:

Data import/export

• getfile() (x) - Provides a quick GUI to grab files for import

• read.genalex() (x) - Reads GenAlEx formatted csv files to a genind object

• genind2genalex() (m) - Converts genind objects to GenAlEx formatted csv files

• genclone2genind() (m) - Removes the @mlg slot from genclone objects

• as.genambig() (m) - Converts genind data to polysat's genambig data structure.

• bootgen2genind() (x) - see aboot() for details)

Data Structures

Data structures "genclone" (based off of adegenet's genind) and "snpclone" (based off of adegenet's genlight for large SNP data sets). Both of these data structures are defined by the presence of an extra MLG slot representing multilocus genotype assignments, which can be a numeric vector or a MLG class object.

• genclone - Handles microsatellite, presence/absence, and small SNP data sets

• snpclone - Designed to handle larger binary SNP data sets.

• MLG - An internal class holding a data frame of multilocus genotype assignments that acts like a vector, allowing the user to easily switch between different MLG definitions.

• bootgen - An internal class used explicitly for aboot() that inherits the gen-class virtual object. It is designed to allow for sampling loci with replacement.

• bruvomat - An internal class designed to handle bootstrapping for Bruvo's distance where blocks of integer loci can be shuffled.

Data manipulation

• as.genclone() (m) - Converts genind objects to genclone objects

• missingno() (m) - Handles missing data

• clonecorrect() (m | s) - Clone-censors at a specified population hierarchy

• informloci() (m) - Detects and removes phylogenetically uninformative loci

• popsub() (m | s) - Subsets genind objects by population

• shufflepop() (m) - Shuffles genotypes at each locus using four different shuffling algorithms

• recode_polyploids() (m | x) - Recodes polyploid data sets with missing alleles imported as "0"

• make_haplotypes() (m | s) - Splits data into pseudo-haplotypes. This is mainly used in AMOVA.

• test_replen() (m) - Tests for inconsistent repeat lengths in microsatellite data. For use in bruvo.dist() functions.

• fix_replen() (m) - Fixes inconsistent repeat lengths. For use in bruvo.dist() functions.

Genetic distances

• bruvo.dist() (m) - Bruvo's distance (see also: fix_replen())

• diss.dist() (m) - Absolute genetic distance (see prevosti.dist())

• nei.dist() (m | x) - Nei's 1978 genetic distance

• rogers.dist() (m | x) - Rogers' euclidean distance

• reynolds.dist() (m | x) - Reynolds' coancestry distance

• edwards.dist() (m | x) - Edwards' angular distance

• prevosti.dist() (m | x) - Prevosti's absolute genetic distance

• bitwise.dist() (s) - Calculates fast pairwise distances for genlight objects.

Bootstrapping

• aboot() (m | s | x) - Creates a bootstrapped dendrogram for any distance measure

• bruvo.boot() (m) - Produces dendrograms with bootstrap support based on Bruvo's distance

• diversity_boot() (x) - Generates boostrap distributions of diversity statistics for multilocus genotypes

• diversity_ci() (m | s | x) - Generates confidence intervals for multilocus genotype diversity.

• resample.ia() (m) - Calculates the index of association over subsets of data.

Multilocus Genotypes

• mlg() (m | s) - Calculates the number of multilocus genotypes

• mll() (m | s) - Displays the current multilocus lineages (genotypes) defined.

• nmll() (m | s) - Same as mlg().

• mlg.crosspop() (m | s) - Finds all multilocus genotypes that cross populations

• mlg.table() (m | s) - Returns a table of populations by multilocus genotypes

• mlg.vector() (m | s) - Returns a vector of a numeric multilocus genotype assignment for each individual

• mlg.id() (m | s) - Finds all individuals associated with a single multilocus genotype

• mlg.filter() (m | s) - Collapses MLGs by genetic distance

• filter_stats() (m | s) - Calculates mlg.filter for all algorithms and plots

• cutoff_predictor() (x) - Predicts cutoff threshold from mlg.filter.

• mll.custom() (m | s) - Allows for the custom definition of multilocus lineages

• mll.levels() (m | s) - Allows the user to change levels of custom MLLs.

• mll.reset() (m | s) - Reset multilocus lineages.

• diversity_stats() (x) - Creates a table of diversity indices for multilocus genotypes.

Index of Association Analysis

Analysis of multilocus linkage disequilibrium.

• ia() (m) - Calculates the index of association

• pair.ia() (m) - Calculates the index of association for all loci pairs.

• win.ia() (s) - Index of association windows for genlight objects.

• samp.ia() (s) - Index of association on random subsets of loci for genlight objects.

Population Genetic Analysis

• poppr.amova() (m | s) - Analysis of Molecular Variance (as implemented in ade4)

• poppr() (m | x) - Returns a diversity table by population

• poppr.all() (m | x) - Returns a diversity table by population for all compatible files specified

• private_alleles() (m) - Tabulates the occurrences of alleles that only occur in one population.

• locus_table() (m) - Creates a table of summary statistics per locus.

• rrmlg() (m | x) - Round-robin multilocus genotype estimates.

• rraf() (m) - Round-robin allele frequency estimates.

• pgen() (m) - Probability of genotypes.

• psex() (m) - Probability of observing a genotype more than once.

• rare_allele_correction (m) - rules for correcting rare alleles for round-robin estimates.

• incomp() (m) - Check data for incomparable samples.

Visualization

• imsn() (m | s) - Interactive construction and visualization of minimum spanning networks

• plot_poppr_msn() (m | s | x) - Plots minimum spanning networks produced in poppr with scale bar and legend

• greycurve() (x) - Helper to determine the appropriate parameters for adjusting the grey level for msn functions

• bruvo.msn() (m) - Produces minimum spanning networks based off Bruvo's distance colored by population

• poppr.msn() (m | s | x) - Produces a minimum spanning network for any pairwise distance matrix related to the data

• info_table() (m) - Creates a heatmap representing missing data or observed ploidy

• genotype_curve() (m | x) - Creates a series of boxplots to demonstrate how many markers are needed to represent the diversity of your data.

Datasets

• Aeut() - (AFLP) Oomycete root rot pathogen Aphanomyces euteiches (Grünwald and Hoheisel, 2006)

• monpop() - (SSR) Peach brown rot pathogen Monilinia fructicola (Everhart and Scherm, 2015)

• partial_clone() - (SSR) partially-clonal data simulated via simuPOP (Peng and Amos, 2008)

• Pinf() - (SSR) Potato late blight pathogen Phytophthora infestans (Goss et. al., 2014)

• Pram() - (SSR) Sudden Oak Death pathogen Phytophthora ramorum (Kamvar et. al., 2015; Goss et. al., 2009)

Author(s)

Zhian N. Kamvar, Jonah C. Brooks, Sydney E. Everhart, Javier F. Tabima, Stacy Krueger-Hadfield, Erik Sotka, Niklaus J. Grünwald

Maintainer: Zhian N. Kamvar

References

——— Papers announcing poppr ———

Kamvar ZN, Tabima JF, Grünwald NJ. (2014) Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2:e281 doi:10.7717/peerj.281

Kamvar ZN, Brooks JC and Grünwald NJ (2015) Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Front. Genet. 6:208. doi:10.3389/fgene.2015.00208

——— Papers referencing data sets ———

Grünwald, NJ and Hoheisel, G.A. 2006. Hierarchical Analysis of Diversity, Selfing, and Genetic Differentiation in Populations of the Oomycete Aphanomyces euteiches. Phytopathology 96:1134-1141 doi: doi:10.1094/PHYTO-96-1134

SE Everhart, H Scherm, (2015) Fine-scale genetic structure of Monilinia fructicola during brown rot epidemics within individual peach tree canopies. Phytopathology 105:542-549 doi: doi:10.1094/PHYTO-03-14-0088-R

Bo Peng and Christopher Amos (2008) Forward-time simulations of nonrandom mating populations using simuPOP. bioinformatics, 24 (11): 1408-1409.

Goss, Erica M., Javier F. Tabima, David EL Cooke, Silvia Restrepo, William E. Fry, Gregory A. Forbes, Valerie J. Fieland, Martha Cardenas, and Niklaus J. Grünwald. (2014) "The Irish potato famine pathogen Phytophthora infestans originated in central Mexico rather than the Andes." Proceedings of the National Academy of Sciences 111:8791-8796. doi: doi:10.1073/pnas.1401884111

Kamvar, Z. N., Larsen, M. M., Kanaskie, A. M., Hansen, E. M., & Grünwald, N. J. (2015). Spatial and temporal analysis of populations of the sudden oak death pathogen in Oregon forests. Phytopathology 105:982-989. doi: doi:10.1094/PHYTO-12-14-0350-FI

Goss, E. M., Larsen, M., Chastagner, G. A., Givens, D. R., and Grünwald, N. J. 2009. Population genetic analysis infers migration pathways of Phytophthora ramorum in US nurseries. PLoS Pathog. 5:e1000583. doi: doi:10.1371/journal.ppat.1000583

Description

Calculate a dendrogram with bootstrap support using any distance applicable to genind or genclone objects.

Usage

aboot(
  x,
  strata = NULL,
  tree = "upgma",
  distance = "nei.dist",
  sample = 100,
  cutoff = 0,
  showtree = TRUE,
  missing = "mean",
  mcutoff = 0,
  quiet = FALSE,
  root = NULL,
  ...
)

Arguments

Details

This function automates the process of bootstrapping genetic data to create a dendrogram with bootstrap support on the nodes. It will randomly sample with replacement the loci of a gen (genind/genpop) object or the columns of a numeric matrix, assuming that all loci/columns are independent. The process of randomly sampling gen objects with replacement is carried out through the use of an internal class called bootgen. This is necessary due to the fact that columns in the genind matrix are defined as alleles and are thus interrelated. This function will specifically bootstrap loci so that results are biologically relevant. With this function, the user can also define a custom distance to be performed on the genind or genclone object. If you have a data frame-like object where all of the columns are independent or pairs of columns are independent, then it may be simpler to use ape::boot.phylo() to calculate your bootstrap support values.

the strata argument

There is an argument called strata. This argument is useful for when you want to bootstrap by populations from a adegenet::genind() object. When you specify strata, the genind object will be converted to adegenet::genpop() with the specified strata.

Value

an object of class ape::phylo().

Note

prevosti.dist() and diss.dist() are exactly the same, but diss.dist() scales better for large numbers of individuals (n > 125) at the cost of required memory.

missing data

Missing data is not allowed by many of the distances. Thus, one of the first steps of this function is to treat missing data by setting it to the average allele frequency in the data set. If you are using a distance that can handle missing data (Prevosti's distance), you can set missing = "ignore" to allow the distance function to handle any missing data. See missingno() for details on missing data.

Bruvo's Distance

While calculation of Bruvo's distance is possible with this function, it is optimized in the function bruvo.boot().

References

Kamvar ZN, Brooks JC and Grünwald NJ (2015) Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Frontiers in Genetics 6:208. doi:10.3389/fgene.2015.00208

See Also

nei.dist() edwards.dist() rogers.dist() reynolds.dist() prevosti.dist() diss.dist() bruvo.boot() ape::boot.phylo() adegenet::dist.genpop() dist() bootgen2genind() bootgen

Examples

data(nancycats)
nan9 <- popsub(nancycats, 9)

set.seed(9999)
# Generate a tree using nei's distance
neinan <- aboot(nan9, dist = nei.dist)

set.seed(9999)
# Generate a tree using custom distance
bindist <- function(x) dist(tab(x), method = "binary")
binnan <- aboot(nan9, dist = bindist)

## Not run: 
# Distances from other packages.
#
# Sometimes, distance functions from other packages will have the constraint
# that the incoming data MUST be genind. Internally, aboot uses the 
# bootgen class ( class?bootgen ) to shuffle loci, and will throw an error
# The function bootgen2genind helps fix that. Here's an example of a function
# that expects a genind class from above
bindist <- function(x){
  stopifnot(is.genind(x))
  dist(tab(x), method = "binary")
}
#
# Fails:
# aboot(nan9, dist = bindist)
## Error: is.genind(x) is not TRUE
#
# Add bootgen2genind to get it working!
# Works:
aboot(nan9, dist = function(x) bootgen2genind(x) %>% bindist)

# AFLP data
data(Aeut)

# Nei's distance
anei <- aboot(Aeut, dist = nei.dist, sample = 1000, cutoff = 50)

# Rogers' distance
arog <- aboot(Aeut, dist = rogers.dist, sample = 1000, cutoff = 50)

# This can also be run on genpop objects
strata(Aeut) <- other(Aeut)$population_hierarchy[-1]
Aeut.gc <- as.genclone(Aeut)
setPop(Aeut.gc) <- ~Pop/Subpop
Aeut.pop <- genind2genpop(Aeut.gc)
set.seed(5000)
aboot(Aeut.pop, sample = 1000) # compare to Grunwald et al. 2006

# You can also use the strata argument to convert to genpop inside the function.
set.seed(5000)
aboot(Aeut.gc, strata = ~Pop/Subpop, sample = 1000)

# And genlight objects 
# From glSim:
## 1,000 non structured SNPs, 100 structured SNPs
x <- glSim(100, 1e3, n.snp.struc=100, ploid=2)
aboot(x, distance = bitwise.dist)

# Utilizing other tree methods

library("ape")

aboot(Aeut.pop, tree = fastme.bal, sample = 1000)

# Utilizing options in other tree methods

myFastME <- function(x) fastme.bal(x, nni = TRUE, spr = FALSE, tbr = TRUE)
aboot(Aeut.pop, tree = myFastME, sample = 1000)


## End(Not run)

Description

The Aeut dataset consists of 187 isolates of the Oomycete root rot pathogen, *Aphanomyces euteiches* collected from two different fields in NW Oregon and W Washington, USA.

Usage

data(Aeut)

Format

a [genind()] object with two populations containing a data frame in the 'other' slot called 'population_hierarchy'. This data frame gives indices of the populations and subpopulations for the data set.

References

Grunwald, NJ and Hoheisel, G.A. 2006. Hierarchical Analysis of Diversity, Selfing, and Genetic Differentiation in Populations of the Oomycete *Aphanomyces euteiches*. Phytopathology 96:1134-1141 doi: doi:10.1094/PHYTO-96-1134

Description

Wrapper for snpclone initializer.

Usage

as.snpclone(x, ..., parallel = FALSE, n.cores = NULL, mlg, mlgclass = TRUE)

Arguments

Author(s)

Zhian N. Kamvar

Examples

(x <- as.snpclone(glSim(100, 1e3, ploid=2)))
## Not run: 
# Without parallel processing
system.time(x <- as.snpclone(glSim(1000, 1e5, ploid=2)))

# With parallel processing... doesn't really save you much time.
system.time(x <- as.snpclone(glSim(1000, 1e5, ploid=2, parallel = TRUE), 
                             parallel = TRUE))

## End(Not run)

Description

This function calculates both dissimilarity and Euclidean distances for genlight or snpclone objects.

Usage

bitwise.dist(
  x,
  percent = TRUE,
  mat = FALSE,
  missing_match = TRUE,
  scale_missing = FALSE,
  euclidean = FALSE,
  differences_only = FALSE,
  threads = 0L
)

Arguments

Details

The default distance calculated here is quite simple and goes by many names depending on its application. The most familiar name might be the Hamming distance, or the number of differences between two strings.

As of poppr version 2.8.0, this function now also calculates Euclidean distance and is considerably faster and more memory-efficient than the standard dist() function.

Value

A dist object containing pairwise distances between samples.

Note

This function is optimized for genlight and snpclone objects. This does not mean that it is a catch-all optimization for SNP data. Three assumptions must be met for this function to work:

1. SNPs are bi-allelic

2. Samples are haploid or diploid

3. All samples have the same ploidy

If the user supplies a genind or genclone object, prevosti.dist() will be used for calculation.

Author(s)

Zhian N. Kamvar, Jonah C. Brooks

See Also

diss.dist(), snpclone, genlight, win.ia(), samp.ia()

Examples

set.seed(999)
x <- glSim(n.ind = 10, n.snp.nonstruc = 5e2, n.snp.struc = 5e2, ploidy = 2)
x
# Assess fraction of different alleles
system.time(xd <- bitwise.dist(x, threads = 1L))
xd

# Calculate Euclidean distance
system.time(xdt <- bitwise.dist(x, euclidean = TRUE, scale_missing = TRUE, threads = 1L))
xdt

## Not run: 

# This function is more efficient in both memory and speed than [dist()] for
# calculating Euclidean distance on genlight objects. For example, we can
# observe a clear speed increase when we attempt a calculation on 100k SNPs
# with 10% missing data:

set.seed(999)
mat <- matrix(sample(c(0:2, NA), 
                     100000 * 50, 
                     replace = TRUE, 
                     prob = c(0.3, 0.3, 0.3, 0.1)),
              nrow = 50)
glite <- new("genlight", mat, ploidy = 2)

# Default Euclidean distance 
system.time(dist(glite))

# Bitwise dist
system.time(bitwise.dist(glite, euclidean = TRUE, scale_missing = TRUE))


## End(Not run)

Description

This function will perform the index of association on a bootstrapped data set multiple times to create a distribution, showing the variation of the index due to repeat observations.

Usage

boot.ia(gid, how = "partial", reps = 999, quiet = FALSE, ...)

Arguments

Value

a data frame with the index of association and standardized index of association in columns. Number of rows represents the number of reps.

Note

This function is experimental. Please do not use this unless you know what you are doing.

See Also

ia(), pair.ia(), psex()

Examples

data(Pinf)
boot.ia(Pinf, reps = 99)

Description

as.genclone will create a genclone object from a genind object OR anything that can be passed to the genind initializer.

Usage

bootgen2genind(bg)

as.genclone(x, ..., mlg, mlgclass = TRUE)

genclone2genind(x)

as.genambig(x)

Arguments

Details

genclone2genind will remove the mlg slot from the genclone object, creating a genind object.

as.genambig will convert a genind or genclone object to a polysat genambig class.

Author(s)

Zhian N. Kamvar

See Also

splitStrata, genclone, read.genalex aboot

Examples

data(Aeut)
Aeut

# Conversion to genclone --------------------------------------------------
Aeut.gc <- as.genclone(Aeut)
Aeut.gc

# Conversion to genind ----------------------------------------------------
Aeut.gi <- genclone2genind(Aeut.gc)
Aeut.gi

# Conversion to polysat's "genambig" class --------------------------------
if (require("polysat")) {
  data(Pinf)
  Pinf.gb <- as.genambig(Pinf)
  summary(Pinf.gb)
}

data(nancycats)

# Conversion to bootgen for random sampling of loci -----------------------
nan.bg  <- new("bootgen", nancycats[pop = 9])
nan.bg

# Conversion back to genind -----------------------------------------------
nan.gid <- bootgen2genind(nan.bg)
nan.gid

Description

Create a tree using Bruvo's Distance with non-parametric bootstrapping.

Usage

bruvo.boot(
  pop,
  replen = 1,
  add = TRUE,
  loss = TRUE,
  sample = 100,
  tree = "upgma",
  showtree = TRUE,
  cutoff = NULL,
  quiet = FALSE,
  root = NULL,
  ...
)

Arguments

Details

This function will calculate a tree based off of Bruvo's distance and then utilize boot.phylo to randomly sample loci with replacement, recalculate the tree, and tally up the bootstrap support (measured in percent success). While this function can take any tree function, it has native support for two algorithms: nj and upgma. If you want to use any other functions, you must load the package before you use them (see examples).

Value

a tree of class phylo with nodelables

Note

Please refer to the documentation for bruvo.dist for details on the algorithm. If the user does not provide a vector of appropriate length for replen , it will be estimated by taking the minimum difference among represented alleles at each locus. IT IS NOT RECOMMENDED TO RELY ON THIS ESTIMATION.

Author(s)

Zhian N. Kamvar, Javier F. Tabima

References

Ruzica Bruvo, Nicolaas K. Michiels, Thomas G. D'Souza, and Hinrich Schulenburg. A simple method for the calculation of microsatellite genotype distances irrespective of ploidy level. Molecular Ecology, 13(7):2101-2106, 2004.

See Also

bruvo.dist, nancycats, upgma, nj, boot.phylo, nodelabels, tab, missingno.

Examples

# Please note that the data presented is assuming that the nancycat dataset 
# contains all dinucleotide repeats, it most likely is not an accurate
# representation of the data.

# Load the nancycats dataset and construct the repeat vector.
data(nancycats)
ssr <- rep(2, 9)

# Analyze the 1st population in nancycats

bruvo.boot(popsub(nancycats, 1), replen = ssr)

## Not run: 

# Always load the library before you specify the function.
library("ape")

# Estimate the tree based off of the BIONJ algorithm.

bruvo.boot(popsub(nancycats, 9), replen = ssr, tree = bionj)

# Utilizing  balanced FastME
bruvo.boot(popsub(nancycats, 9), replen = ssr, tree = fastme.bal)

# To change parameters for the tree, wrap it in a function.
# For example, let's build the tree without utilizing subtree-prune-regraft

myFastME <- function(x) fastme.bal(x, nni = TRUE, spr = FALSE, tbr = TRUE)
bruvo.boot(popsub(nancycats, 9), replen = ssr, tree = myFastME)


## End(Not run)

Description

Calculate the average Bruvo's distance over all loci in a population.

Usage

bruvo.dist(pop, replen = 1, add = TRUE, loss = TRUE, by_locus = FALSE)

bruvo.between(
  query,
  ref,
  replen = 1,
  add = TRUE,
  loss = TRUE,
  by_locus = FALSE
)

Arguments

Details

Bruvo's distance between two alleles is calculated as

$d = 1 - 2^{-\mid x \mid}$

, where x is the number of repeat units between the two alleles (see the Algorithms and Equations vignette for more details). These distances are calculated over all combinations of alleles at a locus and then the minimum average distance between allele combinations is taken as the distance for that locus. All loci are then averaged over to obtain the distance between two samples. Missing data is ignored (in the same fashion as mean(c(1:9, NA), na.rm = TRUE)) if all alleles are missing. See the next section for other cases.

Polyploids

Ploidy is irrelevant with respect to calculation of Bruvo's distance. However, since it makes a comparison between all alleles at a locus, it only makes sense that the two loci need to have the same ploidy level. Unfortunately for polyploids, it's often difficult to fully separate distinct alleles at each locus, so you end up with genotypes that appear to have a lower ploidy level than the organism.

To help deal with these situations, Bruvo has suggested three methods for dealing with these differences in ploidy levels:

• Infinite Model - The simplest way to deal with it is to count all missing alleles as infinitely large so that the distance between it and anything else is 1. Aside from this being computationally simple, it will tend to inflate distances between individuals.

• Genome Addition Model - If it is suspected that the organism has gone through a recent genome expansion, the missing alleles will be replace with all possible combinations of the observed alleles in the shorter genotype. For example, if there is a genotype of [69, 70, 0, 0] where 0 is a missing allele, the possible combinations are: [69, 70, 69, 69], [69, 70, 69, 70], [69, 70, 70, 69], and [69, 70, 70, 70]. The resulting distances are then averaged over the number of comparisons.

• Genome Loss Model - This is similar to the genome addition model, except that it assumes that there was a recent genome reduction event and uses the observed values in the full genotype to fill the missing values in the short genotype. As with the Genome Addition Model, the resulting distances are averaged over the number of comparisons.

• Combination Model - Combine and average the genome addition and loss models.

As mentioned above, the infinite model is biased, but it is not nearly as computationally intensive as either of the other models. The reason for this is that both of the addition and loss models requires replacement of alleles and recalculation of Bruvo's distance. The number of replacements required is equal to n^k where where n is the number of potential replacements and k is the number of alleles to be replaced. To reduce the number of calculations and assumptions otherwise, Bruvo's distance will be calculated using the largest observed ploidy in pairwise comparisons. This means that when comparing [69,70,71,0] and [59,60,0,0], they will be treated as triploids.

Value

an object of class dist or a list of these objects if by_locus = TRUE

Functions

• bruvo.between(): Bruvo's distance between a query and a reference Only diferences between query individuals and reference individuals will be reported All other values are NaN

Note

Do not use missingno with this function.

Missing alleles and Bruvo's distance in poppr versions < 2.5

In earlier versions of poppr, the authors had assumed that, because the calculation of Bruvo's distance does not rely on orderd sets of alleles, the imputation methods in the genome addition and genome loss models would also assume unordered alleles for creating the hypothetical genotypes. This means that the results from this imputation did not consider all possible combinations of alleles, resulting in either an over- or under- estimation of Bruvo's distance between two samples with two or more missing alleles. This version of poppr considers all possible combinations when calculating Bruvo's distance for incomplete genotype with a negligable gain in computation time.

If you want to see the effect of this change on your data, you can use the global poppr option old.bruvo.model. Currently, this option is FALSE and you can set it by using options(old.bruvo.model = TRUE), but make sure to reset it to FALSE afterwards.

Repeat Lengths (replen)

The replen argument is crucial for proper analysis of Bruvo's distance since the calculation relies on the knowledge of the number of steps between alleles. To calculate Bruvo's distance, your raw allele calls are first divided by the repeat lengths and then rounded. This can create a problem with repeat lengths of even size due to the IEC 60559 standard that says rounding at 0.5 is to the nearest even number, meaning that it is possible for two alleles that are one step apart may appear to be exactly the same. This can be fixed by subtracting a tiny number from the repeat length with the function fix_replen. Please consider using this before running Bruvo's distance.

Model Choice

The add and loss arguments modify the model choice accordingly:

• Infinite Model: add = FALSE, loss = FALSE

• Genome Addition Model: add = TRUE, loss = FALSE

• Genome Loss Model: add = FALSE, loss = TRUE

• Combination Model (DEFAULT): add = TRUE, loss = TRUE

Details of each model choice are described in the Details section, above. Additionally, genotypes containing all missing values at a locus will return a value of NA and not contribute to the average across loci.

Repeat Lengths

If the user does not provide a vector of appropriate length for replen , it will be estimated by taking the minimum difference among represented alleles at each locus. IT IS NOT RECOMMENDED TO RELY ON THIS ESTIMATION.

Author(s)

Zhian N. Kamvar

David Folarin

References

Ruzica Bruvo, Nicolaas K. Michiels, Thomas G. D'Souza, and Hinrich Schulenburg. A simple method for the calculation of microsatellite genotype distances irrespective of ploidy level. Molecular Ecology, 13(7):2101-2106, 2004.

See Also

fix_replen, test_replen, bruvo.boot, bruvo.msn

Examples

# Please note that the data presented is assuming that the nancycat dataset 
# contains all dinucleotide repeats, it most likely is not an accurate
# representation of the data.

# Load the nancycats dataset and construct the repeat vector.
data(nancycats)
names(alleles(nancycats)) <- locNames(nancycats) # small bug in this data set
# Assume the alleles are all dinucleotide repeats.
ssr <- rep(2, nLoc(nancycats))
test_replen(nancycats, ssr)         # Are the repeat lengths consistent?
(ssr <- fix_replen(nancycats, ssr)) # Nope. We need to fix them.

# Analyze the first population in nancycats
bruvo.dist(popsub(nancycats, 1), replen = ssr)

## Not run: 

# get the per locus estimates:
bruvo.dist(popsub(nancycats, 1), replen = ssr, by_locus = TRUE)

# View each population as a heatmap.
sapply(popNames(nancycats), function(x) 
heatmap(as.matrix(bruvo.dist(popsub(nancycats, x), replen = ssr)), symm=TRUE))

## End(Not run)

Description

Create minimum spanning network of selected populations using Bruvo's distance.

Usage

bruvo.msn(
  gid,
  replen = 1,
  add = TRUE,
  loss = TRUE,
  mlg.compute = "original",
  palette = topo.colors,
  sublist = "All",
  exclude = NULL,
  blacklist = NULL,
  vertex.label = "MLG",
  gscale = TRUE,
  glim = c(0, 0.8),
  gadj = 3,
  gweight = 1,
  wscale = TRUE,
  showplot = TRUE,
  include.ties = FALSE,
  threshold = NULL,
  clustering.algorithm = NULL,
  ...
)

Arguments

Details

The minimum spanning network generated by this function is generated via igraph's minimum.spanning.tree. The resultant graph produced can be plotted using igraph functions, or the entire object can be plotted using the function plot_poppr_msn, which will give the user a scale bar and the option to layout your data.

node sizes

The area of the nodes are representative of the number of samples. Because igraph scales nodes by radius, the node sizes in the graph are represented as the square root of the number of samples.

mlg.compute

Each node on the graph represents a different multilocus genotype. The edges on the graph represent genetic distances that connect the multilocus genotypes. In genclone objects, it is possible to set the multilocus genotypes to a custom definition. This creates a problem for clone correction, however, as it is very possible to define custom lineages that are not monophyletic. When clone correction is performed on these definitions, information is lost from the graph. To circumvent this, The clone correction will be done via the computed multilocus genotypes, either "original" or "contracted". This is specified in the mlg.compute argument, above.

contracted multilocus genotypes

If your incoming data set is of the class genclone, and it contains contracted multilocus genotypes, this function will retain that information for creating the minimum spanning network. You can use the arguments threshold and clustering.algorithm to change the threshold or clustering algorithm used in the network. For example, if you have a data set that has a threshold of 0.1 and you wish to have a minimum spanning network without a threshold, you can simply add threshold = 0.0, and no clustering will happen.

The threshold and clustering.algorithm arguments can also be used to filter un-contracted data sets.

Value

Note

• Please see the documentation for bruvo.dist for details on the algorithm.

• The edges of these graphs may cross each other if the graph becomes too large.

• The nodes in the graph represent multilocus genotypes. The colors of the nodes are representative of population membership. It is not uncommon to see different populations containing the same multilocus genotype.

Author(s)

Zhian N. Kamvar, Javier F. Tabima

References

Ruzica Bruvo, Nicolaas K. Michiels, Thomas G. D'Souza, and Hinrich Schulenburg. A simple method for the calculation of microsatellite genotype distances irrespective of ploidy level. Molecular Ecology, 13(7):2101-2106, 2004.

See Also

bruvo.dist, nancycats, plot_poppr_msn, mst bruvo.boot, greycurve poppr.msn

Examples

# Load the data set.
data(nancycats)

# View populations 8 and 9 with default colors. 
bruvo.msn(nancycats, replen = rep(2, 9), sublist=8:9, vertex.label="inds", 
          vertex.label.cex=0.7, vertex.label.dist=0.4)
## Not run: 
# View heat colors.
bruvo.msn(nancycats, replen=rep(2, 9), sublist=8:9, vertex.label="inds", 
palette=heat.colors, vertex.label.cex=0.7, vertex.label.dist=0.4)

# View custom colors. Here, we use black and orange.
bruvo.msn(nancycats, replen=rep(2, 9), sublist=8:9, vertex.label="inds", 
palette = colorRampPalette(c("orange", "black")), vertex.label.cex=0.7, 
vertex.label.dist=0.4)

# View with darker shades of grey (setting the upper limit to 1/2 black 1/2 white).
bruvo.msn(nancycats, replen=rep(2, 9), sublist=8:9, vertex.label="inds", 
palette = colorRampPalette(c("orange", "black")), vertex.label.cex=0.7, 
vertex.label.dist=0.4, glim=c(0, 0.5))

# View with no grey scaling.
bruvo.msn(nancycats, replen=rep(2, 9), sublist=8:9, vertex.label="inds", 
palette = colorRampPalette(c("orange", "black")), vertex.label.cex=0.7, 
vertex.label.dist=0.4, gscale=FALSE)

# View with no line widths.
bruvo.msn(nancycats, replen=rep(2, 9), sublist=8:9, vertex.label="inds", 
palette = colorRampPalette(c("orange", "black")), vertex.label.cex=0.7, 
vertex.label.dist=0.4, wscale=FALSE)

# View with no scaling at all.
bruvo.msn(nancycats, replen=rep(2, 9), sublist=8:9, vertex.label="inds", 
palette = colorRampPalette(c("orange", "black")), vertex.label.cex=0.7, 
vertex.label.dist=0.4, gscale=FALSE)

# View the whole population, but without labels.
bruvo.msn(nancycats, replen=rep(2, 9), vertex.label=NA)

## End(Not run)

Description

This function removes any duplicated multilocus genotypes from any specified population strata.

Usage

clonecorrect(pop, strata = 1, combine = FALSE, keep = 1)

Arguments

Details

This function will clone correct based on the stratification provided. To clone correct indiscriminately of population structure, set strata = NA.

Value

a clone corrected genclone, snpclone, or genind object.

Author(s)

Zhian N. Kamvar

Examples

# LOAD A. euteiches data set
data(Aeut)

# Redefine it as a genclone object
Aeut <- as.genclone(Aeut)
strata(Aeut) <- other(Aeut)$population_hierarchy[-1]

# Check the number of multilocus genotypes
mlg(Aeut)
popNames(Aeut)

# Clone correct at the population level.
Aeut.pop <- clonecorrect(Aeut, strata =  ~Pop)
mlg(Aeut.pop)
popNames(Aeut.pop)

## Not run: 
# Clone correct at the subpopulation level with respect to population and
# combine.
Aeut.subpop <- clonecorrect(Aeut, strata = ~Pop/Subpop, combine=TRUE)
mlg(Aeut.subpop)
popNames(Aeut.subpop)

# Do the same, but set to the population level.
Aeut.subpop2 <- clonecorrect(Aeut, strata = ~Pop/Subpop, keep=1)
mlg(Aeut.subpop2)
popNames(Aeut.subpop2)

# LOAD H3N2 dataset
data(H3N2)

strata(H3N2) <- other(H3N2)$x

# Extract only the individuals located in China
country <- clonecorrect(H3N2, strata =  ~country)

# How many isolates did we have from China before clone correction?
sum(strata(H3N2, ~country) == "China") # 155

# How many unique isolates from China after clone correction?
sum(strata(country, ~country) == "China") # 79

# Something a little more complicated. (This could take a few minutes on
# slower computers)

# setting the hierarchy to be Country > Year > Month  
c.y.m <- clonecorrect(H3N2, strata =  ~year/month/country)

# How many isolates in the original data set?
nInd(H3N2) # 1903

# How many after we clone corrected for country, year, and month?
nInd(c.y.m) # 1190

## End(Not run)

Description

Given a series of thresholds for a data set that collapse it into one giant cluster, this will search the top fraction of threshold differences to find the largest difference. The average between the thresholds spanning that difference is the cutoff threshold defining the clonal lineage threshold.

Usage

cutoff_predictor(thresholds, fraction = 0.5)

Arguments

Value

a numeric value representing the threshold at which multilocus lineages should be defined.

Note

This function originally appeared in doi:10.5281/zenodo.17424. This is a bit of a blunt instrument.

Author(s)

Zhian N. Kamvar

References

ZN Kamvar, JC Brooks, and NJ Grünwald. 2015. Supplementary Material for Frontiers Plant Genetics and Genomics 'Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality'. DOI: doi:10.5281/zenodo.17424

Kamvar ZN, Brooks JC and Grünwald NJ (2015) Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Front. Genet. 6:208. doi: doi:10.3389/fgene.2015.00208

See Also

filter_stats mlg.filter

Examples

data(Pinf)
pinfreps <- fix_replen(Pinf, c(2, 2, 6, 2, 2, 2, 2, 2, 3, 3, 2))
pthresh  <- filter_stats(Pinf, distance = bruvo.dist, replen = pinfreps, 
                         plot = TRUE, stats = "THRESHOLD", threads = 1L)

# prediction for farthest neighbor
cutoff_predictor(pthresh$farthest)

# prediction for all algorithms
sapply(pthresh, cutoff_predictor)

Description

diss.dist uses the same discrete dissimilarity matrix utilized by the index of association (see ia for details). By default, it returns a distance reflecting the number of allelic differences between two individuals. When percent = TRUE, it returns a ratio of the number of observed differences by the number of possible differences. Eg. two individuals who share half of the same alleles will have a distance of 0.5. This function can analyze distances for any marker system.

Usage

diss.dist(x, percent = FALSE, mat = FALSE)

Arguments

Details

The distance calculated here is quite simple and goes by many names, depending on its application. The most familiar name might be the Hamming distance, or the number of differences between two strings.

Value

Pairwise distances between individuals present in the genind object.

Note

When percent = TRUE, this is exactly the same as provesti.dist, except that it performs better for large numbers of individuals (n > 125) at the cost of available memory.

Author(s)

Zhian N. Kamvar

See Also

prevosti.dist, bitwise.dist (for SNP data)

Examples

# A simple example. Let's analyze the mean distance among populations of A.
# euteiches.

data(Aeut)
mean(diss.dist(popsub(Aeut, 1)))
## Not run: 
mean(diss.dist(popsub(Aeut, 2)))
mean(diss.dist(Aeut))

## End(Not run)

Description

This function is intended to perform bootstrap statistics on a matrix of multilocus genotype counts in different populations. Results from this function should be interpreted carefully as the default statistics are known to have a downward bias. See the details for more information.

Usage

diversity_boot(
  tab,
  n,
  n.boot = 1L,
  n.rare = NULL,
  H = TRUE,
  G = TRUE,
  lambda = TRUE,
  E5 = TRUE,
  ...
)

Arguments

Details

Bootstrapping is performed in three ways:

• if n.rare is a number greater than zero, then bootstrapping is performed by randomly sampling without replacement n.rare samples from the data.

  \item if `n.boot` is greater than 1, bootstrapping is performed by
  sampling n.boot samples from a multinomial distribution weighted by the
  proportion of each MLG in the data.
  
  \item if `n.boot` is less than 2, bootstrapping is performed by 
  sampling N samples from a multinomial distribution weighted by the
  proportion of each MLG in the data.
  

Downward Bias

When sampling with replacement, the diversity statistics here present a downward bias partially due to the small number of samples in the data. The result is that the mean of the bootstrapped samples will often be much lower than the observed value. Alternatively, you can increase the sample size of the bootstrap by increasing the size of n.boot. Both of these methods should be taken with caution in interpretation. There are several R packages freely available that will calculate and perform bootstrap estimates of Shannon and Simpson diversity metrics (eg. entropart, entropy, simboot, and EntropyEstimation. These packages also offer unbiased estimators of Shannon and Simpson diversity. Please take care when attempting to interpret the results of this function.

Value

a list of objects of class "boot".

Author(s)

Zhian N. Kamvar

See Also

diversity_stats() for basic statistic calculation, diversity_ci() for confidence intervals and plotting, and poppr(). For bootstrap sampling: stats::rmultinom() boot::boot()

Examples

library(poppr)
data(Pinf)
tab <- mlg.table(Pinf, plot = FALSE)
diversity_boot(tab, 10L)
## Not run: 
# This can be done in a parallel fashion (OSX uses "multicore", Windows uses "snow")
system.time(diversity_boot(tab, 10000L, parallel = "multicore", ncpus = 4L))
system.time(diversity_boot(tab, 10000L))

## End(Not run)

Description

This function is for calculating bootstrap statistics and their confidence intervals. It is important to note that the calculation of confidence intervals is not perfect (See Details). Please be cautious when interpreting the results.

Usage

diversity_ci(
  tab,
  n = 1000,
  n.boot = 1L,
  ci = 95,
  total = TRUE,
  rarefy = FALSE,
  n.rare = 10,
  plot = TRUE,
  raw = TRUE,
  center = TRUE,
  ...
)

Arguments

Details

Bootstrapping

For details on the bootstrapping procedures, see diversity_boot(). Default bootstrapping is performed by sampling N samples from a multinomial distribution weighted by the relative multilocus genotype abundance per population where N is equal to the number of samples in the data set. If n.boot > 2, then n.boot samples are taken at each bootstrap replicate. When rarefy = TRUE, then samples are taken at the smallest population size without replacement. This will provide confidence intervals for all but the smallest population.

Confidence intervals

Confidence intervals are derived from the function boot::norm.ci(). This function will attempt to correct for bias between the observed value and the bootstrapped estimate. When center = TRUE (default), the confidence interval is calculated from the bootstrapped distribution and centered around the bias-corrected estimate as prescribed in Marcon (2012). This method can lead to undesirable properties, such as the confidence interval lying outside of the maximum possible value. For rarefaction, the confidence interval is simply determined by calculating the percentiles from the bootstrapped distribution. If you want to calculate your own confidence intervals, you can use the results of the permutations stored in the ⁠$boot⁠ element of the output.

Rarefaction

Rarefaction in the sense of this function is simply sampling a subset of the data at size n.rare. The estimates derived from this method have straightforward interpretations and allow you to compare diversity across populations since you are controlling for sample size.

Plotting

Results are plotted as boxplots with point estimates. If there is no rarefaction applied, confidence intervals are displayed around the point estimates. The boxplots represent the actual values from the bootstrapping and will often appear below the estimates and confidence intervals.

Value

raw = TRUE

• obs a matrix with observed statistics in columns, populations in rows

• est a matrix with estimated statistics in columns, populations in rows

• CI an array of 3 dimensions giving the lower and upper bound, the index measured, and the population.

• boot a list containing the output of boot::boot() for each population.

raw = FALSE

a data frame with the statistic observations, estimates, and confidence intervals in columns, and populations in rows. Note that the confidence intervals are converted to characters and rounded to three decimal places.

Note

Confidence interval calculation

Almost all of the statistics supplied here have a maximum when all genotypes are equally represented. This means that bootstrapping the samples will always be downwardly biased. In many cases, the confidence intervals from the bootstrapped distribution will fall outside of the observed statistic. The reported confidence intervals here are reported by assuming the variance of the bootstrapped distribution is the same as the variance around the observed statistic. As different statistics have different properties, there will not always be one clear method for calculating confidence intervals. A suggestion for correction in Shannon's index is to center the CI around the observed statistic (Marcon, 2012), but there are theoretical limitations to this. For details, see https://stats.stackexchange.com/q/156235/49413.

User-defined functions

While it is possible to use custom functions with this, there are three important things to remember when using these functions:

1. The function must return a single value. 
2. The function must allow for both matrix and vector inputs 
3. The function name cannot match or partially match any arguments 
from [boot::boot()]

Anonymous functions are okay
(e.g. function(x) vegan::rarefy(t(as.matrix(x)), 10)).

Author(s)

Zhian N. Kamvar

References

Marcon, E., Herault, B., Baraloto, C. and Lang, G. (2012). The Decomposition of Shannon’s Entropy and a Confidence Interval for Beta Diversity. Oikos 121(4): 516-522.

See Also

diversity_boot() diversity_stats() poppr() boot::boot() boot::norm.ci() boot::boot.ci()

Examples

library(poppr)
data(Pinf)
diversity_ci(Pinf, n = 100L)
## Not run: 
# With pretty results
diversity_ci(Pinf, n = 100L, raw = FALSE)

# This can be done in a parallel fasion (OSX uses "multicore", Windows uses "snow")
system.time(diversity_ci(Pinf, 10000L, parallel = "multicore", ncpus = 4L))
system.time(diversity_ci(Pinf, 10000L))

# We often get many requests for a clonal fraction statistic. As this is 
# simply the number of observed MLGs over the number of samples, we 
# recommended that people calculate it themselves. With this function, you
# can add it in:

CF <- function(x){
 x <- drop(as.matrix(x))
 if (length(dim(x)) > 1){
   res <- rowSums(x > 0)/rowSums(x)
 } else {
   res <- sum(x > 0)/sum(x)
 }
 return(res)
}
# Show pretty results

diversity_ci(Pinf, 1000L, CF = CF, center = TRUE, raw = FALSE)
diversity_ci(Pinf, 1000L, CF = CF, rarefy = TRUE, raw = FALSE)

## End(Not run)

Description

Calculate diversity statistics on a matrix containing counts of multilocus genotypes per population.

Usage

diversity_stats(z, H = TRUE, G = TRUE, lambda = TRUE, E5 = TRUE, ...)

Arguments

Details

This function will calculate any diversity statistic for counts of multilocus genotypes per population. This does not count allelic diversity. The calculations of H, G, and lambda are all performed by vegan::diversity(). E5 is calculated as

$E_{5} = \frac{(1/\lambda) - 1}{e^{H} - 1}$

.

Value

a numeric matrix giving statistics (columns) for each population (rows).

Author(s)

Zhian N. Kamvar

See Also

diversity_boot() diversity_ci() poppr()

Examples

library(poppr)
data(Pinf)
tab <- mlg.table(Pinf, plot = FALSE)
diversity_stats(tab)
## Not run: 
# Example using the poweRlaw package to calculate the negative slope of the
# Pareto distribution.

library("poweRlaw")
power_law_beta <- function(x){
  xpow <- displ(x[x > 0])                 # Generate the distribution
  xpow$setPars(estimate_pars(xpow))       # Estimate the parameters
  xdat <- plot(xpow, draw = FALSE)        # Extract the data
  xlm <- lm(log(y) ~ log(x), data = xdat) # Run log-log linear model for slope
  return(-coef(xlm)[2])
}

Beta <- function(x){
  x <- drop(as.matrix(x))
  if (length(dim(x)) > 1){
    res <- apply(x, 1, power_law_beta)
  } else {
    res <- power_law_beta(x)
  }
  return(res)
}

diversity_stats(tab, B = Beta)

## End(Not run)

Description

This function is a wrapper to mlg.filter. It will calculate all of the stats for mlg.filter utilizing all of the algorithms.

Usage

filter_stats(
  x,
  distance = bitwise.dist,
  threshold = 1e+06 + .Machine$double.eps^0.5,
  stats = "All",
  missing = "ignore",
  plot = FALSE,
  cols = NULL,
  nclone = NULL,
  hist = "Scott",
  threads = 1L,
  ...
)

Arguments

Value

a list of results from mlg.filter from the three algorithms. (returns invisibly if plot = TRUE)

Note

This function originally appeared in doi:10.5281/zenodo.17424

Author(s)

Zhian N. Kamvar, Jonah C. Brooks

References

ZN Kamvar, JC Brooks, and NJ Grünwald. 2015. Supplementary Material for Frontiers Plant Genetics and Genomics 'Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality'. DOI: doi:10.5281/zenodo.17424

Kamvar ZN, Brooks JC and Grünwald NJ (2015) Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Front. Genet. 6:208. doi: doi:10.3389/fgene.2015.00208

See Also

mlg.filter cutoff_predictor bitwise.dist diss.dist

Examples

# Basic usage example: Bruvo's Distance --------------------------------
data(Pinf)
pinfreps <- fix_replen(Pinf, c(2, 2, 6, 2, 2, 2, 2, 2, 3, 3, 2))
bres <- filter_stats(Pinf, distance = bruvo.dist, replen = pinfreps, plot = TRUE, threads = 1L)
print(bres) # shows all of the statistics

# Use these results with cutoff_filter()
print(thresh <- cutoff_predictor(bres$farthest$THRESHOLDS))
mlg.filter(Pinf, distance = bruvo.dist, replen = pinfreps) <- thresh
Pinf 

# Different distances will give different results -----------------------
nres <- filter_stats(Pinf, distance = nei.dist, plot = TRUE, threads = 1L, missing = "mean")
print(thresh <- cutoff_predictor(nres$farthest$THRESHOLDS))
mlg.filter(Pinf, distance = nei.dist, missing = "mean") <- thresh
Pinf

Description

Attempts to fix inconsistent repeat lengths found by test_replen

Usage

fix_replen(gid, replen, e = 1e-05, fix_some = TRUE)

Arguments

Details

This function is modified from the version used in doi:10.5281/zenodo.13007.
Before being fed into the algorithm to calculate Bruvo's distance, the amplicon length is divided by the repeat unit length. Because of the amplified primer sequence attached to sequence repeat, this division does not always result in an integer and so the resulting numbers are rounded. The rounding also protects against slight mis-calls of alleles. Because we know that

$\frac{(A - e) - (B - e)}{r}$

is equivalent to

$\frac{A - B}{r}$

, we know that the primer sequence will not alter the relationships between the alleles. Unfortunately for nucleotide repeats that have powers of 2, rounding in R is based off of the IEC 60559 standard (see round), that means that any number ending in 5 is rounded to the nearest even digit. This function will attempt to alleviate this problem by adding a very small amount to the repeat length so that division will not result in a 0.5. If this fails, the same amount will be subtracted. If neither of these work, a warning will be issued and it is up to the user to determine if the fault is in the allele calls or the repeat lengths.

Value

a numeric vector of corrected repeat motif lengths.

Author(s)

Zhian N. Kamvar

References

Zhian N. Kamvar, Meg M. Larsen, Alan M. Kanaskie, Everett M. Hansen, & Niklaus J. Grünwald. Sudden_Oak_Death_in_Oregon_Forests: Spatial and temporal population dynamics of the sudden oak death epidemic in Oregon Forests. ZENODO, doi:10.5281/zenodo.13007, 2014.

Kamvar, Z. N., Larsen, M. M., Kanaskie, A. M., Hansen, E. M., & Grünwald, N. J. (2015). Spatial and temporal analysis of populations of the sudden oak death pathogen in Oregon forests. Phytopathology 105:982-989. doi: doi:10.1094/PHYTO-12-14-0350-FI

Ruzica Bruvo, Nicolaas K. Michiels, Thomas G. D'Souza, and Hinrich Schulenburg. A simple method for the calculation of microsatellite genotype distances irrespective of ploidy level. Molecular Ecology, 13(7):2101-2106, 2004.

See Also

test_replen bruvo.dist bruvo.msn bruvo.boot

Examples

data(Pram)
(Pram_replen <- setNames(c(3, 2, 4, 4, 4), locNames(Pram)))
fix_replen(Pram, Pram_replen)
# Let's start with an example of a tetranucleotide repeat motif and imagine
# that there are twenty alleles all 1 step apart:
(x <- 1:20L * 4L)
# These are the true lengths of the different alleles. Now, let's add the
# primer sequence to them. 
(PxP <- x + 21 + 21)
# Now we make sure that x / 4 is equal to 1:20, which we know each have
# 1 difference.
x/4
# Now, we divide the sequence with the primers by 4 and see what happens.
(PxPc <- PxP/4)
(PxPcr <- round(PxPc))
diff(PxPcr) # we expect all 1s

# Let's try that again by subtracting a tiny amount from 4
(PxPc <- PxP/(4 - 1e-5))
(PxPcr <- round(PxPc))
diff(PxPcr)

Description

GENclone is an S4 class that extends the genind object.
SNPclone is an S4 class that extends the genlight object.

They will have all of the same attributes as their parent classes, but they will contain one extra slot to store extra information about multilocus genotypes.

Details

The genclone and snpclone classes will allow for more optimized methods of clone correction.

Previously for genind and genlight objects, multilocus genotypes were not retained after a data set was subset by population. The new mlg slot allows us to assign the multilocus genotypes and retain that information no matter how we subset the data set. This new slot can either contain numeric values for multilocus genotypes OR it can contain a special internal MLG class that allows for custom multilocus genotype definitions and filtering.

Slots

mlg

a vector representing multilocus genotypes for the data set OR an object of class MLG.

Extends

The genclone class extends class "genind", directly.
The snpclone class extends class "genlight", directly.

Note

When calculating multilocus genotypes for genclone objects, a rank function is used, but calculation of multilocus genotypes for snpclone objects is distance-based (via bitwise.dist and mlg.filter). This means that genclone objects are sensitive to missing data, whereas snpclone objects are insensitive.

Author(s)

Zhian N. Kamvar

See Also

as.genclone as.snpclone genind genlight strata setPop MLG mll mlg.filter

Examples

## Not run: 

# genclone objects can be created from genind objects
#
data(partial_clone)
partial_clone
(pc <- as.genclone(partial_clone))

# snpclone objects can be created from genlight objects
#
set.seed(999)
(gl <- glSim(100, 0, n.snp.struc = 1e3, ploidy = 2, parallel = FALSE))
(sc <- as.snpclone(rbind(gl, gl, parallel = FALSE), parallel = FALSE))
# 
# Use mlg.filter to create a distance threshold to define multilocus genotypes.
mlg.filter(sc, threads = 1L) <- 0.25
sc # 82 mlgs


## End(Not run)

Description

genind2genalex will export a genclone or genind object to a \*.csv file formatted for use in genalex.

Usage

genind2genalex(
  gid,
  filename = "",
  overwrite = FALSE,
  quiet = FALSE,
  pop = NULL,
  allstrata = TRUE,
  geo = FALSE,
  geodf = "xy",
  sep = ",",
  sequence = FALSE
)

Arguments

Value

The the file path or connection where the data were written.

Note

If your data set lacks a population structure, it will be coded in the new file as a single population labeled "Pop". Likewise, if you don't have any labels for your individuals, they will be labeled as "ind1" through "indN", with N being the size of your population.

Author(s)

Zhian N. Kamvar

See Also

read.genalex(), clonecorrect(), genclone, genind

Examples

## Not run: 
data(nancycats)
genind2genalex(nancycats, "~/Documents/nancycats.csv", geo=TRUE)

## End(Not run)

Description

Genotype accumulation curves are useful for determining the minimum number of loci necessary to discriminate between individuals in a population. This function will randomly sample loci without replacement and count the number of multilocus genotypes observed.

Usage

genotype_curve(
  gen,
  sample = 100,
  maxloci = 0L,
  quiet = FALSE,
  thresh = 1,
  plot = TRUE,
  drop = TRUE,
  dropna = TRUE
)

Arguments

Details

Internally, this function works by converting the data into a loci object, which represents genotypes as a data frame of factors. Random samples are taken of 1 to n-1 columns of the matrix and the number of unique rows are counted to determine the number of multilocus genotypes in that random sample. This function does not take into account any definitions of MLGs via mlg.filter or mll.custom.

Value

(invisibly by deafuls) a matrix of integers showing the results of each randomization. Columns represent the number of loci sampled and rows represent an independent sample.

Author(s)

Zhian N. Kamvar

Examples

data(nancycats)
nan_geno <- genotype_curve(nancycats)
## Not run: 

# Marker Type Comparison --------------------------------------------------
# With AFLP data, it is often necessary to include more markers for resolution
data(Aeut)
Ageno <- genotype_curve(Aeut)

# Many microsatellite data sets have hypervariable markers
data(microbov)
mgeno <- geotype_curve(microbov)

# Adding a trendline ------------------------------------------------------

# Trendlines: you can add a smoothed trendline with geom_smooth()
library("ggplot2")
p <- last_plot()
p + geom_smooth()

# Producing Figures for Publication ---------------------------------------

# This data set has been pre filtered
data(monpop)
mongeno <- genotype_curve(monpop)

# Here, we add a curve and a title for publication
p <- last_plot()
mytitle <- expression(paste("Genotype Accumulation Curve for ", 
                            italic("M. fructicola")))
p + geom_smooth() + 
  theme_bw() + 
  theme(text = element_text(size = 12, family = "serif")) +
  theme(title = element_text(size = 14)) +
  ggtitle(mytitle)

## End(Not run)

Description

getfile is a convenience function that serves as a wrapper for the functions file.choose(), file.path(), and list.files(). If the user is working in a GUI environment, a window will pop up, allowing the user to choose a specified file regardless of path.

Usage

getfile(multi = FALSE, pattern = NULL, combine = TRUE)

Arguments

Value

Author(s)

Zhian N. Kamvar

Examples

## Not run: 

x <- getfile()
poppr(x$files)

y <- getfile(multi=TRUE, pattern="^.+?dat$") 
#useful for reading in multiple FSTAT formatted files.

yfiles <- poppr.all(y$files)


## End(Not run)

Description

This function has one purpose. It is for deciding the appropriate scaling for a grey palette to be used for edge weights of a minimum spanning network.

Usage

greycurve(
  data = seq(0, 1, length = 1000),
  glim = c(0, 0.8),
  gadj = 3,
  gweight = 1,
  scalebar = FALSE
)

Arguments

Value

A plot displaying a grey gradient from 0.001 to 1 with minimum and maximum values displayed as yellow lines, and an equation for the correction displayed in red.

Author(s)

Zhian N. Kamvar

Examples

# Normal grey curve with an adjustment of 3, an upper limit of 0.8, and
# weighted towards smaller values.
greycurve()
## Not run: 
# 1:1 relationship grey curve.
greycurve(gadj=1, glim=1:0)

# Grey curve weighted towards larger values.
greycurve(gweight=2)

# Same as the first, but the limit is 1.
greycurve(glim=1:0)

# Setting the lower limit to 0.1 and weighting towards larger values.
greycurve(glim=c(0.1,0.8), gweight=2)

## End(Not run)

Description

Calculate the Index of Association and Standardized Index of Association.

Usage

ia(
  gid,
  sample = 0,
  method = 1,
  quiet = FALSE,
  missing = "ignore",
  plot = TRUE,
  hist = TRUE,
  index = "rbarD",
  valuereturn = FALSE
)

pair.ia(
  gid,
  sample = 0L,
  quiet = FALSE,
  plot = TRUE,
  low = "blue",
  high = "red",
  limits = NULL,
  index = "rbarD",
  method = 1L
)

resample.ia(gid, n = NULL, reps = 999, quiet = FALSE, use_psex = FALSE, ...)

jack.ia(gid, n = NULL, reps = 999, quiet = FALSE)

Arguments

Details

• ia() calculates the index of association over all loci in the data set.

• pair.ia() calculates the index of association in a pairwise manner among all loci.

• resample.ia() calculates the index of association on a reduced data set multiple times to create a distribution, showing the variation of values observed at a given sample size (previously jack.ia()).

The index of association was originally developed by A.H.D. Brown analyzing population structure of wild barley (Brown, 1980). It has been widely used as a tool to detect clonal reproduction within populations . Populations whose members are undergoing sexual reproduction, whether it be selfing or out-crossing, will produce gametes via meiosis, and thus have a chance to shuffle alleles in the next generation. Populations whose members are undergoing clonal reproduction, however, generally do so via mitosis. This means that the most likely mechanism for a change in genotype is via mutation. The rate of mutation varies from species to species, but it is rarely sufficiently high to approximate a random shuffling of alleles. The index of association is a calculation based on the ratio of the variance of the raw number of differences between individuals and the sum of those variances over each locus . You can also think of it as the observed variance over the expected variance. If they are the same, then the index is zero after subtracting one (from Maynard-Smith, 1993):

$I_A = \frac{V_O}{V_E}-1$

Since the distance is more or less a binary distance, any sort of marker can be used for this analysis. In the calculation, phase is not considered, and any difference increases the distance between two individuals. Remember that each column represents a different allele and that each entry in the table represents the fraction of the genotype made up by that allele at that locus. Notice also that the sum of the rows all equal one. Poppr uses this to calculate distances by simply taking the sum of the absolute values of the differences between rows.

The calculation for the distance between two individuals at a single locus with a allelic states and a ploidy of k is as follows (except for Presence/Absence data):

$d = \displaystyle \frac{k}{2}\sum_{i=1}^{a} \mid A_{i} - B_{i}\mid$

To find the total number of differences between two individuals over all loci, you just take d over m loci, a value we'll call D:

$D = \displaystyle \sum_{i=1}^{m} d_i$

These values are calculated over all possible combinations of individuals in the data set, ${n \choose 2}$ after which you end up with ${n \choose 2}\cdot{}m$ values of d and ${n \choose 2}$ values of D. Calculating the observed variances is fairly straightforward (modified from Agapow and Burt, 2001):

$V_O = \frac{\displaystyle \sum_{i=1}^{n \choose 2} D_{i}^2 - \frac{(\displaystyle\sum_{i=1}^{n \choose 2} D_{i})^2}{{n \choose 2}}}{{n \choose 2}}$

Calculating the expected variance is the sum of each of the variances of the individual loci. The calculation at a single locus, j is the same as the previous equation, substituting values of D for d:

$var_j = \frac{\displaystyle \sum_{i=1}^{n \choose 2} d_{i}^2 - \frac{(\displaystyle\sum_{i=1}^{n \choose 2} d_i)^2}{{n \choose 2}}}{{n \choose 2}}$

The expected variance is then the sum of all the variances over all m loci:

$V_E = \displaystyle \sum_{j=1}^{m} var_j$

Agapow and Burt showed that $I_A$ increases steadily with the number of loci, so they came up with an approximation that is widely used, $\bar r_d$. For the derivation, see the manual for multilocus.

$\bar r_d = \frac{V_O - V_E} {2\displaystyle \sum_{j=1}^{m}\displaystyle \sum_{k \neq j}^{m}\sqrt{var_j\cdot{}var_k}}$

Value

for pair.ia()

A matrix with two columns and choose(nLoc(gid), 2) rows representing the values for Ia and rbarD per locus pair.

If no sampling has occurred:

A named number vector of length 2 giving the Index of Association, "Ia"; and the Standardized Index of Association, "rbarD"

If there is sampling:

A a named numeric vector of length 4 with the following values:

• Ia - numeric. The index of association.

• p.Ia - A number indicating the p-value resulting from a one-sided permutation test based on the number of samples indicated in the original call.

• rbarD - numeric. The standardized index of association.

• p.rD - A factor indicating the p-value resulting from a one-sided permutation test based on the number of samples indicated in the original call.

If there is sampling and valureturn = TRUE

A list with the following elements:

• index The above vector

• samples A data frame with s by 2 column data frame where s is the number of samples defined. The columns are for the values of Ia and rbarD, respectively.

resample.ia()

a data frame with the index of association and standardized index of association in columns. Number of rows represents the number of reps.

Note

jack.ia() is deprecated as the name was misleading. Please use resample.ia()

Author(s)

Zhian N. Kamvar

References

Paul-Michael Agapow and Austin Burt. Indices of multilocus linkage disequilibrium. Molecular Ecology Notes, 1(1-2):101-102, 2001

A.H.D. Brown, M.W. Feldman, and E. Nevo. Multilocus structure of natural populations of Hordeum spontaneum. Genetics, 96(2):523-536, 1980.

J M Smith, N H Smith, M O'Rourke, and B G Spratt. How clonal are bacteria? Proceedings of the National Academy of Sciences, 90(10):4384-4388, 1993.

See Also

poppr(), missingno(), import2genind(), read.genalex(), clonecorrect(), win.ia(), samp.ia()

Examples

data(nancycats)
ia(nancycats)

# Pairwise over all loci:
data(partial_clone)
res <- pair.ia(partial_clone)
plot(res, low = "black", high = "green", index = "Ia")

# Resampling
data(Pinf)
resample.ia(Pinf, reps = 99)

## Not run: 

# Pairwise IA with p-values (this will take about a minute)
res <- pair.ia(partial_clone, sample = 999)
head(res)

# Plot the results of resampling rbarD. 
library("ggplot2")
Pinf.resamp <- resample.ia(Pinf, reps = 999)
ggplot(Pinf.resamp[2], aes(x = rbarD)) +
  geom_histogram() +
  geom_vline(xintercept = ia(Pinf)[2]) +
  geom_vline(xintercept = ia(clonecorrect(Pinf))[2], linetype = 2) +
  xlab(expression(bar(r)[d]))

# Get the indices back and plot the distributions.
nansamp <- ia(nancycats, sample = 999, valuereturn = TRUE)

plot(nansamp, index = "Ia")
plot(nansamp, index = "rbarD")

# You can also adjust the parameters for how large to display the text
# so that it's easier to export it for publication/presentations.
library("ggplot2")
plot(nansamp, labsize = 5, linesize = 2) +
  theme_bw() +                                      # adding a theme
  theme(text = element_text(size = rel(5))) +       # changing text size
  theme(plot.title = element_text(size = rel(4))) + # changing title size
  ggtitle("Index of Association of nancycats")      # adding a new title

# Get the index for each population.
lapply(seppop(nancycats), ia)
# With sampling
lapply(seppop(nancycats), ia, sample = 999)

# Plot pairwise ia for all populations in a grid with cowplot
# Set up the library and data
library("cowplot")
data(monpop)
splitStrata(monpop) <- ~Tree/Year/Symptom
setPop(monpop)      <- ~Tree

# Need to set up a list in which to store the plots.
plotlist        <- vector(mode = "list", length = nPop(monpop))
names(plotlist) <- popNames(monpop)

# Loop throgh the populations, calculate pairwise ia, plot, and then
# capture the plot in the list
for (i in popNames(monpop)){
  x <- pair.ia(monpop[pop = i], limits = c(-0.15, 1)) # subset, calculate, and plot
  plotlist[[i]] <- ggplot2::last_plot() # save the last plot
}

# Use the plot_grid function to plot.
plot_grid(plotlist = plotlist, labels = paste("Tree", popNames(monpop)))


## End(Not run)

Description

This function will launch an interactive interface that allows you to create, plot, manipulate, and save minimum spanning networks. It runs using the shiny R package.

Usage

imsn()

Details

Creating and plotting MSNs requires three steps:

1. Create a distance matrix from your data

2. Create a minimum spanning network with your data and the matrix

3. Visualize the minimum spanning network

The function plot_poppr_msn is currently the most flexible way of visualizing your minimum spanning network, but with 20 parameters, it can become pretty intimidating trying to find the right display for your MSN.

With this function, all three steps are combined into one interactive interface that will allow you to intuitively modify your minimum spanning network and even save the results to a pdf or png file.

Value

NULL, invisibly

Interface

Buttons

In the left hand panel, there are three buttons to execute the functions. These allow you to run the data set after you manipulate all of the parameters.

• GO! - This button will start the application with the specified parameters

• reData - Use this button when you have changed any parameters under the section Data Parameters. This involves recalculating the distance matrix and msn.

• reGraph - Use this button when you have changed any parameters under the section Graphical Parameters. This involves superficial changes to the display of the minimum spanning network.

Tabs

The right hand panel contains different tabs related to your data set of choice.

• Plot - The minimum spanning network itself

• Data - A display of your data set

• Command - The commands used to create the plot. You can copy and paste this to an R file for reproducibility.

• Save Plot - This provides a tool for you to save the plot to a PDF or PNG image.

• Session Information - displays the result of sessionInfo for reproducibility.

Author(s)

Zhian N. Kamvar

See Also

plot_poppr_msn diss.dist bruvo.dist bruvo.msn poppr.msn nei.dist popsub missingno

Examples

## Not run: 

# Set up some data
library("poppr")
library("magrittr")
data(monpop)
splitStrata(monpop) <- ~Tree/Year/Symptom
summary(monpop)
monpop_ssr <- c(CHMFc4 = 7, CHMFc5 = 2, CHMFc12 = 4, 
                SEA = 4, SED = 4, SEE = 2, SEG = 6, 
                SEI = 3, SEL = 4, SEN = 2, SEP = 4, 
                SEQ = 2, SER = 4)
t26 <- monpop %>% setPop(~Tree) %>% popsub("26") %>% setPop(~Year/Symptom)
t26
if (interactive()) {
  imsn() # select Bruvo's distance and enter "monpop_ssr" into the Repeat Length field.
  
  # It is also possible to run this from github if you are connected to the internet.
  # This allows you to access any bug fixes that may have been updated before a formal
  # release on CRAN

  shiny::runGitHub("grunwaldlab/poppr", subdir = "inst/shiny/msn_explorer")

  # You can also use your own distance matrices, but there's a small catch.
  # in order to do so, you must write a function that will subset the matrix
  # to whatever populations are in your data. Here's an example with the above

  mondist <- bruvo.dist(monpop, replen = monpop_ssr)
  myDist <- function(x, d = mondist){
   dm <- as.matrix(d)          # Convert the dist object to a square matrix
   xi <- indNames(x)           # Grab the sample names that exist
   return(as.dist(dm[xi, xi])) # return only the elements that have the names
                               # in the data set
  }
  # After executing imsn, choose:
  # Distance: custom
  # myDist
  imsn() 
}

## End(Not run)

Description

If two samples share no loci typed in common, they are incomparable and will produce missing data in a distance matrix, which could lead to problems with further analyses. This function finds these samples and returns a matrix of how many other samples these are incomparable with.

Usage

incomp(gid)

Arguments

Value

a square matrix with samples that are incomparable

Examples

data(nancycats)
# These two populations have no samples that are incomparable
incomp(nancycats[pop = c(1, 17)])

# If you reduce the number of loci, we find that there are 
# incomparable samples.
incomp(nancycats[pop = c(1, 17), loc = c(1, 4)])

Description

Create a table summarizing missing data or ploidy information of a genind or genclone object

Usage

info_table(
  gen,
  type = c("missing", "ploidy"),
  percent = TRUE,
  plot = FALSE,
  df = FALSE,
  returnplot = FALSE,
  low = "blue",
  high = "red",
  plotlab = TRUE,
  scaled = TRUE
)

Arguments

Details

Missing data is accounted for on a per-population level.
Ploidy is accounted for on a per-individual level.

For type = 'missing'

This data is potentially useful for identifying areas of systematic missing data. There are a few caveats to be aware of.

• Regarding counts of missing data: Each count represents the number of individuals with missing data at each locus. The last column, "mean" can be thought of as the average number of individuals with missing data per locus.

• Regarding percentage missing data: This percentage is relative to the population and locus, not to the entire data set. The last column, "mean" represents the average percent of the population with missing data per locus.

For type = 'ploidy'

This option is useful for data that has been imported with mixed ploidies. It will summarize the relative levels of ploidy per individual per locus. This is simply based off of observed alleles and does not provide any further estimates.

Value

a matrix, data frame (df = TRUE), or a list (returnplot = TRUE) representing missing data per population (type = 'missing') or ploidy per individual (type = 'ploidy') in a genind or genclone object.

Author(s)

Zhian N. Kamvar

Examples

data(nancycats)
nancy.miss <- info_table(nancycats, plot = TRUE, type = "missing")
data(Pinf)
Pinf.ploid <- info_table(Pinf, plot = TRUE, type = "ploidy")

Description

This function will facilitate in removing phylogenetically uninformative loci from a genclone or genind object. The user has the ability to define what uninformative means by setting a cutoff value for either percentage of differentiating genotypes or minor allele frequency.

Usage

informloci(pop, cutoff = 2/nInd(pop), MAF = 0.01, quiet = FALSE)

Arguments

Details

This function will remove uninformative loci using a traditional MAF cutoff (using isPoly from adegenet) as well as analyzing the number of observed genotypes in a locus. This is important for clonal organisms that can have fixed heterozygous sites not detected by MAF methods.

Value

A genind object with user-defined informative loci.

Note

This will have a few side effects that affect certain analyses. First, the number of multilocus genotypes might be reduced due to the reduced number of markers (if you are only using a genind object). Second, if you plan on using this data for analysis of the index of association, be sure to use the standardized version (rbarD) that corrects for the number of observed loci.

Author(s)

Zhian N. Kamvar

Examples

# We will use a dummy data set to demonstrate how this detects uninformative
# loci using both MAF and a cutoff.

genos <- c("A/A", "A/B", "A/C", "B/B", "B/C", "C/C")

v <- sample(genos, 100, replace = TRUE)
w <- c(rep(genos[2], 99), genos[3])           # found by cutoff
x <- c(rep(genos[1], 98), genos[3], genos[2]) # found by MAF
y <- c(rep(genos[1], 99), genos[2])           # found by both
z <- sample(genos, 100, replace = TRUE)
dat <- df2genind(data.frame(v = v, w = w, x = x, y = y, z = z), sep = "/")

informloci(dat)

## Not run: 
# Ignore MAF
informloci(dat, MAF = 0)

# Ignore cutoff
informloci(dat, cutoff = 0)

# Real data
data(H3N2)
informloci(H3N2)


## End(Not run)

Description

Check for validity of a genclone or snpclone object

Usage

is.snpclone(x)

is.clone(x)

is.genclone(x)

Arguments

Note

a genclone object will always be a valid genind object and a snpclone object will always be a valid genlight object.

Author(s)

Zhian N. Kamvar

Examples

(sc <- as.snpclone(glSim(100, 1e3, ploid=2, parallel = FALSE), 
                   parallel = FALSE, n.cores = 1L))
is.snpclone(sc)
is.clone(sc)
data(nancycats)
nanclone <- as.genclone(nancycats)
is.genclone(nanclone)

Description

Create a table of summary statistics per locus.

Usage

locus_table(
  x,
  index = "simpson",
  lev = "allele",
  population = "ALL",
  information = TRUE
)

Arguments

Value

a table with 4 columns indicating the Number of alleles/genotypes observed, Diversity index chosen, Nei's 1978 gene diversity (expected heterozygosity), and Evenness.

Note

The calculation of Hexp is $(\frac{n}{n-1}) 1 - \sum_{i = 1}^k{p^{2}_{i}}$ where p is the allele frequencies at a given locus and n is the number of observed alleles (Nei, 1978) in each locus and then returning the average. Caution should be exercised in interpreting the results of Hexp with polyploid organisms with ambiguous ploidy. The lack of allelic dosage information will cause rare alleles to be over-represented and artificially inflate the index. This is especially true with small sample sizes.

If lev = "genotype", then all statistics reflect genotypic diversity within each locus. This includes the calculation for Hexp, which turns into the unbiased Simpson's index.

Author(s)

Zhian N. Kamvar

References

Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, and Helene Wagner. vegan: Community Ecology Package, 2012. R package version 2.0-5.

Niklaus J. Gr\"unwald, Stephen B. Goodwin, Michael G. Milgroom, and William E. Fry. Analysis of genotypic diversity data for populations of microorganisms. Phytopathology, 93(6):738-46, 2003

J.A. Ludwig and J.F. Reynolds. Statistical Ecology. A Primer on Methods and Computing. New York USA: John Wiley and Sons, 1988.

E.C. Pielou. Ecological Diversity. Wiley, 1975.

J.A. Stoddart and J.F. Taylor. Genotypic diversity: estimation and prediction in samples. Genetics, 118(4):705-11, 1988.

Masatoshi Nei. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics, 89(3):583-590, 1978.

Claude Elwood Shannon. A mathematical theory of communication. Bell Systems Technical Journal, 27:379-423,623-656, 1948

See Also

vegan::diversity(), poppr()

Examples

data(nancycats)
locus_table(nancycats[pop = 5])
## Not run: 
# Analyze locus statistics for the North American population of P. infestans.
# Note that due to the unknown dosage of alleles, many of these statistics
# will be artificially inflated for polyploids.
data(Pinf)
locus_table(Pinf, population = "North America")

## End(Not run)

Description

Split samples from a genind object into pseudo-haplotypes

Usage

make_haplotypes(gid)

Arguments

Details

Certain analyses, such as amova work best if within-sample variance (error) can be estimated. Practically, this is performed by splitting the genotypes across all loci to create multiple haplotypes. This way, the within-sample distance can be calculated and incorporated into the model. Please note that the haplotypes generated are based on the order of the unphased alleles in the genind object and do not represent true haplotypes.

Haploid data will be returned un-touched.

Value

a haploid genind object with an extra strata column called "Individual".

Note

The other slot will not be copied over to the new genind object.

See Also

poppr.amova() pegas::amova() as.genambig()

Examples

# Diploid data is doubled -------------------------------------------------

data(nancycats)
nan9 <- nancycats[pop = 9]
nan9hap <- make_haplotypes(nan9) 
nan9              # 9 individuals from population 9
nan9hap           # 18 haplotypes
strata(nan9hap)   # strata gains a new column: Individual
indNames(nan9hap) # individuals are renamed sequentially


# Mix ploidy data can be split, but should be treated with caution --------
# 
# For example, the Pinf data set contains 86 tetraploid individuals, 
# but there appear to only be diploids and triploid genotypes. When 
# we convert to haplotypes, those with all missing data are dropped.
data(Pinf)
Pinf
pmiss <- info_table(Pinf, type = "ploidy", plot = TRUE)

# No samples appear to be triploid across all loci. This will cause
# several haplotypes to have a lot of missing data.
p_haps <- make_haplotypes(Pinf)
p_haps
head(genind2df(p_haps), n = 20)

Description

missingno gives the user four options to deal with missing data: remove loci, remove samples, replace with zeroes, or replace with average allele counts.

Usage

missingno(pop, type = "loci", cutoff = 0.05, quiet = FALSE, freq = FALSE)

Arguments

Details

These methods provide a way to deal with systematic missing data and to give a wrapper for adegenet's tab function. ALL OF THESE ARE TO BE USED WITH CAUTION.

Using this function with polyploid data (where missing data is coded as "0") may give spurious results.

Treatment types

• "ignore" - does not remove or replace missing data.

• "loci" - removes all loci containing missing data in the entire data set.

• "genotype" - removes any genotypes/isolates/individuals with missing data.

• "mean" - replaces all NA's with the mean of the alleles for the entire data set.

• "zero" or "0" - replaces all NA's with "0". Introduces more diversity.

Value

a genclone or genind object.

Note

"wild missingno appeared!"

Author(s)

Zhian N. Kamvar

See Also

tab, poppr, poppr.amova, nei.dist, aboot

Examples

data(nancycats)

nancy.locina <- missingno(nancycats, type = "loci")

## Found 617 missing values.
## 2 loci contained missing values greater than 5%.
## Removing 2 loci : fca8 fca45 

nancy.genona <- missingno(nancycats, type = "geno")

## Found 617 missing values.
## 38 genotypes contained missing values greater than 5%.
## Removing 38 genotypes : N215 N216 N188 N189 N190 N191 N192 N302 N304 N310 
## N195 N197 N198 N199 N200 N201 N206 N182 N184 N186 N298 N299 N300 N301 N303 
## N282 N283 N288 N291 N292 N293 N294 N295 N296 N297 N281 N289 N290  

# Replacing all NA with "0" (see tab in the adegenet package).
nancy.0 <- missingno(nancycats, type = "0")

## Replaced 617 missing values 

# Replacing all NA with the mean of each column (see tab in the
# adegenet package).
nancy.mean <- missingno(nancycats, type = "mean")

## Replaced 617 missing values

Description

Create counts, vectors, and matrices of multilocus genotypes.

Usage

mlg(gid, quiet = FALSE)

mlg.table(
  gid,
  strata = NULL,
  sublist = "ALL",
  exclude = NULL,
  blacklist = NULL,
  mlgsub = NULL,
  bar = TRUE,
  plot = TRUE,
  total = FALSE,
  color = FALSE,
  background = FALSE,
  quiet = FALSE
)

mlg.vector(gid, reset = FALSE)

mlg.crosspop(
  gid,
  strata = NULL,
  sublist = "ALL",
  exclude = NULL,
  blacklist = NULL,
  mlgsub = NULL,
  indexreturn = FALSE,
  df = FALSE,
  quiet = FALSE
)

mlg.id(gid)

Arguments

Details

Multilocus genotypes are the unique combination of alleles across all loci. For details of how these are calculated see vignette("mlg", package = "poppr"). In short, for genind and genclone objects, they are calculated by using a rank function on strings of alleles, which is sensitive to missing data. For genlight and snpclone objects, they are calculated with distance methods via bitwise.dist and mlg.filter, which means that these are insensitive to missing data. Three different types of MLGs can be defined in poppr:

• original the default definition of multilocus genotypes as detailed above

• contracted these are multilocus genotypes collapsed into multilocus lineages (mll) with genetic distance via mlg.filter

• custom user-defined multilocus genotypes. These are useful for information such as mycelial compatibility groups

All of the functions documented here will work on any of the MLG types defined in poppr

Value

mlg

an integer describing the number of multilocus genotypes observed.

mlg.table

a matrix with columns indicating unique multilocus genotypes and rows indicating populations. This table can be used with the funciton diversity_stats to calculate the Shannon-Weaver index (H), Stoddart and Taylor's index (aka inverse Simpson's index; G), Simpson's index (lambda), and evenness (E5).

mlg.vector

a numeric vector naming the multilocus genotype of each individual in the dataset.

mlg.crosspop

• default a list where each element contains a named integer vector representing the number of individuals represented from each population in that MLG

• indexreturn = TRUE a vector of integers defining the multilocus genotypes that have individuals crossing populations

• df = TRUE A long form data frame with the columns: MLG, Population, Count. Useful for graphing with ggplot2

mlg.id

a list of multilocus genotypes with the associated individual names per MLG.

Note

The resulting matrix of 'mlg.table' can be used for analysis with the vegan package.

mlg.vector will recalculate the mlg vector for [adegenet::genind] objects and will return the contents of the mlg slot in [genclone][genclone-class] objects. This means that MLGs will be different for subsetted [adegenet::genind] objects.

Author(s)

Zhian N. Kamvar

See Also

vegan::diversity() diversity_stats popsub mll mlg.filter mll.custom

Examples

# Load the data set
data(Aeut)

# Investigate the number of multilocus genotypes.
amlg <- mlg(Aeut)
amlg # 119

# show the multilocus genotype vector 
avec <- mlg.vector(Aeut)
avec 

# Get a table
atab <- mlg.table(Aeut, color = TRUE)
atab

# See where multilocus genotypes cross populations
acrs <- mlg.crosspop(Aeut) # MLG.59: (2 inds) Athena Mt. Vernon

# See which individuals belong to each MLG
aid <- mlg.id(Aeut)
aid["59"] # individuals 159 and 57

## Not run: 

# For the mlg.table, you can also choose to display the number of MLGs across
# populations in the background

mlg.table(Aeut, background = TRUE)
mlg.table(Aeut, background = TRUE, color = TRUE)

# A simple example. 10 individuals, 5 genotypes.
mat1 <- matrix(ncol=5, 25:1)
mat1 <- rbind(mat1, mat1)
mat <- matrix(nrow=10, ncol=5, paste(mat1,mat1,sep="/"))
mat.gid <- df2genind(mat, sep="/")
mlg(mat.gid)
mlg.vector(mat.gid)
mlg.table(mat.gid)

# Now for a more complicated example.
# Data set of 1903 samples of the H3N2 flu virus genotyped at 125 SNP loci.
data(H3N2)
mlg(H3N2, quiet = FALSE)

H.vec <- mlg.vector(H3N2)

# Changing the population vector to indicate the years of each epidemic.
pop(H3N2) <- other(H3N2)$x$country
H.tab <- mlg.table(H3N2, plot = FALSE, total = TRUE)

# Show which genotypes exist accross populations in the entire dataset.
res <- mlg.crosspop(H3N2, quiet = FALSE)

# Let's say we want to visualize the multilocus genotype distribution for the
# USA and Russia
mlg.table(H3N2, sublist = c("USA", "Russia"), bar=TRUE)

# An exercise in subsetting the output of mlg.table and mlg.vector.
# First, get the indices of each MLG duplicated across populations.
inds <- mlg.crosspop(H3N2, quiet = FALSE, indexreturn = TRUE)

# Since the columns of the table from mlg.table are equal to the number of
# MLGs, we can subset with just the columns.
H.sub <- H.tab[, inds]

# We can also do the same by using the mlgsub flag.
H.sub <- mlg.table(H3N2, mlgsub = inds)

# We can subset the original data set using the output of mlg.vector to
# analyze only the MLGs that are duplicated across populations. 
new.H <- H3N2[H.vec %in% inds, ]


## End(Not run)

Description

Multilocus genotypes are initially defined by naive string matching, but this definition does not take into account missing data or genotyping error, casting these as unique genotypes. Defining multilocus genotypes by genetic distance allows you to incorporate genotypes that have missing data o genotyping error into their parent clusters.

Usage

mlg.filter(
  pop,
  threshold = 0,
  missing = "asis",
  memory = FALSE,
  algorithm = "farthest_neighbor",
  distance = "diss.dist",
  threads = 1L,
  stats = "MLGs",
  ...
)

mlg.filter(
  pop,
  missing = "asis",
  memory = FALSE,
  algorithm = "farthest_neighbor",
  distance = "diss.dist",
  threads = 1L,
  ...
) <- value

Arguments

Details

This function will take in any distance matrix or function and collapse multilocus genotypes below a given threshold. If you use this function as the assignment method (mlg.filter(myData, distance = myDist) <- 0.5), the distance function or matrix will be remembered by the object. This means that if you define your own distance matrix or function, you must keep it in memory to further utilize mlg.filter.

Value

Default, a vector of collapsed multilocus genotypes. Otherwise, any combination of the following:

MLGs

a numeric vector defining the multilocus genotype cluster of each individual in the dataset. Each genotype cluster is separated from every other genotype cluster by at least the defined threshold value, as calculated by the selected algorithm.

THRESHOLDS

A numeric vector representing the thresholds beyond which clusters of multilocus genotypes were collapsed.

DISTANCES

A square matrix representing the distances between each cluster.

SIZES

The sizes of the multilocus genotype clusters in order.

Note

mlg.vector makes use of mlg.vector grouping prior to applying the given threshold. Genotype numbers returned by mlg.vector represent the lowest numbered genotype (as returned by mlg.vector) in in each new multilocus genotype. Therefore mlg.filter and mlg.vector return the same vector when threshold is set to 0 or less.

See Also

filter_stats, cutoff_predictor, mll, genclone, snpclone, diss.dist, bruvo.dist

Examples

data(partial_clone)
pc <- as.genclone(partial_clone, threads = 1L) # convert to genclone object

# Basic Use ---------------------------------------------------------------


# Show MLGs at threshold 0.05
mlg.filter(pc, threshold = 0.05, distance = "nei.dist", threads = 1L)
pc # 26 mlgs

# Set MLGs at threshold 0.05
mlg.filter(pc, distance = "nei.dist", threads = 1L) <- 0.05
pc # 25 mlgs

## Not run: 

# The distance definition is persistant
mlg.filter(pc) <- 0.1
pc # 24 mlgs

# But you can still change the definition
mlg.filter(pc, distance = "diss.dist", percent = TRUE) <- 0.1
pc

# Choosing a threshold ----------------------------------------------------


# Thresholds for collapsing multilocus genotypes should not be arbitrary. It
# is important to consider what threshold is suitable. One method of choosing
# a threshold is to find a gap in the distance distribution that represents
# clonal groups. You can look at this by analyzing the distribution of all
# possible thresholds with the function "cutoff_predictor".

# For this example, we'll use Bruvo's distance to predict the cutoff for
# P. infestans.

data(Pinf)
Pinf
# Repeat lengths are necessary for Bruvo's distance
(pinfreps <- fix_replen(Pinf, c(2, 2, 6, 2, 2, 2, 2, 2, 3, 3, 2)))

# Now we can collect information of the thresholds. We can set threshold = 1
# because we know that this will capture the maximum possible distance:
(thresholds <- mlg.filter(Pinf, distance = bruvo.dist, stats = "THRESHOLDS",
                          replen = pinfreps, threshold = 1))
# We can use these thresholds to find an appropriate cutoff
(pcut <- cutoff_predictor(thresholds))
mlg.filter(Pinf, distance = bruvo.dist, replen = pinfreps) <- pcut
Pinf

# This can also be visualized with the "filter_stats" function.

# Special case: threshold = 0 ---------------------------------------------


# It's important to remember that a threshold of 0 is equal to the original
# MLG definition. This example will show a data set that contains genotypes
# with missing data that share all alleles with other genotypes except for 
# the missing one.

data(monpop)
monpop # 264 mlg
mlg.filter(monpop) <- 0
nmll(monpop) # 264 mlg

# In order to merge these genotypes with missing data, we should set the 
# threshold to be slightly higher than 0. We will use the smallest fraction 
# the computer can store.

mlg.filter(monpop) <- .Machine$double.eps ^ 0.5
nmll(monpop) # 236 mlg

# Custom distance ---------------------------------------------------------

# Custom genetic distances can be used either in functions from other
# packages or user-defined functions

data(Pinf)
Pinf
mlg.filter(Pinf, distance = function(x) dist(tab(x))) <- 3
Pinf
mlg.filter(Pinf) <- 4
Pinf

# genlight / snpclone objects ---------------------------------------------


set.seed(999)
gc <- as.snpclone(glSim(100, 0, n.snp.struc = 1e3, ploidy = 2))
gc # 100 mlgs
mlg.filter(gc) <- 0.25
gc # 82 mlgs


## End(Not run)

Description

The following methods allow the user to access and manipulate multilocus lineages in genclone or snpclone objects.

Usage

mll(x, type = NULL)

nmll(x, type = NULL)

mll(x) <- value

Arguments

Details

genclone and snpclone objects have a slot for an internal class of object called MLG. This class allows the storage of flexible mll definitions:

• "original" - naive mlgs defined by string comparison. This is default.

• "contracted" - mlgs defined by a genetic distance threshold.

• "custom" - user-defined MLGs

Value

an object of the same type as x.

Author(s)

Zhian N. Kamvar

See Also

mll.custom mlg.table

Examples

data(partial_clone)
pc <- as.genclone(partial_clone)
mll(pc)
mll(pc) <- "custom"
mll(pc)
mll.levels(pc) <- LETTERS
mll(pc)

Description

This function will allow you to define custom multilocus lineages for your data set.

Usage

mll.custom(x, set = TRUE, value)

mll.custom(x, set = TRUE) <- value

mll.levels(x, set = TRUE, value)

mll.levels(x, set = TRUE) <- value

Arguments

Value

an object of the same type as x

Author(s)

Zhian N. Kamvar

See Also

mll mlg.table

Examples

data(partial_clone)
pc <- as.genclone(partial_clone)
mll.custom(pc) <- LETTERS[mll(pc)]
mll(pc)

# Let's say we had a mistake and the A mlg was actually B. 
mll.levels(pc)[mll.levels(pc) == "A"] <- "B"
mll(pc)

# Set the MLL back to the original definition.
mll(pc) <- "original"
mll(pc)

Description

This function will allow you to reset multilocus lineages for your data set.

Usage

mll.reset(x, value)

Arguments

Value

an object of the same type as x

Note

This method has no assignment method. If "original" is not contained in "value", it is assumed that the "original" definition will be used to reset the MLGs.

Author(s)

Zhian N. Kamvar

See Also

mll mlg.table mll.custom

Examples

# This data set was a subset of a larger data set, so the multilocus
# genotypes are not all sequential
data(Pinf)
mll(Pinf) <- "original"
mll(Pinf)

# If we use mll.reset, then it will become sequential
Pinf.new <- mll.reset(Pinf, TRUE) # reset all
mll(Pinf.new)

## Not run: 

# It is possible to reset only specific mll definitions. For example, let's
# say that we wanted to filter our multilocus genotypes by nei's distance
mlg.filter(Pinf, dist = nei.dist, missing = "mean") <- 0.02

# And we wanted to set those as custom genotypes,
mll.custom(Pinf) <- mll(Pinf, "contracted")
mll.levels(Pinf) <- .genlab("MLG", nmll(Pinf, "custom"))

# We could reset just the original and the filtered if we wanted to and keep
# the custom as it were.

Pinf.new <- mll.reset(Pinf, c("original", "contracted"))

mll(Pinf.new, "original")
mll(Pinf.new, "contracted")
mll(Pinf.new, "custom")

# If "original" is not one of the values, then that is used as a baseline.
Pinf.orig <- mll.reset(Pinf, "contracted")
mll(Pinf.orig, "contracted")
mll(Pinf.new, "contracted")

## End(Not run)