R/qc.R

In chr1swallace/GUESSFM: Using GUESS for Fine Mapping

x2xw <- function(x) {
  w <- x
  x <- as.numeric(names(x))
  return(list(x=x,w=w)) 
}
namedlapply <- function(x,fun) {
  tmp <- lapply(x,fun)
  names(tmp) <- names(x)
  return(tmp)
}
##' Summarise skewness of a posterior distribution on the number of SNPs in a model
##'
##' @title skewness
##' @param x pp.nsnp object
##' @return list of numerics
##' @author Chris Wallace
skewness <- function(x) {
  if(is.list(x))
    return(namedlapply(x,skewness))
  x <- x2xw(x[-1])
  m1 <- sum(x$x * x$w) / sum(x$w)
  m2 <- sum(x$x^2 * x$w) / sum(x$w)
  m3 <- sum(x$x^3 * x$w) / sum(x$w)
  sigma <- sqrt(m2 - m1^2)
  (m3 - 3*m1*sigma^2-m1^3)/(sigma^3)
}
##' Summarise number of modes of a posterior distribution on the number of SNPs in a model
##'
##' @title nmodes
##' @param x pp.nsnp object
##' @return list of numerics
##' @author Chris Wallace
nmodes <- function(x) {
  if(is.list(x))
    return(namedlapply(x,nmodes))
  x <- x2xw(x)
  o <- order(x$w, decreasing=TRUE)
  x$x <- x$x[o]
  mdiff <- sapply(seq_along(x$x)[-1], function(i) min(abs(x$x[i] - x$x[1:(i-1)])))
  sum(mdiff>1)+1
}
##' internal function: summarise snps from highly correlated models.
##'
##' @title snps.from.correlated.models
##' @param ret data.frame generated by qc() on a snpmod
##' @param thr r2 threshold considered unlikely
##' @param nshow number of most frequently occurring SNPs to show
##' @return the nshow top SNPs, if any
##' @author Chris Wallace
snps.from.correlated.models <- function(ret, thr=0.8, nshow=10) {
  bads <- which(ret$maxr2>thr)
  if(!length(bads))
    return(NULL)
  ss <- strsplit(ret$model[bads],"%")
  tt <- table(unlist(ss))/length(ss)
  nbad <- length(bads)
  ntot <- nrow(ret)
  PPbad <- sum(ret$PP[bads])
  message("models containing highly correlated SNPs found: ",nbad,"/",ntot)
  message("accounting for a total PP: ",PPbad)
  message("SNPs found in these models:")
  print(head(sort(tt,decreasing=TRUE),nshow))
  invisible(list(summary=c("nmodels.corr.snps"=nbad,"nmodels.total"=ntot,"sumPP.models.corr.snps"=PPbad),
                 SNPs=head(sort(tt,decreasing=TRUE),nshow)))
}