#=======================================================
#Pacotes usados. Instalacao
#=======================================================
#install.packages("adegenet")
#install.packages("poppr")
#install.packages("dplyr")
#install.packages("hierfstat")
#install.packages("reshape2")
#install.packages("scales")
#install.packages("ape")
#if (!require("BiocManager", quietly = TRUE))
#install.packages("BiocManager")
#BiocManager::install("ggtree")
#install.packages("ggplot2")
#install.packages("ggrepel")
#install.packages("RColorBrewer")
#install.packages("gridExtra")
#=======================================================
#Pacotes usados. Carregamento
#=======================================================
library(adegenet)
library(poppr)
library(dplyr)
library(hierfstat)
library(reshape2)
library(scales)
library(ape)
library(ggtree)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)
library(gridExtra)
#=======================================================
#PASSO 1.Importar arquivo GENPOP
#=======================================================
ESPECIE_gen = import2genind("RadSeq_ESPECIE.gen", ncode = 3, quiet = TRUE)
set.seed(1)
#=======================================================
#PASSO 2.Filtrar os dados: Missing Data - MD
#=======================================================
#Calcular Numero de LOCI
length(ESPECIE_gen$loc.n.all)
#Loci com menos de 75% de MD
filtro1<-missingno(ESPECIE_gen, type="loci", cutoff = 0.75)
#Individuos com menos de 75% de MD
filtro2<-missingno(filtro1, type="geno", cutoff = 0.75)
ESPECIE_gen<-filtro2
rm(filtro1,filtro2)
#Remover individuos/amostras duplicados
dups = mlg.id(ESPECIE_gen)
for (i in dups){ # for each element in the list object
  if (length(dups[i]) > 1){ # if the length is greater than 1
    print(i) # print individuals that are duplicates
  }
}
#Verificar se os loci filtrados sao polimorificos.
isPoly(ESPECIE_gen) %>% summary
#=======================================================
#PASSO 3.Gerar estatisticas descritvas das populacoes
#=======================================================
#Verificar as populacoes
popNames(ESPECIE_gen)
pop(ESPECIE_gen)

#Gerar tabela com as estatisticas
tabela_pop<-bind_cols(as.data.frame(basic.stats(ESPECIE_gen)$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF0027'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF0948'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF1187'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF1274'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF0109'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FA694'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF0711'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF1082'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF1050'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF0192'])$overall),
                      as.data.frame(basic.stats(ESPECIE_gen[ESPECIE_gen$pop=='FF1256'])$overall)
)
colnames(tabela_pop)<-c('GERAL',
                        'ALTO',
                        'MEDIO',
                        'PIRANGA',
                        'CORRENTE',
                        'PIRACICABA',
                        'MATIPO',
                        'MANHUACU',
                        'BAIXO',
                        'GUANDU',
                        'STOANTONIO',
                        'CASCA'
                        )
#Remover estatisticas nao informaticas
tabela_pop<-tabela_pop[ !(row.names(tabela_pop) %in% c("Dst","Htp","Dstp","Fst","Ht","Fstp","Dest")), ]
lab_order = c('GERAL','ALTO','MEDIO','BAIXO','CASCA','CORRENTE','GUANDU','MANHUACU','MATIPO','PIRACICABA','PIRANGA','STOANTONIO') 
tabela_pop = tabela_pop[, lab_order]
grid.table(t(tabela_pop))
#=======================================================
#PASSO 4.FST: Compute pairwise FST (Weir & Cockerham 1984).
#=======================================================
ESPECIE_fst = genet.dist(ESPECIE_gen, method = "WC84") %>% round(digits = 3)
ESPECIE_fst
# Criar um dataframe a partir da matriz de FST
fst.mat=as.matrix(ESPECIE_fst)
colnames(fst.mat)
colnames(fst.mat)<-c('ALTO','MEDIO','PIRANGA','CORRENTE','PIRACICABA','MATIPO','MANHUACU','BAIXO','GUANDU','STOANTONIO','CASCA')  
rownames(fst.mat)  
rownames(fst.mat)<-c('ALTO','MEDIO','PIRANGA','CORRENTE','PIRACICABA','MATIPO','MANHUACU','BAIXO','GUANDU','STOANTONIO','CASCA')  
# Ordenar os dados
lab_order = c('ALTO','BAIXO','CASCA','CORRENTE','GUANDU','MANHUACU','MATIPO','MEDIO','PIRACICABA','PIRANGA','STOANTONIO') 
fst.mat1 = fst.mat[lab_order, ]
fst.mat = fst.mat1[, lab_order]
rm(fst.mat1)
fst.mat
ind = which(upper.tri(fst.mat), arr.ind = TRUE)
ind
?upper.tri()
fst.df = data.frame(Site1 = dimnames(fst.mat)[[2]][ind[,2]],
                    Site2 = dimnames(fst.mat)[[1]][ind[,1]],
                    Fst = fst.mat[ ind ])
