qtl2 (version 0.46)

scan1perm: Permutation test for genome scan with a single-QTL model

Description

Permutation test for a enome scan with a single-QTL model by Haley-Knott regression or a linear mixed model, with possible allowance for covariates.

Usage

scan1perm(
  genoprobs,
  pheno,
  kinship = NULL,
  addcovar = NULL,
  Xcovar = NULL,
  intcovar = NULL,
  weights = NULL,
  reml = TRUE,
  model = c("normal", "binary"),
  n_perm = 1,
  perm_Xsp = FALSE,
  perm_strata = NULL,
  chr_lengths = NULL,
  cores = 1,
  scan_func = NULL,
  ...
)

Value

If perm_Xsp=FALSE, the result is matrix of genome-wide maximum LOD scores, permutation replicates x phenotypes. If perm_Xsp=TRUE, the result is a list of two matrices, one for the autosomes and one for the X chromosome. The object is given class "scan1perm".

Arguments

genoprobs

Genotype probabilities as calculated by calc_genoprob().

pheno

A numeric matrix of phenotypes, individuals x phenotypes.

kinship

Optional kinship matrix, or a list of kinship matrices (one per chromosome), in order to use the LOCO (leave one chromosome out) method.

addcovar

An optional numeric matrix of additive covariates.

Xcovar

An optional numeric matrix with additional additive covariates used for null hypothesis when scanning the X chromosome.

intcovar

An optional numeric matrix of interactive covariates.

weights

An optional numeric vector of positive weights for the individuals. As with the other inputs, it must have names for individual identifiers.

reml

If kinship provided: if reml=TRUE, use REML; otherwise maximum likelihood.

model

Indicates whether to use a normal model (least squares) or binary model (logistic regression) for the phenotype. If model="binary", the phenotypes must have values in \([0, 1]\).

n_perm

Number of permutation replicates.

perm_Xsp

If TRUE, do separate permutations for the autosomes and the X chromosome.

perm_strata

Vector of strata, for a stratified permutation test. Should be named in the same way as the rows of pheno. The unique values define the strata.

chr_lengths

Lengths of the chromosomes; needed only if perm_Xsp=TRUE. See chr_lengths().

cores

Number of CPU cores to use, for parallel calculations. (If 0, use parallel::detectCores().) Alternatively, this can be links to a set of cluster sockets, as produced by parallel::makeCluster().

scan_func

If provided, this function is used for the genome scans. It must take arguments genoprobs, pheno, kinship, addcovar, Xcovar, intcovar, weights, and possibly further arguments.

...

Additional control parameters; see Details.

Details

If kinship is not provided, so that analysis proceeds by Haley-Knott regression, we permute the rows of the phenotype data; the same permutations are also applied to the rows of the covariates (addcovar, Xcovar, and intcovar) are permuted.

If kinship is provided, we instead permute the rows of the genotype data and fit an LMM with the same residual heritability (estimated under the null hypothesis of no QTL).

If Xcovar is provided and perm_strata=NULL, we do a stratified permutation test with the strata defined by the rows of Xcovar. If a simple permutation test is desired, provide perm_strata that is a vector containing a single repeated value.

The ... argument can contain several additional control parameters; suspended for simplicity (or confusion, depending on your point of view). tol is used as a tolerance value for linear regression by QR decomposition (in determining whether columns are linearly dependent on others and should be omitted); default 1e-12. maxit is the maximum number of iterations for converence of the iterative algorithm used when model=binary. bintol is used as a tolerance for converence for the iterative algorithm used when model=binary. eta_max is the maximum value for the "linear predictor" in the case model="binary" (a bit of a technicality to avoid fitted values exactly at 0 or 1).

References

Churchill GA, Doerge RW (1994) Empirical threshold values for quantitative trait mapping. Genetics 138:963--971.

Manichaikul A, Palmer AA, Sen S, Broman KW (2007) Significance thresholds for quantitative trait locus mapping under selective genotyping. Genetics 177:1963--1966.

Haley CS, Knott SA (1992) A simple regression method for mapping quantitative trait loci in line crosses using flanking markers. Heredity 69:315--324.

Kang HM, Zaitlen NA, Wade CM, Kirby A, Heckerman D, Daly MJ, Eskin E (2008) Efficient control of population structure in model organism association mapping. Genetics 178:1709--1723.

See Also

scan1(), chr_lengths(), mat2strata(), scan1gen examples

Examples

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

# 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)

# strata for permutations
perm_strata <- mat2strata(Xcovar)

# permutations with genome scan (just 3 replicates, for illustration)
operm <- scan1perm(probs, pheno, addcovar=covar, Xcovar=Xcovar,
                   n_perm=3, perm_strata=perm_strata)
summary(operm)

# leave-one-chromosome-out kinship matrices
kinship <- calc_kinship(probs, "loco")

# permutations of genome scan with a linear mixed model
operm_lmm <- scan1perm(probs, pheno, kinship, covar, Xcovar, n_perm=3,
                       perm_Xsp=TRUE, perm_strata=perm_strata,
                       chr_lengths=chr_lengths(map))
summary(operm_lmm)


# permutations of scan1gen, fitting glm with probit link
# fitting function
ll_glm <-
   function(pr, pheno, addcovar=NULL, ...)
{
   formula <- ifelse(is.null(pr), "pheno ~ 1", "pheno ~ pr")
   if(!is.null(addcovar)) formula <- paste(formula, "+ addcovar")

    glm_out <- glm(as.formula(formula), family=binomial(link=probit))
    -glm_out$deviance/(2*log(10)) # log10 likelihood
}

# create binary trait
bin_pheno <- setNames(as.numeric(iron$pheno[,1] > median(iron$pheno[,1])),
                      rownames(iron$pheno))

# permutations with glm
operm_glm <- scan1perm(probs, bin_pheno, kinship, covar, Xcovar, n_perm=3,
                       perm_Xsp=TRUE, perm_strata=perm_strata,
                       chr_lengths=chr_lengths(map),
                       scan_func=scan1gen, func=ll_glm)
summary(operm_glm)

if (FALSE) {
# example of scan1perm with scan1snps
file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/main/DOex/DOex.zip")
DOex <- read_cross2(file)
DOex <- DOex[,c("2", "3")] # subset to chr 2 and 3
probs <- calc_genoprob(DOex, error_prob=0.002)

snpdb_file <- system.file("extdata", "cc_variants_small.sqlite", package="qtl2")
queryf <- create_variant_query_func(snpdb_file)

operm <- scan1perm(genoprobs=probs, map=DOex$pmap, pheno=DOex$pheno,
                   scan_func=scan1snps, query_func=queryf, n_perm=3)
}