L1-L3 Lines

Post date: Apr 18, 2018 5:17:5 PM

### Alt lines

#/uufs/chpc.utah.edu/common/home/u6000989/data/cmac_tradeoffs_AR/fullvar

sbatch callvar.sh

#/uufs/chpc.utah.edu/common/home/u6000989/data/cmac_tradeoffs_AR/fullvar/variants

perl vcf2gl.pl 0 calloVars.vcf

#Number of loci: 188; number of individuals 168

perl splitPops.pl cmacLalt.gl

# /uufs/chpc.utah.edu/common/home/u6000989/data/cmac_tradeoffs_AR/likelihood_data

perl wrap_qsub_slurm_popmod.pl cmacseries*

module load gcc gsl hdf5

~/bin/estpost_popmod -p p -s 0 -w 0 -o p_M-14.txt outp_bayes_M-14.hdf5

~/bin/estpost_popmod -p p -s 0 -w 0 -o p_L1-F100.txt outp_bayes_L1-F100.hdf5

~/bin/estpost_popmod -p p -s 0 -w 0 -o p_L2-F87.txt outp_bayes_L2-F87.hdf5

~/bin/estpost_popmod -p p -s 0 -w 0 -o p_L3-F85.txt outp_bayes_L3-F85.hdf5

cut -f 1 -d " " cmacseries_M-14.gl | perl -p -i -e 's/^(\S+)/\1/' > in_pM14.txt

## Delete headers

cut -d "," -f 2 p_M-14.txt | paste in_pM14.txt - > in_pM14-L1.txt

sed -i 's/\t/ /g' in_pM14-L1.txt

cut -d "," -f 2 p_M-14.txt | paste in_pM14.txt - > in_pM14-L2.txt

sed -i 's/\t/ /g' in_pM14-L2.txt

cut -d "," -f 2 p_M-14.txt | paste in_pM14.txt - > in_pM14-L3.txt

sed -i 's/\t/ /g' in_pM14-L3.txt

# Add appropriate headers

# /uufs/chpc.utah.edu/common/home/u0795641/projects/cmac_AR/alt_lines/alt_ne/wfabc

# ~/bin/wfAbcL -i in_pM14-L1.txt -l 198 -n 5000000 -v 0.3 -h 1

## L1 L2 L3 Line analyses

# /uufs/chpc.utah.edu/common/home/u0795641/projects/cmac_AR/alt_lines/alt_ne/wfabc

L1 <- read.table('L1/sel_L1incomp_s.txt', header=FALSE)

L2 <- read.table('L2/sel_L2incomp_s.txt', header=FALSE)

L3 <- read.table('L3/sel_L3incomp_s.txt', header=FALSE)

pm <- read.table('in_pM14.txt',header=FALSE)

p1 <- read.table('L1/in_pM14-L1.txt',header=FALSE,skip=2)

p2 <- read.table('L2/in_pM14-L2.txt',header=FALSE,skip=2)

p3 <- read.table('L3/in_pM14-L3.txt',header=FALSE,skip=2)

L14 <- read.table('~/projects/cmac_AR/L14/alt_ne/wfabc/fsm/sd3/incomp/sel_L14incomp.txt',header=FALSE, sep=' ')

snps<-read.table("~/projects/cmac_AR/L14/alt_ne/popgen/AlleleFreq/pLentilSnps.txt",header=FALSE)

snps[,1] <- sub('(:[[:digit:]]+$)','',snps[,1], perl=TRUE)

ids<-which(as.factor(snps[,1]) %in% as.factor(pm[,1]))

L14pm <- L14[ids,]

## # significant and magnitutde

# L1

mean(abs(L1[(L1[,4] > 0 | L1[,5] <0),1])) #0.1571298

mean(abs(L1[,1])) # 0.06645206

length(L1[(L1[,4] > 0 | L1[,5] <0),1]) #43

# L2

mean(abs(L2[(L2[,4] > 0 | L2[,5] <0),1])) #0.2938663

mean(abs(L2[,1])) # 0.102767

