setwd("/Users/joshuamiller/Documents/My\ Documents/Horns/repeatABLE")
library(GenABEL)
library(RepeatABEL)

convert.snp.tped(tpedfile = "Animals_with_phenos_transposed.tped",tfamfile = "Animals_with_phenos_transposed.tfam","RM_Animals_with_Phenos.raw")
##Make a GenABLE file, note that to do this you have to have a phenotype file with only 1 measure per indivudal
base_genable<-load.gwaa.data(phenofile = "Base_Circ_Typed_IDs_Only_No_Dupes.txt", genofile = "RM_Animals_with_Phenos.raw", force = T)
st_Base<-read.table("Base_Circ_Typed_IDs_Only.txt",header = T)
st_Mass<-read.table("Mass_Typed_IDs_Only.txt", header = T)
st_Length<-read.table("Length_Typed_IDs_Only.txt", header = T)

base_fixed=stbase ~ Age*sex*poly(Day,2)
mass_fixed=stMass ~ Age*sex*poly(Day,2)
length_fixed=stlength ~ Age*sex*poly(Day,2)

##Repeated measures Base circ
##Repeated measure GWAS, using the GenAble file above and a different phenotype file that has repeated measures
Base_Mod1 <- preFitModel(base_fixed, random=~1|id + 1|Year + 1|Yrbirth,
                         genabel.data = base_genable,
                         phenotype.data = st_Base,
                         corStruc=list( id=list("GRM","Ind"), Year=list("Ind"), Yrbirth=list("Ind") ))
Base_GWAS1 <- rGLS(base_fixed, genabel.data = base_genable,
                   phenotype.data = st_Base, V = Base_Mod1$V)
summary(Base_GWAS1)
plot(Base_GWAS1)
##Calculate Genomic inflation factor
estlambda(Base_GWAS1[, "P1df"], plot=TRUE)
text(11, 0, pos= 4,expression(paste(lambda,"=0.91",phantom(.)%+-%phantom(.),"0.0005")))
write.table(Base_GWAS1@results,file = "Base_GWAS.txt", sep = " ")

##Repeated measures Mass
Mass_Mod1 <- preFitModel(mass_fixed, random=~1|id + 1|Year + 1|Yrbirth,
                         genabel.data = base_genable,
                         phenotype.data = st_Mass,
                         corStruc=list( id=list("GRM","Ind"), Year=list("Ind"), Yrbirth=list("Ind") ))
Mass_GWAS1 <- rGLS(mass_fixed, genabel.data = base_genable,
                   phenotype.data = st_Mass, V = Mass_Mod1$V)
summary(Mass_GWAS1)
plot(Mass_GWAS1)
##Calculate Genomic inflation factor
estlambda(Mass_GWAS1[, "P1df"], plot=TRUE)
text(11, 0, pos= 4,expression(paste(lambda,"=1.01",phantom(.)%+-%phantom(.),"0.0004")))
write.table(Mass_GWAS1@results,file = "Mass_GWAS.txt", sep = " ")

##Repeated measures Length
Length_Mod1 <- preFitModel(length_fixed, random=~1|id + 1|Year + 1|Yrbirth,
                         genabel.data = base_genable,
                         phenotype.data = st_Length,
                         corStruc=list( id=list("GRM","Ind"), Year=list("Ind"), Yrbirth=list("Ind") ))
Length_GWAS1 <- rGLS(length_fixed, genabel.data = base_genable,
                   phenotype.data = st_Length, V = Length_Mod1$V)
summary(Length_GWAS1)
plot(Length_GWAS1)
##Calculate Genomic inflation factor
estlambda(Length_GWAS1[, "P1df"], plot=TRUE)
text(11, 0, pos= 4,expression(paste(lambda,"=0.98",phantom(.)%+-%phantom(.),"0.0007")))
length_qt<-qtscore(length_fixed,data = Length_GWAS1, times = 200)
write.table(Length_GWAS1@results,file = "Length_GWAS.txt", sep = " ")

##Manhattan Plots
library(qqman)
R_manhattan<-read.table("RepeatABLE_GWAS_p-values.txt",header = T)
layout(matrix(c(1,1,2,3,3,4,5,5,6), nrow=3,ncol=3,byrow=T))

manhattan(R_manhattan,chr = "CHR", bp="BP",p="Length_p.values",snp = "SNP", 
          suggestiveline = -log10(1/3777), genomewideline = -log10(0.5/3777),ylim = c(0, 5),main = "Horn Length")
estlambda(Length_GWAS1[, "P1df"], plot=TRUE)
text(11, 0, pos= 4,expression(paste(lambda,"=0.98",phantom(.)%+-%phantom(.),"0.0007")))

manhattan(R_manhattan,chr = "CHR", bp="BP",p="Base_p.values",snp = "SNP", 
          suggestiveline = -log10(1/3777), genomewideline = -log10(0.5/3777),ylim = c(0, 5),main = "Horn Base Circumference")
estlambda(Base_GWAS1[, "P1df"], plot=TRUE)
text(11, 0, pos= 4,expression(paste(lambda,"=0.91",phantom(.)%+-%phantom(.),"0.0005")))

manhattan(R_manhattan,chr = "CHR", bp="BP",p="Mass_p.values",snp = "SNP", 
          suggestiveline = -log10(1/3777), genomewideline = -log10(0.5/3777),ylim = c(0, 5),main = "Body Mass")
estlambda(Mass_GWAS1[, "P1df"], plot=TRUE)
text(11, 0, pos= 4,expression(paste(lambda,"=1.01",phantom(.)%+-%phantom(.),"0.0004")))