scan1coef: Calculate QTL effects in scan along one chromosome

Description

Calculate QTL effects in scan along one chromosome with a single-QTL model using Haley-Knott regression or a linear mixed model (the latter to account for a residual polygenic effect), with possible allowance for covariates.

Usage

scan1coef(
  genoprobs,
  pheno,
  kinship = NULL,
  addcovar = NULL,
  nullcovar = NULL,
  intcovar = NULL,
  weights = NULL,
  contrasts = NULL,
  model = c("normal", "binary"),
  zerosum = TRUE,
  se = FALSE,
  hsq = NULL,
  reml = TRUE,
  ...
)

Value

An object of class "scan1coef": a matrix of estimated regression coefficients, of dimension positions x number of effects. The number of effects is n_genotypes + n_addcovar + (n_genotypes-1)*n_intcovar. May also contain the following attributes:

SE - Present if se=TRUE: a matrix of estimated standard errors, of same dimension as coef.
sample_size - Vector of sample sizes used for each phenotype

Arguments

genoprobs: Genotype probabilities as calculated by calc_genoprob().
pheno: A numeric vector of phenotype values (just one phenotype, not a matrix of them)
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.
nullcovar: An optional numeric matrix of additional additive covariates that are used under the null hypothesis (of no QTL) but not under the alternative (with a QTL). This is needed for the X chromosome, where we might need sex as a additive covariate under the null hypothesis, but we wouldn't want to include it under the alternative as it would be collinear with the QTL effects. Only used if kinship is provided but hsq is not, to get estimate of residual heritability.
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.
contrasts: An optional numeric matrix of genotype contrasts, size genotypes x genotypes. For an intercross, you might use cbind(mu=c(1,1,1), a=c(-1, 0, 1), d=c(0, 1, 0)) to get mean, additive effect, and dominance effect. The default is the identity matrix.
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]$ .
zerosum: If TRUE, force the genotype or allele coefficients sum to 0 by subtracting their mean and add another column with the mean. Ignored if contrasts is provided.
se: If TRUE, also calculate the standard errors.
hsq: (Optional) residual heritability; used only if kinship provided.
reml: If kinship provided: if reml=TRUE, use REML; otherwise maximum likelihood.
...: Additional control parameters; see Details;

Details

For each of the inputs, the row names are used as individual identifiers, to align individuals.

If kinship is absent, Haley-Knott regression is performed. If kinship is provided, a linear mixed model is used, with a polygenic effect estimated under the null hypothesis of no (major) QTL, and then taken as fixed as known in the genome scan.

If contrasts is provided, the genotype probability matrix, $P$ , is post-multiplied by the contrasts matrix, $A$ , prior to fitting the model. So we use $P \cdot A$ as the $X$ matrix in the model. One might view the rows of A^-1 as the set of contrasts, as the estimated effects are the estimated genotype effects pre-multiplied by A^-1.

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

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.

Examples

Run this code

# read data
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,c(7,19)] # reduce to chr 7 and 19

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

# calculate genotype probabilities
probs <- calc_genoprob(iron, map, error_prob=0.002)
probs[["7"]] <- probs[["7"]][,,1:5] # reduce to very small number

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

# calculate coefficients for chromosome 7
coef <- scan1coef(probs[,"7"], pheno, addcovar=covar)

# leave-one-chromosome-out kinship matrix for chr 7
kinship7 <- calc_kinship(probs, "loco")[["7"]]

# calculate coefficients for chromosome 7, adjusting for residual polygenic effect
coef_pg <- scan1coef(probs[,"7"], pheno, kinship7, addcovar=covar)

Run the code above in your browser using DataLab