length(L2[(L2[,4] > 0 | L2[,5] <0),1]) #55

# L3

mean(abs(L3[(L3[,4] > 0 | L3[,5] <0),1])) #0.1416834

mean(abs(L3[,1])) # 0.02193066

length(L3[(L3[,4] > 0 | L3[,5] <0),1]) #10

# 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)

#rgb(maxColorValue=255, 215,48,39)

# pdf('selection_plots_alt.pdf', width=21, height=7)

# par(mar=c(6,6,2.5,2), mfrow=c(1,3))

#

# plot(x=L1[,1],y=L2[,1], ylab=('L2 S'),

# xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L1[,4]>0 & L2[,4]>0) | (L1[,5]<0 & L2[,5]<0)),19,

# ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

# ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,1))),

# col=ifelse(

# ((L1[,4]>0 & L2[,4]>0) | (L1[,5]<0 & L2[,5]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L1[,1], L2[,1]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L2", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

# title(main="(a)",adj=0, cex=1.25)

#

# plot(x=L1[,1],y=L3[,1], ylab=('L3 S'),

# xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L1[,4]>0 & L3[,4]>0) | (L1[,5]<0 & L3[,5]<0)),19,

# ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

# ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,1))),

# col=ifelse(

# ((L1[,4]>0 & L3[,4]>0) | (L1[,5]<0 & L3[,5]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L1[,1], L3[,1]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L3", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

# title(main="(b)",adj=0, cex=1.25)

#

# plot(x=L2[,1],y=L3[,1], ylab=('L3 S'),

# xlab=('L2 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L2[,4]>0 & L3[,4]>0) | (L2[,5]<0 & L3[,5]<0)),19,

# ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

# ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,1))),

# col=ifelse(

# ((L2[,4]>0 & L3[,4]>0) | (L2[,5]<0 & L3[,5]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L2[,1], L3[,1]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L3", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

# title(main="(c)",adj=0, cex=1.25)

#

# dev.off()

# pdf('selection_plots_alt_full.pdf', width=25, height=10)

# par(mfrow=c(2,5))

# par(mar=c(4.7,5.5,3,0.5))

# plot(x=L1[,1],y=L2[,1], ylab=('L2 S'),

# xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L1[,4]>0 & L2[,4]>0) | (L1[,5]<0 & L2[,5]<0)),19,

# ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

# ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,1))),

# col=ifelse(

# ((L1[,4]>0 & L2[,4]>0) | (L1[,5]<0 & L2[,5]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L1[,1], L2[,1]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L2", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L1[,1],y=L3[,1], ylab=('L3 S'),

# xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L1[,4]>0 & L3[,4]>0) | (L1[,5]<0 & L3[,5]<0)),19,

# ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

# ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,1))),

# col=ifelse(

# ((L1[,4]>0 & L3[,4]>0) | (L1[,5]<0 & L3[,5]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L1[,1], L3[,1]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L3", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L2[,1],y=L3[,1], ylab=('L3 S'),

# xlab=('L2 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L2[,4]>0 & L3[,4]>0) | (L2[,5]<0 & L3[,5]<0)),19,

# ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

# ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,1))),

# col=ifelse(

# ((L2[,4]>0 & L3[,4]>0) | (L2[,5]<0 & L3[,5]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L2[,1], L3[,1]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L3", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L1[,1],y=L14pm[,6], ylab=('L14A S'),

# xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L1[,4]>0 & L14pm[,9]>0) | (L1[,5]<0 & L14pm[,10]<0)),19,

# ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

# ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), 19,1))),

# col=ifelse(

# ((L1[,4]>0 & L14pm[,9]>0) | (L1[,5]<0 & L14pm[,10]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L1[,1], L14pm[,6]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L14A", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L1[,1],y=L14pm[,11], ylab=('L14B S'),

# xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L1[,4]>0 & L14pm[,14]>0) | (L1[,5]<0 & L14pm[,15]<0)),19,

# ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

# ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), 19,1))),

