### re done on August 31st 2017 ########
##### SURVIVAL ##########
#getting the files with all quantile values and converting them to a data frame to get xfold values
## gla ms
g_enrich_all_ms = as.data.frame(all_glams_surv_s)$glams_xfold
g_enrich_east_ms = as.data.frame(east_glams_surv_s)$glams_xfold
g_enrich_west_ms = as.data.frame(west_glams_surv_s)$glams_xfold
## gla ac
g_enrich_all_ac = as.data.frame(all_glaac_surv_s)$glaac_xfold
g_enrich_east_ac = as.data.frame(east_glaac_surv_s)$glaac_xfold
g_enrich_west_ac = as.data.frame(west_glaac_surv_s)$glaac_xfold
## gla ms pvalues
g_enrich_all_ms_p = as.data.frame(all_glams_surv_s)$glams_pvals
g_enrich_east_ms_p= as.data.frame(east_glams_surv_s)$glams_pvals
g_enrich_west_ms_p = as.data.frame(west_glams_surv_s)$glams_pvals
## gla ac pvalues
g_enrich_all_ac_p = as.data.frame(all_glaac_surv_s)$glaac_pvals
g_enrich_east_ac_p = as.data.frame(east_glaac_surv_s)$glaac_pvals
g_enrich_west_ac_p = as.data.frame(west_glaac_surv_s)$glaac_pvals
g_ptsallms = g_enrich_all_ms_p
g_ptseastms = g_enrich_east_ms_p
g_ptswestms = g_enrich_west_ms_p
## for pch
g_ptsallms[g_ptsallms > 0.05] <- 1
g_ptsallms[g_ptsallms <= 0.05] <- 2
g_ptseastms[g_ptseastms > 0.05] <- 1
g_ptseastms[g_ptseastms <= 0.05]<- 2
g_ptswestms[g_ptswestms > 0.05] <- 1
g_ptswestms[g_ptswestms <= 0.05] <- 2
##for pch ac
g_ptsallac = g_enrich_all_ac_p
g_ptseastac = g_enrich_east_ac_p
g_ptswestac = g_enrich_west_ac_p
## for pch
g_ptsallac[g_ptsallac > 0.05] <- 1
g_ptsallac[g_ptsallac <= 0.05] <- 2
g_ptseastac[g_ptseastac > 0.05] <- 1
g_ptseastac[g_ptseastac <= 0.05]<- 2
g_ptswestac[g_ptswestac > 0.05] <- 1
g_ptswestac[g_ptswestac <= 0.05] <- 2
setEPS()
postscript("glaSurv_enrichments.eps", width=10, height=8)
pdf('glaSurv_enrichments.pdf',width=24, height=12)
par(mfrow = c(1,2))
#plot ms
par(mar=c(4,5.5,3,1.5))
plot(g_enrich_all_ms,lwd=2, col="blue", type = "o", pch = c(16,1)[g_ptsallms], ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
lines(g_enrich_east_ms, type="o",pch=c(15,0)[g_ptseastms], lty=2, col="red", lwd=2)
lines(g_enrich_west_ms, type="o", pch=c(18,5)[g_ptswestms], lty=3, col="black", lwd=2)
g_range = range(0, g_enrich_all_ms, g_enrich_east_ms, g_enrich_west_ms)
legend(6,2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3, lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="Survival - GLA on Medicago", cex.main=1.5,adj=0.1)
#plot ac
plot(g_enrich_all_ac, type="o", lwd=2,col="blue",pch = c(16,1)[g_ptsallac], ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
lines(g_enrich_east_ac, type="o",pch=c(15,0)[g_ptseastac],lty=2, col="red",lwd=2)
lines(g_enrich_west_ac, type="o", pch=c(18,5)[g_ptswestac],lty=3, col="black",lwd=2)
g_range = range(0, g_enrich_all_ac, g_enrich_east_ac, g_enrich_west_ac)
legend(6, 2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3,lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="Survival - GLA on Astragalus", cex.main=1.5,adj=0.1)
dev.off()
## SLA #########
#getting the files with all quantile values and converting them to a data frame to get xfold values
## sla ms
s_enrich_all_ms = as.data.frame(all_slams_surv_s)$slams_xfold
s_enrich_east_ms = as.data.frame(east_slams_surv_s)$slams_xfold
s_enrich_west_ms = as.data.frame(west_slams_surv_s)$slams_xfold
## sla ac
s_enrich_all_ac = as.data.frame(all_slaac_surv_s)$slaac_xfold
s_enrich_east_ac = as.data.frame(east_slaac_surv_s)$slaac_xfold
s_enrich_west_ac = as.data.frame(west_slaac_surv_s)$slaac_xfold
## sla ms pvalues
s_enrich_all_ms_p = as.data.frame(all_slams_surv_s)$slams_pvals
s_enrich_east_ms_p= as.data.frame(east_slams_surv_s)$slams_pvals
s_enrich_west_ms_p = as.data.frame(west_slams_surv_s)$slams_pvals
## sla ac pvalues
s_enrich_all_ac_p = as.data.frame(all_slaac_surv_s)$slaac_pvals
s_enrich_east_ac_p = as.data.frame(east_slaac_surv_s)$slaac_pvals
s_enrich_west_ac_p = as.data.frame(west_slaac_surv_s)$slaac_pvals
s_ptsallms = s_enrich_all_ms_p
s_ptseastms = s_enrich_east_ms_p
s_ptswestms = s_enrich_west_ms_p
## for pch
s_ptsallms[s_ptsallms > 0.05] <- 1
s_ptsallms[s_ptsallms <= 0.05] <- 2
s_ptseastms[s_ptseastms > 0.05] <- 1
s_ptseastms[s_ptseastms <= 0.05]<- 2
s_ptswestms[s_ptswestms > 0.05] <- 1
s_ptswestms[s_ptswestms <= 0.05] <- 2
##for pch ac
s_ptsallac = s_enrich_all_ac_p
s_ptseastac = s_enrich_east_ac_p
s_ptswestac = s_enrich_west_ac_p
## for pch
s_ptsallac[s_ptsallac > 0.05] <- 1
s_ptsallac[s_ptsallac <= 0.05] <- 2
s_ptseastac[s_ptseastac > 0.05] <- 1
s_ptseastac[s_ptseastac <= 0.05]<- 2
s_ptswestac[s_ptswestac > 0.05] <- 1
s_ptswestac[s_ptswestac <= 0.05] <- 2
setEPS()
postscript("slaSurv_enrichments.eps", width=10, height=8)
pdf('slaSurv_enrichments.pdf',width=24, height=12)
par(mfrow = c(1,2))
#plot ms
par(mar=c(4,5.5,3,1.5))
plot(s_enrich_all_ms,lwd=2, col="blue", type = "o", pch = c(16,1)[s_ptsallms], ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
lines(s_enrich_east_ms, type="o",pch=c(15,0)[s_ptseastms], lty=2, col="red", lwd=2)
lines(s_enrich_west_ms, type="o", pch=c(18,5)[s_ptswestms], lty=3, col="black", lwd=2)
s_range = range(0, s_enrich_all_ms, s_enrich_east_ms, s_enrich_west_ms)
legend(6,2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3, lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="Survival - SLA on Medicago", cex.main=1.5,adj=0.1)
#plot ac
plot(s_enrich_all_ac, type="o", lwd=2,col="blue",pch = c(16,1)[s_ptsallac], ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
lines(s_enrich_east_ac, type="o",pch=c(15,0)[s_ptseastac],lty=2, col="red",lwd=2)
lines(s_enrich_west_ac, type="o", pch=c(18,5)[s_ptswestac],lty=3, col="black",lwd=2)
s_range = range(0, s_enrich_all_ac, s_enrich_east_ac, s_enrich_west_ac)
legend(6, 2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3,lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="Survival - SLA on Astragalus", cex.main=1.5,adj=0.1)
dev.off()
## plotting all 4 together
pdf('expCor_surv_enrichments.pdf',width=24, height=12)
par(mfrow = c(2,2))
#### gla #########
#plot ms
par(mar=c(5,5.5,3,1.5))
plot(g_enrich_all_ms,lwd=2.5, col="red2", type = "o", pch = c(1,16)[g_ptsallms],cex=1.5, ylim = c(0,2.3), xaxt="n", xlab = " ", ylab="Xfold enrichments", cex.lab=2, cex.axis=1.5)
lines(g_enrich_east_ms, type="o",pch=c(1,16)[g_ptseastms], lty=1, col="limegreen", lwd=2.5)
lines(g_enrich_west_ms, type="o", pch=c(1,16)[g_ptswestms], lty=1, col="dodgerblue3", lwd=2.5)
g_range = range(0, g_enrich_all_ms, g_enrich_east_ms, g_enrich_west_ms)
legend(8,2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("red2","limegreen","dodgerblue3"), lty=1, lwd=2.5)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1), cex.axis=1.5)
title(main="Survival - GLA on Medicago", cex.main=2,adj=0)
#plot ac
plot(g_enrich_all_ac, type="o", lwd=2,col="red2",pch = c(1,16)[g_ptsallac], ylim = c(0,2.3), xaxt="n", xlab = " ", ylab = " ", cex.lab=2, cex.axis=1.5)
lines(g_enrich_east_ac, type="o",pch=c(1,16)[g_ptseastac],lty=1, col="limegreen",lwd=2.5)
lines(g_enrich_west_ac, type="o", pch=c(1,16)[g_ptswestac],lty=1, col="dodgerblue3",lwd=2.5)
g_range = range(0, g_enrich_all_ac, g_enrich_east_ac, g_enrich_west_ac)
legend(8, 2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("red2","limegreen","dodgerblue3"), lty=1,lwd=2.5)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1), cex.axis=1.5)
title(main="Survival - GLA on Astragalus", cex.main=2,adj=0)
###### sla ################
#plot ms
par(mar=c(5,5.5,3,1.5))
plot(s_enrich_all_ms,lwd=2, col="red2", type = "o", pch = c(1,16)[s_ptsallms], ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2, cex.axis=1.5)
lines(s_enrich_east_ms, type="o",pch=c(1,16)[s_ptseastms], lty=1, col="limegreen", lwd=2.5)
lines(s_enrich_west_ms, type="o", pch=c(1,16)[s_ptswestms], lty=1, col="dodgerblue3", lwd=2.5)
s_range = range(0, s_enrich_all_ms, s_enrich_east_ms, s_enrich_west_ms)
legend(8,2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("red2","limegreen","dodgerblue3"),lty=1, lwd=2.5)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1),cex.axis=1.5)
title(main="Survival - SLA on Medicago", cex.main=2,adj=0)
#plot ac
plot(s_enrich_all_ac, type="o", lwd=2,col="red2",pch = c(1,16)[s_ptsallac], ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab=" ", cex.lab=2, cex.axis=1.5)
lines(s_enrich_east_ac, type="o",pch=c(1,16)[s_ptseastac],lty=1, col="limegreen",lwd=2.5)
lines(s_enrich_west_ac, type="o", pch=c(1,16)[s_ptswestac],lty=1, col="dodgerblue3",lwd=2.5)
s_range = range(0, s_enrich_all_ac, s_enrich_east_ac, s_enrich_west_ac)
legend(8, 2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("red2","limegreen","dodgerblue3"), lty=1,lwd=2.5)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1),cex.axis=1.5)
title(main="Survival - SLA on Astragalus", cex.main=2,adj=0)
dev.off()
##plots
#enrichments for ms
enrich_all_ms = c(1.121,1.213,1.16,1.31,1.34,1.33,1.29,1.25,1.22,1.19)
enrich_east_ms = c(2.24, 1.65, 1.69, 1.56, 1.5, 1.4, 1.33, 1.28, 1.25, 1.21)
enrich_west_ms = c(2.107, 1.53, 1.43, 1.50, 1.45, 1.45, 1.42, 1.34, 1.26, 1.25)
#enrichment for ac
enrich_all_ac = c(1.82,1.43,1.5,1.47,1.46,1.43,1.33,1.30,1.26,1.27)
enrich_east_ac = c(1.82, 1.43, 1.49, 1.41, 1.40, 1.36, 1.31, 1.26, 1.23, 1.21)
enrich_west_ac = c(2.07, 1.28,1.39, 1.45, 1.44, 1.33, 1.3, 1.27,1.24,1.21)
setEPS()
postscript("glaSurv_enrichments.eps", width=10, height=8)
pdf('glaSurv_enrichments.pdf',width=24, height=12)
par(mfrow = c(1,2))
#plot ms
par(mar=c(4,5.5,3,1.5))
plot(enrich_all_ms, type="o",lwd=2, col="blue", ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
lines(enrich_east_ms, type="o",pch=22, lty=2, col="red", lwd=2)
lines(enrich_west_ms, type="o", pch=23, lty=3, col="black", lwd=2)
g_range = range(0, enrich_all_ms, enrich_east_ms, enrich_west_ms)
legend(6,2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3, lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="GLA on Medicago", cex.main=1.5,adj=0.1)
#plot ac
plot(enrich_all_ac, type="o", lwd=2,col="blue", ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
lines(enrich_east_ac, type="o",pch=22, lty=2, col="red",lwd=2)
lines(enrich_west_ac, type="o", pch=23, lty=3, col="black",lwd=2)
g_range = range(0, enrich_all_ac, enrich_east_ac, enrich_west_ac)
legend(6, 2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3,lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="GLA on Astragalus", cex.main=1.5,adj=0.1)
dev.off()
setEPS()
postscript("allms_Exp_enrichments.eps", width=12, height=9)
par(mar=c(7,8,8,7), mgp=c(4,1,0))
barplot(enrich_ms, xlab = "top SNPs %", ylab = "xfold enrichment", ylim = c(0,1.5), col = "mistyrose4", cex.lab=2, cex.axis=1.5, names.arg=c("1%","2%","3%","4%","5%","6%","7%","8%","9%" ,"10%"), cex.names=0.8)
dev.off()
pdf('glaSurv_enrichments.pdf',width=8, height=7)
#par(mfrow = c(1,2))
#plot ms
par(mar=c(4,5.5,3,1.5))
plot(enrich_east_ms, type="o",lwd=2, col="blue", ylim = c(0.9,2.3), xaxt="n", xlab = "Top SNPs %", ylab="Xfold enrichments", cex.lab=2)
lines(enrich_west_ms, type="o",pch=22, lty=2, col="red", lwd=2)
lines(enrich_east_ac, type="o",pch=22, lty=3, col="darkgreen",lwd=2)
lines(enrich_west_ac, type="o", pch=23, lty=4, col="black",lwd=2)
g_range = range(0, enrich_all_ms, enrich_east_ms, enrich_west_ms)
legend(6,2.3, c("Melissa-east-ms", "Melissa-west-ms", "Melissa-east-ac","Melissa-west-ac"), cex=1.3, col=c("blue","red","darkgreen","black"),pch=21:23, lty=1:3, lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
dev.off()
title(main="GLA on Medicago", cex.main=1.5,adj=0.1)
#plot ac
plot(enrich_all_ac, type="o", lwd=2,col="blue", ylim = c(0,2.3), xaxt="n", xlab = "Top SNPs %", ylab="xfold enrichments", cex.lab=2)
g_range = range(0, enrich_all_ac, enrich_east_ac, enrich_west_ac)
legend(6, 2.3, c("Melissa-all", "Melissa-east", "Melissa-west"), cex=1.3, col=c("blue","red","black"),pch=21:23, lty=1:3,lwd=2)
axis(1, at=1:10, labels=c(0.01, 0.02, 0.03, 0.04,0.05,0.06,0.07,0.08,0.09,0.1))
title(main="GLA on Astragalus", cex.main=1.5,adj=0.1)
dev.off()
setEPS()
postscript("Exp_results.eps", width =10, height = 10)
pdf('Exp_results.pdf',width=10, height=10)
#par(mfrow = c(2,1))
par(mar=c(6,8,8,6),mgp=c(4,1,0))
#survival
survcounts = c(17,17,13,9,18,15,5,8)
par(mar=c(7,8,2,8) + 0.1,mgp=c(4,1,0))
barplot(survcounts, col = c("darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue"), xlab = "Population and host plants", ylab = "Number of SNPs observed", cex.lab = 1.5, ylim = c(0,20))
#title(main="(A) ",cex.main = 1.5, adj=0)
legend(x=14.5,y=19,pch=c(22,22,22,22),pt.bg=c("darkgreen","darkolivegreen1","slateblue1","skyblue"),legend=c("GLA - Alfalfa","GLA - Astragalus","SLA - Alfalfa","SLA - Astragalus"),bty="o",cex=1)
dev.off()
#weight
wgtcounts = c(12,16,7,11,12,13,9,18,13,14,12,8)
par(mar=c(7,8,2,6) + 0.1,mgp=c(4,1,0))
barplot(survcounts, col = c("darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue"), xlab = "Population and host plants", ylab = "Number observed", cex.lab = 1.5, ylim = c(0,20))
title(main="(B) ",cex.main = 1.5, adj=0)
dev.off()