clines

##multipanel plot

pdf("multipanel.pdf",width=15,height=10, bg="white")

layout(matrix(c(1,3,4,2,3,4),nrow=3, ncol=2))

#hybrid index plot

h<-read.table("bgcout_pct3_3.hdf5_hybrid.txt",sep=",", header=F)

par(mar=c(4,6,3,4))

hist(h[,1],xlim=c(0,1),col="lavender",xlab="Hybrid index",ylab="Frequency",main="",cex.lab=2.5, cex.axis=1.5)

title(main="(A) ",cex.main = 2, adj=0)

#genomic clines plot

par(mar=c(4,6,3,4))

h<-seq(0,1,0.01)

phi<-calcCline(h,A3est[1,1],B3est[1,1])

plot(h,phi,type='l',xlab="Hybrid index",ylab=expression("Probs. of ancestry"),col="darkgray",lwd=.7,cex.lab=2.3, cex.axis=1.5)

for(i in subNeu){

phi<-calcCline(h,A3est[1,i],B3est[1,i])

lines(h,phi,col="darkgray",lwd=2)

}

for(i in ohiA){

phi<-calcCline(h,A3est[1,i],B3est[1,i])

lb<-which(phi == 1)[1]

ub<-length(phi)

phi[c(lb:ub)]<-1

lines(h,phi,col="darkgreen",lwd=2)

}

for(i in olowA){

phi<-calcCline(h,A3est[1,i],B3est[1,i])

lines(h,phi,col="darkolivegreen3",lwd=2)

}

for(i in ohiB){

phi<-calcCline(h,A3est[1,i],B3est[1,i])

lines(h,phi,col="darkorchid",lwd=3)

}

abline(a=0,b=1,lty=2,lwd=2)

title(main="(B) ",cex.main = 2, adj=0)

##scatterplot

#alpha

par(mar=c(4,6,3,4))

plot(alpha_mp_f[,3][1:1126],col=alpha_mp_f[,1],pch=21,cex=1,bg=alpha_mp_f[,1], xaxt="n", ylab="Probs. of ancestry",cex.lab=2.5, cex.axis=1.5, xlab="")

points(x=xsex, y=x[,3], col="red1", pch=20,cex=2)

axis(side=1, at=cticks, labels=c("Z","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22"),cex.axis=1.5)

abline(v=225, lwd=4)

abline(h=1, lwd=4)

title(main="(C) Alpha",cex.main = 2, adj=0)

#beta

#palette(mycolsb)

par(mar=c(5,6,3,4))

plot(beta_mp_f[,3][1:1126],col=beta_mp_f[,1],pch=21,cex=1,bg=beta_mp_f[,1], xaxt="n", ylab="Probs. of ancestry",cex.lab=2.5, cex.axis=1.5, xlab="")

points(x=xsexb, y=xb[,3], col="red1", pch=20,cex=2)

axis(side=1, at=cticks, labels=c("Z","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22"),cex.axis=2)

abline(v=225, lwd=4)

abline(h=1, lwd=4)

mtext("Chromosome",side=1,line=-1.5,outer=TRUE,cex=1.5)

title(main="(D) Beta",cex.main = 2, adj=0)

dev.off()