Post date: Jun 03, 2014 11:56:5 PM
Note, the LD data here are r2 estimates from bcftools that only consider loci with maf > 5%.
## read in LD
ldGla<-read.table("glaLd.txt",header=FALSE)
ldSla<-read.table("slaLd.txt",header=FALSE)
## calulcate distances in bp
dbpGla<-ldGla$V4-ldGla$V2
dbpSla<-ldSla$V4-ldSla$V2
## calculate averages, 5th quantile, and 95th quantile for different defined bins
bins<-c(seq(0,1000,100),seq(1000,100000,1000),seq(110000,1250000,10000))
N<-length(bins)-1
mnsGla<-rep(NA,N)
mnsSla<-rep(NA,N)
qsGla<-matrix(NA,nrow=N,ncol=2)
qsSla<-matrix(NA,nrow=N,ncol=2)
for (i in 1:N){
qsGla[i,]<-quantile(ldGla$V5[dbpGla > bins[i] & dbpGla <= bins[i+1]],probs=c(0.05,0.95))
qsSla[i,]<-quantile(ldSla$V5[dbpSla > bins[i] & dbpSla <= bins[i+1]],probs=c(0.05,0.95))
mnsGla[i]<-mean(ldGla$V5[dbpGla > bins[i] & dbpGla <= bins[i+1]])
mnsSla[i]<-mean(ldSla$V5[dbpSla > bins[i] & dbpSla <= bins[i+1]])
}
## allele frequency distn.
L<-206047
gaf<-scan("glaAfs.txt",n=L)
saf<-scan("slaAfs.txt",n=L)
mafGaf<-gaf
mafGaf[gaf > 0.5]<-1-gaf[gaf > 0.5]
mafSaf<-saf
mafSaf[saf > 0.5]<-1-saf[saf > 0.5]
outSaf<-hist(mafSaf[mafSaf > 0])
outGaf<-hist(mafGaf[mafGaf > 0])
## remove rare variants
x<-which((gaf < 0.01 | gaf > 0.99) | (saf < 0.01 | saf > 0.99))
## plots
pdf("varFig.pdf",width=8,height=12)
layout(matrix(c(1,2,3,3,4,4),nrow=3,ncol=2,byrow=T),widths=c(4,4),heights=c(4,4,4))
par(mar=c(5,5.5,2.5,2))
barplot(as.matrix(rbind(outGaf$counts/sum(outGaf$counts),outSaf$counts/sum(outSaf$counts))),beside=T,names=outGaf$mids,cex.names=0.7,legend.text=c("GLA","SLA"),xlab="minor allele frequency",ylab="proportion of loci",cex.lab=1.5,cex.axis=1.2)
mtext("(a) site allele frequency distribution",side=3,line=1,cex=1.2,adj=0)
plot(gaf[-x],saf[-x],cex=0.4,col="lightblue",xlab="GLA allele frequency",ylab="SLA allele frequency",cex.lab=1.5,cex.axis=1.2)
mtext("(b) joint allele frequency distribution",side=3,line=1,cex=1.2,adj=0)
plot(dbpGla,ldGla$V5,cex=0.4,col="lightblue",xlab="physical distance (bp)",ylab=expression(paste("linkage disequilibrium (",r^2,")")),cex.lab=1.5,cex.axis=1.2)
lines(bins[-1],mnsGla)
#lines(bins[-1],qsGla[,1],col="darkgray")
lines(bins[-1],qsGla[,2],col="darkgray")
mtext("(c) linkage disequilibrium in GLA",side=3,adj=0,line=1,cex=1.2)
plot(dbpSla,ldSla$V5,cex=0.4,col="lightblue",xlab="physical distance (bp)",ylab=expression(paste("linkage disequilibrium (",r^2,")")),cex.lab=1.5,cex.axis=1.2)
lines(bins[-1],mnsSla)
#lines(bins[-1],qsSla[,1],col="darkgray")
lines(bins[-1],qsSla[,2],col="darkgray")
title(main="(c) linkage disequilibrium in SLA",side=3,adj=0,line=1,cex=1.2)
dev.off()