Post date: Mar 13, 2017 7:15:37 PM
#### SELECTION
####
#### WF ABC
## /uufs/chpc.utah.edu/common/home/u0795641/projects/cmac_AR/L14/alt_ne/wfabc/fsm/sd3/incomp
## sbatch abc.sh
## mv out* sd3
# /uufs/chpc.utah.edu/common/home/u0795641/projects/cmac_AR/wfabc/fsm/sd3/incomp
#/uufs/chpc.utah.edu/common/home/u0795641/bin/wfabc-dyn -i pLentilSnps.txt -l 198 -n 5000000 -v 0.3 -h 1 > out_incomp_snp198.txt
## Order selection plots by LD
library(RColorBrewer)
library(dendextend)
library(fields)
library(gplots)
N<-48
L<-21342
snps<-read.table("~/projects/cmac_AR/L14/alt_ne/ld/locusids.txt",header=FALSE)
cl<-1.6
ca<-1.2
soi<-read.table("~/projects/cmac_AR/L14/alt_ne/ld/pLentilSnps.txt",header=FALSE,skip=0)
ids<-which(as.character(snps[,1]) %in% as.character(soi[,1]))
GF1<-t(matrix(scan("~/projects/cmac_AR/L14/alt_ne/ld/L14-F1_p.txt",n=N*L,sep=","),nrow=N,ncol=L,byrow=T))
ldF1<-cor(t(GF1[ids,]))^2
dend <- as.dendrogram(hclust(dist(ldF1)))
dend1 <- color_branches(dend,k=7)
col_labels <- get_leaves_branches_col(dend1)
cols <- unique(col_labels)
cols <- c("#CC476B22", "#A86B0022", "#5D840022", "#00944F22", "#0093A922", "#4473D722", "#C03FBE22")
tree <- cutree(dend,k=7)
otree <- soi[order(tree),]
# length(which(tree==1)) # 31
#
# length(which(tree==2)) # 15
#
# length(which(tree==3)) # 45
#
# length(which(tree==4)) # 41
#
# length(which(tree==5)) # 23
#
# length(which(tree==6)) # 21
#
# length(which(tree==7)) # 22
rand <- order(tree)
# vector of x-axis points to change color:
X1 <- c(0, 31.5, 46.5, 91.5, 132.5, 155.5, 176.5)
X2 <- c(31.5, 46.5, 91.5, 132.5, 155.5, 176.5, 199)
# ## Individual selection plots
# V<-rep(1:198, each=1)
#
# set.seed(001) # just to make it reproducible
# rand <- sample(V)
L <- read.table("sel_L14incomp.txt",header=FALSE,sep=" ") ##for sd =.3
## manhattan plots
pdf("sel_manhat_L14.pdf",width=14,height=14)
par(mar=c(2.5, 1, 2, 1),oma=c(5,6,2,0), mfrow=c(3,1))
sym<-19
plot(1:198,L[rand, 1],ylim=c(-1.0,1.0),cex.axis=1.7,pch=sym,cex=1,
type="n", ylab='', xlab='')
segments(1:198,L[rand,4],1:198,L[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L[rand,1],pch=sym,cex=1, col=ifelse(L[rand,4] > 0, "red",
ifelse(L[rand,5]<0,"red","black")))
abline(h=0,lty=2,lwd=1.1)
lim <- par("usr")
rect(X1, lim[1]-1 , X2, lim[2]+1, border = NA, col = cols)
title(main="(a) F0 to F4",adj=0,cex.main=2.5, line=.4)
plot(1:198,L[rand,6],ylim=c(-1.0,1.0),cex.axis=1.7,pch=sym,cex=0.6,
type="n", ylab='', xlab='')
segments(1:198,L[rand,9],1:198,L[rand,10],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L[rand,6],pch=sym,cex=1, col=ifelse(L[rand,9] > 0, "red",
ifelse(L[rand,10]<0,"red","black")))
abline(h=0,lty=2,lwd=1.1)
lim <- par("usr")
rect(X1, lim[1]-1 , X2, lim[2]+1, border = NA, col = cols)
title(main="(b) F4 to F16A",adj=0,cex.main=2.5, line=.4)
plot(1:198,L[rand,11],ylim=c(-1.0,1.0),cex.axis=1.7,
pch=sym,cex=0.6,type="n", ylab='', xlab='')
segments(1:198,L[rand,14],1:198,L[rand,15],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L[rand,11],pch=sym,cex=1, col=ifelse(L[rand,14] > 0, "red",
ifelse(L[rand,15]<0,"red","black")))
abline(h=0,lty=2,lwd=1.1)
lim <- par("usr")
rect(X1, lim[1]-1 , X2, lim[2]+1, border = NA, col = cols)
title(main="(c) F4 to F16B",adj=0,cex.main=2.5, line=.4)
mtext(text='locus (SNP)', side=1, line=2.5, outer=TRUE, cex=3)
mtext(text='selection on lentil', side=2, line=2.5, outer=TRUE, cex=3)
dev.off()
## # significant and magnitutde
# P to F4
mean(abs(L[(L[,4] > 0 | L[,5] <0),1])) #0.7493587
mean(abs(L[,1])) # 0.356571
length(L[(L[,4] > 0 | L[,5] <0),1]) #45
# F4 to F16A
mean(abs(L[(L[,9] > 0 | L[,10] <0),6])) #0.3691194
mean(abs(L[,6])) # 0.1967997
length(L[(L[,9] > 0 | L[,10] <0),6]) #16
# F4 to F16B
mean(abs(L[(L[,14] > 0 | L[,15] <0),11])) #0.4621066
mean(abs(L[,11])) # 0.2044151
length(L[(L[,14] > 0 | L[,15] <0),11]) #13
## how many stay the same/flip for A
same <- which(L[,4]>0 & L[,9]>0 | L[,5]<0 & L[,10]<0)
flip <- which(L[,4]>0 & L[,10]<0 | L[,5]<0 & L[,9]>0)
## how many stay the same/flip for B
same <- which(L[,4]>0 & L[,9]>0 | L[,5]<0 & L[,10]<0)
flip <- which(L[,4]>0 & L[,10]<0 | L[,5]<0 & L[,9]>0)
## how many stay the same/flip for A B
same <- which((L[,9]>0 & L[,14]>0) | (L[,10]<0 & L[,15]<0))
flip <- which(L[,4]>0 & L[,15]<0 | L[,5]<0 & L[,14]>0)
(L[,4]>0 & L[,9]>0) | (L[,5]<0 & L[,10]<0)
## pairwise plots
# if same significant selection for both lines for both epochs, red (215,48,39)
# if same selection for epoch 1 and epoch 2a, orange (253,174,97)
# if same selection for epoch 1 and epoch2b, blue (5,113,176)
pdf('selection_plots1.pdf', height=21, width=7)
par(mar=c(6,6,2.5,2), mfrow=c(3,1))
plot(x=L[,1], y=L[,6], ylab=('S from F4 to F16A'), cex=1.5,
xlab=('S from F0 to F4'), xlim=c(-1,1), ylim=c(-1,1), cex.lab=2, cex.axis=1.7,
pch=ifelse(
((L[,4]>0 & L[,9]>0) | (L[,5]<0 & L[,10]<0)),19,
ifelse(((L[,4]>0) | (L[,5]<0)), 19,
ifelse(((L[,9]>0) | (L[,10]<0)), 19,1))),
col=ifelse(
((L[,4]>0 & L[,9]>0) | (L[,5]<0 & L[,10]<0)),rgb(maxColorValue=255, 215,48,39),
ifelse(((L[,4]>0) | (L[,5]<0)), rgb(maxColorValue=255, 253,174,97),
ifelse(((L[,9]>0) | (L[,10]<0)), rgb(maxColorValue=255, 5,113,176),"darkgrey"))))
legend('bottomright',legend=(paste('cor = 0.437')), bty='n',
cex=1.5)
legend('topleft', legend=(c("F0 to F4 Significance", "F4 to F16A Significance",
"Significant in Both",'Neither')), col=c(rgb(maxColorValue=255, 253,174,97),
rgb(maxColorValue=255, 5,113,176), rgb(maxColorValue=255, 215,48,39),'darkgrey'),
pch=c(19,19,19,1), bty='n', cex=1.5)
title(main="(a)",adj=0, cex=2.5, line=.4)
abline(h=0, v=0, lty=2)
plot(x=L[,1], y=L[,11], ylab=('S from F4 to F16B'), cex=1.5,
xlab=('S from F0 to F4'), xlim=c(-1,1), ylim=c(-1,1), cex.lab=2, cex.axis=1.7,
pch=ifelse(
((L[,4]>0 & L[,14]>0) | (L[,5]<0 & L[,15]<0)),19,
ifelse(((L[,4]>0) | (L[,5]<0)), 19,
ifelse(((L[,14]>0) | (L[,15]<0)), 19,1))),
col=ifelse(
((L[,4]>0 & L[,14]>0) | (L[,5]<0 & L[,15]<0)),rgb(maxColorValue=255, 215,48,39),
ifelse(((L[,4]>0) | (L[,5]<0)), rgb(maxColorValue=255, 253,174,97),
ifelse(((L[,14]>0) | (L[,15]<0)), rgb(maxColorValue=255, 5,113,176),"darkgrey"))))
legend('bottomright',legend=(paste('cor = 0.438')), bty='n',
cex=1.5)
legend('topleft', legend=(c("F0 to F4 Significance", "F4 to F16B Significance",
"Significant in Both","Neither")), col=c(rgb(maxColorValue=255, 253,174,97),
rgb(maxColorValue=255, 5,113,176), rgb(maxColorValue=255, 215,48,39),'darkgrey'),
pch=c(19,19,19,1), bty='n', cex=1.5)
title(main="(b)",adj=0,cex=2, line=.4)
abline(h=0, v=0, lty=2)
plot(y=L[,11], x=L[,6], xlab='S from F4 to F16A', cex=1.5,
ylab='S from F4 to F16B', cex.lab=2, xlim=c(-1,1), ylim=c(-1,1), cex.axis=1.7,
pch=ifelse(
((L[,9]>0) | (L[,10]<0)) &
((L[,14]>0) | (L[,15]<0)),19,
ifelse(((L[,14]>0) | (L[,15]<0)), 19,
ifelse(((L[,9]>0) | (L[,10]<0)), 19,1))),
col=ifelse(
((L[,9]>0) | (L[,10]<0)) &
((L[,14]>0) | (L[,15]<0)),rgb(maxColorValue=255,215,48,39),
ifelse(((L[,14]>0) | (L[,15]<0)), rgb(maxColorValue=255,5,113,176),
ifelse(((L[,9]>0) | (L[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))
legend('bottomright',legend=(paste('cor = 0.865')), bty='n',
cex=1.5)
legend('topleft', legend=(c("F4 to F16A Significance", "F4 to F16B Significance",
"Significant in Both","Neither")), col=c(rgb(maxColorValue=255,253,174,97),
rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,215,48,39),'darkgrey'),
pch=c(19,19,19,1), bty='n', cex=1.5)
title(main="(c)",adj=0, cex=2, line=.4)
abline(h=0, v=0, lty=2)
dev.off()
## Summary of H's
H <- read.table("het_L14incomp.txt",header=FALSE,sep=" ") ##for sd =.3
L <- read.table("sel_L14incomp.txt",header=FALSE,sep=" ") ##for sd =.3
## Individual selection plots
V<-rep(1:198, each=1)
set.seed(001) # just to make it reproducible
rand <- sample(V)
pdf("hL.pdf",width=14,height=14)
par(mfrow=c(3,1))
par(mar=c(2.5, 1, 2, 1), oma=c(5,5.5,2,0))
sym<-19
plot(1:198,H[rand,1],ylim=c(0,1.0),cex.lab=2,cex.axis=1.7,xlab="",ylab="",pch=sym,type="n",)
segments(1:198,H[rand,4],1:198,H[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,H[rand,1],pch=sym,cex=1)
abline(h=.5,lty=2,lwd=1.1)
title(main="(a) F0 to F4",adj=0,cex.main=2)
#mean(H[,1]) ## .5346348
plot(1:198,H[,6],ylim=c(0,1.0),cex.lab=2,cex.axis=1.7,xlab="",ylab="",pch=sym,type="n")
segments(1:198,H[,9],1:198,H[,10],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,H[,6],pch=sym,cex=1)
abline(h=.5,lty=2,lwd=1.1)
title(main="(b) F4 to F16A",adj=0,cex.main=2)
#mean(H[,6]) ## .476356
plot(1:198,H[,11],ylim=c(0,1.0),cex.lab=2,cex.axis=1.7,xlab="",ylab="",pch=sym,type="n")
segments(1:198,H[,14],1:198,H[,15],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,H[,11],pch=sym,cex=1)
abline(h=0.5,lty=2,lwd=1.1)
title(main="(c) F4 to F16B",adj=0,cex.main=2)
#mean(H[,11]) ## .4723527
mtext(text="locus (SNP)", side=1, line=1.5, outer=TRUE,cex=2)
mtext(text="heterozygous effect", side=2, line=2.5, outer=TRUE,cex=2.5)
dev.off()
pdf('scatter_sh.pdf', width=7, height=21)
par(mar=c(2.5,2,2.5,2), mfrow=c(3,1),oma=c(5,5.5,2,0))
plot(L[,1],H[,1], xlab='', ylab='', cex=1.75, cex.axis=2, ylim=c(0,1), xlim=c(-1,1))
abline(v=0,h=.5,lty=2)
legend('bottomright',legend=(paste('cor =', round(cor(L[,1], H[,1]),3))), bty='n', cex=1.5)
title(main="(a) P to F4",adj=0,cex.main=2, line=.4)
plot(L[,6],H[,6], xlab='', ylab='', cex=1.75, cex.axis=2, ylim=c(0,1), xlim=c(-1,1))
abline(v=0,h=.5, lty=2)
legend('bottomright',legend=(paste('cor =', round(cor(L[,6], H[,6]),3))), bty='n', cex=1.5)
title(main="(b) F4 to F16A",adj=0,cex.main=2, line=.4)
plot(L[,11],H[,11], xlab='',ylab='',cex=2, cex.axis=2, ylim=c(0,1), xlim=c(-1,1))
abline(v=0,h=.5,lty=2)
legend('bottomright',legend=(paste('cor =', round(cor(L[,11], H[,11]),3))), bty='n',cex=1.5)
title(main="(c) F4 to F16B",adj=0,cex.main=2, line=.4)
mtext(text='Selection Coefficients', side=1, line=2, outer=TRUE, cex=3)
mtext(text='Heterozygous Effect', side=2, line=2, outer=TRUE, cex=3)
dev.off()
hL14 <- read.table("correlations/Corrs_L14_Het.txt")
hL14A <- read.table("correlations/Corrs_L14A_Het.txt")
hL14B <- read.table("correlations/Corrs_L14B_Het.txt")
sL14 <- read.table('correlations/Corrs_L14.txt')
sL14A <- read.table('correlations/Corrs_L14A.txt')
sL14B <- read.table('correlations/Corrs_L14B.txt')
## Early
corrL14HetL14 <- matrix(NA, ncol=1, nrow=1000)
for(i in 1:1000){
corrL14HetL14[i,] <- cor(x=t(hL14[i,]),y=t(sL14[i,]))
}
mean(corrL14HetL14)
quantile(corrL14HetL14, c(.25,.5,.75))
## A
corrL14HetL14A <- matrix(NA, ncol=1, nrow=1000)
for(i in 1:1000){
corrL14HetL14A[i,] <- cor(x=t(hL14A[i,]),y=t(sL14A[i,]))
}
mean(corrL14HetL14A)
quantile(corrL14HetL14A, c(.25,.5,.75))
## B
corrL14HetL14B <- matrix(NA, ncol=1, nrow=1000)
for(i in 1:1000){
corrL14HetL14B[i,] <- cor(x=t(hL14B[i,]),y=t(sL14B[i,]))
}
mean(corrL14HetL14B)
quantile(corrL14HetL14B, c(.25,.5,.75))
listcors <- list(corrL14HetL14,corrL14HetL14A,corrL14HetL14B)
cors <- matrix(NA,ncol=4,nrow=3)
for(i in 1:3){
cors[i,] <- c(mean(listcors[[i]]),quantile(listcors[[i]], c(.25, 0.5, 0.75)))
}
xtable(cors, caption='To be filled', label='To be filled', digits=3)
## posterior probabilities of each Selection model
# 0 - same selection for both early and late
# 1 - two estimates of selection, early and late
# 2 - three estimates of selection, early and late(A & B)
S <- read.table("pp_incomp.txt",header=FALSE,sep=" ")
# S <- read.table("pp_incomp.txt",header=FALSE,sep=" ")
## which selection model 'won' for each SNP
model <- apply(S,1,function(x) which(x==max(x)))
length(which(model==1)) ## 53
length(which(model==2)) ## 93
length(which(model==3)) ## 105
pdf('pp_box.pdf')
par(mar=c(5,5,2,2))
boxplot(S, xlab='scenario', ylab='posterior probability', col='grey', xaxt='n', cex.lab=2, cex.main=2)
axis(1, at=1:3, labels=c(1,2,3))
dev.off()