calcFst: Estimate Population Differentiation Statistics

In polysat: Tools for Polyploid Microsatellite Analysis

Estimate Population Differentiation Statistics

Description

Given a data frame of allele frequencies and population sizes, calcPopDiff calculates a matrix of pairwise FST, GST, Jost's D, or RST values, or a single global value for any of these four statistics. calcFst is a wrapper for calcPopDiff to allow backwards compatibility with previous versions of polysat.

Usage

calcPopDiff(freqs, metric, pops = row.names(freqs), 
            loci = unique(gsub("\\..*$", "", names(freqs))), global = FALSE,
            bootstrap = FALSE, n.bootstraps = 1000, object = NULL)

calcFst(freqs, pops = row.names(freqs),
        loci = unique(gsub("\\..*$", "", names(freqs))), ...)

Arguments

freqs	A data frame of allele frequencies and population sizes such as that produced by simpleFreq or deSilvaFreq. Each population is in one row, and a column called Genomes (or multiple columns containing the locus names and “Genomes” seperated by a period) contains the relative size of each population. All other columns contain allele frequencies. The names of these columns are the locus name and allele name, separated by a period.
metric	The population differentiation statistic to estimate. Can be “Fst”, “Gst”, “Jost's D”, or “Rst”.
pops	A character vector. Populations to analyze, which should be a subset of row.names(freqs).
loci	A character vector indicating which loci to analyze. These should be a subset of the locus names as used in the column names of freqs.
global	Boolean. If TRUE, a single global statistic will be estimated across all populations. If FALSE, pairwise statistics will be estimated between populations.
bootstrap	Boolean. If TRUE, a set of replicates bootstrapped across loci will be returned. If FALSE, a single value will be returned for each pair of populations (if global = FALSE) or for the whole set (if global = TRUE).
n.bootstraps	Integer. The number of bootstrap replicates to perform. Ignored if bootstrap = FALSE.
object	A "genambig" or "genbinary" object with the Usatnts slot filled in. Required for metric = "Rst" only.
...	Additional arguments to be passed to calcPopDiff (i.e. global, bootstrap, and/or n.bootstraps).

Details

For metric = "Fst" or calcFst:

HS and HT are estimate directly from allele frequencies for each locus for each pair of populations, then averaged across loci. Wright's FST is then calculated for each pair of populations as (HT-HS)/HT.

H values (expected heterozygosities for populations and combined populations) are calculated as one minus the sum of all squared allele frequencies at a locus. To calculte HT, allele frequencies between two populations are averaged before the calculation. To calculate HS, H values are averaged after the calculation. In both cases, the averages are weighted by the relative sizes of the two populations (as indicated by freqs$Genomes).

For metric = "Gst":

This metric is similar to FST, but mean allele frequencies and mean expected heterozygosities are not weighted by population size. Additionally, unbiased estimators of HS and HT are used according to Nei and Chesser (1983; equations 15 and 16, reproduced also in Jost (2008)). Instead of using twice the harmonic mean of the number of individuals in the two subpopulations as described by Nei and Chesser (1983), the harmonic mean of the number of allele copies in the two subpopulations (taken from freq$Genomes) is used, in order to allow for polyploidy. GST is estimated for each locus and then averaged across loci.

For metric = "Jost's D":

The unbiased estimators of HS and HT and calculated as with GST, without weighting by population size. They are then used to estimate D at each locus according to Jost (2008; equation 12):

2 * (HT - HS)/(1 - HS)

Values of D are then averaged across loci.

For metric = "Rst":

RST is estimated on a per-locus basis according to Slatkin (1995), but with populations weighted equally regardless of size. Values are then averaged across loci.

For each locus:

S_w = 1/d * ∑_j^d (∑_i ∑_{i' < i} p_ij * p_i'j * n_j^2 * (a_i - a_i')^2)/(n_j * (n_j - 1)))

S = (∑_i ∑_{i' < i} p_i * p_i' * n^2 * (a_i - a_i')^2)/(n * (n - 1))

RST = (S - S_w)/S

where d is the number of populations, j is an individual population, i and i' are individual alleles, p_ij is the frequency of an allele in a population, n_j is the number of allele copies in a population, a_i is the size of an allele expressed in number of repeat units, p_i is an allele frequency averaged across populations (with populations weighted equally), and n is the total number of allele copies across all populations.

Value

If global = FALSE and bootstrap = FALSE, a square matrix containing FST, GST, D, or RST values. The rows and columns of the matrix are both named by population.

If global = TRUE and bootstrap = FALSE, a single value indicating the FST, GST, D, or RST value.

If global = TRUE and bootstrap = TRUE, a vector of bootstrapped replicates of FST, GST, D, or RST.

If global = FALSE and bootstrap = TRUE, a square two-dimensional list, with rows and columns named by population. Each item is a vector of bootstrapped values for FST, GST, D, or RST for that pair of populations.

Author(s)

Lindsay V. Clark

References

Nei, M. (1973) Analysis of gene diversity in subdivided populations. Proceedings of the National Academy of Sciences of the United States of America 70, 3321–3323.

Nei, M. and Chesser, R. (1983) Estimation of fixation indices and gene diversities. Annals of Human Genetics 47, 253–259.

Jost, L. (2008) GST and its relatives to not measure differentiation. Molecular Ecology 17, 4015–4026.

Slatkin, M. (1995) A measure of population subdivision based on microsatellite allele frequencies. Genetics 139, 457–462.

See Also

simpleFreq, deSilvaFreq

Examples

# create a data set (typically done by reading files)
mygenotypes <- new("genambig", samples = paste("ind", 1:6, sep=""),
                   loci = c("loc1", "loc2"))
Genotypes(mygenotypes, loci = "loc1") <- list(c(206), c(208,210),
                                              c(204,206,210),
    c(196,198,202,208), c(196,200), c(198,200,202,204))
Genotypes(mygenotypes, loci = "loc2") <- list(c(130,134), c(138,140),
                                              c(130,136,140),
    c(138), c(136,140), c(130,132,136))
PopInfo(mygenotypes) <- c(1,1,1,2,2,2)
mygenotypes <- reformatPloidies(mygenotypes, output="sample")
Ploidies(mygenotypes) <- c(2,2,4,4,2,4)
Usatnts(mygenotypes) <- c(2,2)

# calculate allele frequencies
myfreq <- simpleFreq(mygenotypes)

# calculate pairwise differentiation statistics
myfst <- calcPopDiff(myfreq, metric = "Fst")
mygst <- calcPopDiff(myfreq, metric = "Gst")
myD <- calcPopDiff(myfreq, metric = "Jost's D")
myrst <- calcPopDiff(myfreq, metric = "Rst", object = mygenotypes)

# examine the results
myfst
mygst
myD
myrst

# get global statistics
calcPopDiff(myfreq, metric = "Fst", global = TRUE)
calcPopDiff(myfreq, metric = "Gst", global = TRUE)
calcPopDiff(myfreq, metric = "Jost's D", global = TRUE)
calcPopDiff(myfreq, metric = "Rst", global = TRUE, object = mygenotypes)

polysat documentation built on Aug. 23, 2022, 5:07 p.m.