fst.df
#Manter a ordem da matriz ao plotar os dados
fst.df$Site1 = factor(fst.df$Site1, levels = unique(fst.df$Site1))
fst.df$Site2 = factor(fst.df$Site2, levels = unique(fst.df$Site2))
#Converter valores negativos em zero
fst.df$Fst[fst.df$Fst < 0] = 0
# Print data.frame summary
fst.df %>% str
# Fst italic label
fst.label = expression(italic("F")[ST])
# Extract middle Fst value for gradient argument
mid = max(fst.df$Fst) / 2

# Plotar o HeatMap da matriz de FST
ggplot(data = fst.df, aes(x = Site1, y = Site2, fill = Fst))+
  geom_tile(colour = "black")+
  geom_text(aes(label = Fst), color="black", size = 4)+
  scale_fill_gradient2(low = "#46a0f5", mid = "white", high = "#f22300", midpoint = mid, name = fst.label, limits = c(0, max(fst.df$Fst)), breaks = c(0,0.05,0.10,0.15,0.20))+
  scale_x_discrete(expand = c(0,0))+
  scale_y_discrete(expand = c(0,0), position = "right")+
  theme(axis.text = element_text(colour = "black", size = 10),
        axis.title = element_blank(),
        panel.grid = element_blank(),
        panel.background = element_blank(),
        legend.position = "right",
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 10),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
  )
#=======================================================
#PASSO 5.DAPC - discriminant analysis of principal components
#=======================================================
set.seed(123)
#Substituir os missing data pela média das populações
x = tab(ESPECIE_gen, NA.method = "zero")
#Fazer um teste de cross validation
crossval = xvalDapc(x, ESPECIE_gen$pop, result = "groupMean", xval.plot = TRUE)
# Number of PCs with best stats (lower score = better)
crossval$`Root Mean Squared Error by Number of PCs of PCA`
crossval$`Number of PCs Achieving Highest Mean Success`
crossval$`Number of PCs Achieving Lowest MSE`
numPCs = as.numeric(crossval$`Number of PCs Achieving Lowest MSE`)
# Run a DAPC using site IDs as priors
dapc1 = dapc(ESPECIE_gen, ESPECIE_gen$pop, n.pca = numPCs, n.da = 3)
# Analyse how much percent of genetic variance is explained by each axis
percent = dapc1$eig/sum(dapc1$eig)*100
barplot(percent, ylab = "Genetic variance explained by eigenvectors (%)", ylim = c(0,100),
        names.arg = round(percent, 1))
# Create a data.frame containing individual coordinates
ind_coords = as.data.frame(dapc1$ind.coord)
# Rename columns of dataframe
colnames(ind_coords) = c("Axis1","Axis2","Axis3")
# Add a column containing individuals
ind_coords$Ind = indNames(ESPECIE_gen)
# Add a column with the site IDs
ind_coords$Site = ESPECIE_gen$pop
# Calculate centroid (average) position for each population
centroid = aggregate(cbind(Axis1, Axis2, Axis3) ~ Site, data = ind_coords, FUN = mean)
# Add centroid coordinates to ind_coords dataframe
ind_coords = left_join(ind_coords, centroid, by = "Site", suffix = c("",".cen"))
# Custom x and y labels
xlab = paste("Axis 1 (", format(round(percent[1], 1), nsmall=1)," %)", sep="")
ylab = paste("Axis 2 (", format(round(percent[2], 1), nsmall=1)," %)", sep="")
#Ajustar os IDs das amostras
ind_coords$Site
ind_coords$Site <- as.character(ind_coords$Site)
ind_coords['Site'][ind_coords['Site']=='FF0027']<-'ALTO'
ind_coords['Site'][ind_coords['Site']=='FF0948']<-'MEDIO'
ind_coords['Site'][ind_coords['Site']=='FF1187']<-'PIRANGA'
ind_coords['Site'][ind_coords['Site']=='FF1274']<-'CORRENTE'
ind_coords['Site'][ind_coords['Site']=='FF0109']<-'PIRACICABA'
ind_coords['Site'][ind_coords['Site']=='FA694']<-'MATIPO'
ind_coords['Site'][ind_coords['Site']=='FF0711']<-'MANHUACU'
ind_coords['Site'][ind_coords['Site']=='FF1082']<-'BAIXO'
ind_coords['Site'][ind_coords['Site']=='FF1050']<-'GUANDU'
ind_coords['Site'][ind_coords['Site']=='FF0192']<-'STOANTONIO'
ind_coords['Site'][ind_coords['Site']=='FF1256']<-'CASCA'
ind_coords$Site <- as.factor(ind_coords$Site)
ind_coords$Site

