Post date: Apr 18, 2018 5:20:1 PM
## ABC Check of Selection
## check of ABC for L14A/B
## original simulation, 0.5,0.25,0.75,0.025,0.975
check <- read.table('check_L1incomp_s.txt', header=FALSE)
check <- check[which(!is.na(check[,1])),]
pdf('abc_check.pdf')
par(mar=c(4,4,3,1))
plot(x=check[,1],y=check[,2],ylab='', xlab='', cex.lab=2)
title(main="true vs estimated selection L14")
legend('bottomright',legend=(paste('cor =', round(cor(check[,1], check[,2]),3))), bty='n', cex=1.5)
mtext(text='True S', side=1, line=2, cex=1.5)
mtext(text='Estimated S', side=2, line=2, cex=1.5)
dev.off()
## RMSE
sqrt(mean((check[,1]-check[,2])^2)) #0.1954163
cor(check[,1], check[,2]) #0.39816
1-length(which((check[,1] < check[,5]) | (check[,1] > check[,6])))/2510 ## 0.9466135
mean(abs(check[,1]) < abs(check[,2])) # 0.1816733
mean((check[,1] > 0 & check[,2] > 0) | (check[,1] < 0 & check[,2] < 0)) # 0.6055777
### CHECKING ALT SD'S
L1 <- read.table("sd1/sel_L14incomp.txt",header=FALSE,sep=" ")
L2 <- read.table("sd2/sel_L14incomp.txt",header=FALSE,sep=" ")
L3 <- read.table("sd3/incomp/sel_L14incomp.txt",header=FALSE,sep=" ")
L4 <- read.table("sd4/sel_L14incomp.txt",header=FALSE,sep=" ")
L5 <- read.table("sd5/sel_L14incomp.txt",header=FALSE,sep=" ")
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)))
cols <- c("#CC476B22", "#A86B0022", "#5D840022", "#00944F22", "#0093A922", "#4473D722", "#C03FBE22")
tree <- cutree(dend,k=7)
otree <- soi[order(tree),]
rand <- order(tree)
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)
# V<-rep(1:251, each=1)
#
# set.seed(001) # just to make it reproducible
# rand <- sample(V)
pdf("sd_check_manhat.pdf",width=12,height=18)
par(mar=c(2,2,2,1), mfrow=c(5,1), oma=c(4,5,1,1))
sym<-19
plot(1:198,L1[rand, 1],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:198,L1[rand,4],1:198,L1[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L1[rand,1],pch=sym,cex=1, col=ifelse(L1[rand,4] > 0, "red", ifelse(L1[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) SD = 0.1",adj=0,cex.main=2)
plot(1:198,L2[rand, 1],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:198,L2[rand,4],1:198,L2[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L2[rand,1],pch=sym,cex=1, col=ifelse(L2[rand,4] > 0, "red", ifelse(L2[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="(b) SD = 0.2",adj=0,cex.main=2)
plot(1:198,L3[rand, 1],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:198,L3[rand,4],1:198,L3[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L3[rand,1],pch=sym,cex=1, col=ifelse(L3[rand,4] > 0, "red", ifelse(L3[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="(c) SD = 0.3",adj=0,cex.main=2)
plot(1:198,L4[rand, 1],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:198,L4[rand,4],1:198,L4[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L4[rand,1],pch=sym,cex=1, col=ifelse(L4[rand,4] > 0, "red", ifelse(L4[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="(d) SD = 0.4",adj=0,cex.main=2)
plot(1:198,L5[rand, 1],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:198,L5[rand,4],1:198,L5[rand,5],lwd=1,lty=3,col="black")
points(lwd=2.5, 1:198,L5[rand,1],pch=sym,cex=1, col=ifelse(L5[rand,4] > 0, "red", ifelse(L5[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="(e) SD = 0.5",adj=0,cex.main=2)
mtext(text='locus (SNP)', side=1, line=2, outer=TRUE, cex=2)
mtext(text='Selection on Lentil', side=2, line=2, outer=TRUE, cex=2)
dev.off()
mean(abs(L1[(L1[,4] > 0 | L1[,5] <0),1]), na.rm=TRUE) #.7087367
mean(abs(L1[,1]), na.rm=TRUE) # 3543269
length(L1[(L1[,4] > 0 | L1[,5] <0),1], na.rm=TRUE) #49
mean(abs(L2[(L2[,4] > 0 | L2[,5] <0),1])) #.7087367
mean(abs(L2[,1])) # 3543269
length(L2[(L2[,4] > 0 | L2[,5] <0),1]) #49
mean(abs(L3[(L3[,4] > 0 | L3[,5] <0),1])) #.7087367
mean(abs(L3[,1])) # 3543269
length(L3[(L3[,4] > 0 | L3[,5] <0),1]) #49
mean(abs(L4[(L4[,4] > 0 | L4[,5] <0),1])) #.7087367
mean(abs(L4[,1])) # 3543269
length(L4[(L4[,4] > 0 | L4[,5] <0),1]) #49
mean(abs(L5[(L5[,4] > 0 | L5[,5] <0),1])) #.7087367
mean(abs(L5[,1])) # 3543269
length(L5[(L5[,4] > 0 | L5[,5] <0),1]) #49
## UNCONSTRAINED MODEL Check
L <- read.table("sel_L14incomp_unc.txt",header=FALSE,sep=" ") ##for sd =.3
V<-rep(1:251, each=1)
set.seed(001) # just to make it reproducible
rand <- sample(V)
pdf("sL14_unc_check.pdf",width=12,height=12)
par(mar=c(2,2,2,1), mfrow=c(3,1), oma=c(4,5,1,1))
par(mar=c(4.7,5.5,3,0.5))
sym<-19
plot(1:251,L[rand,1],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:251,L[rand,4],1:251,L[rand,5],lwd=0.35,lty=3,col="black")
points(lwd=2.5, 1:251,L[rand,1],pch=sym,cex=0.25, col=ifelse(L[rand,4] > 0, "red", ifelse(L[rand,5]<0,"red","black")))
abline(h=0,lty=2,lwd=1.1)
title(main="(a) F0 to F4",adj=0,cex.main=2)
plot(1:251,L[rand,6],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:251,L[rand,9],1:251,L[rand,10],lwd=0.35,lty=3,col="black")
points(lwd=2.5, 1:251,L[rand,6],pch=sym,cex=0.25, col=ifelse(L[rand,9] > 0, "red", ifelse(L[rand,10]<0,"red","black")))
abline(h=0,lty=2,lwd=1.1)
title(main="(b) F4 to F16A",adj=0,cex.main=2)
plot(1:251,L[rand,11],ylim=c(-1.0,1.0),cex.lab=2,cex.axis=1.4,xlab="",ylab="",pch=sym,cex=0.6,type="n")
segments(1:251,L[rand,14],1:251,L[rand,15],lwd=0.35,lty=3,col="black")
points(lwd=2.5, 1:251,L[rand,11],pch=sym,cex=0.25, col=ifelse(L[rand,14] > 0, "red", ifelse(L[rand,15]<0,"red","black")))
abline(h=0,lty=2,lwd=1.1)
title(main="(c) F4 to F16B",adj=0,cex.main=2)
mtext(text='locus (SNP)', side=1, line=2, outer=TRUE, cex=2)
mtext(text='Selection on Lentil', side=2, line=2, outer=TRUE, cex=2)
dev.off()
## # significant and magnitutde
# P to F4
mean(abs(L[(L[,4] > 0 | L[,5] <0),1])) # 0.6080787
mean(abs(L[,1])) # 0.3538437
length(L[(L[,4] > 0 | L[,5] <0),1]) # 65
# F4 to F16A
mean(abs(L[(L[,9] > 0 | L[,10] <0),6])) # 0.3786581
mean(abs(L[,6])) #0.1379851
length(L[(L[,9] > 0 | L[,10] <0),6]) #46
# F4 to F16B
mean(abs(L[(L[,14] > 0 | L[,15] <0),11])) # 0.401417
mean(abs(L[,11])) # 0.1543956
length(L[(L[,14] > 0 | L[,15] <0),11]) # 51
### CORRELATIONS
sL14_unc <- read.table('Corrs_L1_unc.txt')
sL14A_unc <- read.table('Corrs_L2_unc.txt')
sL14B_unc <- read.table('Corrs_L3_unc.txt')
corrL14AL14B <- matrix(NA, ncol=1, nrow=300)
for(i in 1:300){
corrL14AL14B[i,] <- cor(x=t(sL14A_unc[i,]),y=t(sL14B_unc[i,]))
}
corrL14L14A <- matrix(NA, ncol=1, nrow=300)
for(i in 1:300){
corrL14L14A[i,] <- cor(x=t(sL14_unc[i,]),y=t(sL14A_unc[i,]))
}
corrL14L14B <- matrix(NA, ncol=1, nrow=300)
for(i in 1:300){
corrL14L14B[i,] <- cor(x=t(sL14_unc[i,]),y=t(sL14B_unc[i,]))
}
listcors <- list(corrL14AL14B, corrL14L14A, corrL14L14B)
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)
#### Tolerance of 1000 simulations
sL14_unc <- read.table('Corrs_L14_unc_tol.txt')
sL14A_unc <- read.table('Corrs_L14A_unc_tol.txt')
sL14B_unc <- read.table('Corrs_L14B_unc_tol.txt')
corrL14AL14B <- matrix(NA, ncol=1, nrow=1000)
for(i in 1:1000){
corrL14AL14B[i,] <- cor(x=t(sL14A_unc[i,]),y=t(sL14B_unc[i,]))
}
corrL14L14A <- matrix(NA, ncol=1, nrow=1000)
for(i in 1:1000){
corrL14L14A[i,] <- cor(x=t(sL14_unc[i,]),y=t(sL14A_unc[i,]))
}
corrL14L14B <- matrix(NA, ncol=1, nrow=1000)
for(i in 1:1000){
corrL14L14B[i,] <- cor(x=t(sL14_unc[i,]),y=t(sL14B_unc[i,]))
}
listcors <- list(corrL14AL14B, corrL14L14A, corrL14L14B)
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)