Post date: Nov 14, 2019 4:37:19 PM
I summarized my bgc runs. I used these rather than what Maria had as hers used a MAF of 0.2 instead of an allele frequency difference of 0.2 to define ancestry informative markers.
This is in /uufs/chpc.utah.edu/common/home/u6000989/projects/frogs/clines/ and the script is bgcPlots.R (see below). Quantitative summaries in Maria's paper come from here.
## hybrid index
hi<-read.table("h.txt",sep=",",header=F)
hi_est<-hi[,1]
## alpha
f_a<-list.files(pattern="^alpha")
a_list<-vector("list",5)
for(i in 1:5){
a_list[[i]]<-read.table(f_a[i],sep=",")
}
a_mcmc<-cbind(a_list[[1]],a_list[[2]],a_list[[3]],a_list[[4]],a_list[[5]])
a_est<-apply(a_mcmc,1,quantile,probs=c(.5,0.025,0.975))
a_cred<-a_est[2,] > 0 | a_est[3,] < 0
## beta
f_b<-list.files(pattern="^beta")
b_list<-vector("list",5)
for(i in 1:5){
b_list[[i]]<-read.table(f_b[i],sep=",")
}
b_mcmc<-cbind(b_list[[1]],b_list[[2]],b_list[[3]],b_list[[4]],b_list[[5]])
b_est<-apply(b_mcmc,1,quantile,probs=c(.5,0.025,0.975))
sum(b_est[2,] > 0 | b_est[3,] < 0)
## 0, thus drop it
a_pp<-apply(a_mcmc>0,1,mean)
calcCline<-function(h=NULL, a=NULL, b=NULL){
z <- 2 * h * (1-h) * (a + (2 * h - 1) * b)
p <- h + z
p[p<0]<-0
p[p>1]<-1
return(p)
}
save(list=ls(),file="cline.rdat")
## plots
cm<-1.4
cl<-1.4
ca<-1.1
library(scales)
pdf("frogBgc.pdf",width=8,height=8)
layout(matrix(c(1,2,3,3),nrow=2,ncol=2,byrow=TRUE),widths=c(4,4),heights=c(4,4))
## hybrid index
par(mar=c(4.2,5,2.5,1))
hist(hi_est,col="purple",xlab="Hybrid Index",ylab="Frequency",main="",cex.lab=cl,cex.axis=ca,xlim=c(0,1))
title("A) Hybrid indexes",cex.main=cm)
## clines
h<-seq(0,1,0.01)
phi<-calcCline(h,a_est[1,1],b_est[1,1])
plot(h,phi,type='l',xlab="Hybrid Index",ylab="Ancestry Prob.",col="darkgray",lwd=.4,cex.lab=cl)
cls<-sample(2:7000,1000,replace=FALSE)
for(i in cls){
if(a_est[2,i] > 0){acol<-"dodgerblue4"}
else if(a_est[3,i] < 0){acol<-"firebrick2"}
else {acol<-"darkgray"}
phi<-calcCline(h,a_est[1,i],b_est[1,i])
lines(h,phi,col=acol,lwd=.4)
}
abline(a=0,b=1,lty=2)
title("B) Genomic clines",cex.main=cm)
plot(a_est[1,],pch=20,cex=c(0.1,0.25)[a_cred+1],col="darkgray",ylim=c(-5.5,5.5),xlab="SNP",ylab=expression(paste("Cline parameter ",alpha)),cex.lab=cl,cex.axis=ca)
x<-1:8002
segments(x[a_cred],a_est[2,a_cred],x[a_cred],a_est[3,a_cred],col=alpha("cadetblue",0.2),lwd=.2,lty=1)
abline(h=0,lty=2)
title("C) Cline parameter estiamtes",cex.main=cm)
dev.off()
phis<-matrix(NA,nrow=8002,ncol=length(h))
for(i in 1:8002){
phis[i,]<-calcCline(h,a_est[1,i],b_est[1,i])
}