Post date: Jul 29, 2015 9:32:31 PM
We are estimating Nm for the non T. cristinae populations from the whole genome data using the same ABC approach (and priors) described here. However, here I am only using 5000 random SNPs as the stripe and selection candidate windows are not relevant. Files described are in /labs/evolution/projects/timema_speciation/demog. The R script subSetLociRadiation.R was used to grab the 5000 random SNPs and make the infiles, it is here:
## grab random set of 5000 SNPs (same set) for radiation pops
L<-5074942
X<-sample(1:L,5000,replace=FALSE)
## population pairs
p1<-read.table("p_timemaBCBOG_C.txt",header=FALSE)
p2<-read.table("p_timemaBCBOG_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanBCBOG_CxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaBCTUR_C.txt",header=FALSE)
p2<-read.table("p_timemaBCTUR_P.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanBCTUR_CxP.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaBMCG3_IC.txt",header=FALSE)
p2<-read.table("p_timemaBMCG3_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanBMCG3_ICxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaBMT_C.txt",header=FALSE)
p2<-read.table("p_timemaBMT_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanBMT_CxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaBS_C.txt",header=FALSE)
p2<-read.table("p_timemaBS_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanBS_CxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaCR_C.txt",header=FALSE)
p2<-read.table("p_timemaCR_CY.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanCR_CxCY.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaHFRS_M.txt",header=FALSE)
p2<-read.table("p_timemaHFRS_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanHFRS_MxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaLP_DF.txt",header=FALSE)
p2<-read.table("p_timemaLP_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanLP_DFxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaSM_Q.txt",header=FALSE)
p2<-read.table("p_timemaSM_RW.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanSM_QxRW.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
p1<-read.table("p_timemaVP_C.txt",header=FALSE)
p2<-read.table("p_timemaVP_Q.txt",header=FALSE)
out<-cbind(p1[X,3],p2[X,3])
write.table(out,"pRanVP_CxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
The main ABC runs involve 20 reps of 25000 simulations per population pair using the average (rounded down to nearest integer) sample size (the program uses a single same size). Here is an example:
cd /local/scratch/
sleep 19
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanVP_CxQ.txt -l 5000 -x 25000 -n 20 -a 1 > abcSim19_pRanVP_CxQ.txt
rsync -avz abcSim19_pRanVP_CxQ.txt /labs/evolution/projects/timema_speciation/demog/radsims/
rm abcSim19_pRanVP_CxQ.txt
I then calculated the observed summary statistics as,
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanBCBOG_CxQ.txt -l 5000 -n 19 -a 0 > obsSS_pRanBCBOG_CxQ.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanBCTUR_CxP.txt -l 5000 -n 16 -a 0 > obsSS_pRanBCTUR_CxP.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanBMCG3_ICxQ.txt -l 5000 -n 20 -a 0 > obsSS_pRanBMCG3_ICxQ.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanBMT_CxQ.txt -l 5000 -n 19 -a 0 > obsSS_pRanBMT_CxQ.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanBS_CxQ.txt -l 5000 -n 20 -a 0 > obsSS_pRanBS_CxQ.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanCR_CxCY.txt -l 5000 -n 20 -a 0 > obsSS_pRanCR_CxCY.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanHFRS_MxQ.txt -l 5000 -n 19 -a 0 > obsSS_pRanHFRS_MxQ.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanLP_DFxQ.txt -l 5000 -n 18 -a 0 > obsSS_pRanLP_DFxQ.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanSM_QxRW.txt -l 5000 -n 19 -a 0 > obsSS_pRanSM_QxRW.txt
~/bin/wfabc -i /labs/evolution/projects/timema_speciation/demog/pRanVP_CxQ.txt -l 5000 -n 20 -a 0 > obsSS_pRanVP_CxQ.txt