R/quantifySNPs.R

In stnava/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Simple quantifySNPs function.
#'
#' quantifySNPs converts trinary snps to frequency data
#'
#' @param snps input matrix
#' @param freqthresh remove snps below this frequency
#' @param shiftit shift the snps to smooth the estimate
#' @param replaceWithF replaces snps with frequency values
#' @param traitvecin map snps to trait vector
#' @param trainvec defines training data
#' @return matrix is output
#' @author Avants BB
#' @examples
#'
#' mat <- matrix(c(0, 1, 2, 0, 0, 1, 2, 2, 2), ncol = 3)
#' wmat <- quantifySNPs(mat, freqthresh = 0)
#'
#' @export quantifySNPs
quantifySNPs <- function(
    snps, freqthresh = 0.1,
    shiftit = FALSE, replaceWithF = TRUE, traitvecin = NA, trainvec = NA) {
  if (nargs() == 0) {
    print("Usage:  x_b<-quantifySNPs( x ) ")
    return(1)
  }
  qsnps <- snps
  okrow <- rep(FALSE, ncol(qsnps))
  traitvec <- as.numeric(traitvecin)
  NN <- ncol(qsnps)
  progress <- txtProgressBar(min = 0, max = NN, style = 3)
  for (y in 1:NN) # or dd for ADNI_SNPS
  {
    temp <- snps[, y]
    if (sum(is.na(as.numeric(temp))) == 0) okrow[y] <- TRUE
    if (okrow[y]) {
      t1 <- which(temp == 0)
      t2 <- which(temp == 1)
      t3 <- which(temp == 2)
      f1 <- length(t1) / length(temp)
      f2 <- length(t2) / length(temp)
      f3 <- length(t3) / length(temp)
      if (f1 < freqthresh) okrow[y] <- FALSE
      if (f2 < freqthresh) okrow[y] <- FALSE
      if (f3 < freqthresh) okrow[y] <- FALSE
      if (replaceWithF & is.na(traitvec)) {
        qsnps[t1, y] <- f1
        qsnps[t2, y] <- f2
        qsnps[t3, y] <- f3
      }
      if (replaceWithF & !is.na(traitvec) & is.na(trainvec)) {
        qsnps[t1, y] <- mean(traitvec[t1])
        qsnps[t2, y] <- mean(traitvec[t2])
        qsnps[t3, y] <- mean(traitvec[t3])
      }
      if (replaceWithF & !is.na(traitvec) & !is.na(trainvec)) {
        t1tr <- which(temp == 0 & trainvec)
        t2tr <- which(temp == 1 & trainvec)
        t3tr <- which(temp == 2 & trainvec)
        qsnps[t1, y] <- mean(traitvec[t1tr])
        qsnps[t2, y] <- mean(traitvec[t2tr])
        qsnps[t3, y] <- mean(traitvec[t3tr])
      }
    }
    if (y %% 100 == 0 || y == NN) {
      setTxtProgressBar(progress, y)
    }
  }
  on.exit({
    close(progress)
  })
  qsnps <- qsnps[, okrow]
  if (shiftit & usePkg("magic")) {
    qsnps <- qsnps + magic::ashift(qsnps, c(0, 1)) + magic::ashift(qsnps, c(0, -1))
    qsnps <- qsnps + magic::ashift(qsnps, c(0, 1)) + magic::ashift(qsnps, c(0, -1))
  }
  return(qsnps)
}