On June 18th 2017 (for evolution):
I redid the PCA for the map and treemix. Here is the link to Zach's notes on this:
https://sites.google.com/site/gompertlabnotes/home/researcher-pages/zach-gompert/lycaeides/lycaeides-host-associated-selection/entropyresultspca
Here is the folder where the R workspace for this plot is:
/uufs/chpc.utah.edu/common/home/u6007910/projects/lmelissa_hostAdaptation/entropy/results
Here is where the plot is:
/uufs/chpc.utah.edu/common/home/u6007910/projects/lmelissa_hostAdaptation/entropy/figs/pcplotPlnts_sam.pdf
Here is my code for redoing the plot:
pdf("../figs/pcplotPlnts_sam.pdf",width=8,heigh=8)
par(pty='s')
par(mar=c(5,5,0.5,0.5))
plot(pcgcov$x[,1],pcgcov$x[,2],col="gray",xlab=paste("PC1 (",pc1,"%)"),ylab=paste("PC2 (",pc2,"%)"),cex.lab=1.3,cex.axis=1.6,type='n')
for(j in 1:dim(plnt)[1]){
X<-which(ids==as.character(plnt[j,1]))
if (is.na(plnt[j,3])==TRUE){
myc<-"gray"
}
else{
if (as.character(plnt[j,3])=='M'){
myc<-"palegreen2"
}
else if (as.character(plnt[j,3])=='A'){
myc<-"plum"
}
else if (as.character(plnt[j,3])=='L'){
myc<-"royalblue"
}
else if (as.character(plnt[j,3])=='G'){
myc<-"sandybrown"
}
}
points(pcgcov$x[X,1],pcgcov$x[X,2],col=myc, pch=19)
}
text(popmns1,popmns2,unique(ids))
#legend(-0.5,-0.25,legend=c("alfalfa","Astragalus","Lupinus","Glycyrrhiza","Native sp."),fill=c("forestgreen","purple","blue","orange","gray"),bty='n')
dev.off()
On May 11th 2017.
This is on the cluster in the folder: lmelissa_hostAdapt/Lmelissa_treemix/april_final2017
This is to create the final tree. For this I unzipped the out_treemix1.treeout.gz file. We chose the one migration event tree using gunzip out_treemix1.treeout.gz > out_treemix1.treeout.
Go to R
#plot the pve figure
pve = read.table("pveMigrationEvents.txt", header=FALSE)
plot(pve$V1, pve$V2, xlab = "Number of migration events", ylab = "Proportion of variation explained")
setEPS()
postscript("pve_migrationevents.eps",width=12,height = 9)
par(mar=c(6,8,8,6),mgp=c(4,1,0))
plot(pve$V1, pve$V2, xlab = "Number of migration events", ylab = "Proportion of variation explained", type ="h", cex.lab=2, cex.axis = 1.5)
dev.off()
#using ape to get final tree
library(ape)
mytree <- read.tree("out_treemix1.treeout")
#check the node numbers to rotate the tree
nodelabels()
#to rotate to get the root of the tree in the bottom
rot.phy = rotate(mytree, node=29)
plot(rot.phy0)
#plot like this
mycol = c("green4", "green4", "mediumblue", "green4", "mediumblue", "mediumblue", "green4", "mediumblue", "mediumblue", "mediumblue", "mediumblue", "mediumblue", "mediumblue", "green4", "green4", "green4", "green4","green4", "green4", "green4", "green4", "mediumblue", "mediumblue", "green4", "green4", "green4","black","black")
plot(rot.phy, type = "phylogram", use.edge.length = TRUE, node.pos = NULL, show.tip.label = TRUE, show.node.label = FALSE,edge.color = "black", edge.width = 1, edge.lty = 1, font = 4,cex = par("cex"), adj = NULL, srt = 0, no.margin = FALSE,root.edge = FALSE, label.offset = 0.001, underscore = FALSE,x.lim = NULL, y.lim = NULL, direction = "rightwards", lab4ut = NULL, tip.color = mycol, plot = TRUE,rotate.tree = 0, open.angle = 0, node.depth = 1,align.tip.label = FALSE)
setEPS()
postscript("final_tree.eps",width=12,height = 12)
par(mar=c(6,8,8,6),mgp=c(4,1,0))
plot(rot.phy, type = "phylogram", use.edge.length = TRUE,
node.pos = NULL, show.tip.label = TRUE, show.node.label = FALSE,
edge.color = "black", edge.width = 1, edge.lty = 1, font = 3,
cex = par("cex"), adj = NULL, srt = 0, no.margin = FALSE,
root.edge = FALSE, label.offset = 0.001, underscore = FALSE,
x.lim = NULL, y.lim = NULL, direction = "rightwards",
lab4ut = NULL, tip.color = mycol, plot = TRUE,
rotate.tree = 0, open.angle = 0, node.depth = 1,
align.tip.label = FALSE)
dev.off()
Good resource for ape : http://schmitzlab.info/phylo.html
On May 9th 2017.
This saved on my desktop melissa_hostAdapt/final_figures/manuscript_figures
EXPERIMENTAL RESULTS BARPLOTS
pdf("Exp_barplots.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(14,11,9,3,13,14,5,9,13,16,11,9)
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 observed", cex.lab = 1.5, ylim = c(0,20))
title(main="(A) Survival ",cex.main = 1.5, adj=0)
legend(x=14.5,y=21,pch=c(22,22,22,22),pt.bg=c("darkgreen","darkolivegreen1","slateblue1","skyblue"),legend=c("GLA-Medicago","GLA-Astragalus","SLA-Medicago","SLA-Astragalus"),bty="o",cex=1, xpd=TRUE)
#weight
wgtcounts = c(11,11,10,6,16,13,19,12,8,10,4,17)
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) Weight ",cex.main = 1.5, adj=0)
dev.off()
********************************************************
ALL POPULATIONS SIMPLE SCATTERPLOT
setEPS()
postscript("allpop_scatter.eps", width = 12, height = 9)
par(mar=c(6,8,8,6),mgp=c(4,1,0))
bf = all$bfmeans
colors = rep("gray", length(bf))
colors[bf > 10 ] <- "slateblue1"
plot(all$position, bf, col = colors, pch=20, xlab = "Scaffold positions", ylab="Bayes Factors", cex.lab=2, cex.axis = 1.5)
dev.off()
EAST WEST PARALLEL SCATTERPLOT
pdf("eastWestParallel_scatter.pdf", width = 12, height = 9)
par(mar=c(6,8,8,6),mgp=c(4,1,0))
plot(wcomb, ecomb, pch=19, cex.lab=2, cex.axis = 1.5, xlab = "Bayesfactors for Western Populations", ylab = "Bayesfactors for Eastern Populations", main = "")
title(main="Environmental correlations between Eastern and Western populations",cex.main=1.6, adj=0.5)
rect(5,-2, 55, 5, density = NULL, angle = 45, col = "darkorange", border = NULL, lty = par("lty"), lwd = par("lwd"))
rect(-2, 5, 5, 34, density = NULL, angle = 45, col = "plum2", border = NULL, lty = par("lty"), lwd = par("lwd"))
rect(-20, -20, 5, 5, density = NULL, angle = 45, col = "darkgray", border = NULL, lty = par("lty"), lwd = par("lwd"))
dev.off()
EAST WEST HISTOGRAM
setEPS()
postscript("eastWest_hist.eps", width =9, height = 9)
par(mar=c(6,8,8,6),mgp=c(4,1,0))
hist(null, col="slateblue1",xlab="Number of shared SNPs",ylab ="Density",cex.lab=2,cex.axis=1.5, main =" ")
abline(v=58, lwd=4, lty=5)
dev.off()
SEX CHROMOSOMES
null<-rbinom(n=10000,size=58,prob=0.0237)
setEPS()
postscript("sexchromsnps_null.eps", width = 9, height =9)
pdf("sexchromsnps_null.pdf", width = 9, height =9)
par(mar=c(6,8,8,6),mgp=c(4,1,0))
hist(null, xlim=c(-0.5, 6), breaks=5,xlab ="Number of shared SNPs", ylab="Density", main=" ",cex.lab=2, cex.axis = 1.5, col="slateblue1")
abline(v=6,lwd=4, lty=5)
#title(main="Histogram of null",cex.main=2,adj=0.5)
dev.off()
STRUCTURAL ANNOTATION
FUNCTIONAL ANNOTATION