R/susieR.r

In MRCIEU/gwasglue: GWAS summary data sources connected to analytical tools

#' Perform fine mapping pipeline using susieR
#'
#' Clumps data, then maps those to LD representative regions. Within each detected LD representative region, fine mapping is performed
#'
#' @param vcffile Path to vcf file
#' @param bfile Path to ld reference panel
#' @param plink_bin Path to plink
#' @param pop EUR, ASN or AFR
#' @param clump_kb <what param does>
#' @param clump_r2 <what param does>
#' @param clump_p <what param does>
#' @param ... Optional arguments to be passed to susie_rss
#'
#' @export
#' @return List
susieR_pipeline <- function(vcffile, bfile, plink_bin, pop, threads=1, clump_kb=1000, clump_r2=0.001, clump_p=5e-8, ...)
{
	message("Performing clumping")
	clumped <- clump_gwasvcf(vcffile, plink=plink_bin, bfile=bfile)

	message("Map clumps to regions")
	regions <- map_variants_to_regions(clumped$chrpos, pop=pop)
	regions$rsid <- clumped$rsid[match(regions$variant, clumped$chrpos)]

	message("Obtain LD matrices for each region")
	m <- gwasvcf_to_finemapr(region=regions$region, vcf=vcffile, bfile=ldref, plink_bin=plink, threads=threads)

	message("Perform susieR finemapping in each region")
	m2 <- parallel::mclapply(1:length(m), function(i) {
		message(i)
		res <- susieR::susie_rss(
			m[[i]]$z$zscore, 
			m[[i]]$ld, 
			...
		)
		res$fmset <- sapply(m[[i]]$susieR$sets$cs, function(x){
			m[[i]]$z$snp[x[which.max(m[[i]]$susieR$pip[x])]]
		})
		return(res)
	}, mc.cores=threads)

	for(i in 1:length(m))
	{
		m[[i]]$susieR <- m2[[i]]
	}
	out <- list(clumped=clumped, regions=regions, res=m)
	return(out)
}