#plot results from phase1 and verify with phase2 allele freq patterns
library(scattermore)
library(gridExtra)
library(VennDiagram)
## Loading required package: grid
## Loading required package: futile.logger
library(qvalue)
## Warning: package 'qvalue' was built under R version 3.5.2
library(outliers)
library(e1071)
## Warning: package 'e1071' was built under R version 3.5.2
#read in all sig values
all.sigs<-read.csv("~/Downloads/all.sigs.csv")
head(all.sigs)
## Chromosome Position lfmm.prec lfmm.hum lfmm.temp lfmm.corr lfmm.rand
## 1 2R 7075 0.9999652 0.9890107 0.9979580 0.9495994 0.9998185
## 2 2R 7376 0.9827711 0.9396630 0.9999695 0.9340263 0.9995228
## 3 2R 8031 0.8614518 0.8521034 0.9999695 0.9923323 0.9998978
## 4 2R 8146 0.9999652 0.9992322 0.9999695 0.9969829 0.9998185
## 5 2R 8753 0.9999652 0.9992322 0.9999695 0.9999976 0.9940057
## 6 2R 8855 0.9999652 0.9995471 0.9999695 0.9999976 0.9998185
## baye.prec baye.hum baye.temp baye.corr baye.rand vim.prec vim.hum
## 1 0.925163 0.913844 0.888063 0.935426 0.894285 -0.0003265053 0
## 2 0.948070 0.948620 0.951980 0.968750 0.950960 0.0000000000 0
## 3 0.928480 0.888700 0.937270 0.959680 0.958050 0.0000000000 0
## 4 0.966990 0.951600 0.955110 0.961830 0.953680 0.0000000000 0
## 5 0.762070 0.781880 0.914980 0.942730 0.849330 -0.0004747510 0
## 6 0.878700 0.902380 0.908420 0.899330 0.872020 0.0000000000 0
## vim.temp vim.corr vim.rand cors.prec cors.hum cors.temp
## 1 -0.015040239 0 -0.0001561354 0.0009345581 0.6290288 0.04542161
## 2 0.000000000 0 0.0000000000 0.7568032010 0.2930638 0.06728453
## 3 -0.005494650 0 -0.0690642596 0.0001598981 0.1069027 0.37926075
## 4 -0.007923008 0 0.0000000000 0.0802181561 0.4060133 0.02754343
## 5 -0.261914239 0 0.0000000000 0.2055432673 0.2021491 0.25747780
## 6 -0.018258204 0 -0.0100515393 0.1680477375 0.1342927 0.25407633
## cors.corr cors.rand
## 1 0.27584719 0.24754816
## 2 0.46752964 0.01687357
## 3 0.30538410 0.61262120
## 4 0.08819597 0.14306009
## 5 0.20319027 0.72081667
## 6 0.14479965 0.56095031
#bring in phase1 allele frequencies
phase1.all.freqs<-read.csv(file = "~/Desktop/anoph.phase2/phase1.all.freqs.csv")
#make column for id
phase1.all.freqs$id<-paste(phase1.all.freqs$chrom,phase1.all.freqs$POS)
head(phase1.all.freqs)
## chrom POS X.8.821 X.3.862 X.3.635 X0.384 X0.77 X4.341
## 1 X 247 0.1666667 0.4230769 0.6129032 0.05357143 0.59708738 0.18811881
## 2 X 2370 0.0000000 0.4230769 0.6129032 0.05357143 0.58737864 0.19306931
## 3 X 2895 0.0000000 0.4230769 0.6129032 0.05357143 0.59708738 0.18811881
## 4 X 3508 0.1666667 0.0000000 0.0000000 0.00000000 0.05339806 0.05940594
## 5 X 4144 0.6666667 0.5769231 0.3870968 0.00000000 0.03398058 0.08910891
## 6 X 4273 0.6666667 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
## X4.777 X5.747 X8.5 X9.25 X11.15 X11.233 X11.235
## 1 0.16265060 0.17582418 0.06666667 0.12500 0.12244898 0.02727273 0.04347826
## 2 0.16867470 0.18681319 0.06666667 0.12500 0.12244898 0.02727273 0.06521739
## 3 0.16265060 0.18131868 0.06666667 0.12500 0.12244898 0.02727273 0.06521739
## 4 0.05421687 0.04945055 0.10000000 0.03125 0.02040816 0.02727273 0.05434783
## 5 0.06024096 0.06593407 0.13333333 0.06250 0.28571429 0.62727273 0.47826087
## 6 0.00000000 0.00000000 0.00000000 0.00000 0.00000000 0.00000000 0.00000000
## X11.891 id
## 1 0.00000000 X 247
## 2 0.00000000 X 2370
## 3 0.00000000 X 2895
## 4 0.09782609 X 3508
## 5 0.33695652 X 4144
## 6 0.00000000 X 4273
#bring in phase2 allele frequencies
phase2.all.freqs<-read.csv(file = "~/Desktop/anoph.phase2/phase2.all.freqs.csv")
#make id column
phase2.all.freqs$id<-paste(phase2.all.freqs$chrom,phase2.all.freqs$POS)
#bring in dataframe with locality and environmental data for each individual phase1
phase1.alldata<-read.csv(file="~/Downloads/phase1_all_vairables_last.csv")
#subset to only the variables we need (lat,long,hum,temp,precip)
meta<-phase1.alldata[,c("latitude","longitude","p_annual","t_mean","h_mean","autocorrelated","random")]
#subset only the unique lat/longs
metapops<-unique(meta)
#sort by latitude to match the order of the allele frequency files
sort.pops.phase1 <-metapops[order(metapops$latitude),]
#bring in dataframe with locality and environmental data for each individual phase2
phase2.alldata<-read.csv(file = "~/Downloads/phase2_all_vairables_last.csv")
#subset to only the variables we need (lat,long,hum,temp,precip)
meta<-phase2.alldata[,c("latitude","longitude","p_annual","t_mean","h_mean","autocorrelated","random")]
#subset only the unique lat/longs
metapops<-unique(meta)
#sort by latitude to match the order of the allele frequency files
sort.pops.phase2<-metapops[order(metapops$latitude),]
sort.pops.phase2<-sort.pops.phase2[c(3,8,10:28),]
make dataframes where each row contains allele frequency data for the given SNP and the environmental values used to generate significant GEA test for that SNP (for phase 1 data)
#pull all significant outlier SNP IDs for LFMM humidity
out<- all.sigs[all.sigs$lfmm.hum < .01,c(1,2)]
#add the allele frequency at all 14 sites for each outlier SNP to the dataframe
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
#initialize empty dataframe
z<-data.frame()
#make matching dataframe with the value of the environmental variable used for testing for each SNP
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$h_mean)}
#fix column names to be consistent
colnames(z)<-c(sort.pops.phase1$latitude)
#combine into a single data.frame
lfmm.hum.outliers<-cbind(x,z,test=rep("lfmm.hum", times=nrow(x)))
#each row contains the identifier for an outlier SNP, the allele frequency of that SNP across phase1,
#the environmental values used for testing that SNP, and the name of the test
#lfmm precipitation
out<- all.sigs[all.sigs$lfmm.prec < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$p_annual)}
colnames(z)<-c(sort.pops.phase1$latitude)
lfmm.prec.outliers<-cbind(x,z,test=rep("lfmm.prec", times=nrow(x)))
#lfmm temperature
out<- all.sigs[all.sigs$lfmm.temp < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$t_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
lfmm.temp.outliers<-cbind(x,z,test=rep("lfmm.temp", times=nrow(x)))
#lfmm spatially autocorrelated variable
out<- all.sigs[all.sigs$lfmm.corr < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$autocorrelated)}
colnames(z)<-c(sort.pops.phase1$latitude)
lfmm.corr.outliers<-cbind(x,z,test=rep("lfmm.corr", times=nrow(x)))
#lfmm random variable
out<- all.sigs[all.sigs$lfmm.rand < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$random)}
colnames(z)<-c(sort.pops.phase1$latitude)
lfmm.rand.outliers<-cbind(x,z,test=rep("lfmm.rand", times=nrow(x)))
#bayescenv prec
out<- all.sigs[all.sigs$baye.prec < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$p_annual)}
colnames(z)<-c(sort.pops.phase1$latitude)
baye.prec.outliers<-cbind(x,z,test=rep("baye.prec", times=nrow(x)))
#bayescenv hum
out<- all.sigs[all.sigs$baye.hum < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$h_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
baye.hum.outliers<-cbind(x,z,test=rep("baye.hum", times=nrow(x)))
#bayescenv temp
out<- all.sigs[all.sigs$baye.temp < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$t_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
baye.temp.outliers<-cbind(x,z,test=rep("baye.temp", times=nrow(x)))
#bayescenv corr
out<- all.sigs[all.sigs$baye.corr < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$autocorrelated)}
colnames(z)<-c(sort.pops.phase1$latitude)
baye.corr.outliers<-cbind(x,z,test=rep("baye.corr", times=nrow(x)))
#bayescenv rand
out<- all.sigs[all.sigs$baye.rand < .01,c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$random)}
colnames(z)<-c(sort.pops.phase1$latitude)
baye.rand.outliers<-cbind(x,z,test=rep("baye.rand", times=nrow(x)))
#vim prec
out<- all.sigs[all.sigs$vim.prec > 1, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$p_annual)}
colnames(z)<-c(sort.pops.phase1$latitude)
vim.prec.outliers<-cbind(x,z,test=rep("vim.prec", times=nrow(x)))
#vim hum
out<- all.sigs[all.sigs$vim.hum > 1, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$h_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
vim.hum.outliers<-cbind(x,z,test=rep("vim.hum", times=nrow(x)))
#vim temp
out<- all.sigs[all.sigs$vim.temp > 1, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$t_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
vim.temp.outliers<-cbind(x,z,test=rep("vim.temp", times=nrow(x)))
#vim corr
out<- all.sigs[all.sigs$vim.corr > 1, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$autocorrelated)}
colnames(z)<-c(sort.pops.phase1$latitude)
vim.corr.outliers<-cbind(x,z,test=rep("vim.corr", times=nrow(x)))
#vim rand
out<- all.sigs[all.sigs$vim.rand > 1, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$random)}
colnames(z)<-c(sort.pops.phase1$latitude)
vim.rand.outliers<-cbind(x,z,test=rep("vim.rand", times=nrow(x)))
#cors prec
out<- all.sigs[all.sigs$cors.prec < .0001, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$p_annual)}
colnames(z)<-c(sort.pops.phase1$latitude)
cors.prec.outliers<-cbind(x,z,test=rep("cors.prec", times=nrow(x)))
#cors hum
out<- all.sigs[all.sigs$cors.hum < .0001, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$h_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
cors.hum.outliers<-cbind(x,z,test=rep("cors.hum", times=nrow(x)))
#cors temp
out<- all.sigs[all.sigs$cors.temp < .0001, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$t_mean)}
colnames(z)<-c(sort.pops.phase1$latitude)
cors.temp.outliers<-cbind(x,z,test=rep("cors.temp", times=nrow(x)))
#cors corr
out<- all.sigs[all.sigs$cors.corr < .0001, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$autocorrelated)}
colnames(z)<-c(sort.pops.phase1$latitude)
cors.corr.outliers<-cbind(x,z,test=rep("cors.corr", times=nrow(x)))
#cors rand
out<- all.sigs[all.sigs$cors.rand < .0001, c(1,2)]
x <- phase1.all.freqs[phase1.all.freqs$id %in% paste(out[,1],out[,2]),c(1:16)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase1$random)}
colnames(z)<-c(sort.pops.phase1$latitude)
cors.rand.outliers<-cbind(x,z,test=rep("cors.rand", times=nrow(x)))
#combine all the dataframes
full<-rbind(lfmm.hum.outliers,lfmm.prec.outliers,lfmm.temp.outliers,lfmm.rand.outliers,lfmm.corr.outliers,
baye.hum.outliers,baye.prec.outliers,baye.temp.outliers,baye.rand.outliers,baye.corr.outliers,
vim.hum.outliers,vim.prec.outliers,vim.temp.outliers,vim.rand.outliers,vim.corr.outliers,
cors.hum.outliers,cors.prec.outliers,cors.temp.outliers,cors.rand.outliers,cors.corr.outliers)
#build correlation dataframe for phase1 data
store<-c() #init vector
for (i in 1:nrow(full)){
store[i]<-cor(as.numeric(as.vector(full[i,17:30])), as.numeric(as.vector(full[i,3:16]))) #calculate correlation coefficient for each SNP
}
#add correlation coefficient to df
full$cor<-store
make dataframes where each row contains allele frequency data for the given SNP and the environmental values used to generate significant GEA test for that SNP (for phase 2 data)
#####step2 repeat this process with phase2 data
####
###
##
#
#pull all significant outlier SNP IDs for LFMM humidity
out<- all.sigs[all.sigs$lfmm.hum < .01,c(1,2)]
#add the allele frequency at all 14 sites for each outlier SNP to the dataframe
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
#initialize empty dataframe
z<-data.frame()
#make matching dataframe with the value of the environmental variable used for testing for each SNP
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$h_mean)}
#fix column names to be consistent
colnames(z)<-c(sort.pops.phase2$latitude)
#combine into a single data.frame
lfmm.hum.outliers<-cbind(x,z,test=rep("lfmm.hum", times=nrow(x)))
#each row contains the identifier for an outlier SNP, the allele frequency of that SNP across phase1,
#the environmental values used for testing that SNP, and the name of the test
#lfmm precipitation
out<- all.sigs[all.sigs$lfmm.prec < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$p_annual)}
colnames(z)<-c(sort.pops.phase2$latitude)
lfmm.prec.outliers<-cbind(x,z,test=rep("lfmm.prec", times=nrow(x)))
#lfmm temperature
out<- all.sigs[all.sigs$lfmm.temp < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$t_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
lfmm.temp.outliers<-cbind(x,z,test=rep("lfmm.temp", times=nrow(x)))
#lfmm spatially autocorrelated variable
out<- all.sigs[all.sigs$lfmm.corr < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$autocorrelated)}
colnames(z)<-c(sort.pops.phase2$latitude)
lfmm.corr.outliers<-cbind(x,z,test=rep("lfmm.corr", times=nrow(x)))
#lfmm random variable
out<- all.sigs[all.sigs$lfmm.rand < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$random)}
colnames(z)<-c(sort.pops.phase2$latitude)
lfmm.rand.outliers<-cbind(x,z,test=rep("lfmm.rand", times=nrow(x)))
#bayescenv prec
out<- all.sigs[all.sigs$baye.prec < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$p_annual)}
colnames(z)<-c(sort.pops.phase2$latitude)
baye.prec.outliers<-cbind(x,z,test=rep("baye.prec", times=nrow(x)))
#bayescenv hum
out<- all.sigs[all.sigs$baye.hum < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$h_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
baye.hum.outliers<-cbind(x,z,test=rep("baye.hum", times=nrow(x)))
#bayescenv temp
out<- all.sigs[all.sigs$baye.temp < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$t_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
baye.temp.outliers<-cbind(x,z,test=rep("baye.temp", times=nrow(x)))
#bayescenv corr
out<- all.sigs[all.sigs$baye.corr < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$autocorrelated)}
colnames(z)<-c(sort.pops.phase2$latitude)
baye.corr.outliers<-cbind(x,z,test=rep("baye.corr", times=nrow(x)))
#bayescenv rand
out<- all.sigs[all.sigs$baye.rand < .01,c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$random)}
colnames(z)<-c(sort.pops.phase2$latitude)
baye.rand.outliers<-cbind(x,z,test=rep("baye.rand", times=nrow(x)))
#vim prec
out<- all.sigs[all.sigs$vim.prec > 1, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$p_annual)}
colnames(z)<-c(sort.pops.phase2$latitude)
vim.prec.outliers<-cbind(x,z,test=rep("vim.prec", times=nrow(x)))
#vim hum
out<- all.sigs[all.sigs$vim.hum > 1, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$h_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
vim.hum.outliers<-cbind(x,z,test=rep("vim.hum", times=nrow(x)))
#vim temp
out<- all.sigs[all.sigs$vim.temp > 1, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$t_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
vim.temp.outliers<-cbind(x,z,test=rep("vim.temp", times=nrow(x)))
#vim corr
out<- all.sigs[all.sigs$vim.corr > 1, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$autocorrelated)}
colnames(z)<-c(sort.pops.phase2$latitude)
vim.corr.outliers<-cbind(x,z,test=rep("vim.corr", times=nrow(x)))
#vim rand
out<- all.sigs[all.sigs$vim.rand > 1, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$random)}
colnames(z)<-c(sort.pops.phase2$latitude)
vim.rand.outliers<-cbind(x,z,test=rep("vim.rand", times=nrow(x)))
#cors prec
out<- all.sigs[all.sigs$cors.prec < .0001, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$p_annual)}
colnames(z)<-c(sort.pops.phase2$latitude)
cors.prec.outliers<-cbind(x,z,test=rep("cors.prec", times=nrow(x)))
#cors hum
out<- all.sigs[all.sigs$cors.hum < .0001, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$h_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
cors.hum.outliers<-cbind(x,z,test=rep("cors.hum", times=nrow(x)))
#cors temp
out<- all.sigs[all.sigs$cors.temp < .0001, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$t_mean)}
colnames(z)<-c(sort.pops.phase2$latitude)
cors.temp.outliers<-cbind(x,z,test=rep("cors.temp", times=nrow(x)))
#cors corr
out<- all.sigs[all.sigs$cors.corr < .0001, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$autocorrelated)}
colnames(z)<-c(sort.pops.phase2$latitude)
cors.corr.outliers<-cbind(x,z,test=rep("cors.corr", times=nrow(x)))
#cors rand
out<- all.sigs[all.sigs$cors.rand < .0001, c(1,2)]
x <- phase2.all.freqs[phase2.all.freqs$id %in% paste(out[,1],out[,2]),c(1:23)]
z<-data.frame()
for (i in 1:nrow(x)){z<-rbind(z,sort.pops.phase2$random)}
colnames(z)<-c(sort.pops.phase2$latitude)
cors.rand.outliers<-cbind(x,z,test=rep("cors.rand", times=nrow(x)))
#combine all the dataframes
full.2<-rbind(lfmm.hum.outliers,lfmm.prec.outliers,lfmm.temp.outliers,lfmm.rand.outliers,lfmm.corr.outliers,
baye.hum.outliers,baye.prec.outliers,baye.temp.outliers,baye.rand.outliers,baye.corr.outliers,
vim.hum.outliers,vim.prec.outliers,vim.temp.outliers,vim.rand.outliers,vim.corr.outliers,
cors.hum.outliers,cors.prec.outliers,cors.temp.outliers,cors.rand.outliers,cors.corr.outliers)
calculate correlation coefficient for phase 2 data
#build correlation dataframe for phase2 data
store<-c() #init vector
for (i in 1:nrow(full.2)){
store[i]<-cor(as.numeric(as.vector(full.2[i,24:44])), as.numeric(as.vector(full.2[i,3:23]))) #perform linear reg
}
#add as column in full.2
full.2$cor<-store
Test each SNP for matching correlation coefficient and < .5 difference in magnitude
#subset corframe 1 SNPs found in cor.frame 2
hom.full <- full[paste(full[,1],full[,2],full[,31]) %in% paste(full.2[,1],full.2[,2],full.2[,45]),]
#check to make sure all cells are homologous between these data frames
table(paste(hom.full$chrom,hom.full$pos,hom.full$test) == paste(full.2$chrom,full.2$pos,full.2$test))
##
## TRUE
## 41622
#absolute value difference in correlation coefficient
hist(abs(hom.full$cor-full.2$cor))
# test how many SNPs have mismatched directional relationships between phase 1 and 2
table(hom.full$cor > 0 & full.2$cor < 0 | hom.full$cor < 0 & full.2$cor > 0)
##
## FALSE TRUE
## 37299 4323
#make dataframe containing info and vetting results for each SNP
vetted<-data.frame(chrom=hom.full$chrom,
pos=hom.full$POS,
test=hom.full$test,
pass.direction.test=hom.full$cor < 0 & full.2$cor < 0 | hom.full$cor > 0 & full.2$cor > 0,
pass.strength.test=abs(hom.full$cor-full.2$cor) < .5)
estimate false positive rate overall and for each candidate SNP dataset
#calculate overall false positive rate
1-sum(vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/nrow(vetted) #total tested SNPs
## [1] 0.253736
#because we expect our random variable (random between phases 1 and 2) to inflate the overall false positive rate,
#calculate false positive rate for only the real environmental variables
vet<-vetted[vetted$test != "lfmm.rand" & vetted$test != "baye.rand" & vetted$test != "vim.rand"
& vetted$test != "lfmm.corr" & vetted$test != "baye.corr" & vetted$test != "vim.corr"
& vetted$test != "cors.rand" & vetted$test != "cors.corr",]
1-sum(vet$pass.direction.test == TRUE & vet$pass.strength.test == TRUE)/nrow(vet) #total tested SNPs
## [1] 0.2479263
#lfmm calculate number of successfully vetted SNPs and the false positive rate for each test
#lfmm hum
sum(vetted$test == "lfmm.hum") #total SNPs
## [1] 1183
sum(vetted$test == "lfmm.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 535
1-sum(vetted$test == "lfmm.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "lfmm.hum") #false positive rate
## [1] 0.5477599
#lfmm prec
sum(vetted$test == "lfmm.prec") #total SNPs
## [1] 306
sum(vetted$test == "lfmm.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 259
1-sum(vetted$test == "lfmm.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "lfmm.prec") #false positive rate
## [1] 0.1535948
#lfmm temp
sum(vetted$test == "lfmm.temp") #total SNPs
## [1] 5013
sum(vetted$test == "lfmm.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 4418
1-sum(vetted$test == "lfmm.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "lfmm.temp") #false positive rate
## [1] 0.1186914
#lfmm rand
sum(vetted$test == "lfmm.rand") #total SNPs
## [1] 254
sum(vetted$test == "lfmm.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 105
1-sum(vetted$test == "lfmm.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "lfmm.rand") #false positive rate
## [1] 0.5866142
#lfmm corr
sum(vetted$test == "lfmm.corr") #total SNPs
## [1] 2553
sum(vetted$test == "lfmm.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 1930
1-sum(vetted$test == "lfmm.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "lfmm.corr") #false positive rate
## [1] 0.2440266
#baye calculate number of successfully vetted SNPs and the false positive rate for each test
#baye hum
sum(vetted$test == "baye.hum") #total SNPs
## [1] 1954
sum(vetted$test == "baye.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 1105
1-sum(vetted$test == "baye.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "baye.hum") #false positive rate
## [1] 0.4344933
#baye prec
sum(vetted$test == "baye.prec") #total SNPs
## [1] 1730
sum(vetted$test == "baye.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 1592
1-sum(vetted$test == "baye.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "baye.prec") #false positive rate
## [1] 0.07976879
#baye temp
sum(vetted$test == "baye.temp") #total SNPs
## [1] 2453
sum(vetted$test == "baye.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 2390
1-sum(vetted$test == "baye.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "baye.temp") #false positive rate
## [1] 0.02568284
#baye rand
sum(vetted$test == "baye.rand") #total SNPs
## [1] 684
sum(vetted$test == "baye.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 431
1-sum(vetted$test == "baye.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "baye.rand") #false positive rate
## [1] 0.369883
#baye corr
sum(vetted$test == "baye.corr") #total SNPs
## [1] 808
sum(vetted$test == "baye.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 651
1-sum(vetted$test == "baye.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "baye.corr") #false positive rate
## [1] 0.1943069
#vim calculate number of successfully vetted SNPs and the false positive rate for each test
#vim hum
sum(vetted$test == "vim.hum") #total SNPs
## [1] 384
sum(vetted$test == "vim.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 298
1-sum(vetted$test == "vim.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "vim.hum") #false positive rate
## [1] 0.2239583
#vim prec
sum(vetted$test == "vim.prec") #total SNPs
## [1] 193
sum(vetted$test == "vim.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 144
1-sum(vetted$test == "vim.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "vim.prec") #false positive rate
## [1] 0.253886
#vim temp
sum(vetted$test == "vim.temp") #total SNPs
## [1] 176
sum(vetted$test == "vim.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 139
1-sum(vetted$test == "vim.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "vim.temp") #false positive rate
## [1] 0.2102273
#vim rand
sum(vetted$test == "vim.rand") #total SNPs
## [1] 388
sum(vetted$test == "vim.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 297
1-sum(vetted$test == "vim.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "vim.rand") #false positive rate
## [1] 0.2345361
#vim corr
sum(vetted$test == "vim.corr") #total SNPs
## [1] 332
sum(vetted$test == "vim.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 250
1-sum(vetted$test == "vim.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "vim.corr") #false positive rate
## [1] 0.246988
#cors calculate number of successfully vetted SNPs and the false positive rate for each test
#cors hum
sum(vetted$test == "cors.hum") #total SNPs
## [1] 239
sum(vetted$test == "cors.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 131
1-sum(vetted$test == "cors.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "cors.hum") #false positive rate
## [1] 0.4518828
#cors prec
sum(vetted$test == "cors.prec") #total SNPs
## [1] 12245
sum(vetted$test == "cors.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 8442
1-sum(vetted$test == "cors.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "cors.prec") #false positive rate
## [1] 0.3105757
#cors temp
sum(vetted$test == "cors.temp") #total SNPs
## [1] 164
sum(vetted$test == "cors.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 131
1-sum(vetted$test == "cors.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "cors.temp") #false positive rate
## [1] 0.2012195
#cors rand
sum(vetted$test == "cors.rand") #total SNPs
## [1] 68
sum(vetted$test == "cors.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 18
1-sum(vetted$test == "cors.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "cors.rand") #false positive rate
## [1] 0.7352941
#cors corr
sum(vetted$test == "cors.corr") #total SNPs
## [1] 10495
sum(vetted$test == "cors.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE) #nSNPs
## [1] 7795
1-sum(vetted$test == "cors.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE)/sum(vetted$test == "cors.corr") #false positive rate
## [1] 0.2572654
#add column saying whether passed both
vetted$vet<-vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE
vetted$vet[vetted$vet == TRUE]<-"pass"
vetted$vet[vetted$vet == FALSE]<-"fail"
#write out each vetted SNP dataset
write.csv(vetted$test == "lfmm.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/lfmm.hum.vetted.snps.csv" , quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "lfmm.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/lfmm.prec.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "lfmm.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/lfmm.temp.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "lfmm.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/lfmm.rand.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "lfmm.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/lfmm.corr.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "baye.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/baye.hum.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "baye.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/baye.prec.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "baye.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/baye.temp.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "baye.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/baye.rand.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "baye.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/baye.corr.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "vim.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/vim.hum.vetted.snps.csv", quote=F, row.names = F)#nSNPs
write.csv(vetted$test == "vim.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/vim.prec.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "vim.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/vim.temp.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "vim.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/vim.rand.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "vim.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/vim.corr.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "cors.hum" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/cors.hum.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "cors.prec" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/cors.prec.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "cors.temp" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/cors.temp.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "cors.rand" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/cors.rand.vetted.snps.csv", quote=F, row.names = F) #nSNPs
write.csv(vetted$test == "cors.corr" & vetted$pass.direction.test == TRUE & vetted$pass.strength.test == TRUE, "~/Downloads/vet.output/cors.corr.vetted.snps.csv", quote=F, row.names = F) #nSNPs
plot 100 randomly selected points to assure visually that this vetting scheme is doing a satisfactory job in identifying SNPs where the putatively significant relationship identified by the GEA is not supported by the semi-independent phase 2 sampling (whether the given SNP passed vetting printed as the x axis label)
#plot 100 random points to visually confirm that the classification worked
set.seed(62399)
plot.samples<-vetted[sample(nrow(vetted),100),]
par(mfrow=c(2,2))
for (i in paste(plot.samples$chrom,plot.samples$pos,plot.samples$test)){
plot(as.numeric(as.vector(full[paste(full$chrom,full$POS,full$test) == i,c(17:30)])),
as.numeric(as.vector(full[paste(full$chrom,full$POS,full$test) == i,c(3:16)])), main = i, ylim =c(0,1), xlab=plot.samples$vet[paste(plot.samples$chrom,plot.samples$pos,plot.samples$test) == i],ylab="allele frequency")
abline(lm(as.numeric(as.vector(full[paste(full$chrom,full$POS,full$test) == i,c(3:16)]))~
as.numeric(as.vector(full[paste(full$chrom,full$POS,full$test) == i,c(17:30)]))))
points(as.numeric(as.vector(full.2[paste(full.2$chrom,full.2$POS,full.2$test) == i,c(24:44)])),
as.numeric(as.vector(full.2[paste(full.2$chrom,full.2$POS,full.2$test) == i,c(3:23)])), col="red")
abline(lm(as.numeric(as.vector(full.2[paste(full.2$chrom,full.2$POS,full.2$test) == i,c(3:23)]))~
as.numeric(as.vector(full.2[paste(full.2$chrom,full.2$POS,full.2$test) == i,c(24:44)]))), col="red")
}
