R/LD-prune.R

In mattwarkentin/genetools: What the Package Does (One Line, Title Case)

#' @title LD Prune
#'
#' @description Prune a set of variants based on linkage-disequilibrium.
#'
#' @param ... Additional options passed to PLINK.
#' @param type Type of LD pruning.
#' @param window Window width (in kb).
#' @param step Variant count to shift the window.
#' @param threshold R-squared threshold.
#' @param load Logical. Should results be loaded into \code{R}?
#' @param recode Recode text fileset.
#' @param fam not sure yet...
#'
#' @export

prune <- function(...,
                  type = c("indep", "indep-pairwise", "indep-pairphase")[2],
                  window = 200,
                  step = 50,
                  threshold = ifelse(type=="indep", 2, 0.9),
                  load = FALSE,
                  recode = NULL,
                  fam = NULL) {

  options <- list(...)
  if(!is.null(options$cmd)) stop("Don't specify cmd")

  save.out <- options$out
  options.out <- NULL
  options$cmd <- paste0("--", type)
  options$cmd <- paste(options$cmd, window, step, threshold)
  out <- do.call("plink", options)

  ### We should be able to over-write the extract arguments, since only those
  ### SNPs selected are in .prune.in
  options$extract <- paste0(out, ".prune.in")
  options$cmd <- "--make-bed"

  if(!load) {
    options$out <- save.out
    return(do.call("plink", options))
  }
  else {
    out <- do.call("plink", options)
    options$out <- save.out
    options$cmd <- NULL
    options$recode <- recode
    options$fam <- fam
    options$bfile <- out
    options$keep <-
      options$remove <-
      options$extract <-
      options$exclude <-
      options$chr <- NULL
    return(do.call("load.to.R", options))
  }
}