ld and imputation demo

ld.r

n <- 1000
nsnp <- 10

r <- matrix(0.8, nsnp, nsnp)
diag(r) <- 1

x1 <- rmvnorm(n, rep(0, nsnp), sigma=r)
x2 <- rmvnorm(n, rep(0, nsnp), sigma=r)

x1[1:500, nsnp] <- x2[1:500, 1]

for(i in 1:nsnp)
{
	x2[,i] <- x2[,i] + rnorm(n, sd=i/(nsnp*2))
	x1[,i] <- x1[,i] + rnorm(n, sd=(nsnp-i+1)/(nsnp*2))
}

x <- cbind(x1, x2)

cv <- 10
g1 <- 5
g2 <- 15
melted_correlation_matrix <- melt(cor(x))
dot_dat <- data.frame(Var1=c(cv, g1, g2), Var2=c(cv, g1, g2), value=1, val=c("Causal variant", "Genotyped", "Genotyped"))
# dot_dat <- data.frame(Var1=nsnp, Var2=nsnp, value=1)
ggplot(melted_correlation_matrix, aes(x=Var1, y=Var2, fill=value)) + 
	geom_tile() + 
	scale_fill_gradient(low="white", high="red") +
	theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
	geom_point(data=dot_dat, aes(colour=val), size=3) +
	labs(x="SNP number", y="SNP number", fill="R-square", colour="")

y <- x[,10] * sqrt(0.1) + rnorm(n, sd=sqrt(0.9))
cor(y, x[,10])^2

summary(lm(y ~ x[,10]))
summary(lm(y ~ x[,5]))
summary(lm(y ~ x[,15]))
summary(lm(y ~ x[,5] + x[,15]))

pval <- rep(0, ncol(x))
for(i in 1:ncol(x))
{
	pval[i] <- summary(lm(y ~ x[,i]))$coefficients[2,4]
}
dat <- data.frame(pval=-log10(pval), snp=1:ncol(x), lab="Untyped SNP", stringsAsFactors=FALSE)
dat$lab[cv] <- "Causal variant"
dat$lab[g1] <- "Genotyped variant"
dat$lab[g2] <- "Genotyped variant"

ggplot(dat, aes(x=snp, y=pval)) + geom_point(aes(colour=lab))

plot(-log10(pval), ylim=c(0,max(-log10(pval))))


y <- x[,5] * sqrt(0.1) + rnorm(n, sd=sqrt(0.9))
cor(y, x[,5])^2

summary(lm(y ~ x[,5]))
summary(lm(y ~ x[,4]))
summary(lm(y ~ x[,15]))
summary(lm(y ~ x[,4] + x[,15]))

pval <- rep(0, ncol(x))
for(i in 1:ncol(x))
{
	pval[i] <- summary(lm(y ~ x[,i]))$coefficients[2,4]
}

plot(-log10(pval), ylim=c(0,max(-log10(pval))))