bayes_int: Calculate Bayes credible intervals

Description

Calculate Bayes credible intervals for a single LOD curve on a single chromosome, with the ability to identify intervals for multiple LOD peaks.

Usage

bayes_int(
  scan1_output,
  map,
  chr = NULL,
  lodcolumn = 1,
  threshold = 0,
  peakdrop = Inf,
  prob = 0.95,
  expand2markers = TRUE
)

Value

A matrix with three columns:

ci_lo - lower bound of interval
pos - peak position
ci_hi - upper bound of interval

Each row corresponds to a different peak.

Arguments

scan1_output: An object of class "scan1" as returned by scan1().
map: A list of vectors of marker positions, as produced by insert_pseudomarkers().
chr: Chromosome ID to consider (must be a single value).
lodcolumn: LOD score column to consider (must be a single value).
threshold: Minimum LOD score for a peak.
peakdrop: Amount that the LOD score must drop between peaks, if multiple peaks are to be defined on a chromosome.
prob: Nominal coverage for the interval.
expand2markers: If TRUE, QTL intervals are expanded so that their endpoints are at genetic markers.

Details

We identify a set of peaks defined as local maxima that exceed the specified threshold, with the requirement that the LOD score must have dropped by at least peakdrop below the lowest of any two adjacent peaks.

At a given peak, if there are ties, with multiple positions jointly achieving the maximum LOD score, we take the average of these positions as the location of the peak.

The default is to use threshold=0, peakdrop=Inf, and prob=0.95. We then return results a single peak, no matter the maximum LOD score, and give a 95% Bayes credible interval.

Examples

Run this code

# read data
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,7]

# insert pseudomarkers into map
map <- insert_pseudomarkers(iron$gmap, step=1)

# calculate genotype probabilities
probs <- calc_genoprob(iron, map, error_prob=0.002)

# grab phenotypes and covariates; ensure that covariates have names attribute
pheno <- iron$pheno
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)

# perform genome scan
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)

# 95% Bayes credible interval for QTL on chr 7, first phenotype
bayes_int(out, map, chr=7, lodcolum=1)