##################### LGN0313 - Melhoramento Genético #####################

##################### Aula Prática 12 - Análises do R #####################

################################
#carregando arquivos da Turma 1
################################

#Lembre-se de antes de tudo informar qual a pasta do seu computador que o R buscará seus arquivos.
#Para isso aperte CTRL+SHIFT+H e escolha a pasta em que os .txt com os dados estão.

turma1<-read.table("T1.txt", h=T)
turma1$Bloco<-as.factor(turma1$Bloco) #informar que a coluna Bloco do objeto "turma1" deve ser tratada como fator
str(turma1)

###############################################
# anova e herdabilidade - NF e AP
###############################################

#agrupar os gráficos que serão gerados
par(mfrow=c(1,2))

# NF
tapply(turma1$NF, turma1[,1:2], mean) # genótipos por bloco
interaction.plot(turma1$Bloco, turma1$Hibrido, turma1$NF, legend = T, trace.label ="Hibrido", ylab="NF T1", xlab = "Bloco", col = 2)

anova.t1.NF<-aov(NF~Bloco+Hibrido, data = turma1)
anova(anova.t1.NF)

(R<-max(as.numeric(turma1$Bloco))) # Buscar número de repetições
(Vg.t1.NF<-(anova(anova.t1.NF)[[3]][2]-anova(anova.t1.NF)[[3]][3])/R)
(Ve.t1.NF<-anova(anova.t1.NF)[[3]][3])
(h.t1.NF<-Vg.t1.NF/(Vg.t1.NF+Ve.t1.NF/R)) # herdabilidade
(CV.t1.NF<-sqrt(anova(anova.t1.NF)[[3]][3])/mean(turma1$NF)*100) #CV%

# AP
tapply(turma1$AP, turma1[,1:2], mean) # genótipos por bloco
interaction.plot(turma1$Bloco, turma1$Hibrido, turma1$AP, legend = T, trace.label ="Hibrido", ylab="AP T1", xlab = "Bloco", col = 3)

anova.t1.AP<-aov(AP~Bloco+Hibrido, data = turma1)
anova(anova.t1.AP)

(Vg.t1.AP<-(anova(anova.t1.AP)[[3]][2]-anova(anova.t1.AP)[[3]][3])/R)
(Ve.t1.AP<-anova(anova.t1.AP)[[3]][3])
(h.t1.AP<-Vg.t1.AP/(Vg.t1.AP+Ve.t1.AP/R))
(CV.t1.AP<-sqrt(anova(anova.t1.AP)[[3]][3])/mean(turma1$AP)*100) #CV%

##### Pacote para construção de gráficos #### 
install.packages("ggplot2") # Instalar o pacote, se você já intalou antes não precisa reinstalar
library("ggplot2") #"ativando" o pacote
attach(turma1)
ggplot() +
  aes(x = Bloco, y = NF, color = Hibrido) +
  geom_line(aes(group = Hibrido)) +
  geom_point()
##############################################

########################
# teste de médias
########################

#agrupar os gráficos que serão gerados
par(mfrow=c(1,2))

tukey.t1.NF<-TukeyHSD(x = anova.t1.NF, which = "Hibrido", ordered = T, conf.level = 0.95)
tukey.t1.NF
plot(tukey.t1.NF)

tukey.t1.AP<-TukeyHSD(x = anova.t1.AP, which = "Hibrido", ordered = T, conf.level = 0.95)
tukey.t1.AP
plot(tukey.t1.AP)

########## Outra forma de realizar TUKEY ###########
app<-lm(AP~Bloco+Hibrido,turma1)
anova(app)
install.packages("emmeans") #mesma coisa, se já instalou este pacote não precisa reinstalar
library(emmeans)
(fm_1<-emmeans(app,"Hibrido"))
plot(fm_1)
cld(fm_1,adjust="tukey",Letters=letters)
##################################

##################
#Ganho de seleção
##################

# NF 
(Xo.t1.NF<-mean(turma1$NF)) # média populacional
(Xs.t1.NF<-mean(sort(tapply(turma1$NF, turma1$Hibrido, mean), decreasing = T)[1:2])) #média dos selecionados
(DS.t1.NF<-Xs.t1.NF-Xo.t1.NF) # diferencial de seleção
(GS.t1.NF<-DS.t1.NF*h.t1.NF) #ganho de seleção
(GS.t1.NF.per<-GS.t1.NF/Xo.t1.NF*100) #ganho de seleção em percentagem

