est_maxLnL: Numerically Estimating the Maximum LnL parameter values from...
In tylerhether/HMMancestry: R package using the Forward-Backward algorithm to infer genotypes, recombination hotspots, and gene conversion tracts from low-coverage next-generation sequence data

Description Usage Arguments Details Value Author(s) References See Also Examples

Infers the genome-wide recombination rate (Morgans / bp) and the assignment probabilty directly from the data using maximum likelihood estimates.

1
2
3

est_maxLnL(dat, ploidy = "haploid", initial_p_assign = "NULL",
  initial_scale = "NULL", tolerance = 0.001, n_coarse_steps = 5,
  n_iterations = 30, plot = FALSE, n_cores = 1)

`dat`	a data.frame with five columns: `Ind` The individual ID `Chr` The chromosome ID `Snp` The (sorted) snp location (in bps) `p0` The number of reads that mapped to parent 0 `p1` The number of reads that mapped to parent 1 Note that the names of the columns can be different and there can be additional columns to the right (`est_maxLnL` only uses the first 5 columns of the data.frame).
`ploidy`	The ploidy to analyze. If `ploidy="diploid"` (default) then `fb_diploid` is used to infer states and estimate likelihood. If `ploidy="haploid"` then `fb_haploid` will be used.
`initial_p_assign`	The initial assignment probability (see details). If "NULL" (default), then a coarse search will be performed with `n_coarse_steps` equally spaced from 0.95 to 0.999.
`initial_scale`	The initial genome-wide recombination rate (Morgans / bp). If "NULL" (default), then a coarse search will be performed `n_coarse_steps` equally spaced from 1e-06 to 1e-04.
`tolerance`	The tolerance used to stop the search.
`n_coarse_steps`	The size of the 1D (if either `initial_p_assign` or `initial_scale` is "NULL") or 2D grid (if both are "NULL"). Increasing `n_coarse_steps` can greatly increase computation time.
`n_iterations`	The number of iterations during the fine scale parameter estimate serach (see details).
`n_cores`	Integer specifying the number of processors to use

est_maxLnL has a coarse and fine scale method to estimate the genome-wide recombination rate, \hat{c}, and assignement probability, \hat{p}. The coarse scale uses either fb_haploid or fb_diploid output and estimates the natural log-likelihood, LnL, of the data across a coarse (in parameter space) grid of varying p or c values and the fine scale uses a two-variable Newton-Raphson method to hone in on a closer estimate. For further details see Hether et al. (in prep).

A data.frame containing the following columns:

p_assign_hat An estimate of the assignment probability
scale_hat An estimate of mean recombination rate
n_iterations The number of iterations carried out

Tyler D. Hether

Hether, T.D., C. G. Wiench1, and P.A. Hohenlohe (in review). 2015. Novel molecular and analytical tools for efficient estimation of rates of meiotic crossover, non-crossover and gene conversion

P. Deuflhard, Newton Methods for Nonlinear Problems. Affine Invariance and Adaptive Algorithms. Springer Series in Computational Mathematics, Vol. 35. Springer, Berlin, 2004.

fb_haploid, fb_haploid

# Simulating 30 haploid recombinants
set.seed(1234567)        # For reproducibility
n_spores <- 30           # number of recombinants
l <- 1000                # number of snps to simulate
c <- 1e-06               # recombination rate between snps (Morgan/bp)
snps <- c(1:l)*1e2       # snps are evenly spaced 100bps apart
p_a <- 0.985             # assignment probability
coverage <- 1            # mean coverage
# Now simulate
sim <- sim_en_masse(n.spores=n_spores, scale=c, snps=snps,
	p.assign=p_a, mu.rate=0, f.cross=1, f.convert=0,
	length.conversion=0, coverage=coverage)
# Now estmate params
res <- est_maxLnL(sim, ploidy="haploid", plot=TRUE, n_cores=5)
res