R/polygenicQTL.R

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

#' @title Run rrBLUP polygenic inheritance on an R/qtl model
#'
#' @description
#' \code{polygenicQTL} Parses and pushes qtl model data int rrBLUP's mixed model
#' genomic selection pipe. Outputs polygenic inheritance stats.
#'
#' @param cross A qtl cross object containing conditional genotype probabilites calculated with
#' stepwidth = "fixed".
#' @param pheno.col Character or numeric vector of length 1,
#' indicating the phenotype to be tested.
#' @param qtl A QTL model, generated by  makeqtl (etc.) .
#' @param formula The formula for the qtl, names of columns in the covariate dataframe
#' and altnames of the QTL need to match exactly, or will be ignored.
#' @param covar the covariate specified in the model / formula. See fitqtl.
#' @param rrBLUP.method should reml or ml be run?
#' @param verbose Should a normal or binary model be fit?
#' @param ... additional arguments passed on to rrBLUP::mixed.solve
#'
#' @details See rrBLUP::mixed.solve and rrBLUP::kin.blup
#'
#' @return A dataframe with 3 columns:
#' \itemize{
#'  \item{"Vg"}{The genetic variance attributable to the kinship (polygenic)}
#'  \item{"Ve"}{The amount of environmental (residual) varaince}
#'  \item{"polyH2"}{the proportion of Vtotal in Vg. Equivalent to the amount of
#'  heritability attributable to kinship (polygenic) variance following incoporpation of
#'  any fixed effects}
#' }
#' @examples
#' \dontrun{
#' library(qtlTools)
#' library(rrBLUP)
#' data(fake.f2)
#' cross<-fake.f2
#' cross<-calc.genoprob(cross, step = 12, stepwidth = "fixed")
#' qtl<-makeqtl(cross, chr = 13, pos = 28, what = "prob")
#' formula = " ~ Q1 + Q1*sex + sex"
#' covar = data.frame(sex = pull.pheno(cross, pheno.col = "sex"))
#' polygenicQTL(cross = cross, pheno.col = 1)
#' polygenicQTL(cross = cross, covar = covar, formula = "~ sex", pheno.col = 1)
#' polygenicQTL(cross = cross, qtl = qtl, formula = "~ Q1", pheno.col = 1)
#' polygenicQTL(cross = cross, qtl = qtl,covar = covar, formula = "~ Q1 + Q1*sex + sex", pheno.col = 1)
#' }
#'
#' @import qtl
#' @export
polygenicQTL<-function(cross,
                       formula = NULL, qtl = NULL,
                       pheno.col = 1, covar=NULL,
                       rrBLUP.method = "REML",
                       verbose = TRUE,
                       ...){

  y <- pull.pheno(cross, pheno.col = pheno.col)
  if(!requireNamespace("rrBLUP", quietly = TRUE)){
    stop("install the rrBLUP package to run polygenicQTL\n")
  }else{
    requireNamespace("rrBLUP", quietly = TRUE)
  }

  if("prob" %in% names(cross$geno[[1]])){
    atr<-attributes(cross$geno[[1]]$prob)
    genotypes<-attr(cross$geno[[1]]$prob,"dimnames")[[3]]
  }else{
    if("draws" %in% names(cross$geno[[1]])){
      atr<-attributes(cross$geno[[1]]$draws)
      tmp<-calc.genoprob(cross)
      genotypes<-attr(tmp$geno[[1]]$prob,"dimnames")[[3]]
    }else{
      stop("run either calc.genoprob or sim.geno first.\n")
    }
  }

  if(atr$step == 0) atr$step = 1
  ag<-argmax.geno(cross, step = atr$step, error.prob = atr$error.prob,
                  off.end = atr$off.end, map.function = atr$map.function,
                  stepwidth = "fixed")
  gp<-pull.argmaxgeno(ag, include.pos.info=F)

  gpn<-gp
  for(i in 1:length(genotypes)) gpn[gpn == i]<-genotypes[i]
  if(!is.null(covar)){
    gpn<-cbind(covar, gp)
  }

  genos<-unique(as.numeric(gp))[order(unique(as.numeric(gp)))]
  ngeno<-length(genos)
  if(ngeno == 2){
    gp[gp == genos[1]]<- (-1)
    gp[gp == genos[2]]<- 1
  }else{
    if(ngeno == 3){
      gp[gp == genos[1]]<- (-1)
      gp[gp == genos[2]]<- 0
      gp[gp == genos[3]]<- 1
    }else{
      if(ngeno == 4){
        gp1<-gp
        gp2<-gp
        gp1[gp1 == genos[1]]<- (-1)
        gp1[gp1 == genos[2]]<- 1
        gp2[gp2 == genos[3]]<- (-1)
        gp2[gp2 == genos[4]]<- 1
        gp<-cbind(gp1,gp2)
      }else{
        stop(cat("don't know how to deal with cross type",class(cross)[1]))
      }
    }
  }


  if(verbose) cat("calculating kinship matrix \n")
  kinship.matrix<-rrBLUP::A.mat(gp)

  if(verbose) cat("running polygenic selection \n")
  if(is.null(covar) & is.null(qtl)){
    test<-rrBLUP::mixed.solve(y=y,
                              K = kinship.matrix,
                              method = rrBLUP.method, ...)
  }else{
    if(is.null(formula)){
      warning("no formula specified, assuming an additive QTL model")
      if(is.null(covar)){
        formula <- paste0("~ ", paste0(qtl$altname, collapse = " + "))
      }else{
        if(is.null(qtl)){
          formula <- paste0("~ ", paste0(colnames(covar), collapse = " + "))
        }else{
          formula <- paste0("~ ", paste0(qtl$altname, collapse = " + "), " + ",
                            paste0(colnames(covar), collapse = " + "))
        }
      }
    }
    if(!is.null(qtl)){
      mars<-find.marker(cross, chr = qtl$chr, pos = qtl$pos)
      for(i in 1:nqtl(qtl)) formula<-gsub(qtl$altname[i],mars[i],formula)
    }
    fixed.matrix<-model.matrix(as.formula(formula), data = data.frame(gpn))
    test<-rrBLUP::mixed.solve(y=y,
                              K = kinship.matrix,
                              X = fixed.matrix,
                              method = rrBLUP.method, ...)
  }
  return(data.frame(Vg = test$Vu,
                    Ve = test$Ve,
                    polyH2 = test$Vu/(test$Vu+test$Ve)))
}