# col=ifelse(

# ((L1[,4]>0 & L14pm[,14]>0) | (L1[,5]<0 & L14pm[,15]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L1[,1], L14pm[,11]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L14B", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L2[,1],y=L14pm[,6], ylab=('L14A S'),

# xlab=('L2 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L2[,4]>0 & L14pm[,9]>0) | (L2[,5]<0 & L14pm[,10]<0)),19,

# ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

# ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), 19,1))),

# col=ifelse(

# ((L2[,4]>0 & L14pm[,9]>0) | (L2[,5]<0 & L14pm[,10]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L2[,1], L14pm[,6]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L14A", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L2[,1],y=L14pm[,11], ylab=('L14B S'),

# xlab=('L2 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L2[,4]>0 & L14pm[,14]>0) | (L2[,5]<0 & L14pm[,15]<0)),19,

# ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

# ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), 19,1))),

# col=ifelse(

# ((L2[,4]>0 & L14pm[,14]>0) | (L2[,5]<0 & L14pm[,15]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L2[,1], L14pm[,11]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L14B", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L3[,1],y=L14pm[,6], ylab=('L14A S'),

# xlab=('L3 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L3[,4]>0 & L14pm[,9]>0) | (L3[,5]<0 & L14pm[,10]<0)),19,

# ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,

# ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), 19,1))),

# col=ifelse(

# ((L3[,4]>0 & L14pm[,9]>0) | (L3[,5]<0 & L14pm[,10]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L3[,1], L14pm[,6]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L3", "L14A", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(x=L3[,1],y=L14pm[,11], ylab=('L14B S'),

# xlab=('L3 S'), cex.lab=2, cex.axis=1.7,, ylim=c(-1,1), xlim=c(-1,1),pch=ifelse(

# ((L3[,4]>0 & L14pm[,14]>0) | (L3[,5]<0 & L14pm[,15]<0)),19,

# ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,

# ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), 19,1))),

# col=ifelse(

# ((L3[,4]>0 & L14pm[,14]>0) | (L3[,5]<0 & L14pm[,15]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,5,113,176),

# ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L3[,1], L14pm[,11]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("L3", "L14B", "Significant in Both", "Neither")),

# col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

# rgb(maxColorValue=255,215,48,39), "darkgrey"), pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

#

# plot(y=L14[,11], x=L14[,6], xlab='S from F4 to F16A',

# ylab='S from F4 to F16B', cex.lab=2, cex.axis=1.7,, xlim=c(-1,1), ylim=c(-1,1),

# pch=ifelse(

# ((L14[,9]>0) | (L14[,10]<0)) &

# ((L14[,14]>0) | (L14[,15]<0)),19,

# ifelse(((L14[,14]>0) | (L14[,15]<0)), 19,

# ifelse(((L14[,9]>0) | (L14[,10]<0)), 19,1))),

# col=ifelse(

# ((L14[,9]>0) | (L14[,10]<0)) &

# ((L14[,14]>0) | (L14[,15]<0)),rgb(maxColorValue=255,215,48,39),

# ifelse(((L14[,14]>0) | (L14[,15]<0)), rgb(maxColorValue=255,253,174,97),

# ifelse(((L14[,9]>0) | (L14[,10]<0)), rgb(maxColorValue=255,5,113,176),"darkgrey"))))

# legend('bottomright', cex=1.5,legend=(paste('cor =', round(cor(L14[,6], L14[,11]),3))), bty='n')

# legend('topleft', bty='n', cex=1.5, legend=(c("F4 to F16A Significance", "F4 to F16B Significance",

# "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

# rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255,215,48,39), "darkgrey"),

# pch=c(19,19,19,1), bty='n')

# abline(h=0, v=0, lty=2)

# dev.off()

#### CORRELATIONS of uncertainty

sL1 <- read.table('L1/tol/Corrs_L1.txt',header=FALSE)

sL2 <- read.table('L2/tol/Corrs_L2.txt',header=FALSE)

sL3 <- read.table('L3/tol/Corrs_L3.txt',header=FALSE)

