R/qtlStats.R

In jtlovell/qtlTools: Data Processing and Plotting in Association with R/qtl

#' @title Extract statistics from a QTL model
#'
#' @description
#' \code{qtlStats} Pushes a QTL model through fitqtl and calculates confidence intervals.
#'
#' @param cross The qtl cross object with marker names that need to be changed.
#' @param pheno.col Character or numeric vector of length 1, indicating the phenotype to be tested
#' @param mod A QTL model, generated by either makeqtl (etc.) or stepwiseqtl.
#' @param form If a generic (not stepwiseqtl) qtl model is supplied and does
#' not contain a formula, this is required. The formula must match the QTL in the
#' fitted QTL (mod) object.
#' @param covar the covariate specified in the model / formula. See fitqtl.
#' @param calcConfint Should confidence intervals be calculated? If so, see ... Also
#' ensure that lodprofiles are included in the model.
#' @param estEffects Should allelic effects be calculated?
#' @param model Should a normal or binary model be fit?
#' @param ... additional arguments passed on to calcCis.
#'
#' @details A simple method to produce stats from a fitted model. Calls several R/qtl
#' functions depending on the user specifications. These include fitqtl, summary.fitqtl,
#' bayesint, dropint, and qtlTools::calcCis. The method of QTL inference is extracted from the QTL model.
#'
#' @examples
#' library(qtlTools)
#' data(fake.bc)
#' cross<-fake.bc
#' cross <- calc.genoprob(cross, step=2.5)
#' mod <- makeqtl(cross, chr = c(2,5), pos = c(40,25), what = "prob")
#' nform <- "y ~ Q1 + Q2 + Q1*sex + sex"
#' sex <- data.frame(sex = cross$phe$sex)
#' mod <- refineqtl(cross, mod, pheno.col = "pheno1", qtl = mod,
#'    formula = nform, covar = sex, method="hk")
#' qtlStats(cross, pheno.col = "pheno1",form = nform, mod = mod, covar=sex)
#'
#' data(fake.f2)
#' cross<-fake.f2
#' cross <- calc.genoprob(cross, step=2.5)
#' mod <- makeqtl(cross, chr = c(1,13), pos = c(27.5,27.5), what = "prob")
#' nform <- "y ~ Q1 + Q2 + Q1*sex + sex"
#' sex <- data.frame(sex = cross$phe$sex)
#' mod <- refineqtl(cross, mod, pheno.col = 1, qtl = mod,
#'    formula = nform, covar = sex, method="hk")
#' qtlStats(cross, pheno.col = 1,form = nform, mod = mod, covar=sex)
#'
#' @return A dataframe of statistics.
#'
#' @import qtl
#' @export

qtlStats<-function(cross,
                   mod,
                   pheno.col,
                   form=NULL,
                   covar=NULL,
                   calcConfint=TRUE,
                   estEffects = TRUE,
                   model="normal",
                   ...){

  if(nqtl(mod)<1) stop("supplied model must contain at least 1 QTL\n")

  if(length(pheno.col)>1){
    cat(">1 Phenotype supplied, using the first \n")
    pheno.col<-pheno.col[1]
  }

  if(is.null(form) & is.null(attr(mod, "formula")))
    stop("formula must either be supplied, or included in mod (QTL model)\n")

  if(is.null(form))
    form <- formula(mod)

  if("prob" %in% names(mod)){
    method.qtl = "hk"
  }else{
    method.qtl = "imp"
  }

  # 1. parse the formula
  form1<-form
  form<-as.formula(form)
  terms<-attr(terms(form), "term.labels")
  out<-data.frame(qtlnames = mod$name, altnames = mod$altname, stringsAsFactors=F)
  terms2<-terms
  for(i in 1:nrow(out)){
    terms2<-gsub(out$altnames[i], out$qtlnames[i],terms2, fixed=T)
  }
  form.out<-data.frame(altnames = terms, terms = terms2, stringsAsFactors=F)
  out<-merge(form.out,out , by = "altnames", all=T)
  out$phenotype = pheno.col
  out$formula <- form1
  covar.name<-colnames(covar)

  # 2. Get the dropone stats and merge
  f<-summary(fitqtl(cross, qtl = mod, formula = form, pheno.col = pheno.col,
                    covar = covar, method = method.qtl, model=model,
            dropone = TRUE, get.ests = estEffects))
  if("result.drop" %in% names(f)){
    out.drop<-data.frame(f$result.drop)
    out.drop$terms = rownames(out.drop)
  }else{
    out.drop<-with(data.frame(f$result.full),
                   data.frame(df = df[1], Type.III.SS = SS[1], LOD = LOD[1],
                              X.var = X.var[1], F.value = NA,
                              Pvalue.Chi2.=Pvalue.Chi2.[1],
                              Pvalue.F. = Pvalue.F.[1], terms = terms2))
  }

  out<-merge(out, out.drop, by = "terms")
  out$modelLod<-f$lod
  out$modelPercVar<-f$result.full[1,5]

  # 3. Get the estimates and merge
  if(estEffects){
    out.est<-data.frame(f$ests)
    out.est<-out.est[-which(rownames(out.est)=="Intercept"),]
    out.est$terms = rownames(out.est)
    if(nrow(out.est) == nrow(out)){
      out<-merge(out, out.est, by = "terms")
    }else{
      qtype<-out.est$terms
      for(i in mod$name) qtype<-gsub(i,"q",qtype, fixed=T)
      qtype<-gsub(paste0(":",covar.name),"",qtype, fixed = T)
      qtype<-gsub(paste0(covar.name,":"),"",qtype, fixed = T)
      oed<-out.est[qtype == "qd",]
      oed$terms<-rownames(oed)
      names(oed)[1:3]<-paste(names(oed)[1:3],"dom",sep="_")
      oea<-out.est[qtype == "qa" | out.est$terms %in% out$terms,]
      oea$terms<-rownames(oea)
      for(i in out$terms){
        oea$terms<-gsub(paste0(i,"a"),i,oea$terms)
        oed$terms<-gsub(paste0(i,"d"),i,oed$terms)
      }
      out<-merge(out, oea, by = "terms", all.x=T)
      out<-merge(out, oed, by = "terms", all.x=T)
    }
  }

  # 4. Get Cis and merge
  if(calcConfint){
    cis<-calcCis(cross=cross,
                 mod=mod,
                 qtlnames = NULL, ...)
    cis<-cis[-which(colnames(cis)=="pheno")]
    colnames(cis)[1]<-"terms"
    out<-merge(out, cis, by = "terms", all=T)
  }

  return(out)
}