Post date: Aug 01, 2015 1:39:20 AM
While the results from the 10 population pairs mostly make sense, the migration rates seem low (compared to T. cristinae and to equilibrium island-model estimates). I am (jointly) examining two possible causes of this. First, unlike the T. cristinae results, variants here were called for multiple species together (all 10 pairs of pops.), which mean many 'variants' will be invariant in any given pair. This likely affects the results. Second, the prior on divergence time had a low upper bound (20,000 generations) for some or many of the pairs. To overcome these issues I am re-running the analysis with only variants that have mean MAF > 5% for a given pair (averaged between the two populations), and I increased the upper bound of the prior on T to 1 million. Here are the new priors:
// default values for priors
param.Nlb = 200; // log unfiorm prior on Nanc, N0, N1
param.Nub = 20000;
param.Tlb = 10; // log uniform prior on split time
param.Tub = 1000000;
param.Mlb = 0; // uniform prior on migration rate
param.Mub = 0.2;
Here is the revised subSetLociRadiation.R script:
## grab random set of 5000 SNPs (same set) for radiation pops
L<-5074942
X<-sample(1:L,L,replace=FALSE)## this randomizes order
#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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5BCBOG_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5BCTUR_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5BMCG3_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5BMT_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5BS_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5CR_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5HFRS_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5LP_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5SM_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])
a<-apply(out,1,mean)
b<-which(a > 0.05 & a < 0.95)
x<-b[1:5000]
write.table(out[x,],"pMaf5VP_CxQ.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
The new ABC result files will have 'Maf5' in their names.