sL14 <- read.table('../../../L14/alt_ne/wfabc/fsm/sd3/incomp/correlations/Corrs_L14.txt')

sL14A <- read.table('../../../L14/alt_ne/wfabc/fsm/sd3/incomp/correlations/Corrs_L14A.txt')

sL14B <- read.table('../../../L14/alt_ne/wfabc/fsm/sd3/incomp/correlations/Corrs_L14B.txt')

pm <- read.table('in_pM14.txt',header=FALSE)

snps<-read.table("~/projects/cmac_AR/L14/alt_ne/popgen/AlleleFreq/pLentilSnps.txt",header=FALSE)

snps[,1] <- sub('(:[[:digit:]]+$)','',snps[,1], perl=TRUE)

ids<-which(as.factor(snps[,1]) %in% as.factor(pm[,1]))

sL14Apm <- sL14A[,ids]

sL14Bpm <- sL14B[,ids]

cors <- matrix(NA,ncol=4,nrow=12)

corrL1L2 <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL1L2[i,] <- cor(x=t(sL1[i,]),y=t(sL2[i,]))

}

mean(corrL1L2) # 0.1169224

##25%, 50%, 75% quantiles

quantile(corrL1L2, c(.25,.5,.75)) # 0.06679478 & 0.09597502 & 0.13027982

corrL1L3 <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL1L3[i,] <- cor(x=t(sL1[i,]),y=t(sL3[i,]))

}

mean(corrL1L3) # 0.1224706

quantile(corrL1L3, c(.25,.5,.75)) # 0.08234315 & 0.11991164 & 0.16207567

corrL2L3 <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL2L3[i,] <- cor(x=t(sL2[i,]),y=t(sL3[i,]))

}

mean(corrL2L3) # 0.2721624

quantile(corrL2L3, c(.25,.5,.75)) # 0.2388623 & 0.2703653 & 0.3053912

#### A

corrL1L14A <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL1L14A[i,] <- cor(x=t(sL1[i,]),y=t(sL14Apm[i,]))

}

mean(corrL1L14A) # 0.04317622

quantile(corrL1L14A, c(.25,.5,.75)) # 0.009543956 & 0.040645554 & 0.074248739

corrL2L14A <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL2L14A[i,] <- cor(x=t(sL2[i,]),y=t(sL14Apm[i,]))

}

mean(corrL2L14A) # -0.01664427

quantile(corrL2L14A, c(.25,.5,.75)) # -0.05068051 & -0.01170325 & 0.02199597

corrL3L14A <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL3L14A[i,] <- cor(x=t(sL3[i,]),y=t(sL14Apm[i,]))

}

mean(corrL3L14A) # 0.01548586

quantile(corrL3L14A, c(.25,.5,.75)) # -0.01751444 0.01721027 0.05112255

#### B

corrL1L14B <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL1L14B[i,] <- cor(x=t(sL1[i,]),y=t(sL14Bpm[i,]))

}

mean(corrL1L14B) # 0.04598432

quantile(corrL1L14B, c(.25,.5,.75)) # 0.01289220 & 0.04151493 & 0.07711225

corrL2L14B <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL2L14B[i,] <- cor(x=t(sL2[i,]),y=t(sL14Bpm[i,]))

}

mean(corrL2L14B) # -0.008915894

quantile(corrL2L14B, c(.25,.5,.75)) # -0.044136752 & -0.003007554 & 0.029787808

corrL3L14B <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL3L14B[i,] <- cor(x=t(sL3[i,]),y=t(sL14Bpm[i,]))

}

mean(corrL3L14B) # 0.02058023

quantile(corrL3L14B, c(.25,.5,.75)) # -0.01421123 & 0.02174151 & 0.06090165

### A and B

corrL14AL14B <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL14AL14B[i,] <- cor(x=t(sL14A[i,]),y=t(sL14B[i,]))

}

mean(corrL14AL14B) # 0.8081001

