#Code for adult gonads 
D = read.csv(file = "FILE NAME")

# load necessary libraries

library(car)
library(compute.es)
library(effects)
library(ggplot2)
library(multcomp)
library(grid)
library(gridExtra)
###########################################################################



D$temptrt = interaction(D$Temp,D$Treatment) # combine vars into
#a single interaction var

#account for pseudoreplication 

D = aggregate(D,by=list(D$temptrt,D$Rep),data=D,FUN=mean)
D$temptrt <- D$Group.1

#two way anova with tukey's
m1 = aov(HSD17 ~ temptrt, data = D)
summary.aov(m1)
summary(glht(m1,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(HSD17~temptrt, D, mean)

### Replace gene response variable (HSD17) with others (e.g., GPR30) ###

m2 = aov(HSD3 ~ temptrt, data = D)
summary.aov(m2)
summary(glht(m2,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(HSD3~temptrt, D, mean)

m3 = aov(Arx ~ temptrt, data = D)
summary.aov(m3)
summary(glht(m3,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(Arx ~ temptrt, D, mean)

m4 = aov(CYP19a ~ temptrt, data = D)
summary.aov(m4)
summary(glht(m4,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(CYP19a~temptrt, D, mean)

m5 = aov(ESR1 ~ temptrt, data = D)
summary.aov(m5)
summary(glht(m5,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(ESR1~temptrt, D, mean)

m6 = aov(ESR2 ~ temptrt, data = D)
summary.aov(m6)
summary(glht(m6,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(ESR2~temptrt, D, mean)

m7 = aov(ESR3 ~ temptrt, data = D)
summary.aov(m7)
summary(glht(m7,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(ESR3~temptrt, D, mean)

m8 = aov(FSHR ~ temptrt, data = D)
summary.aov(m8)
summary(glht(m8,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(FSHR~temptrt, D, mean)

m9 = aov(HSP90 ~ temptrt, data = D)
summary.aov(m9)
summary(glht(m9,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr')) 

aggregate(HSP90~temptrt, D, mean)

m10 = aov(INHA ~ temptrt, data = D)
summary.aov(m10)
summary(glht(m10,linfct=mcp(temptrt='Tukey')),adjusted(type='fdr'))

aggregate(INHA ~ temptrt, D, mean)

#end code for adult gonads

#code for larval samples

D = read.csv(file = "File name")

# loads necessary libraries
library(car)
library(compute.es)
library(effects)
library(ggplot2)
library(multcomp)
library(grid)
library(gridExtra)

D$temptrt1 = interaction(D$Treatment, D$Temp) # combine two vars into
#a single interaction var

# account for psuedoreplication
D = aggregate(D,by=list(D$temptrt1,D$Rep),data=D,FUN=mean)
D$temptrt1 <- D$Group.1

#ANOVAs
m1 = aov(HSD17~temptrt1, data=D)
summary.aov(m1)
summary(glht(m1,linfct=mcp(temptrt1='Tukey'))) 

#averages 
aggregate(HSD17~temptrt1, D, mean)

m2 = aov(HSD3~temptrt1, data=D)
summary.aov(m2)
summary(glht(m2,linfct=mcp(temptrt1='Tukey'))) 

aggregate(HSD3~temptrt1, D, mean)

m3 = aov(GPR30~temptrt1, data=D)
summary.aov(m3)
summary(glht(m3,linfct=mcp(temptrt1='Tukey'))) 

aggregate(GPR30~temptrt1, D, mean)

m4 = aov(ESR1~temptrt1, data=D)
summary.aov(m4)
summary(glht(m4,linfct=mcp(temptrt1='Tukey')))

aggregate(ESR1~temptrt1, D, mean)

m5 = aov(ESR2~temptrt1, data=D)
summary.aov(m5)
summary(glht(m5,linfct=mcp(temptrt1='Tukey')))

aggregate(ESR2~temptrt1, D, mean)

m6 = aov(ESR3~temptrt1, data=D)
summary.aov(m6)
summary(glht(m6,linfct=mcp(temptrt1='Tukey')))

aggregate(ESR3~temptrt1, D, mean)

m7 = aov(GnRHR~temptrt1, data=D)
summary.aov(m7)
summary(glht(m7,linfct=mcp(temptrt1='Tukey')))

aggregate(GnRHR~temptrt1, D, mean)

m8 = aov(IGF2~temptrt1, data=D)
summary.aov(m8)
summary(glht(m8,linfct=mcp(temptrt1='Tukey')))

aggregate(IGF2~temptrt1, D, mean)

m9 = aov(CYP19b~temptrt1, data=D)
summary.aov(m9)
summary(glht(m9,linfct=mcp(temptrt1='Tukey')))

aggregate(CYP19b~temptrt1, D, mean)

m10 = aov(FSHR~temptrt1, data=D)
summary.aov(m10)
summary(glht(m10,linfct=mcp(temptrt1='Tukey')))

aggregate(FSHR~temptrt1, D, mean)