R/top_snps_pval_ld.R

In jean997/sherlockAsh: CAUSE: Causal Analysis Using Summary Effect Estimates

#'@title Rank SNPs by pval, prune for LD
#'@description Selects a set of top SNP based on p-value that are in minimal LD
#'@param ld list of LD matrices
#'@param pvals Vector of pvalues (optional)
#'@param pval_thresh Maximum pvalue
#'@param max_snps Maximum number of SNPs to retain
#'@param ld_thresh Maximum pairwise correlation between SNPs
#'@export
top_snps_pval_ld <- function(ld, pvals, pval_thresh = Inf, max_snps = Inf, ld_thresh=0.1){
  if(missing(ld)) stop("List of LD matrices is required\n")

  lens <- sapply(ld, function(x){nrow(x)})
  blocks <- rep(1:length(ld), lens)
  p <- length(blocks)

  if(is.finite(pval_thresh) & missing(pvals)){
    stop("pval_thresh is given but pvals is missing.\n")
  }

  if(missing(pvals)){
    o <- sample(1:p, size=p, replace=FALSE)
  }else{
    stopifnot(length(pvals)==p)
    o <- order(pvals, decreasing=F)
    max_ix <- sum(pvals < pval_thresh)
    o <- o[1:max_ix]
  }
  keep <- c()
  while(length(o) > 0 & length(keep) < max_snps){
    keep <- c(keep, o[1])
    ix <- o[1]
    block <- blocks[o[1]]
    blockix <- which(blocks==block)
    myix <- which(blockix == o[1])
    remove_ix <- blockix[ which(ld[[block]][,myix] > ld_thresh)]
    o <- o[!o %in% remove_ix]
  }
  return(keep)
}