quantile(corrL14AL14B, c(.25,.5,.75)) # 0.7771013 0.8158660 0.8482890

## Through time

corrL14L14A <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL14L14A[i,] <- cor(x=t(sL14[i,]),y=t(sL14A[i,]))

}

mean(corrL14L14A) # 0.3465065

quantile(corrL14L14A, c(.25,.5,.75)) # 0.2960399 0.3541672 0.4037436

corrL14L14B <- matrix(NA, ncol=1, nrow=1000)

for(i in 1:1000){

corrL14L14B[i,] <- cor(x=t(sL14[i,]),y=t(sL14B[i,]))

}

mean(corrL14L14B) # 0.3415268

quantile(corrL14L14B, c(.25,.5,.75)) # 0.2925448 0.3489446 0.3973663

listcors <- list(corrL1L2,corrL1L3,corrL2L3,corrL1L14A,corrL2L14A,corrL3L14A,

corrL1L14B,corrL2L14B,corrL3L14B, corrL14AL14B, corrL14L14A, corrL14L14B)

cors <- matrix(NA,ncol=4,nrow=12)

for(i in 1:12){

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)

pdf('selection_plots_alt_full_2.pdf', width=14, height=35)

par(mfrow=c(5,2))

par(mar=c(4.7,5.5,3,0.5))

plot(x=L1[,1],y=L2[,1], ylab=('L2 S'), cex=1.5,

xlab=('L1 S'), cex.lab=2, cex.axis=1.7,,pch=ifelse(

((L1[,4]>0 & L2[,4]>0) | (L1[,5]<0 & L2[,5]<0)),19,

ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,1))),

col=ifelse(

((L1[,4]>0 & L2[,4]>0) | (L1[,5]<0 & L2[,5]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.096')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L2", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(a)",adj=0, cex.main=2.5, line=.4)

plot(x=L1[,1],y=L3[,1], ylab=('L3 S'), cex=1.5,

xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L1[,4]>0 & L3[,4]>0) | (L1[,5]<0 & L3[,5]<0)),19,

ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,1))),

col=ifelse(

((L1[,4]>0 & L3[,4]>0) | (L1[,5]<0 & L3[,5]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.120')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L3", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(b)",adj=0, cex.main=2.5, line=.4)

plot(x=L2[,1],y=L3[,1], ylab=('L3 S'), cex=1.5,

xlab=('L2 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L2[,4]>0 & L3[,4]>0) | (L2[,5]<0 & L3[,5]<0)),19,

ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,1))),

col=ifelse(

((L2[,4]>0 & L3[,4]>0) | (L2[,5]<0 & L3[,5]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.270')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L3", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(c)",adj=0, cex.main=2.5, line=.4)

plot(x=L1[,1],y=L14pm[,6], ylab=('L14A S'), cex=1.5,

xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L1[,4]>0 & L14pm[,9]>0) | (L1[,5]<0 & L14pm[,10]<0)),19,

ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), 19,1))),

