estimateSH: Estimate Selection Coefficient and Dominance

In ThomasTaus/poolSeq: Simulate and Analyze Pool-seq Data

Description

Based on a set of allele frequency trajectories and the effective population size both the selection coefficient (s) and dominance (h) can be estimated.

Usage

estimateSH(traj, t, Ne, haploid = FALSE, h = NA, N.ctraj = 0, simulate.p.value = FALSE, N.pval = 1000, cov = NA, approximate = TRUE, method=c("LLS", "NLS", "automatic"))

## S3 method for class 'estsh'
print(x, ...)

## S3 method for class 'estsh'
coef(object, ...)

## S3 method for class 'estsh'
confint(object, parm, level = 0.95, N.ci = 1000, warn = FALSE, ...)

Arguments

x, object	object of class "estsh" as returned by estimateSH.
traj	numeric matrix of allele frequency trajectories, where columns and rows correspond to measurement time points and replicates, respectively.
t	numeric vector indicating the number of generations for each column in traj.
Ne	numeric specifying the effective population size
haploid	logical indicating wether a population of haploid (haploid = TRUE) or diploid (default) individuals should be simulated.
h	numeric defining the dominance that should be assumed. If h = NA (default) dominance will be estimated alongside with the selection coefficient.
N.ctraj	numeric determining the number of simulations that should be performed to correct for bias in the consensus trajectory if both p0 and Ne are small. By default, no bias correction is performed (sim = 0). If a bias correction should be performed, a value of sim = 1000 is recommended.
simulate.p.value	locical indicating whether the p-value should be estimated by simulations.
N.pval	integer, the number of simulations performed to estimate the p-value (by default N.pval = 1000).
cov	numeric matrix of sequence coverages, where columns and rows correspond to measurement time points and replicates, respectively
approximate	logical. If approximate = true a continuous time approximation is used to estimate s (and h). Alternatively, forward in time simulations are performed. The differences between continuous and discrete models increase with s.
method	character string defining the method that should be used to estimate s (and p0 or h).
parm	a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered.
level	the confidence level required
N.ci	integer specifying the number of simulations that should be performed to get the bootstrap confidence interval.
warn	logical specifying if a warning should be sent in case the confidence interval may be inaccurate. Otherwise (default) NA will be returned in such cases.
...	currently not used.

Value

A list with class "estsh" containing the following components:

traj	matrix with observed allele frequency trajectories.
t	generations of each measurement time point.
Ne	effective population size.
haploid	logical for diploid or haploid individuals.
ctraj	consensus trajectory used to estimate s and h.
N.ctraj	number of simulations to correct 'ctraj' from bias.
s	estimated selection coefficient.
h	if h == NA, estimated dominance.
h.given	dominance defined in h.
cov	matrix with observed sequence coverages.
p.value	p-value estimated by simulations under neutrality.
N.pval	number of simulations to estimate the p-value.

Author(s)

Thomas Taus

See Also

consensus.traj and compareSH

Examples

# simulate allele frequency trajectories for a given selection coefficient and dominance
tp <- seq(0, 100, by=10)
Ne <- 1000
traj <- wf.traj(p0=rep(0.05, times=10), Ne=Ne, t=tp, s=0.1, h=0.5, haploid=FALSE)

# estimate only selection coefficient, assuming that dominance is 0.5
estimateSH(traj, t=tp, Ne=Ne, haploid=FALSE, h=0.5)

# estimate both selection coefficient and dominance for simulated data
estimateSH(traj, t=tp, Ne=Ne, haploid=FALSE, method="NLS")

# compute p-value for estimated selection coefficient
estimateSH(traj, t=tp, Ne=Ne, haploid=FALSE, h=0.5, simulate.p.value=TRUE) ##### takes about 15 sec #####

# estimate confidence interval for both s and h
res <- estimateSH(traj, t=tp, Ne=Ne, haploid=FALSE, h=0.5)
confint(res)  ##### takes about 15 sec #####

ThomasTaus/poolSeq documentation built on Feb. 17, 2020, 1:52 p.m.