#Genetic distances between all individuals across all loci
#2/26/18
#updated 2/28/18
#updated 3/30/18
#set working directory
setwd("/uufs/chpc.utah.edu/common/home/u6000989/data/grasses/variants/")
#read genotype file, it should have the dimension of 21581 rows for SNPs and 617 columns for individuals
genotype.matrix<-read.table("PSSP_BayGen_metadata")
#genotype file has treatment and quad ids at the end, want to get rid of last two rows
#transpose genotypes to have individuals as rows and SNPs as columns
genotypes<-genotype.matrix[1:21581,1:617]
#I had to write a file and then read it to make everything numeric
g<-read.table("PSSP_gen.txt",sep=",",header=F)
write.table(All,file="PSSP_gen.txt",row.names=F,col.names=F,sep=",")
setwd("/uufs/chpc.utah.edu/common/home/u1065430/Data/")
#read in genotype file
genotypes<-read.table("PSSP_gen.txt",header=F,sep=",")
#transpose data to have individuals as rows and get rid of removals
genotypes2<-t(genotypes[-c(232:307)])
#function to calculate the genetic distance
distance<-function(x,y){
output<-mean(abs(x-y))/2
return(output)
}
dist.matrix<-matrix(NA, dim(genotypes2)[1],dim(genotypes2)[1])
for(i in 1:nrow(genotypes2)){
for(j in 1:nrow(genotypes2)){
dist.matrix[i,j]<-distance(genotypes2[i,],genotypes2[j,])
}
}
write.table(dist.matrix,file="PSSP_GenDistMatrix.txt",row.names=F,col.names=F,sep=",")
#########################################################
#Ploting densities
# import distance matrix
D_mat <-read.table("PSSP_GenDistMatrix.txt",header=F,sep=",")
setwd("/uufs/chpc.utah.edu/common/home/u1065430/variants/")
# import treatment data
IDs <- read.csv("pops_split_trt.txt",header=F)
names(IDs)[1] <- "fullname"
IDs$fullname <- as.character((IDs$fullname))
IDs$trt <- sapply(strsplit(IDs$fullname," "),"[[",2)
IDs <- subset(IDs,trt !="GrassRemoval" & trt!="ShrubRemoval")
# create treatment matrix
trt_mat <- matrix("c",NROW(D_mat),NCOL(D_mat))
for(i in 1:NROW(trt_mat)){
for(j in 1:NCOL(trt_mat)){
trt_mat[i,j] <- paste0(IDs$trt[i],IDs$trt[j])
}
}
# get rid of diagonals and lower.tri in matrices
D_mat[lower.tri(D_mat,diag=T)] <- NA
trt_mat[lower.tri(trt_mat,diag=T)] <- "duplicate"
# summarize treatment comparison labels
table(trt_mat)
# Looks like we have ControlDrought but no DroughtControl, due to the ordering
# of the plots by treatment. This makes things a little easier.
# store treatment comparisons
CC <- D_mat[trt_mat=="ControlControl"]
CD <- D_mat[trt_mat=="ControlDrought"]
CI <- D_mat[trt_mat=="ControlIrrigation"]
DI <- D_mat[trt_mat=="DroughtIrrigation"]
DD <- D_mat[trt_mat=="DroughtDrought"]
II <- D_mat[trt_mat=="IrrigationIrrigation"]
# store elevation comparisons
LeLe <- D_mat[trt_mat=="LELE"]
HeHe <- D_mat[trt_mat=="HEHE"]
CLe <- D_mat[trt_mat=="ControlLE"]
CHe <- D_mat[trt_mat=="ControlHE"]
HeLe <- D_mat[trt_mat=="HELE"]
# plot figure--experiment
#put plots on the same page
par(mfrow=c(1,2))
plot(density(CC),col="black",ylim=c(0,10000),main="Treatment")
lines(density(CD),col="red",lty="dashed")
lines(density(CI),col="blue",lty="dashed")
lines(density(DI),col="magenta",lty="dashed")
lines(density(DD),col="red",lty="solid")
lines(density(II),col="blue",lty="solid")
legend("topleft",c("C-C","C-D","C-I","D-I","D-D","I-I"),
col=c("black","red","blue","magenta","red","blue"),
lty=c("solid","dashed","dashed","dashed","solid","solid"),bty="n")
# plot figure--elevation
plot(density(CC),col="black",ylim=c(0,30000),main="Elevation")
lines(density(CLe),col="red",lty="dashed")table(trt_mat)
lines(density(CHe),col="blue",lty="dashed")
lines(density(HeLe),col="green",lty="dashed")
lines(density(LeLe),col="red",lty="solid")
lines(density(HeHe),col="blue",lty="solid")
legend("topleft",c("C-C","C-Low","C-High","High-Low","Low-Low","High-High"),
col=c("black","red","blue","green","red","blue"),
trt.summary
Mean Sd Median X95..Quantile
CC 0.07626518 0.02136499 0.07914129 0.10349169
CD 0.07656189 0.02244092 0.08051883 0.10447506
CI 0.07320773 0.02184020 0.07670120 0.10189244
DI 0.07118397 0.02156030 0.07626953 0.09757835
DD 0.07438466 0.02275348 0.07938610 0.10341004
II 0.06714136 0.01957670 0.07342020 0.08821582
elev.summary
Mean Sd Median X95..Quantile
CC 0.07626518 0.021364994 0.07914129 0.10349169
CLe 0.07908616 0.016216632 0.08090534 0.10034185
CHe 0.08005260 0.009594479 0.08041281 0.09162712
HeLe 0.08029773 0.007821255 0.07934597 0.09119724
LeLe 0.07960833 0.013049962 0.07992496 0.09734773
HeHe 0.07538988 0.004537135 0.07473375 0.08384693
lty=c("solid","dashed","dashed","dashed","solid","solid"),bty="n")