Heterozygosity

#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=",")