col=ifelse(

((L1[,4]>0 & L14pm[,9]>0) | (L1[,5]<0 & L14pm[,10]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.041')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L14A", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(d)",adj=0, cex.main=2.5, line=.4)

plot(x=L1[,1],y=L14pm[,11], ylab=('L14B S'), cex=1.5,

xlab=('L1 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L1[,4]>0 & L14pm[,14]>0) | (L1[,5]<0 & L14pm[,15]<0)),19,

ifelse(((L1[,4]>0) | (L1[,5]<0)), 19,

ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), 19,1))),

col=ifelse(

((L1[,4]>0 & L14pm[,14]>0) | (L1[,5]<0 & L14pm[,15]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L1[,4]>0) | (L1[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.042')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L1", "L14B", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(e)",adj=0, cex.main=2.5, line=.4)

plot(x=L2[,1],y=L14pm[,6], ylab=('L14A S'), cex=1.5,

xlab=('L2 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L2[,4]>0 & L14pm[,9]>0) | (L2[,5]<0 & L14pm[,10]<0)),19,

ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), 19,1))),

col=ifelse(

((L2[,4]>0 & L14pm[,9]>0) | (L2[,5]<0 & L14pm[,10]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = -0.012')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L14A", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(f)",adj=0, cex.main=2.5, line=.4)

plot(x=L2[,1],y=L14pm[,11], ylab=('L14B S'), cex=1.5,

xlab=('L2 S'), cex.lab=2, cex.axis=1.7,,pch=ifelse(

((L2[,4]>0 & L14pm[,14]>0) | (L2[,5]<0 & L14pm[,15]<0)),19,

ifelse(((L2[,4]>0) | (L2[,5]<0)), 19,

ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), 19,1))),

col=ifelse(

((L2[,4]>0 & L14pm[,14]>0) | (L2[,5]<0 & L14pm[,15]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L2[,4]>0) | (L2[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = -0.003')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L2", "L14B", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(g)",adj=0, cex.main=2.5, line=.4)

plot(x=L3[,1],y=L14pm[,6], ylab=('L14A S'), cex=1.5,

xlab=('L3 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L3[,4]>0 & L14pm[,9]>0) | (L3[,5]<0 & L14pm[,10]<0)),19,

ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,

ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), 19,1))),

col=ifelse(

((L3[,4]>0 & L14pm[,9]>0) | (L3[,5]<0 & L14pm[,10]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L14pm[,9]>0) | (L14pm[,10]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.017')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L3", "L14A", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(h)",adj=0, cex.main=2.5, line=.4)

plot(x=L3[,1],y=L14pm[,11], ylab=('L14B S'), cex=1.5,

xlab=('L3 S'), cex.lab=2, cex.axis=1.7,, pch=ifelse(

((L3[,4]>0 & L14pm[,14]>0) | (L3[,5]<0 & L14pm[,15]<0)),19,

ifelse(((L3[,4]>0) | (L3[,5]<0)), 19,

ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), 19,1))),

col=ifelse(

((L3[,4]>0 & L14pm[,14]>0) | (L3[,5]<0 & L14pm[,15]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L3[,4]>0) | (L3[,5]<0)), rgb(maxColorValue=255,5,113,176),

ifelse(((L14pm[,14]>0) | (L14pm[,15]<0)), rgb(maxColorValue=255,253,174,97),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.022')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("L3", "L14B", "Significant in Both", "Neither")), col=c(rgb(maxColorValue=255,5,113,176),

rgb(maxColorValue=255,253,174,97), rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(i)",adj=0, cex.main=2.5, line=.4)

plot(y=L14[,11], x=L14[,6], xlab='S from F4 to F16A', cex=1.5,

ylab='S from F4 to F16B', cex.lab=2, cex.axis=1.7,,

pch=ifelse(

((L14[,9]>0) | (L14[,10]<0)) &

((L14[,14]>0) | (L14[,15]<0)),19,

ifelse(((L14[,14]>0) | (L14[,15]<0)), 19,

ifelse(((L14[,9]>0) | (L14[,10]<0)), 19,1))),

col=ifelse(

((L14[,9]>0) | (L14[,10]<0)) &

((L14[,14]>0) | (L14[,15]<0)),rgb(maxColorValue=255, 215,48,39),

ifelse(((L14[,14]>0) | (L14[,15]<0)), rgb(maxColorValue=255,253,174,97),

ifelse(((L14[,9]>0) | (L14[,10]<0)), rgb(maxColorValue=255,5,113,176),"darkgrey"))))

legend('bottomright', cex=1.5,legend=(paste('cor = 0.865')), bty='n')

legend('topleft', bty='n', cex=1.5, legend=(c("F4 to F16A Significance", "F4 to F16B Significance", "Significant in Both", "Neither")),

col=c(rgb(maxColorValue=255,5,113,176), rgb(maxColorValue=255,253,174,97),

rgb(maxColorValue=255, 215,48,39),"darkgrey"), pch=c(19,19,19,1))

abline(h=0, v=0, lty=2)

title(main="(j)",adj=0, cex.main=2.5, line=.4)

dev.off()