QTLseqr like G prime analysis — runGprimeAnalysis

Identify QTL using a smoothed G statistic

Usage

runGprimeAnalysis_local(
  SNPset,
  windowSize = 1e+06,
  outlierFilter = "deltaSNP",
  filterThreshold = 0.1,
  ...
)

Arguments

SNPset: Data frame SNP set containing previously filtered SNPs
windowSize: the window size (in base pairs) bracketing each SNP for which to calculate the statitics. Magwene et. al recommend a window size of ~25 cM, but also recommend optionally trying several window sizes to test if peaks are over- or undersmoothed. Default: 1e+06
outlierFilter: one of either "deltaSNP" or "Hampel". Method for filtering outlier (ie QTL) regions for p-value estimation. Default: 'deltaSNP'
filterThreshold: The absolute delta SNP index to use to filter out putative QTL (default = 0.1)
...: Other arguments passed to locfit and subsequently to locfit.raw() (or the lfproc). Usefull in cases where you get "out of vertex space warnings"; Set the maxk higher than the default 100. See locfit.raw(). But if you are getting that warning you should seriously consider increasing your window size. Default: 0.1

Value

The supplied SNP set tibble after G' analysis. Includes five new columns:

Gprime - The tricube smoothed G statistic based on the supplied window size
pvalue - the pvalue at each SNP calculatd by non-parametric estimation
negLog10Pval - the -Log10(pvalue) supplied for quick plotting
qvalue - the Benajamini-Hochberg adjusted p-value

Details

A wrapper for all the functions that perform the full G prime analysis to identify QTL. The following steps are performed:
1) Genome-wide G statistics are calculated by getG.
G is defined by the equation: $$G = 2*\sum_{i=1}^{4} n_{i}*ln\frac{obs(n_i)}{exp(n_i)}$$ Where for each SNP, $n_i$ from i = 1 to 4 corresponds to the reference and alternate allele depths for each bulk, as described in the following table:

Allele	High Bulk	Low Bulk
Reference	$n_1$	$n_2$
Alternate	$n_3$	$n_4$

...and $obs(n_i)$ are the observed allele depths as described in the data frame. getG calculates the G statistic using expected values assuming read depth is equal for all alleles in both bulks: $$exp(n_1) = ((n_1 + n_2)*(n_1 + n_3))/(n_1 + n_2 + n_3 + n_4)$$ $$exp(n_2) = ((n_2 + n_1)*(n_2 + n_4))/(n_1 + n_2 + n_3 + n_4)$$ $$exp(n_3) = ((n_3 + n_1)*(n_3 + n_4))/(n_1 + n_2 + n_3 + n_4)$$ $$exp(n_4) = ((n_4 + n_2)*(n_4 + n_3))/(n_1 + n_2 + n_3 + n_4)$$
2) G'

A tricube-smoothed G statistic is predicted by local regression within each chromosome using tricubeStat. This works as a weighted average across neighboring SNPs that accounts for Linkage disequilibrium (LD) while minizing noise attributed to SNP calling errors. G values for neighboring SNPs within the window are weighted by physical distance from the focal SNP.

3) P-values are estimated based using the non-parametric method described by Magwene et al. 2011 with the function getPvals. Breifly, using the natural log of Gprime a median absolute deviation (MAD) is calculated. The Gprime set is trimmed to exclude outlier regions (i.e. QTL) based on Hampel's rule. An alternate method for filtering out QTL is proposed using absolute delta SNP indeces greater than 0.1 to filter out potential QTL. An estimation of the mode of the trimmed set is calculated using the mlv function from the package modeest. Finally, the mean and variance of the set are estimated using the median and mode and p-values are estimated from a log normal distribution.

4) Negative Log10- and Benjamini-Hochberg adjusted p-values are calculated using p.adjust