directory: /uufs/chpc.utah.edu/common/home/u6007910/projects/lmelissa_hostAdaptation/bayPass
## Did this in the bayPass folder. The plots are saved as shared_enrichment.eps and shared_enrichment. pdf.
##This is the code for creating the enrichment barplot. This is to see if how the xfold enrichments change with percentages
Here are the results:
##to get the enrichments
bndsw = quantile(wcomb,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndse = quantile(ecomb,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
wtop1 = which(wcomb > quantile(wcomb, probs=0.99,na.rm=T))
etop1 = which(ecomb > quantile(ecomb, probs=0.99,na.rm=T))
for (i in 1:10000){
x<-round(0.01*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop1 %in% etop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
wtop2 = which(wcomb > quantile(wcomb, probs=0.98,na.rm=T))
etop2 = which(ecomb > quantile(ecomb, probs=0.98,na.rm=T))
for (i in 1:10000){
x<-round(0.02*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop2 %in% etop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop3 = which(wcomb > quantile(wcomb, probs=0.97,na.rm=T))
etop3 = which(ecomb > quantile(ecomb, probs=0.97,na.rm=T))
for (i in 1:10000){
x<-round(0.03*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop3 %in% etop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop4 = which(wcomb > quantile(wcomb, probs=0.96,na.rm=T))
etop4 = which(ecomb > quantile(ecomb, probs=0.96,na.rm=T))
for (i in 1:10000){
x<-round(0.04*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop4 %in% etop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop5 = which(wcomb > quantile(wcomb, probs=0.95,na.rm=T))
etop5 = which(ecomb > quantile(ecomb, probs=0.95,na.rm=T))
for (i in 1:10000){
x<-round(0.05*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop5 %in% etop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop6 = which(wcomb > quantile(wcomb, probs=0.94,na.rm=T))
etop6 = which(ecomb > quantile(ecomb, probs=0.94,na.rm=T))
for (i in 1:10000){
x<-round(0.06*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop6 %in% etop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop7 = which(wcomb > quantile(wcomb, probs=0.93,na.rm=T))
etop7 = which(ecomb > quantile(ecomb, probs=0.93,na.rm=T))
for (i in 1:10000){
x<-round(0.07*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop7 %in% etop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop8 = which(wcomb > quantile(wcomb, probs=0.92,na.rm=T))
etop8 = which(ecomb > quantile(ecomb, probs=0.92,na.rm=T))
for (i in 1:10000){
x<-round(0.08*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop8 %in% etop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop9 = which(wcomb > quantile(wcomb, probs=0.91,na.rm=T))
etop9 = which(ecomb > quantile(ecomb, probs=0.91,na.rm=T))
for (i in 1:10000){
x<-round(0.09*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop9 %in% etop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
wtop10 = which(wcomb > quantile(wcomb, probs=0.9,na.rm=T))
etop10 = which(ecomb > quantile(ecomb, probs=0.9,na.rm=T))
for (i in 1:10000){
x<-round(0.1*206028,0)
wnull = sample(1:206028,x,replace=FALSE)
enull = sample(1:206028,x,replace=FALSE)
null[i] = sum(wnull %in% enull)
}
obs = sum(wtop10 %in% etop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
##the plot
setEPS()
postscript("shared_enrichments.eps", width=12, height=9)
par(mar=c(7,8,8,7), mgp=c(4,1,0))
barplot(enrich, xlab = "top SNPs %", ylab = "xfold enrichment", ylim = c(0,3), col = "yellow3", cex.lab=2, cex.axis=1.5, names.arg=c("0.01", "0.02", "0.03", "0.04","0.05","0.06",0.07","0.08","0.09","0.1"), cex.names=0.8)
dev.off()