#R code to calculate heterozygosity across all quads
#2/5/18
#Read in all text files of allele frequencines that includes the mean, median, CI LB, and CI UP as the rows
#Set working directory to Zach's sub-directory freqs
setwd("/uufs/chpc.utah.edu/common/home/6000989/data/grasses/freqs/")
X1.all<-read.table("p_o_PSSP_X1.txt",sep=",")
X2.all<-read.table("p_o_PSSP_X2.txt",sep=",")
X5.all<-read.table("p_o_PSSP_X5.txt",sep=",")
X6.all<-read.table("p_o_PSSP_X6.txt",sep=",")
X7.all<-read.table("p_o_PSSP_X7.txt",sep=",")
X8.all<-read.table("p_o_PSSP_X8.txt",sep=",")
X11.all<-read.table("p_o_PSSP_X11.txt",sep=",")
X12.all<-read.table("p_o_PSSP_X12.txt",sep=",")
X13.all<-read.table("p_o_PSSP_X13.txt",sep=",")
X14.all<-read.table("p_o_PSSP_X14.txt",sep=",")
X15.all<-read.table("p_o_PSSP_X15.txt",sep=",")
X17.all<-read.table("p_o_PSSP_X17.txt",sep=",")
X18.all<-read.table("p_o_PSSP_X18.txt",sep=",")
X19.all<-read.table("p_o_PSSP_X19.txt",sep=",")
X20.all<-read.table("p_o_PSSP_X20.txt",sep=",")
LE.all<-read.table("p_o_PSSP_LE.txt",sep=",")
HE.all<-read.table("p_o_PSSP_HE.txt",sep=",")
E1Q1.all<-read.table("p_o_PSSP_E1Q1.txt",sep=",")
E1Q2.all<-read.table("p_o_PSSP_E1Q2.txt",sep=",")
E1Q3.all<-read.table("p_o_PSSP_E1Q3.txt",sep=",")
E1Q4.all<-read.table("p_o_PSSP_E1Q4.txt",sep=",")
E1Q5.all<-read.table("p_o_PSSP_E1Q5.txt",sep=",")
E1Q6.all<-read.table("p_o_PSSP_E1Q6.txt",sep=",")
P10E1Q1.all<-read.table("p_o_PSSP_P10E1Q1.txt",sep=",")
P10E1Q2.all<-read.table("p_o_PSSP_P10E1Q2.txt",sep=",")
P10E1Q3.all<-read.table("p_o_PSSP_P10E1Q3.txt",sep=",")
P10E1Q4.all<-read.table("p_o_PSSP_P10E1Q4.txt",sep=",")
P7E1Q1.all<-read.table("p_o_PSSP_P7E1Q1.txt",sep=",")
P7E1Q2.all<-read.table("p_o_PSSP_P7E1Q2.txt",sep=",")
P7E1Q3.all<-read.table("p_o_PSSP_P7E1Q3.txt",sep=",")
P7E1Q4.all<-read.table("p_o_PSSP_P7E1Q4.txt",sep=",")
#Get rid of median and confidence interavals
X1<-t(X1.all[,1])
X2<-t(X2.all[,1])
X5<-t(X5.all[,1])
X6<-t(X6.all[,1])
X7<-t(X7.all[,1])
X8<-t(X8.all[,1])
X11<-t(X11.all[,1])
X12<-t(X12.all[,1])
X13<-t(X13.all[,1])
X14<-t(X14.all[,1])
X17<-t(X17.all[,1])
X18<-t(X18.all[,1])
X19<-t(X19.all[,1])
X20<-t(X20.all[,1])
LE<-t(LE.all[,1])
HE<-t(HE.all[,1])
E1Q1<-t(E1Q1.all[,1])
E1Q2<-t(E1Q2.all[,1])
E1Q3<-t(E1Q3.all[,1])
E1Q4<-t(E1Q4.all[,1])
E1Q5<-t(E1Q5.all[,1])
E1Q6<-t(E1Q6.all[,1])
P10E1Q1<-t(P10E1Q1.all[,1])
P10E1Q2<-t(P10E1Q2.all[,1])
P10E1Q3<-t(P10E1Q3.all[,1])
P10E1Q4<-t(P10E1Q4.all[,1])
P7E1Q1<-t(P7E1Q1.all[,1])
P7E1Q2<-t(P7E1Q2.all[,1])
P7E1Q3<-t(P7E1Q3.all[,1])
P7E1Q4<-t(P7E1Q4.all[,1])
#combine quds based on treatment
Drought<-rbind(X2,X5,X8,X11,X13,X17,X20)
Irrigation<-rbind(X1,X6,X7,X12,X14,X18,X19)
Natural<-rbind(LE,HE)
BigExclosure<-rbind(E1Q1,E1Q2,E1Q3,E1Q4,E1Q5,E1Q6)
P10Exclosure<-rbind(P10E1Q1,P10E1Q2,P10E1Q3,P10E1Q4)
P7Exclosure<-rbind(P7E1Q1,P7E1Q2,P7E1Q3,P7E1Q4)
All<-rbind(Drought,Irrigation,Natural,BigExclosure,P10Exclosure,P7Exclosure)
#Function to calculate heterozygosity
He<-function(p){
He<-2*p*(1-p)
return(He)
}
#Calcualte heterozyogisty
Drought.He<-He(Drought)
Irrigation.He<-He(Irrigation)
Natural.He<-He(Natural)
BigExclosure.He<-He(BigExclosure)
P10Exclosure.He<-He(P10Exclosure)
P7Exclosure.He<-He(P7Exclosure)
All.He<-He(All)
#write a text file that has heterozygote values to Jacque's directory
setwd("/uufs/chpc.utah.edu/common/home/u1065430/")
write.table(All.He,file="PSSP_Heterozygosity.txt",row.names=F,col.names=F,sep=",")
#Heterozygosity across all quads
X1.He<-He(X1)
X2.He<-He(X2)
X5.He<-He(X5)
X6.He<-He(X6)
X7.He<-He(X7)
X8.He<-He(X8)
X11.He<-He(X11)
X12.He<-He(X12)
X13.He<-He(X13)
X14.He<-He(X14)
X17.He<-He(X17)
X18.He<-He(X18)
X19.He<-He(X19)
X20.He<-He(X20)
E1Q1.He<-He(E1Q1)
E1Q2.He<-He(E1Q2)
E1Q3.He<-He(E1Q3)
E1Q4.He<-He(E1Q4)
E1Q5.He<-He(E1Q5)
E1Q6.He<-He(E1Q6)
P10E1Q1.He<-He(P10E1Q1)
P10E1Q2.He<-He(P10E1Q2)
P10E1Q3.He<-He(P10E1Q3)
P10E1Q4.He<-He(P10E1Q4)
P7E1Q1.He<-He(P7E1Q1)
P7E1Q2.He<-He(P7E1Q2)
P7E1Q3.He<-He(P7E1Q3)
P7E1Q4.He<-He(P7E1Q4)
X1.HeMean<-as.matrix(mean(X1.He))
X2.HeMean<-as.matrix(mean(X2.He))
X5.HeMean<-as.matrix(mean(X5.He))
X6.HeMean<-as.matrix(mean(X6.He))
X7.HeMean<-as.matrix(mean(X7.He))
X8.HeMean<-as.matrix(mean(X8.He))
X11.HeMean<-as.matrix(mean(X11.He))
X12.HeMean<-as.matrix(mean(X12.He))
X13.HeMean<-as.matrix(mean(X13.He))
X14.HeMean<-as.matrix(mean(X14.He))
X17.HeMean<-as.matrix(mean(X17.He))
X18.HeMean<-as.matrix(mean(X18.He))
X19.HeMean<-as.matrix(mean(X19.He))
X20.HeMean<-as.matrix(mean(X20.He))
E1Q1.HeMean<-as.matrix(mean(E1Q1.He))
E1Q2.HeMean<-as.matrix(mean(E1Q2.He))
E1Q3.HeMean<-as.matrix(mean(E1Q3.He))
E1Q4.HeMean<-as.matrix(mean(E1Q4.He))
E1Q5.HeMean<-as.matrix(mean(E1Q5.He))
E1Q6.HeMean<-as.matrix(mean(E1Q6.He))
P10E1Q1.HeMean<-as.matrix(mean(P10E1Q1.He))
P10E1Q2.HeMean<-as.matrix(mean(P10E1Q2.He))
P10E1Q3.HeMean<-as.matrix(mean(P10E1Q3.He))
P10E1Q4.HeMean<-as.matrix(mean(P10E1Q4.He))
P7E1Q1.HeMean<-as.matrix(mean(P7E1Q1.He))
P7E1Q2.HeMean<-as.matrix(mean(P7E1Q2.He))
P7E1Q3.HeMean<-as.matrix(mean(P7E1Q3.He))
P7E1Q4.HeMean<-as.matrix(mean(P7E1Q4.He))
Drought.HeMean<-rbind(X2.HeMean,X5.HeMean,X8.HeMean,X11.HeMean,X13.HeMean,X17.He.Mean,X20.HeMean)
Irrigation.HeMean<-rbind(X1.HeMean,X6.HeMean,X7.HeMean,X12.HeMean,X14.HeMean,X18.HeMean,X19.HeMean)
Natural.HeMean<-rbind(LE.HeMean,HE.HeMean)
BigExclosure.HeMean<-rbind(E1Q1.HeMean,E1Q2.HeMean,E1Q3.HeMean,E1Q4.HeMean,E1Q5.HeMean,E1Q6.HeMean)
P10Exclosure.HeMean<-rbind(P10E1Q1.HeMean,P10E1Q2.HeMean,P10E1Q3.HeMean,P10E1Q4.HeMean)
P7Exclosure.HeMean<-rbind(P7E1Q1.HeMean,P7E1Q2.HeMean,P7E1Q3.HeMean,P7E1Q4.HeMean)
All.HeMean<-rbind(Drought.HeMean,Irrigation.HeMean,Natural.HeMean,BigExclosure.HeMean,P10Exclosure.HeMean,P7Exclosure.HeMean)
Treatment.HeMean<-rbind(Drought.HeMean,Irrigation.HeMean)
write.table(Treatment.HeMean,file="PSSP_Treatment_HeMean.txt",row.names=F,col.names=F,sep=",")
write.table(Natural.HeMean,file="PSSP_Natural_HeMean.txt",row.names=F,col.names=F,sep=",")
write.table(BigExclosure.HeMean,file="PSSP_BigExclosure_HeMean.txt",row.names=F,col.names=F,sep=",")
write.table(P10Exclosure.HeMean,file="PSSP_P10Exlosure_HeMean.txt",row.names=F,col.names=F,sep=",")
write.table(P7Exclosure.HeMean,file="PSSP_P7Exclosure_HeMean.txt",row.names=F,col.names=F,sep=",")