Post date: Oct 02, 2019 5:23:59 PM
The directory for the effective migration surface analyses is /uufs/chpc.utah.edu/common/home/gompert-group2/projects/lmel_dimensions/mig_surface/.
I ran analyses with rare vs. common variants, and with a deme number of 200 or 400. The analysis was run with:
sbatch RunEems.sh
#!/bin/sh
#SBATCH --time=72:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=12
#SBATCH --account=gompert-kp
#SBATCH --partition=gompert-kp
#SBATCH --job-name=eems
#SBATCH --mail-type=FAIL
#SBATCH --mail-user=zach.gompert@usu.edu
cd /uufs/chpc.utah.edu/common/home/gompert-group2/projects/lmel_dimensions/mig_surface
perl RunEemsFork.pl params*ini
This runs 3 chains for each of the sets of conditions. Each chain is set up in its own *ini file. Here is the chain 1 file for each.
datapath = ./mel_common
mcmcpath = ./out/mel_common-400-chain1
nIndiv = 541
nSites = 20449
nDemes = 400
diploid = true
numMCMCIter = 4000000
numBurnIter = 2000000
numThinIter = 9999
mrateMuProposalS2 = 0.001
datapath = ./mel_common
mcmcpath = ./out/mel_common-chain1
nIndiv = 541
nSites = 20449
nDemes = 200
diploid = true
numMCMCIter = 4000000
numBurnIter = 2000000
numThinIter = 9999
mrateMuProposalS2 = 0.001
datapath = ./mel_rare
mcmcpath = ./out/mel_rare-400-chain1
nIndiv = 541
nSites = 47470
nDemes = 400
diploid = true
numMCMCIter = 4000000
numBurnIter = 2000000
numThinIter = 9999
mrateMuProposalS2 = 0.001
datapath = ./mel_rare
mcmcpath = ./out/mel_rare-chain1
nIndiv = 541
nSites = 47470
nDemes = 200
diploid = true
numMCMCIter = 4000000
numBurnIter = 2000000
numThinIter = 9999
mrateMuProposalS2 = 0.001
Then, I made some standard/summary plots (using specialized functions from the R package that comes with eems), and extracted effective migration estimates from the MCMC output/parameter estimates. I used averages across chains within a set of conditions. This is all in commands.R:
library("rEEMSplots")
library("rgdal")
library("rworldmap")
library("rworldxtra")
eems.plots(mcmcpath="mel_common-chain1",plotpath="o_mel_common-c1",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_common-chain2",plotpath="o_mel_common-c2",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_common-chain3",plotpath="o_mel_common-c3",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_rare-chain1",plotpath="o_mel_rare-c1",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_rare-chain2",plotpath="o_mel_rare-c2",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_rare-chain3",plotpath="o_mel_rare-c3",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_common-400-chain1",plotpath="o_mel_common-400-c1",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_common-400-chain2",plotpath="o_mel_common-400-c2",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_common-400-chain3",plotpath="o_mel_common-400-c3",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_rare-400-chain1",plotpath="o_mel_rare-400-c1",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_rare-400-chain2",plotpath="o_mel_rare-400-c2",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
eems.plots(mcmcpath="mel_rare-400-chain3",plotpath="o_mel_rare-400-c3",longlat=FALSE,add.grid=TRUE,add.demes=TRUE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84")
dimS<-read.tabel("dimSites.txt",header=FALSE)
Np<-dim(dimS)[1]
m<-matrix(NA,nrow=Np,ncol=4)
rf<-list.files(pattern="RData")
for(i in 1:4){
mx<-matrix(NA,nrow=Np,ncol=3)
for(j in 1:3){
n<-i+(j-1)*4
load(rf[n])
for(k in 1:Np){
a<-which.min(abs(dimS[k,2]-xym.values[,1])+abs(dimS[k,3]-xym.values[,2]))
mx[k,j]<-xym.values[a,3]
}
}
m[,i]<-apply(mx,1,mean)
}
row.names(m)<-dimS[,1]
write.table(m,file="effMigTab.txt",row.names=T,col.names=T,quote=FALSE)
## playing with adding dimension sites
coords_merc<-sp::spTransform(SpatialPoints(as.matrix(dimS[,2:3]),CRS("+proj=longlat +datum=WGS84")),CRS("+proj=merc +datum=WGS84"))
coord_merc<-coords_merc@coords
eems.plots(mcmcpath="mel_common-chain1",plotpath="o_mel_common-c1",longlat=FALSE,add.grid=TRUE,add.demes=FALSE,add.outline=TRUE,out.png=FALSE,projection.in="+proj=longlat +datum=WGS84",projection.out="+proj=merc +datum=WGS84",m.plot.xy={text(coord_merc);})
Matt added this to the path analysis and it has a non-zero coefficient (positive effect of migration on melissa). I want to make a nice plot for rare 400 (the one we used based on ridge regression) that has the sites for genetic data and the dimensions sites. The code at the end is a start towards this, but isn't working yet.