bin read depths of SNPs into categories having at least S observations

Description

bin read depths of SNPs into categories having at least S observations

Usage

bin_depths(D, S)

Arguments

Value

This returns a list with two components. dp_bins is a matrix of the same shape as D with the bin categories (as 1, 2, ...) and -1 for this cells corresponding to missing genotypes. num_cats is the number of depth bins. tidy_bins is a long format description of the bins. bin_stats is a tibble giving summary information about the read depth bins which is useful for plotting things, etc.

Examples

# get a matrix of read depths and make it an integer matrix
depths <- vcfR::extract.gt(lobster_buz_2000, element = "DP")
storage.mode(depths) <- "integer"

# get a character matrix of genotypes, so we can figure out which
# are missing and mask those from depths
genos <- vcfR::extract.gt(lobster_buz_2000, element = "GT")

# make missing in depths if missing in genos
depths[is.na(genos)] <- NA

# bin the read depths into bins with at least 1000 observations in each bin
bins <- bin_depths(depths, 1000)

Estimate heterozygote miscall rate for different read depth categories (no nulls)

Description

To see how this Rcpp function is applied, see the code in infer_m.

Usage

estimate_m_rd(Y, R, init_m, num_cats, p_prior, m_prior, num_reps)

Arguments

Computed expected and observed genotype frequencies from a ‘vcfR’ object

Description

Under the assumption of Hardy-Weinberg equilibrium, this function uses the estimated allele frequencies from the data set in v to compute expected genotype frequencies, and then reports these along with the observed genotype frequencies. Loci come out named as CHROM–POS.

Usage

exp_and_obs_geno_freqs(
  v = NULL,
  d012 = NULL,
  prop_indv_required = 0.5,
  prop_loci_required = 0.5
)

Arguments

Value

Returns a tibble with the following columns: snp = the locus name as CHROM–POS; p = The frequency of the alternate (ALT) allele; ntot = the total number of individuals with no missing data at the locus; geno = column telling which genotype (0, 1, or 2) is referred to; p_exp = expected frequency of the genotype; p_obs = observed frequency of genotype; n_exp = expected number of such genotypes; n_obs = observed number of such genotypes; z_score = simple statistic giving how far the observed genotype frequency is from that expected under Hardy-Weinberg equilibrium.

Examples

eao <- exp_and_obs_geno_freqs(v = lobster_buz_2000)

# if you wanted to run that on an 012 matrix,
# it would be like this:
eao012 <- exp_and_obs_geno_freqs(d012 = lobster_buz_2000_as_012_matrix)

function to return a tibble with the min/max values of genotype freqs possible

Description

These mins and maxes occur because the genotypes are used to estimate the allele frequencies. The output of this function is useful for putting boundaries on plots.

Usage

geno_freq_boundaries()

Value

a tibble

Examples

gfb <- geno_freq_boundaries()

from on 012 file compute expected (assuming HW equilbrium) and observed genotype counts

Description

This is an internal function.

Usage

geno_freq_calc_single(g012)

Arguments

Value

See information for return value of exp_and_obs_geno_freqs.

Examples

# get an 012 matrix from the lobster data
tmp <- t(vcfR::extract.gt(lobster_buz_2000, element = "GT"))
locnames <- colnames(tmp)
g <- make_it_012(tmp)
colnames(g) <- locnames # put these back on since make_it_012 removes them
gf <- geno_freq_calc_single(g)

return a ‘ggplot2’ plot object of observed and expected genotype freqs

Description

return a ‘ggplot2’ plot object of observed and expected genotype freqs

Usage

geno_freqs_scatter(gfc, alpha = 0.2, max_plot_loci = 500)

Arguments

Examples

# get the expected and observed geno freqs
gfreqs <- exp_and_obs_geno_freqs(lobster_buz_2000)
g <- geno_freqs_scatter(gfreqs)

# now g is a 'ggplot2' object.

simulate a new miscall rate for each read depth category given X and Y

Description

This just writes new values into M as if it were an output variable. Internal function used within estimate_m_rd.

Usage

gibbsM(M, num_cats, X, Y, R, pri)

simulate new reference allele frequencies from their beta full conditional

Description

This just writes new values into P as if it were an output variable. Internal function used within estimate_m_rd.

Usage

gibbsP(p, X, pri)

compute full conditional for each X given Y, p, R, and m, and then sample from it

Description

This just writes new values into X as if it were an output variable. Internal function used within estimate_m_rd.

Usage

gibbsX(X, Y, R, p, M)

get posterior estimates for m from different read depth categories

Description

This function calls internal C++ routines that perform Markov chain Monte Carlo (MCMC) to sample from the posterior distribution of the heterozygote miscall rate for each read depth category.

Usage

infer_m(
  v,
  minBin,
  indivs = NULL,
  init_m = 0.1,
  num_sweeps = 500,
  burn_in = 100
)

Arguments

Details

