proportion of variance due to rare variants

rare_variant_prop.R

calc.proprare <- function(nsnp, shape1, eff.func, nid)
{
	# Generate allele frequency distribution
	mafs <- rbeta(nsnp, shape1, 1) / 2
	# Limit min maf to be function of hypothetical sample size
	mafs <- pmax(mafs, 1/(nid*2))
	# Create distribution of effect sizes
	eff <- eff.func(mafs)
	# Calculate variance explained by each SNP
	vexp <- 2*eff^2*mafs*(1-mafs)
	# Return proportion of variance explained due to SNPs with MAF < 0.01
	return(sum(vexp[mafs < 0.01]) / sum(vexp))
}

# 
eff.func1 <- function(x) rep(1, length(x))
eff.func2 <- function(x) abs(log(x))
eff.func3 <- function(x) abs(log(x)^2)


params <- expand.grid(
	nsnp=100000,
	nid=10000, 
	shape1=c(0.3, 0.1, 0.05), 
	eff.func=c("eff.func1", "eff.func2", "eff.func3"), 
	proprare=NA
)

for(i in 1:nrow(params))
{
	cat(i, "\n")
	params$proprare[i] <- calc.proprare(params$nsnp[i], params$shape1[i], get(as.character(params$eff.func[i])), params$nid[i])
}
params



shape1 <- c(0.3, 0.1, 0.05)
nid <- 10000
nsnp <- 100000

par(mfrow=c(3,3))
hist(rbeta(nsnp, 0.3, 1) / 2, breaks=100)
hist(rbeta(nsnp, 0.1, 1) / 2, breaks=100)
hist(rbeta(nsnp, 0.05, 1) / 2, breaks=100)

mafs <- pmax(rbeta(nsnp, shape1[1], 1) / 2, 1/(nid*2))
eff <- eff.func1(mafs)
plot(x=mafs, y=eff)

mafs <- pmax(rbeta(nsnp, shape1[1], 1) / 2, 1/(nid*2))
eff <- eff.func2(mafs)
plot(x=mafs, y=eff)

mafs <- pmax(rbeta(nsnp, shape1[1], 1) / 2, 1/(nid*2))
eff <- eff.func3(mafs)
plot(x=mafs, y=eff)

mafs <- pmax(rbeta(nsnp, shape1[1], 1) / 2, 1/(nid*2))
eff <- eff.func1(mafs)
vexp <- 2*eff^2*mafs*(1-mafs)
plot(x=mafs, y=vexp)

mafs <- pmax(rbeta(nsnp, shape1[1], 1) / 2, 1/(nid*2))
eff <- eff.func2(mafs)
vexp <- 2*eff^2*mafs*(1-mafs)
plot(x=mafs, y=vexp)

mafs <- pmax(rbeta(nsnp, shape1[1], 1) / 2, 1/(nid*2))
eff <- eff.func3(mafs)
vexp <- 2*eff^2*mafs*(1-mafs)
plot(x=mafs, y=vexp)

levels(params$eff.func) <- c("Neutral", "log(MAF)", expression(log(MAF)^2))
params$shape1 <- as.factor(params$shape1)
levels(params$shape1) <- c("~ Beta(0.3, 1)", "~ Beta(0.1, 1)", "~ Beta(0.05, 1)")

library(ggplot2)
ggplot(params, aes(x=shape1, y=proprare)) +
geom_bar(stat="identity", aes(fill=eff.func), position="dodge") +
scale_fill_brewer() +
theme_bw()