Map marginal or response QTLs for feature-SNP pairs

Source: R/core.R

This function takes as input individual phenotype and genotype data and tests for marginal and response QTLs.

map_qtl(
  candidate,
  fdr = NULL,
  fn,
  geno,
  pheno,
  anno,
  covar = NULL,
  scale = FALSE,
  rint = FALSE,
  ranef = FALSE,
  kinship = NULL,
  filter.geno = TRUE,
  excluded = FALSE,
  type = "response",
  seed = 1,
  epsilon = 0.01,
  control = NULL,
  gls = FALSE
)

Arguments

candidate

A data frame containing feature IDs and SNP IDs. The first column must contain feature IDs and the second column must contain SNP IDs.

fdr

A scalar specifying a false discovery rate (FDR) threshold.

fn

A character string specifying the function. This must be one of "nonlinear" and "linear", corresponding to nonlinear and linear models, respectively.

geno

A data frame formatted as input for MatrixEQTL. The first column must contain SNP IDs. The rest of the columns must contain genotypes, which can be dosage or coded as (0, 1, 2). Each column corresponds to each sample.

pheno

A data frame formatted as input for MatrixEQTL. The first column must contain feature IDs (e.g., genes and ATAC-peak regions). The rest of the columns must contain pre-processed molecular phenotypes. Each column corresponds to each sample.

anno

A data frame containing sample IDs, subject IDs, and treatment conditions. The columns must be named as "sample", "subject", and "condition".

covar

A data frame containing fixed effect covariates. The first column must contain sample IDs. The rest of the columns must contain covariate values. Each row corresponds to each sample. Each column corresponds to each covariate.

scale

A logical scalar as to whether to standardize the phenotype, y. This can be set to TRUE only if fn is set to "linear".

rint

A logical scalar as to whether the phenotypes need to be RINT-transformed. This cannot be set to TRUE if fn is set to "nonlinear".

ranef

A logical scalar as to whether to include random effect.

kinship

A matrix containing pairwise genetic relatedness between subjects. The row and column names must match the set of unique elements of the "subject" column in "anno" in the corresponding order (i.e., unique(anno$subject)). This is set to NULL by default, in which case the identity matrix is used.

filter.geno

A logical scalar as to whether to ensure that homozygotes for both alleles have at least one observation in both control and treated conditions.

excluded

A logical scalar as to whether to return a data frame containing feature-SNP pairs that have been excluded from the analysis. If TRUE, the output is a list of data frames.

type

A character string specifying the type of QTLs. This must be one of "marginal" and "response".

seed

An integer specifying a seed for RNG.

epsilon

A small positive scalar. The upper bound of the \(h^2\) search interval is set to \(1 - \epsilon\). Default is 1e-2. Ignored when ranef is FALSE.

control

A list of control parameters passed to optim for the inner BFGS step. If NULL, the defaults are used. Ignored when gls = TRUE.

gls

A logical scalar. If TRUE, use the closed-form generalized least squares (GLS) estimator instead of BFGS. Only valid when fn is "linear". Default is FALSE.

Value

A data frame containing some of the following columns:

  • feat_id - A character string specifying the feature ID.

  • snp_id - A character string specifying the SNP ID.

  • beta0 - A scalar value representing the intercept estimate.

  • beta_g - A scalar value representing the genotype effect estimate based on the marginal model (when type = "marginal") or the interaction model (when type = "response").

  • beta_t - A scalar value representing the treatment effect estimate based on the interaction model.

  • beta_gxt - A scalar value representing the interaction effect estimate based on the interaction model.

  • sigma - A scalar value representing the estimate of the residual error standard deviation.

  • sigma_u - A scalar value representing the estimate of the random effect standard deviation.

  • p_value - A scalar value representing the nominal p-value from LRT.

  • return_code - An integer: 0 if all optimizations converged, NA otherwise.

Examples

anno <- data.frame(
    sample = paste0("s", 1:20),
    subject = rep(paste0("d", 1:10), each = 2),
    condition = rep(c(0, 1), 10)
)
dat <- make_sim_data(
    num = c(0, 0, 0, 0, 0, 0, 0, 1),
    anno = anno, fn = "nonlinear",
    coef = c(0.5, 0.3, 0.2)
)
result <- map_qtl(
    candidate = dat$candidate[1, ],
    fn = "nonlinear",
    geno = dat$geno,
    pheno = dat$pheno,
    anno = anno
)