The read depth bins are determined by passing to the function minBin—the minimum number of observations desired for each read depth bin. The function then breaks the observations into bins so that each read depth bin has at least minBin observations.

Note that if you want to estimate the heterozygote miscall rate overall (i.e., not conditioning each estimate on a read depth bin), then simply give a very large number (larger than the number of markers times the number of individuals) for minBin. For example, you could use a number like 1e15 for minBin. As a consequence, all the genotypes will be put into a single read depth bin.

Value

A list with six components:

m_posteriors

A tibble with 6 columns: bin = the index of the read depth bin; mean = the posterior mean estimate of the the heterozygote miscall rate in that bin; lo95 = the low endpoint of the 95 95 and mean_dp = the mean read depth in the bin.

m_traces

A tibble with all the values visited for m for every read depth bin. This tibble has three columns: bin = the index of the read depth bin; sweep = the sweep number, value = the value of m for that read depth bin in that particular sweep.

dp_summary

A tibble summarizing how many genotypes of different read depths appear in each bin.

bin_stats

A tibble with a different summary of the read-depth bins.

num_sweeps

Number of MCMC sweeps used.

burn_in

Number of sweeps discarded as burn in.

Examples

# Shorter run than recommended (for quick example...)
im <- infer_m(lobster_buz_2000, minBin = 1000, num_sweeps = 100, burn_in = 20)

Restriction-associated digest (RAD) sequence data from 36 lobsters at 2000 single nucleotide polymorphisms (SNPs)

Description

A data set from the study by Benestan et al. (2016).

Usage

lobster_buz_2000

Format

A ‘vcfR’ object in which the original data set has been reduced to just the 36 lobsters from the BUZ population and a randomly sampled 2000 SNPs from the 10,156 originally available.

This is the sort of object obtained after calling vcfR::read.vcfR() on a variant call format (VCF) file.

Source

https://datadryad.org/stash/dataset/doi:10.5061/dryad.q771r

An 012 matrix from (RAD) sequence data from 36 lobsters at 2000 single nucleotide polymorphisms (SNPs)

Description

A data set from the study by Benestan et al. (2016).

Usage

lobster_buz_2000_as_012_matrix

Format

This is an integer matrix with positions in rows and samples in columns. It is an 012 matrix that corresponds to lobster_buz_2000. It is useful as an example of the necessary format for the d012 argument to exp_and_obs_geno_freqs.

Source

https://datadryad.org/stash/dataset/doi:10.5061/dryad.q771r

just a quick function for making an 012 matrix from a character matrix

Description

The standard way within R of pulling values out of a named vector really bogs down on large data sets. So I will do this instead.

Usage

make_it_012(M)

Arguments

Value

An integer matrix of values which are 0, 1, 2, or -1.

Examples

# get an 012 matrix from the lobster data
tmp <- t(vcfR::extract.gt(lobster_buz_2000, element = "GT"))
dimnames(tmp) <- NULL
g <- make_it_012(tmp)

return a ‘ggplot2’ plot object of the posterior esimates for heterozyote miscall rates

Description

This just returns a ‘ggplot2’ plot object that plots the read depth bins on the x-axis and the posterior mean m estimates (and credible intervals) on the y-axis, and depicts the number of genotypes in each read depth bin using color.

Usage

posteriors_plot(P)

Arguments

Value

a ‘ggplot2’ plot object.

Examples

# get something to plot (short run for example)
im <- infer_m(lobster_buz_2000, minBin = 1000, num_sweeps = 100, burn_in = 20)

# then plot it
g <- posteriors_plot(im$m_posteriors)

# now g is a 'ggplot2' object

convert a VCF into an 012,-1 matrix and read_depth bin matrix for estimation

Description

convert a VCF into an 012,-1 matrix and read_depth bin matrix for estimation

Usage

prep_vcf_for_est_m_rd(v, DF = "DP", minBin)

Arguments

Value

This sends back a list with mat012: the 012,-1 matrix of genotypes, and dp_bins_list: the list returned by bin_depths().

Examples

pv <- prep_vcf_for_est_m_rd(lobster_buz_2000, minBin = 1000)

tidy up the estimate_m_rd output into something you can plot

Description

Internal function. See use in infer_m.

Usage

tidy_m_ests(E, S, TB, burn = 50)

Arguments

whoa: Evaluation of genotyping error in genotype by sequencing data

Description

The name is found in the capitals here: Where's my Heterozygotes?! Observations of genotyping Accuracy.

the whoa main high-level functions

The main function in the package whoa is infer_m. This function infers the heterozygote miscall rate (the rate at which true heterozygotes have been miscalled as homozygotes in genotype-by-sequencing data) for calls made upon genotypes falling within different read depth categories.

The output from infer_m is easily plotted by passing the m_posteriors component of the output list from infer_m into posteriors_plot.

example data

The package comes with a small data set, lobster_buz_2000, which was read in from a variant call format (VCF) file and is now stored in the package as a ‘vcfR’ object.