R/calcCis.R

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

#' @title Calculate all confidence intervals from a QTL model
#'
#' @description
#' \code{calcCis} Takes a multiple QTL model or scanone and
#' permutations output and returns a dataframe of confidence intervals.
#'
#' @param cross The QTL cross object.
#' @param mod A QTL model containing LOD profiles. Can be generated by
#' stepwiseqtl (with keepLodProfile = TRUE), or by refineqtl.
#' @param s1.output Output from a scanone call. Not needed if qtl models are supplied to mod.
#' If mod is not supplied, it is necessary to provide both s1.output and perm.output.
#' @param perm.output Output from a scanone call with n.perm > 1. Needs to have the identical
#' parameters as s1.output.
#' @param qtlnames An optional vector of user-specified QTL identifiers. Used only if mod is
#' supplied.
#' @param pheno.col If passing a QTL model, specify the name of the phenotype.
#' @param lodint Should lodint (default) or bayesint be used to estimate confidence intervals?
#' @param drop If lodint = TRUE, what drop threshold should be used?
#' @param prob If lodint = FALSE, what significance should be used for bayesint?
#' @param expandtomarkers Logical. Should the intervals be expanded to the nearest real marker?
#' @param ... Additional arguments passed on to pullSigQTL, such as the alpha threshold for a qtl.
#' @details This function iterates through the qtl in the model, returning confidence intervals.
#' @return A dataframe of confidence intervals, including the maximum LOD score,
#' bounding markers and positions.
#'
#' @examples
#' library(qtlTools)
#' data(fake.bc)
#' cross<-fake.bc
#' cross <- calc.genoprob(cross, step=2.5)
#'
#' # Make a QTL object and formula
#' mod <- makeqtl(cross, chr = c(2,5), pos = c(40,25), what = "prob")
#' sex <- data.frame(sex = cross$phe$sex)
#' nform <- "y ~ Q1 + Q2 + Q1*sex + sex"
#'
#' #Calculate lodprofiles for confidence interval estimation
#' mod <- refineqtl(cross, mod, pheno.col = "pheno1",
#'   qtl = mod, formula = nform, covar = sex, method="hk")
#' calcCis(cross = cross, mod=mod)
#' \dontrun{
#' s1<-scanone(cross, method="hk", pheno.col=c("pheno1", "pheno2"))
#' perm<-scanone(cross, n.perm=100, method="hk",pheno.col=c("pheno1", "pheno2"), verbose=FALSE)
#' calcCis(cross,s1.output=s1, perm.output=perm)
#' }
#' @import qtl
#' @export

calcCis<-function(cross,
                  mod = NULL,
                  s1.output = NULL,
                  pheno.col = NA,
                  perm.output = NULL,
                  qtlnames = NULL,
                  lodint = TRUE,
                  drop=1.5,
                  prob=0.95,
                  expandtomarkers = FALSE,
                  chr = NULL,
                  ...){
  if(is.null(mod) & is.null(s1.output))
    stop("either scanone peaks, or qtl models must be supplied\n")
  if(!is.null(s1.output) & is.null(perm.output))
    stop("if scanones are provided, permutation results must also be included\n")
  if(!is.null(mod)){
    if(is.null(attr(mod, "lodprofile")))
      stop("LOD profiles must be included in the QTL model\n")
    if(is.null(qtlnames))  qtlnames <- mod$name
    if(is.null(chr)){
      chr = chrnames(cross)
    }
    out<-lapply(1:nqtl(mod), function(j){
      if(lodint){
        ciout<-lodint(mod,
                      chr,
                      qtl.index=j,
                      drop=drop,
                      expandtomarkers=expandtomarkers)
      }else{
        ciout<-bayesint(mod,
                        chr,
                        qtl.index=j,
                        prob=prob,
                        expandtomarkers=expandtomarkers)
      }
      lowmarker<-rownames(ciout)[1]
      highmarker<-rownames(ciout)[3]
      lowposition<-as.numeric(ciout[1,2])
      highposition<-as.numeric(ciout[3,2])
      maxLod<-as.numeric(sapply(attr(mod, "lodprofile"), function(x)
        max(x$lod))[[j]])

      chr=mod$chr[j]
      pos=mod$pos[j]
      return(data.frame(pheno = pheno.col,
                        qtlname = qtlnames[j],
                        chr, pos, maxLod,
                        lowmarker,highmarker,
                        lowposition,highposition))
    })
  }else{
    phes<-colnames(s1.output)[-c(1:2)]
    qtl.peaks <- pullSigQTL(cross, pheno.col=phes,
                            s1.output = s1.output,
                            perm.output = perm.output,
                            returnQTLModel = FALSE, ...)
    out<-lapply(1:nrow(qtl.peaks), function(j){
      dat<-qtl.peaks[j,]
      phe<-dat$pheno
      lodcolumn <- which(names(s1.output) == phe)-2
      chr<-dat$chr
      pos<-dat$pos
      if(lodint){
        ciout<-lodint(s1.output,
                      chr=chr,
                      lodcolumn=lodcolumn,
                      drop=drop,
                      expandtomarkers=expandtomarkers)
      }else{
        ciout<-bayesint(s1.output,
                        chr=chr,
                        lodcolumn=lodcolumn,
                        prob=prob,
                        expandtomarkers=expandtomarkers)
      }
      lowmarker<-rownames(ciout)[1]
      highmarker<-rownames(ciout)[3]
      lowposition<-as.numeric(ciout[1,2])
      highposition<-as.numeric(ciout[3,2])
      maxLod<-dat$lod1

      return(data.frame(pheno = phe,
                        qtlname = NA,
                        chr, pos, maxLod,
                        lowmarker,highmarker,
                        lowposition,highposition))
    })
  }
  return(data.frame(do.call(rbind, out)))
}