centroid$Site
centroid$Site <- as.character(centroid$Site)
centroid['Site'][centroid['Site']=='FF0027']<-'ALTO'
centroid['Site'][centroid['Site']=='FF0948']<-'MEDIO'
centroid['Site'][centroid['Site']=='FF1187']<-'PIRANGA'
centroid['Site'][centroid['Site']=='FF1274']<-'CORRENTE'
centroid['Site'][centroid['Site']=='FF0109']<-'PIRACICABA'
centroid['Site'][centroid['Site']=='FA694']<-'MATIPO'
centroid['Site'][centroid['Site']=='FF0711']<-'MANHUACU'
centroid['Site'][centroid['Site']=='FF1082']<-'BAIXO'
centroid['Site'][centroid['Site']=='FF1050']<-'GUANDU'
centroid['Site'][centroid['Site']=='FF0192']<-'STOANTONIO'
centroid['Site'][centroid['Site']=='FF1256']<-'CASCA'
centroid$Site <- as.factor(centroid$Site)
centroid$Site
# Definir a paleta de cores do grafico de DAPC
#cols = brewer.pal(nPop(ESPECIE_gen), "Set3")
centroid$Site
ind_coords$Site

Piranga="#5c1c17"
Casca="#b155b3"
Matipo="#5ccc0e"
Guandu="#427d15"
Alto="#394cb8"
Suacui="#ffa600"
Piracicaba="#cef542"
StoAntonio="#b55720"
Manhuacu="#53ac13"
Medio="#63bff0"
Caratinga="#681d6d"
Baixo="#a7d5ed"
Corrente="#ff90ff"

cols=c(Alto,Baixo,Casca,Corrente,Guandu,Manhuacu,Matipo,Medio,Piracicaba,Piranga,StoAntonio)

# Plotar o gráfico de DAPC: Scatter plot axis 1 vs. 2
ggplot(data = ind_coords, aes(x = Axis1, y = Axis2))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  # spider segments
  geom_segment(aes(xend = Axis1.cen, yend = Axis2.cen, colour = Site), show.legend = FALSE)+
  # points
  geom_point(aes(fill = Site), shape = 21, size = 3, show.legend = FALSE)+
  # centroids
  geom_label_repel(data = centroid, aes(label = Site, fill = Site), size = 4, show.legend = FALSE)+
    # colouring
  scale_fill_manual(values = cols)+
  scale_colour_manual(values = cols)+
  # custom labels
  labs(x = xlab, y = ylab)+
  ggtitle("ESPECIE DAPC")+theme( plot.title = element_text(hjust = 0.5))
#=======================================================
#PASSO 6. Agrupamento dos SNPs
#=======================================================
distancia_SNPs<-dist(x, method = "euclidean")
clustering_SNPs<-hclust(distancia_SNPs, method="ward.D2")

tree<-as.phylo(clustering_SNPs)

matriz_plot <- ggtree(tree, layout="fan", branch.length = "none") + geom_tiplab(size =4)
meta_data <- read.table("coleta.tab", sep="\t", header=TRUE,check.names=FALSE, stringsAsFactor=F, row.names = 1)
hm <- gheatmap(matriz_plot,meta_data, offset = 4, width=0.15, font.size=14, colnames = FALSE) + scale_fill_manual(values=c(
  "calha"="#005cd4",
  "Alto_rio_Doce"="#394cb8",
  "Médio_rio_Doce"="#63bff0",
  "Baixo_rio_Doce"="#a7d5ed",
  "tributario"="#d40000",
  "Tributário_-_Caratinga"="#681d6d",
  "Tributário_-_Casca"="#b155b3",
  "Tributário_-_Corrente"="#ff90ff",
  "Tributário_-_Guandu"="#427d15",
  "Tributário_-_Manhuaçu"="#53ac13",
  "Tributário_-_Matipó"="#5ccc0e",
  "Tributário_-_Piracicaba"="#cef542",
  "Tributário_-_Piranga"="#5c1c17",
  "Tributário_-_Sto_Antônio"="#b55720",
  "Tributário_-_Suaçuí"="#ffa600",
  "A"="#807680",
  "B"="#d4c8b8",
  "C"="#505050")
)
hm<-hm+ggtitle("SNP clustering: distance=euclidean; method=Ward.D2")+theme( plot.title = element_text(hjust = 0.5))
plot(hm)
#=======================================================