# AP
(Xo.t1.AP<-mean(turma1$AP)) # média populacional
(Xs.t1.AP<-mean(sort(tapply(turma1$AP, turma1$Hibrido, mean))[1:2])) #média dos selecionados
(DS.t1.AP<-Xs.t1.AP-Xo.t1.AP) # diferencial de seleção
(GS.t1.AP<-DS.t1.AP*h.t1.AP) #ganho de seleção
(GS.t1.AP.per<-GS.t1.AP/Xo.t1.AP*100) #ganho de seleção em percentagem

#################
# Interação GxA
#################

#carregando arquivo da outra turma
turma2<-read.table("T2.txt", h=T)
turma2$Bloco<-as.factor(turma2$Bloco)
str(turma2)

datagxa<-rbind(turma1, turma2) #unindo os dados das duas turmas
datagxa$Turma<-as.factor(c(rep(1,20), rep(2, 20))) #criando o fator turma (ambiente)
str(datagxa)
head(datagxa); tail(datagxa) #Só para ver o começo e o fim da nova tabela de dados

#agrupar os gráficos que serão gerados
par(mfrow=c(1,2))

# Anova conjunta com efeito de bloco dentro de turma

#NF
tapply(datagxa$NF, datagxa[,c(1,5)], mean) # genótipos por turma
gxa.NF<-aov(NF~Turma:Bloco+Turma+Hibrido+Turma*Hibrido, data=datagxa)
anova(gxa.NF)  
interaction.plot(datagxa$Turma, datagxa$Hibrido, datagxa$NF, legend = T, trace.label ="Hibrido", ylab="NF", xlab = "Turma", col = 10)

#AP
tapply(datagxa$AP, datagxa[,c(1,5)], mean) # genótipos por turma
gxa.AP<-aov(AP~Turma:Bloco+Turma+Hibrido+Turma*Hibrido, data=datagxa)
anova(gxa.AP)
interaction.plot(datagxa$Turma, datagxa$Hibrido, datagxa$AP, legend = T, trace.label ="Hibrido", ylab="AP", xlab = "Turma", col = 11)

###### Estimativa dos componentes de variância #####
(A<-max(as.numeric(datagxa$Turma)))
(Vg.conj.AP<-(anova(gxa.AP)[[3]][2]-anova(gxa.AP)[[3]][5])/(R*A))
(Vha.conj.AP<-(anova(gxa.AP)[[3]][4]-anova(gxa.AP)[[3]][5])/R)
(Ve.conj.AP<-anova(gxa.AP)[[3]][5])
(h.conj.AP<-Vg.conj.AP/(Vg.conj.AP+(Vha.conj.AP/A)+Ve.conj.AP/(R*A)))
###################################################

######################
# correlação genética 
######################

turma1$APNF<-as.numeric(turma1$AP+turma1$NF) #criando variável acessória
head(turma1)

t1.APNF<-aov(APNF~Bloco+Hibrido, data = turma1)
anova(t1.APNF)

#covariância fenotípica
(COVF.t1<-(anova(t1.APNF)[[3]][2]-anova(anova.t1.NF)[[3]][2]-anova(anova.t1.AP)[[3]][2])/2)

#covariância ambiental
(COVe.t1<-(anova(t1.APNF)[[3]][3]-anova(anova.t1.NF)[[3]][3]-anova(anova.t1.AP)[[3]][3])/2)

#covariância genética
(COVg.t1<-(COVF.t1-COVe.t1)/R)

#correlação genética
(corg.t1<-COVg.t1/sqrt(Vg.t1.NF*Vg.t1.AP))


#######################################
# ganho indireto e direto com a seleção
#######################################

# função para estimar o i da intensidade de seleção 
i<-function(persel){
  estimate<-dnorm(qnorm(1-(persel/100), 0, 1), 0, 1)/(persel/100)
  return(estimate)
}

(isel<-i(40)) # selecionando os dois melhores dos cinco híbridos (40% de 5)

# Ganho indireto de seleção em selecionando NF e esperando ganho em AP
(GSi.t1_NF_AP<-isel*corg.t1*sqrt(h.t1.NF)*sqrt(Vg.t1.AP))

# Ganho indireto de seleção em selecionando AP e esperando ganho em NF
(GSi.t1_AP_NF<-isel*corg.t1*sqrt(h.t1.AP)*sqrt(Vg.t1.NF))

# Ganhos diretos de seleção
(GSd.t1.NF<-isel*1.00*sqrt(h.t1.NF)*sqrt(Vg.t1.NF))
(GSd.t1.AP<-isel*1.00*sqrt(h.t1.AP)*sqrt(Vg.t1.AP))

##################################### FIM ##################################