Last updated: 2025-07-28
Checks: 6 1
Knit directory:
Machine-learning-e-redes-neurais-artificiais-no-melhoramento-genetico-do-cafeeiro/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of
the R Markdown file created these results, you’ll want to first commit
it to the Git repo. If you’re still working on the analysis, you can
ignore this warning. When you’re finished, you can run
wflow_publish to commit the R Markdown file and build the
HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20250709) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version ce1b0cf. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rproj.user/
Ignored: analysis/figure/
Ignored: analysis/modelos_mistos_dados_reais.Rmd
Ignored: data/dados_reais/Fenótipos 2014_2015_2016_C.arabica (1).xlsx
Ignored: data/dados_reais/Fenótipos 2014_2015_2016_C.arabica.xlsx
Ignored: data/dados_reais/FiltStep1_minQ10.0_minDP3_DPrange15-maxMeanDepth750_miss0.4_maf0.01_mac1_EMB_291701_SNPs_Raw.vcf.012.012
Ignored: data/dados_reais/FiltStep1_minQ10.0_minDP3_DPrange15-maxMeanDepth750_miss0.4_maf0.01_mac1_EMB_291701_SNPs_Raw.vcf.012.012.indv
Ignored: data/dados_reais/FiltStep1_minQ10.0_minDP3_DPrange15-maxMeanDepth750_miss0.4_maf0.01_mac1_EMB_291701_SNPs_Raw.vcf.012.012.pos
Ignored: data/dados_reais/FiltStep1_minQ10_minDP3_DPrange15-750_miss0.4_maf0.01_mac1_EMB_291701_filt1_snps_RAW.012.indv.txt
Ignored: data/dados_reais/FiltStep1_minQ10_minDP3_DPrange15-750_miss0.4_maf0.01_mac1_EMB_291701_filt1_snps_RAW.012.pos.txt
Ignored: data/dados_reais/FiltStep1_minQ10_minDP3_DPrange15-750_miss0.4_maf0.01_mac1_EMB_291701_filt1_snps_RAW.012.txt
Ignored: data/dados_reais/fpls-09-01934.pdf
Ignored: data/dados_reais/plants-13-01876-v2.pdf
Ignored: output/BLUPS_par_mmer.Rdata
Ignored: output/blups_all_wide.csv
Ignored: output/gwas_combined_results.rds
Ignored: output/gwas_results.rds
Ignored: output/importance_ML.rds
Ignored: output/importance_ML_mean.rds
Ignored: output/manhattan_BM.tiff
Ignored: output/manhattan_Cer.tiff
Ignored: output/manhattan_Prod.Cor.tiff
Ignored: output/mean_pheno.csv
Ignored: output/mod.RDS
Ignored: output/plot_lddecay.Rda
Ignored: output/pred_mod.RDS
Ignored: output/real_lddecay.tiff
Ignored: output/real_results_consolidated_10r_3f.xlsx
Ignored: output/real_results_consolidated_5r_5f.xlsx
Ignored: output/res.RDS
Ignored: output/result_sommer_pl_geracao_random.RDS
Ignored: output/result_sommer_random.RDS
Ignored: output/results_cv_GBLUP_a.rds
Ignored: output/results_cv_GBLUP_ad.rds
Ignored: output/results_cv_GBLUP_ade.rds
Ignored: output/simulated_results_consolidated.xlsx
Unstaged changes:
Modified: analysis/gws_dados_reais.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/gws_dados_reais.Rmd) and
HTML (docs/gws_dados_reais.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | ce1b0cf | WevertonGomesCosta | 2025-07-28 | update gwas_dados_reais .html |
| Rmd | 76b5e77 | WevertonGomesCosta | 2025-07-28 | add scripts .Rmd to gws and gwas for real and simulated data |
| Rmd | 82b3f68 | WevertonGomesCosta | 2025-07-24 | update save RDS pheno, geno, pos geno |
| Rmd | f8be55a | WevertonGomesCosta | 2025-07-24 | add gws and gwas_dados_reais.Rmd |
Este documento apresenta a implementação de modelos mistos para seleção genômica ampla (GWS) utilizando dados reais de Coffea arabica.
library(cvTools) # Para validação cruzada
library(kableExtra) # Para tabelas
library(parallel) # Para paralelização
library(doParallel) # Para paralelização
library(foreach) # Para paralelização
library(sommer) # Para modelos mistos
library(earth) # Para MARS
library(rpart) # Para árvores de decisão
library(randomForest) # Para Random Forest
library(gbm) # Para Boosting
library(neuralnet) # Para Redes Neurais
library(ANN2) # Para Redes Neurais
library(caret) # Para importância de variaveis
library(tidyverse) # Para manipulação de dadosOs dados fenotípicos foram obtidos a partir de um arquivo Excel contendo informações de ID e BLUPS das características fenotípicas, além de alguns marcadores.
# Carregar os dados fenotípicos
dados <- readxl::read_xlsx("data/dados_reais/dados_cafe.xlsx")
head(dados[1:4])# A tibble: 6 × 4
ID Prod.Cor Cer BM
<dbl> <dbl> <dbl> <dbl>
1 78 15.0 2.84 2.48
2 7 14.6 1.73 1.98
3 94 13.1 2.71 2.44
4 22 12.7 2.64 1.49
5 172 12.7 2.64 2.99
6 139 12.2 3.73 2.35Os rótulos dos indivíduos foram extraídos de um arquivo Excel, contendo informações sobre os códigos dos indivíduos e suas respectivas progenies. Esses rótulos foram utilizados para alinhar os dados fenotípicos com os genotípicos.
O conjunto de dados originais foi dividido em dois bancos de dados, cada um contendo informações sobre os indivíduos e seus respectivos genótipos.
Em todos os dois conjuntos os dados foram filtrados para maf 0.01 e mac 1, e os SNPs foram removidos se apresentassem mais de 40% de missing.
Para o primeiro banco de dados, os rótulos dos indivíduos foram extraídos de um arquivo VCF e um arquivo Excel contendo as informações dos progenitores. Os códigos dos indivíduos foram limpos e organizados para facilitar a análise.
Foi necessário remover um indivíduo que apresentava inconsistências nos dados. Após a limpeza, os rótulos foram unidos com os dados dos indivíduos, criando uma coluna de genótipo para cada progenie.
# Banco de dados 1
# Carregar os rótulos dos indivíduos
ind_geno1 <- read.table(
"data/dados_reais/FiltStep1_minQ10.0_minDP3_DPrange15-maxMeanDepth750_miss0.4_maf0.01_mac1_EMB_291701_SNPs_Raw.vcf.012.012.indv"
) %>%
mutate(codigos_limpos = gsub("EMB_291701_|_", "", V1))
# Criando uma coluna de ID
ind_geno1$Obs <- rownames(ind_geno1)
ind_geno1$Obs <- as.numeric(ind_geno1$Obs)
ind_geno1 <- ind_geno1 %>%
arrange(Obs)
head(ind_geno1) V1 codigos_limpos Obs
1 EMB_291701_P01_WH09 P01WH09 1
2 EMB_291701_P02_WB07 P02WB07 2
3 EMB_291701_P01_WH10 P01WH10 3
4 EMB_291701_P02_WE10 P02WE10 4
5 EMB_291701_P02_WA06 P02WA06 5
6 EMB_291701_P02_WG11 P02WG11 6# Carregar os rótulos dos indivíduos
rotulo_geno1 <- readxl::read_xls("data/dados_reais/Rótulos Novos SNPs.xls",
sheet = 2,
col_names = FALSE)
colnames(rotulo_geno1) <- c("Obs", "codigos_limpos", "progenie", "progenie_limpos")
# Limpeza dos códigos dos indivíduos
rotulo_geno1 <- rotulo_geno1 %>%
mutate(
codigos_limpos = gsub("P03", "P01", codigos_limpos),
codigos_limpos = gsub("P04", "P02", codigos_limpos),
codigos_limpos = gsub("P05", "P03", codigos_limpos),
ID = as.numeric(str_extract(progenie, "(?<=\\()[0-9]+(?=\\))"))
) %>%
arrange(Obs)
head(rotulo_geno1)# A tibble: 6 × 5
Obs codigos_limpos progenie progenie_limpos ID
<dbl> <chr> <chr> <chr> <dbl>
1 1 P01WH09 (95) 7-5 7-5 95
2 2 P02WB07 (171) 12-6 12-6 171
3 3 P01WH10 (96) 7-6 7-6 96
4 4 P02WE10 (74) 5-14 5-14 74
5 5 P02WA06 (104) 7-14 7-14 104
6 6 P02WG11 (36) 3-6 3-6 36# Verificação de consistência entre os códigos dos indivíduos
ind_geno1$codigos_limpos == rotulo_geno1$codigos_limpos [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE# Removendo individuo
ind_geno1 <- ind_geno1 %>%
slice(-37)
# Verificação de consistência entre os códigos dos indivíduos
all(unique(ind_geno1$codigos_limpos == rotulo_geno1$codigos_limpos))[1] TRUE# Criando uma coluna de ID
ind_geno1$Obs <- rownames(ind_geno1)
ind_geno1$Obs <- as.numeric(ind_geno1$Obs)
ind_geno1 <- ind_geno1 %>%
arrange(Obs)
# Unindo os dados dos indivíduos com os rótulos
ind_geno1 <- ind_geno1 %>%
full_join(rotulo_geno1) %>%
mutate(genotipo = paste("PROG", gsub("-", "_", progenie_limpos)))
head(ind_geno1) V1 codigos_limpos Obs progenie progenie_limpos ID
1 EMB_291701_P01_WH09 P01WH09 1 (95) 7-5 7-5 95
2 EMB_291701_P02_WB07 P02WB07 2 (171) 12-6 12-6 171
3 EMB_291701_P01_WH10 P01WH10 3 (96) 7-6 7-6 96
4 EMB_291701_P02_WE10 P02WE10 4 (74) 5-14 5-14 74
5 EMB_291701_P02_WA06 P02WA06 5 (104) 7-14 7-14 104
6 EMB_291701_P02_WG11 P02WG11 6 (36) 3-6 3-6 36
genotipo
1 PROG 7_5
2 PROG 12_6
3 PROG 7_6
4 PROG 5_14
5 PROG 7_14
6 PROG 3_6Para o segundo banco de dados, os rótulos dos indivíduos foram extraídos de um arquivo de texto contendo informações sobre os indivíduos e seus respectivos genótipos. Os códigos dos indivíduos foram limpos e organizados, e uma coluna de ID foi criada para facilitar a análise.
Os rótulos dos indivíduos foram então unidos com os dados dos indivíduos, criando uma coluna de genótipo para cada progenie. A consistência entre os códigos dos indivíduos foi verificada, garantindo que os dados estivessem corretos.
# Banco de dados 2
# Carregar os rótulos dos indivíduos
ind_geno2 <- read.table(
"data/dados_reais/FiltStep1_minQ10_minDP3_DPrange15-750_miss0.4_maf0.01_mac1_EMB_291701_filt1_snps_RAW.012.indv.txt"
) %>%
mutate(
codigos_limpos = gsub("EMB_291701_|_", "", V1),
codigos_limpos2 = substr(codigos_limpos, 4, nchar(codigos_limpos))
)
# Criando uma coluna de ID
ind_geno2$Obs <- rownames(ind_geno2)
ind_geno2$Obs <- as.numeric(ind_geno2$Obs)
ind_geno2 <- ind_geno2 %>%
arrange(Obs)
head(ind_geno2) V1 codigos_limpos codigos_limpos2 Obs
1 EMB_291701_P03_WD11 P03WD11 WD11 1
2 EMB_291701_P03_WF06 P03WF06 WF06 2
3 EMB_291701_P03_WC12 P03WC12 WC12 3
4 EMB_291701_P03_WE08 P03WE08 WE08 4
5 EMB_291701_P03_WE03 P03WE03 WE03 5
6 EMB_291701_P03_WC01 P03WC01 WC01 6# Carregar os rótulos dos indivíduos
rotulo_geno2 <- readxl::read_xlsx("data/dados_reais/Rótulos C. arabica_Proj_Piloto.xlsx",
col_names = FALSE)
colnames(rotulo_geno2) <- c("Obs", "codigos_limpos2", "progenie", "progenie_limpos", "info")
# Limpeza dos códigos dos indivíduos
rotulo_geno2 <- rotulo_geno2 %>%
mutate(ID = as.numeric(str_extract(progenie, "(?<=\\()[0-9]+(?=\\))"))) %>%
arrange(Obs)
head(rotulo_geno2)# A tibble: 6 × 6
Obs codigos_limpos2 progenie progenie_limpos info ID
<dbl> <chr> <chr> <chr> <chr> <dbl>
1 1 WD11 (184) 13-4 13-4 F2_1 184
2 2 WF06 (47) 4-2 4-2 RTs_2 47
3 3 WC12 (168) 12-3 12-3 F2_3 168
4 4 WE08 (79) 6-4 6-4 RTs_4 79
5 5 WE03 (63) 5-3 5-3 RTr_5 63
6 6 WC01 (20) 2-5 2-5 RTs_6 20# Verificação de consistência entre os códigos dos indivíduos
all(unique(ind_geno2$codigos_limpos2 == rotulo_geno2$codigos_limpos2))[1] TRUE# Unindo os dados dos indivíduos com os rótulos
ind_geno2 <- ind_geno2 %>%
full_join(rotulo_geno2) %>%
mutate(genotipo = paste("PROG", gsub("-", "_", progenie_limpos)))
head(ind_geno2) V1 codigos_limpos codigos_limpos2 Obs progenie
1 EMB_291701_P03_WD11 P03WD11 WD11 1 (184) 13-4
2 EMB_291701_P03_WF06 P03WF06 WF06 2 (47) 4-2
3 EMB_291701_P03_WC12 P03WC12 WC12 3 (168) 12-3
4 EMB_291701_P03_WE08 P03WE08 WE08 4 (79) 6-4
5 EMB_291701_P03_WE03 P03WE03 WE03 5 (63) 5-3
6 EMB_291701_P03_WC01 P03WC01 WC01 6 (20) 2-5
progenie_limpos info ID genotipo
1 13-4 F2_1 184 PROG 13_4
2 4-2 RTs_2 47 PROG 4_2
3 12-3 F2_3 168 PROG 12_3
4 6-4 RTs_4 79 PROG 6_4
5 5-3 RTr_5 63 PROG 5_3
6 2-5 RTs_6 20 PROG 2_5Os dois conjuntos de dados de indivíduos foram unidos, e os indivíduos que estão presentes no conjunto de dados fenotípicos foram filtrados. A coluna de genótipo foi ajustada para refletir os genótipos dos indivíduos presentes no conjunto de dados fenotípicos.
# Unindo os dois conjuntos de dados de indivíduos
ind_geno <- ind_geno1 %>%
full_join(ind_geno2)
# Filtrando os indivíduos que estão presentes no conjunto de dados fenotípicos
ind_geno <- ind_geno %>%
mutate(Obs = 1:nrow(ind_geno)) %>%
filter(ID %in% dados$ID)
head(ind_geno) V1 codigos_limpos Obs progenie progenie_limpos ID
1 EMB_291701_P01_WH09 P01WH09 1 (95) 7-5 7-5 95
2 EMB_291701_P02_WB07 P02WB07 2 (171) 12-6 12-6 171
3 EMB_291701_P01_WH10 P01WH10 3 (96) 7-6 7-6 96
4 EMB_291701_P02_WE10 P02WE10 4 (74) 5-14 5-14 74
5 EMB_291701_P02_WA06 P02WA06 5 (104) 7-14 7-14 104
6 EMB_291701_P02_WG11 P02WG11 6 (36) 3-6 3-6 36
genotipo codigos_limpos2 info
1 PROG 7_5 <NA> <NA>
2 PROG 12_6 <NA> <NA>
3 PROG 7_6 <NA> <NA>
4 PROG 5_14 <NA> <NA>
5 PROG 7_14 <NA> <NA>
6 PROG 3_6 <NA> <NA># Ajustando os nomes dos genótipos
dados <- dados %>%
left_join(ind_geno[,c("ID", "genotipo")])
head(dados[1:4])# A tibble: 6 × 4
ID Prod.Cor Cer BM
<dbl> <dbl> <dbl> <dbl>
1 78 15.0 2.84 2.48
2 7 14.6 1.73 1.98
3 94 13.1 2.71 2.44
4 22 12.7 2.64 1.49
5 172 12.7 2.64 2.99
6 139 12.2 3.73 2.35Os dados genotípicos foram filtrados para remover SNPs específicos e garantir a consistência entre os dois conjuntos de dados. As colunas correspondentes aos SNPs indesejados foram removidas, e os genótipos foram ajustados para refletir apenas os indivíduos presentes no conjunto de dados fenotípicos.
Para os dados genotípicos do banco de dados 1 foi realizada a remoção do índivíduo sem complementação, como feito no item 1.2.1. Após isso, os dados foram carregados e as posições dos SNPs foram extraídas. Os SNPs indesejados foram removidos, e os genótipos foram ajustados para refletir apenas os indivíduos presentes no conjunto de dados fenotípicos.
# Banco de dados 1
# Passo 1: carregar os dados genotípicos e remover individuos indesejados
geno1 <- read.table("data/dados_reais/FiltStep1_minQ10.0_minDP3_DPrange15-maxMeanDepth750_miss0.4_maf0.01_mac1_EMB_291701_SNPs_Raw.vcf.012.012")[-1]
geno1 <- geno1 %>%
slice(-37)
# Passo 2: carregar as posições dos SNPs
pos_geno1 <- read.table("data/dados_reais/FiltStep1_minQ10.0_minDP3_DPrange15-maxMeanDepth750_miss0.4_maf0.01_mac1_EMB_291701_SNPs_Raw.vcf.012.012.pos")
pos_geno1 <- pos_geno1 %>%
mutate(pos = paste0(V1, ":", V2))
pos_geno1$n_markers = 1:nrow(pos_geno1)
colnames(geno1) <- pos_geno1$pos
# Passo 3: remover SNPs indesejados
markers1 <- pos_geno1 %>%
filter(pos == "chr0:177480018" | pos == "chr0:177480063" | pos == "chr0:177480100")
geno1 <- geno1 %>%
select(-all_of(markers1$n_markers))Para os dados genotípicos do banco de dados 2, os SNPs indesejados foram removidos, e as posições dos SNPs foram extraídas. Os genótipos foram ajustados para refletir apenas os indivíduos presentes no conjunto de dados fenotípicos.
# Banco de dados 2
# Passo 1: carregar os dados genotípicos
geno2 <- read.table("data/dados_reais/FiltStep1_minQ10_minDP3_DPrange15-750_miss0.4_maf0.01_mac1_EMB_291701_filt1_snps_RAW.012.txt")[-1]
pos_geno2 <- read.table("data/dados_reais/FiltStep1_minQ10_minDP3_DPrange15-750_miss0.4_maf0.01_mac1_EMB_291701_filt1_snps_RAW.012.pos.txt")
# Passo 2: carregar as posições dos SNPs
pos_geno2 <- pos_geno2 %>%
mutate(pos = paste0(V1, ":", V2))
pos_geno2$n_markers = 1:nrow(pos_geno2)
colnames(geno2) <- pos_geno2$pos
# Passo 3: remover SNPs indesejados
markers2 <- pos_geno2 %>%
filter(pos == "chr0:177480018" | pos == "chr0:177480063" | pos == "chr0:177480100")
geno2 <- geno2 %>%
select(-all_of(markers2$n_markers))Os dados genotípicos dos dois bancos foram unidos, e as colunas foram nomeadas de acordo com os genótipos dos indivíduos. Os genótipos foram filtrados para garantir que apenas os indivíduos presentes no conjunto de dados fenotípicos fossem incluídos. As colunas comuns entre os dois conjuntos de dados foram alinhadas, e os valores dos genótipos foram ajustados para refletir a presença de SNPs.
# Passo 4: Nomear as colunas
geno1$genotipo1 <- ind_geno1$genotipo
geno2$genotipo1 <- ind_geno2$genotipo
# Passo 5: Filtrar os genótipos que estão presentes no conjunto de dados fenotípicos
geno1 <- geno1 %>%
filter(genotipo1 %in% dados$genotipo)
geno2 <- geno2 %>%
filter(genotipo1 %in% dados$genotipo)
# Passo 6: Alinhar os genótipos
colunas_comuns <- intersect(names(geno1), names(geno2))
# Passo 7: Unir os genótipos
geno <- bind_rows(
select(geno1, all_of(colunas_comuns)),
select(geno2, all_of(colunas_comuns))
)
# Passo 8: Ajustar os valores dos genótipos
geno[geno == -1] <- 0
# Passo 9: Ajustar os nomes das colunas
rownames(geno) <- geno$genotipo1
# Passo 10: Remover a coluna de genotipo1
geno <- geno %>%
select(-genotipo1)
# Passo 11: Padronizar os valores dos genótipos para -1, 0 e 1
geno <- geno - 1
geno[1:5,1:5] CAonly_Contig15745:814-1134:165 CAonly_Contig2698:0-267:74
PROG 7_5 -1 0
PROG 12_6 -1 0
PROG 7_6 0 0
PROG 5_14 -1 0
PROG 7_14 -1 0
CAonly_Contig2698:0-267:88 CAonly_Contig2698:0-267:110
PROG 7_5 0 0
PROG 12_6 0 0
PROG 7_6 0 0
PROG 5_14 0 0
PROG 7_14 0 0
CAonly_Contig2698:0-267:125
PROG 7_5 0
PROG 12_6 -1
PROG 7_6 0
PROG 5_14 -1
PROG 7_14 -1# Passo 12: Ajustar o arquivo de fenotipos
pheno <- dados[c(1:4, ncol(dados))]
rownames(pheno) <- pheno$genotipo
pheno[1:5,]# A tibble: 5 × 5
ID Prod.Cor Cer BM genotipo
<dbl> <dbl> <dbl> <dbl> <chr>
1 78 15.0 2.84 2.48 PROG 6_3
2 7 14.6 1.73 1.98 PROG 1_7
3 94 13.1 2.71 2.44 PROG 7_4
4 22 12.7 2.64 1.49 PROG 2_7
5 172 12.7 2.64 2.99 PROG 12_7# Passo 13: Agrupar os arquivos de posição dos SNPs
pos_geno <- bind_rows(pos_geno1, pos_geno2)
pos_geno[1:5,] V1 V2 pos n_markers
1 CAonly_Contig15745:814-1134 165 CAonly_Contig15745:814-1134:165 1
2 CAonly_Contig2698:0-267 74 CAonly_Contig2698:0-267:74 2
3 CAonly_Contig2698:0-267 88 CAonly_Contig2698:0-267:88 3
4 CAonly_Contig2698:0-267 110 CAonly_Contig2698:0-267:110 4
5 CAonly_Contig2698:0-267 125 CAonly_Contig2698:0-267:125 5Para validação cruzada, definimos o número de dobras (nfolds) e repetições (nrept). Em seguida, alinhamos os IDs dos genótipos entre os dados fenotípicos e genotípicos, garantindo que correspondam corretamente. Por fim, geramos uma lista de dobras para validação cruzada.
O modelo G-BLUP (Genomic Best Linear Unbiased Prediction) é uma abordagem amplamente utilizada para predição genômica, que considera tanto os efeitos aditivos quanto os efeitos de dominância dos marcadores genéticos. Neste estudo, implementamos três variantes do G-BLUP: Aditivo, Aditivo-Dominante e Aditivo-Dominante-Epistástico.
O modelo é implementado utilizando a função mmes do
pacote sommer, que permite ajustar modelos mistos com
estruturas de parentesco complexas. A validação cruzada é realizada para
avaliar a acurácia preditiva do modelo.
Para o modelo G-BLUP aditivo, utilizamos a matriz de parentesco aditivo (K) calculada a partir dos dados genotípicos. Este modelo assume que os efeitos dos marcadores são puramente aditivos, sem considerar interações de dominância ou epistasia.
# Construção da matriz de parentesco aditivo (K)
K <- A.mat(as.matrix(geno))
colnames(K) <- rownames(K) <- rownames(pheno) <- pheno$genotipoO G-BLUP aditivo é dada pelo modelo:
\[ \text{y} = \mu + \text{X}\beta + \text{Z}u + e \]
onde: - \(\text{y}\) é o vetor de fenótipos, - \(\mu\) é a média geral, - \(\text{X}\) é a matriz de efeitos fixos (neste caso, intercepto), - \(\beta\) é o vetor de efeitos fixos, - \(\text{Z}\) é a matriz de efeitos aleatórios (neste caso, a matriz de parentesco aditivo \(K\)), - \(u\) é o vetor de efeitos aleatórios (GEBVs), - \(e\) é o vetor de erros aleatórios.
# Configuração do cluster para paralelização
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para G-BLUP aditivo
results_cv_GBLUP_a <- foreach(
idx = seq_len(nrow(combos)),
.packages = c("sommer", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
idx=3
# Extração dos parâmetros da combinação atual
tr <- combos$trait[idx]
r <- combos$rep[idx]
f <- combos$fold[idx]
# Preparação dos dados para o modelo
ph <- pheno[,c("genotipo", tr)]
# Extração da matriz de parentesco aditivo
K1 <- K[ph$genotipo, ph$genotipo]
# Ajuste dos nomes das colunas e linhas da matriz K1
colnames(K1) <- rownames(K1) <- ph$genotipo
# Criação do dataframe para o modelo
td <- ph
# Definição dos índices de validação cruzada
val_idx <- fold.list[[r]][[f]]
td[val_idx, tr] <- NA
# Ajuste do modelo G-BLUP aditivo
fit <- mmes(
fixed = as.formula(paste(tr, "~ 1")),
random = ~ vsm(ism(genotipo), Gu = K1),
data = td,
rcov = ~ units
)
# Extração dos preditos e verdadeiros valores
pred <- fit$u
truth <- ph[[tr]]
# Predição e métricas
pred <- data.frame(fit$u, ph)
# Cálculo dos GEBVs para treino e validação
pred_T <- pred[-val_idx, "fit.u"]
truth_T <- pred[[tr]][-val_idx]
pred_V <- pred[val_idx, "fit.u"]
truth_V <- pred[[tr]][val_idx]
# Criação do dataframe de resultados
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, truth_T), 4),
R2_Valid = round(cor(pred_V, truth_V), 4),
RMSE_Train = round(sqrt(mean((pred_T - truth_T)^2)), 4),
RMSE_Valid = round(sqrt(mean((pred_V - truth_V)^2)), 4),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V)
)
}
# Finalização do cluster
stopCluster(cl)
# Salvando os resultados da validação cruzada
saveRDS(results_cv_GBLUP_a, file = "output/real_results_cv_GBLUP_a.rds")O modelo G-BLUP aditivo-dominante estende o G-BLUP aditivo ao incluir efeitos de dominância, permitindo capturar interações entre alelos homozigotos e heterozigotos. A matriz de parentesco aditivo-dominante (D) é utilizada para modelar esses efeitos.
# Construção das matrizes de parentesco aditivo (K) e aditivo-dominante (D)
K <- A.mat(as.matrix(geno))
D <- D.mat(as.matrix(geno))
colnames(K) <- rownames(K) <- rownames(pheno) <- pheno$genotipo
colnames(D) <- rownames(D) <- rownames(pheno) O G-BLUP aditivo-dominante é dada pelo modelo:
\[ \text{y} = \mu + \text{X}\beta + \text{Z}_a u_a + \text{Z}_d u_d + e \]
onde: - \(\text{y}\) é o vetor de fenótipos, - \(\mu\) é a média geral, - \(\text{X}\) é a matriz de efeitos fixos (neste caso, intercepto), - \(\beta\) é o vetor de efeitos fixos, - \(\text{Z}_a\) é a matriz de efeitos aleatórios aditivos (matriz K), - \(u_a\) é o vetor de efeitos aleatórios aditivos (GEBVs aditivos), - \(\text{Z}_d\) é a matriz de efeitos aleatórios dominantes (matriz D), - \(u_d\) é o vetor de efeitos aleatórios dominantes (GEBVs dominantes), - \(e\) é o vetor de erros aleatórios.
# Configuração do cluster para paralelização
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para G-BLUP aditivo-dominante
results_cv_GBLUP_ad <- foreach(
idx = seq_len(nrow(combos)),
.packages = c("sommer", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
# Extração dos parâmetros da combinação atual
tr <- combos$trait[idx]
r <- combos$rep[idx]
f <- combos$fold[idx]
# Preparação dos dados para o modelo
ph <- pheno[,c("genotipo", tr)]
K1 <- K[ph$genotipo, ph$genotipo]
D1 <- D[ph$genotipo, ph$genotipo]
# Ajuste dos nomes das colunas e linhas das matrizes K1 e D1
colnames(K1) <- rownames(K1) <- ph$genotipo
colnames(D1) <- rownames(D1) <- ph$genotipo
# Criação do dataframe para o modelo
td <- ph
td$genotipod <- td$genotipo
# Definição dos índices de validação cruzada
val_idx <- fold.list[[r]][[f]]
td[val_idx, tr] <- NA
# Ajuste do modelo G-BLUP aditivo-dominante
fit <- mmes(
fixed = as.formula(paste(tr, "~ 1")),
random = ~ vsm(ism(genotipo), Gu = K1) + vsm(ism(genotipod), Gu = D1),
data = td,
rcov = ~ units,
tolParInv = ifelse(idx == 29 | idx == 45, 5e-3, 1e-6) # Ajuste de tolerância para a matriz inversa,
)
# Extração dos preditos e verdadeiros valores
pred_a <- fit$uList$`vsm(ism(genotipo), Gu = K1`
pred_d <- fit$uList$`vsm(ism(genotipod), Gu = D1`
# Soma dos preditos aditivos e dominantes
pred <- pred_a + pred_d
truth <- ph[[tr]]
# Predição e métricas
pred <- data.frame(pred, ph)
# Cálculo dos GEBVs para treino e validação
pred_T <- pred[-val_idx, "mu"]
truth_T <- pred[[tr]][-val_idx]
pred_V <- pred[val_idx, "mu"]
truth_V <- pred[[tr]][val_idx]
# Criação do dataframe de resultados
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, truth_T), 4),
R2_Valid = round(cor(pred_V, truth_V), 4),
RMSE_Train = round(sqrt(mean((pred_T - truth_T)^2)), 4),
RMSE_Valid = round(sqrt(mean((pred_V - truth_V)^2)), 4),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V)
)
}
# Finalização do cluster
stopCluster(cl)
# Salvando os resultados da validação cruzada
saveRDS(results_cv_GBLUP_ad, file = "output/real_results_cv_GBLUP_ad.rds")O modelo G-BLUP aditivo-dominante-epistástico estende o G-BLUP aditivo-dominante ao incluir efeitos epistáticos, permitindo capturar interações complexas entre alelos. A matriz de parentesco epistático (E) é utilizada para modelar esses efeitos.
# Construção das matrizes de parentesco aditivo (K), aditivo-dominante (D) e aditivo-dominante-epistástico (E)
K <- A.mat(as.matrix(geno))
D <- D.mat(as.matrix(geno))
E <- E.mat(as.matrix(geno))
colnames(K) <- rownames(K) <- rownames(pheno) <- pheno$genotipo
colnames(D) <- rownames(D) <- rownames(pheno)
colnames(E) <- rownames(E) <- rownames(pheno)O G-BLUP aditivo-dominante-epistástico é dada pelo modelo:
\[ \text{y} = \mu + \text{X}\beta + \text{Z}_a u_a + \text{Z}_d u_d + \text{Z}_e u_e + e \]
onde: - \(\text{y}\) é o vetor de fenótipos, - \(\mu\) é a média geral, - \(\text{X}\) é a matriz de efeitos fixos (neste caso, intercepto), - \(\beta\) é o vetor de efeitos fixos, - \(\text{Z}_a\) é a matriz de efeitos aleatórios aditivos (matriz K), - \(u_a\) é o vetor de efeitos aleatórios aditivos (GEBVs aditivos), - \(\text{Z}_d\) é a matriz de efeitos aleatórios dominantes (matriz D), - \(u_d\) é o vetor de efeitos aleatórios dominantes (GEBVs dominantes), - \(\text{Z}_e\) é a matriz de efeitos aleatórios epistáticos (matriz E), - \(u_e\) é o vetor de efeitos aleatórios epistáticos (GEBVs epistáticos), - \(e\) é o vetor de erros aleatórios.
# Configuração do cluster para paralelização
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para G-BLUP aditivo-dominante-epistástico
results_cv_GBLUP_ade <- foreach(
idx = seq_len(nrow(combos)),
.packages = c("sommer", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
# Extração dos parâmetros da combinação atual
tr <- combos$trait[idx]
r <- combos$rep[idx]
f <- combos$fold[idx]
# Preparação dos dados para o modelo
ph <- pheno[,c("genotipo", tr)]
K1 <- K[ph$genotipo, ph$genotipo]
D1 <- D[ph$genotipo, ph$genotipo]
E1 <- E[ph$genotipo, ph$genotipo]
# Ajuste dos nomes das colunas e linhas das matrizes K1, D1 e E1
colnames(K1) <- rownames(K1) <- ph$genotipo
colnames(D1) <- rownames(D1) <- ph$genotipo
colnames(E1) <- rownames(E1) <- ph$genotipo
# Criação do dataframe para o modelo
td <- ph
td$genotipod <- td$genotipo
td$genotipoe <- td$genotipo
# Definição dos índices de validação cruzada
val_idx <- fold.list[[r]][[f]]
td[val_idx, tr] <- NA
# Ajuste do modelo G-BLUP aditivo-dominante-epistástico
fit <- mmes(
fixed = as.formula(paste(tr, "~ 1")),
random = ~ vsm(ism(genotipo), Gu = K1) + vsm(ism(genotipod), Gu = D1) + vsm(ism(genotipoe), Gu = E1) ,
data = td,
rcov = ~ units,
tolParInv = ifelse(idx == 29 | idx == 45, 5e-3, 1e-6) # Ajuste de tolerância para a matriz inversa,
)
# Extração dos preditos e verdadeiros valores
pred_a <- fit$uList$`vsm(ism(genotipo), Gu = K1`
pred_d <- fit$uList$`vsm(ism(genotipod), Gu = D1`
pred_e <- fit$uList$`vsm(ism(genotipoe), Gu = E1`
# Soma dos preditos aditivos, dominantes e epistáticos
pred <- pred_a + pred_d + pred_e
truth <- ph[[tr]]
# Predição e métricas
pred <- data.frame(fit$u, ph)
# Cálculo dos GEBVs para treino e validação
pred_T <- pred[-val_idx, "mu"]
truth_T <- pred[[tr]][-val_idx]
pred_V <- pred[val_idx, "mu"]
truth_V <- pred[[tr]][val_idx]
# Criação do dataframe de resultados
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, truth_T), 4),
R2_Valid = round(cor(pred_V, truth_V), 4),
RMSE_Train = round(sqrt(mean((pred_T - truth_T)^2)), 4),
RMSE_Valid = round(sqrt(mean((pred_V - truth_V)^2)), 4),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V)
)
}
# Finalização do cluster
stopCluster(cl)
# Salvando os resultados da validação cruzada
saveRDS(results_cv_GBLUP_ade, file = "output/real_results_cv_GBLUP_ade.rds")A Multivariate Adaptative Regression Splines (MARS) é uma técnica de
modelagem não paramétrica que pode capturar relações complexas entre
variáveis preditoras e resposta. Neste estudo, implementamos três
variantes do MARS: MARS 1 (grau 1) e MARS 2 (grau 2) e MARS 3 (grau 3) .
Ambas as variantes são ajustadas utilizando a função earth
do pacote earth, que é otimizada para lidar com grandes
conjuntos de dados.
O modelo MARS 1 é uma versão simplificada que utiliza splines lineares para modelar a relação entre as variáveis preditoras e a resposta. Este modelo é adequado para capturar relações lineares e interações simples entre os marcadores genéticos.
# Precisamos dividir nosso espaço de memoria para não estourar a memoria
max_workers <- floor(32 / 6.5) # = 4 workers
cores_para_usar <- min(max_workers, detectCores() - 1)
cl <- makeCluster(cores_para_usar)
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos_mars <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para MARS 1
results_cv_MARS_1 <- foreach(
idx = seq_len(nrow(combos_mars)),
.packages = c("earth", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
# Extração dos parâmetros da combinação atual
tr <- combos_mars$trait[idx]
r <- combos_mars$rep[idx]
f <- combos_mars$fold[idx]
# Preparação dos dados para o modelo
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
# Criando dados de treinamento e validação
yT <- pheno[train_idx, tr][[1]]
xT <- as.matrix(geno[train_idx, ])
yV <- pheno[valid_idx, tr][[1]]
xV <- as.matrix(geno[valid_idx, ])
# Ajuste do modelo MARS 1
fit <- earth(x = xT, y = yT, degree = 1)
# Predição nos conjuntos de treinamento e validação
pred_T <- predict(fit, xT)
pred_V <- predict(fit, xV)
# Extrair importância das variáveis
imp <- evimp(fit, trim = FALSE)
# Seleciona colunas relevantes e converte para data.frame
imp_df <- as.data.frame(imp[, c("rss", "nsubsets")]) # rss: Overall importance
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
# Ajusta nomes das colunas e remove sufixo indesejado dos marcadores
imp_df <- imp_df %>%
select(Overall = rss, marker, Rep, Fold) %>%
mutate(marker = stringr::str_replace_all(marker, "-unused", ""))
# Limpeza de memória
gc()
# Criação do dataframe de resultados
tibble(
Trait = tr,
Degree = 1,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 4),
R2_Valid = round(cor(pred_V, yV), 4),
RMSE_Train = round(sqrt(mean(((pred_T - yT)^2))), 4),
RMSE_Valid = round(sqrt(mean(((pred_V - yV)^2))), 4),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df))
}
# Finalização do cluster
stopCluster(cl)
# Salvando os resultados da validação cruzada
saveRDS(results_cv_MARS_1, file = "output/real_results_cv_MARS_1.rds")O modelo MARS 2 é uma versão mais complexa que utiliza splines quadráticos para modelar a relação entre as variáveis preditoras e a resposta. Este modelo é adequado para capturar relações não lineares mais complexas entre os marcadores genéticos.
# Precisamos dividir nosso espaço de memoria para não estourar a memoria
max_workers <- floor(32 / 6.5) # = 4 workers
cores_para_usar <- min(max_workers, detectCores() - 1)
cl <- makeCluster(cores_para_usar)
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos_mars <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para MARS 2
results_cv_MARS_2 <- foreach(
idx = seq_len(nrow(combos_mars)),
.packages = c("earth", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
# Extração dos parâmetros da combinação atual
tr <- combos_mars$trait[idx]
r <- combos_mars$rep[idx]
f <- combos_mars$fold[idx]
# Preparação dos dados para o modelo
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
# Criando dados de treinamento e validação
yT <- pheno[train_idx, tr][[1]]
xT <- as.matrix(geno[train_idx, ])
yV <- pheno[valid_idx, tr][[1]]
xV <- as.matrix(geno[valid_idx, ])
# Ajuste do modelo MARS 2
fit <- earth(x = xT, y = yT, degree = 2)
# Predição nos conjuntos de treinamento e validação
pred_T <- predict(fit, xT)
pred_V <- predict(fit, xV)
# Extrair importância das variáveis
imp <- evimp(fit, trim = FALSE)
# Seleciona colunas relevantes e converte para data.frame
imp_df <- as.data.frame(imp[, c("rss", "nsubsets")]) # rss: Overall importance
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
# Ajusta nomes das colunas e remove sufixo indesejado dos marcadores
imp_df <- imp_df %>%
select(Overall = rss, marker, Rep, Fold) %>%
mutate(marker = stringr::str_replace_all(marker, "-unused", ""))
# Limpeza de memória
gc()
# Criação do dataframe de resultados
tibble(
Trait = tr,
Degree = 2,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 4),
R2_Valid = round(cor(pred_V, yV), 4),
RMSE_Train = round(sqrt(mean(((pred_T - yT)^2))), 4),
RMSE_Valid = round(sqrt(mean(((pred_V - yV)^2))), 4),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df)
)
}
# Finalização do cluster
stopCluster(cl)
# Salvando os resultados da validação cruzada
saveRDS(results_cv_MARS_2, file = "output/real_results_cv_MARS_2.rds")O modelo MARS 3 é uma versão ainda mais complexa que utiliza splines cúbicos para modelar a relação entre as variáveis preditoras e a resposta. Este modelo é adequado para capturar relações não lineares muito complexas entre os marcadores genéticos.
# Precisamos dividir nosso espaço de memoria para não estourar a memoria
max_workers <- floor(32 / 6.5) # = 4 workers
cores_para_usar <- min(max_workers, detectCores() - 1)
cl <- makeCluster(cores_para_usar)
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos_mars <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para MARS 3
results_cv_MARS_3 <- foreach(
idx = seq_len(nrow(combos_mars)),
.packages = c("earth", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
# Extração dos parâmetros da combinação atual
tr <- combos_mars$trait[idx]
r <- combos_mars$rep[idx]
f <- combos_mars$fold[idx]
# Preparação dos dados para o modelo
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
# Criando dados de treinamento e validação
yT <- pheno[train_idx, tr][[1]]
xT <- as.matrix(geno[train_idx, ])
yV <- pheno[valid_idx, tr][[1]]
xV <- as.matrix(geno[valid_idx, ])
# Ajuste do modelo MARS 3
fit <- earth(x = xT, y = yT, degree = 3)
# Predição nos conjuntos de treinamento e validação
pred_T <- predict(fit, xT)
pred_V <- predict(fit, xV)
# Extrair importância das variáveis
imp <- evimp(fit, trim = FALSE)
# Seleciona colunas relevantes e converte para data.frame
imp_df <- as.data.frame(imp[, c("rss", "nsubsets")]) # rss: Overall importance
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
# Ajusta nomes das colunas e remove sufixo indesejado dos marcadores
imp_df <- imp_df %>%
select(Overall = rss, marker, Rep, Fold) %>%
mutate(marker = stringr::str_replace_all(marker, "-unused", ""))
# Limpeza de memória
gc()
# Criação do dataframe de resultados
tibble(
Trait = tr,
Degree = 3,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 4),
R2_Valid = round(cor(pred_V, yV), 4),
RMSE_Train = round(sqrt(mean(((pred_T - yT)^2))), 4),
RMSE_Valid = round(sqrt(mean(((pred_V - yV)^2))), 4),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df)
)
}
# Finalização do cluster
stopCluster(cl)
# Salvando os resultados da validação cruzada
saveRDS(results_cv_MARS_3, file = "output/real_results_cv_MARS_3.rds")A árvore de decisão é um modelo de aprendizado de máquina que utiliza
uma estrutura hierárquica para tomar decisões baseadas em condições
lógicas. Neste estudo, utilizamos a função rpart do pacote
rpart para construir árvores de decisão, que são adequadas
para capturar relações não lineares entre os marcadores genéticos e os
fenótipos.
# Configuração do cluster para paralelização
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
# Combinações de traits, repetições e dobras
combos_dt <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
# Execução da validação cruzada para Decision Tree
results_cv_DT <- foreach(
idx = seq_len(nrow(combos_dt)),
.packages = c("rpart", "dplyr", "tibble", "caret"),
.combine = bind_rows
) %dopar% {
# Extração dos parâmetros da combinação atual
tr <- combos_dt$trait[idx]
r <- combos_dt$rep[idx]
f <- combos_dt$fold[idx]
# Preparação dos dados para o modelo
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
# Criando dados de treinamento e validação
yT <- pheno[train_idx, tr][[1]]
xT <- as.matrix(geno[train_idx, ])
yV <- pheno[valid_idx, tr][[1]]
xV <- as.matrix(geno[valid_idx, ])
# Ajuste do modelo Decision Tree
y <- pheno[[tr]]
x <- geno
# Ajuste do modelo Decision Tree
fit <- rpart(y ~ as.matrix(x), subset = train_idx)
# Predição nos conjuntos de treinamento e validação
pred_T <- predict(fit)
pred_V <- predict(fit, newdata = x)[valid_idx]
# Importância das variáveis
imp_df <- data.frame(varImp(fit, scale = TRUE, value = "rss"))
rownames(imp_df) <- gsub("^as\\.matrix\\(x\\)", "", rownames(imp_df), perl = TRUE)
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
# Cria dataframe com a importância das variáveis para os marcadores não selecionados na arvore
imp_df_setdiff <- data.frame(
Overall = 0,
marker = setdiff(colnames(x), imp_df$marker),
Rep = r,
Fold = f
)
imp_df <- imp_df %>%
full_join(imp_df_setdiff)
# Limpeza de memória
gc()
# Criação do dataframe de resultados
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 3),
R2_Valid = round(cor(pred_V, yV), 3),
RMSE_Train = round(sqrt(mean((pred_T - yT )^2)), 3),
RMSE_Valid = round(sqrt(mean((pred_V - yV)^2)), 3),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df)
)
}
stopCluster(cl)
saveRDS(results_cv_DT, file = "output/real_results_cv_DT.rds")O Bagging (Bootstrap Aggregating) é uma técnica de ensemble que
combina múltiplos modelos treinados em subconjuntos aleatórios dos dados
para melhorar a robustez e a precisão das previsões. Neste estudo,
utilizamos o pacote randomForest para implementar o
Bagging, ajustando modelos de floresta aleatória com todos os marcadores
genéticos.
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
combos_bag <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
results_cv_Bag <- foreach(
idx = seq_len(nrow(combos_bag)),
.packages = c("randomForest", "dplyr", "tibble", "caret"),
.combine = bind_rows
) %dopar% {
tr <- combos_bag$trait[idx]
r <- combos_bag$rep[idx]
f <- combos_bag$fold[idx]
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
yT <- pheno[[tr]][train_idx] # vetor
xT <- geno[train_idx, ]
yV <- pheno[[tr]][valid_idx]
xV <- geno[valid_idx, ]
fit <- randomForest(x = xT, y = yT, mtry = ncol(xT))
pred_T <- predict(fit, xT)
pred_V <- predict(fit, xV)
imp_df <- data.frame(varImp(fit, scale = TRUE, value = "rss"))
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
gc()
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 3),
R2_Valid = round(cor(pred_V, yV), 3),
RMSE_Train = round(sqrt(mean((pred_T - yT)^2)), 3),
RMSE_Valid = round(sqrt(mean((pred_V - yV)^2)), 3),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df)
)
}
stopCluster(cl)
saveRDS(results_cv_Bag, file = "output/real_results_cv_Bag.rds")O Random Forest é um método de aprendizado de máquina que combina
múltiplas árvores de decisão para melhorar a precisão e reduzir o
overfitting. Neste estudo, utilizamos o pacote randomForest
para implementar o Random Forest, ajustando modelos com um número
reduzido de marcadores genéticos.
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
combos_rf <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
results_cv_RF <- foreach(
idx = seq_len(nrow(combos_rf)),
.packages = c("randomForest", "dplyr", "tibble", "caret"),
.combine = bind_rows
) %dopar% {
tr <- combos_rf$trait[idx]
r <- combos_rf$rep[idx]
f <- combos_rf$fold[idx]
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
yT <- pheno[[tr]][train_idx] # vetor
xT <- geno[train_idx, ]
yV <- pheno[[tr]][valid_idx]
xV <- geno[valid_idx, ]
fit <- randomForest(xT, yT, mtry = floor(ncol(xT) / 3), ntree = 500)
pred_T <- predict(fit, xT)
pred_V <- predict(fit, xV)
imp_df <- data.frame(varImp(fit, scale = TRUE, value = "rss"))
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
gc()
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 3),
R2_Valid = round(cor(pred_V, yV), 3),
RMSE_Train = round(sqrt(mean((pred_T - yT)^2)), 3),
RMSE_Valid = round(sqrt(mean((pred_V - yV)^2)), 3),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df)
)
}
stopCluster(cl)
saveRDS(results_cv_RF, file = "output/real_results_cv_RF.rds")O Boosting é uma técnica de aprendizado de máquina que combina
múltiplos modelos fracos para criar um modelo forte. Neste estudo,
utilizamos o pacote gbm para implementar o Boosting,
ajustando modelos com todos os marcadores genéticos.
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
combos_gbm <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
stringsAsFactors = FALSE
)
results_cv_GBM <- foreach(
idx = seq_len(nrow(combos_gbm)),
.packages = c("gbm", "dplyr", "tibble", "caret"),
.combine = bind_rows
) %dopar% {
tr <- combos_gbm$trait[idx]
r <- combos_gbm$rep[idx]
f <- combos_gbm$fold[idx]
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
yT <- pheno[[tr]][train_idx] # vetor
xT <- geno[train_idx, ]
yV <- pheno[[tr]][valid_idx]
xV <- geno[valid_idx, ]
y <- pheno[[tr]]
x <- geno
fit <-
gbm.fit(
x = xT,
y = yT,
distribution = "gaussian",
n.trees = 500,
interaction.depth = 2
)
pred_T <- predict(fit, xT, n.trees = 500)
pred_V <- predict(fit, xV, n.trees = 500)
imp_df <- data.frame(varImp(fit, scale = TRUE, value = "rss" , numTrees = 500))
# Adiciona informações de repetição e fold, mantendo nomes dos marcadores
imp_df$marker <- rownames(imp_df)
imp_df$Rep <- r
imp_df$Fold <- f
gc()
tibble(
Trait = tr,
Rep = r,
Fold = f,
R2_Train = round(cor(pred_T, yT), 3),
R2_Valid = round(cor(pred_V, yV), 3),
RMSE_Train = round(sqrt(mean((pred_T - yT)^2)), 3),
RMSE_Valid = round(sqrt(mean((pred_V - yV)^2)), 3),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V),
imp = list(imp_df)
)
}
stopCluster(cl)
saveRDS(results_cv_GBM, file = "output/real_results_cv_GBM.rds")O Perceptron Multicamadas (PMC) é uma rede neural feedforward que
pode capturar relações complexas entre os marcadores genéticos e os
fenótipos. Neste estudo, utilizamos o pacote neuralnet para
implementar o PMC, ajustando modelos com diferentes configurações de
camadas ocultas.
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
combos_nn <- expand.grid(
trait = traits,
rep = seq_len(nrept),
fold = seq_len(nfolds),
h1 = c(1, 2, 4, 8, 16, 32),
h2 = c(0, 1, 2, 4, 8, 16, 32),
stringsAsFactors = FALSE
)
results_cv_NN <- foreach(
idx = seq_len(nrow(combos_nn)),
.packages = c("ANN2", "neuralnet", "dplyr", "tibble"),
.combine = bind_rows
) %dopar% {
tr <- combos_nn$trait[idx]
r <- combos_nn$rep[idx]
f <- combos_nn$fold[idx]
h1 <- combos_nn$h1[idx]
h2 <- combos_nn$h2[idx]
hidden.layers <- case_when(h2 != 0 ~ c(h1, h2),
TRUE ~ h1)
valid_idx <- fold.list[[r]][[f]]
train_idx <- setdiff(seq_len(nrow(pheno)), valid_idx)
yT <- pheno[[tr]][train_idx] # vetor
xT <- geno[train_idx, ]
yV <- pheno[[tr]][valid_idx]
xV <- geno[valid_idx, ]
# Construção e Treinamento da Rede Neural
NN <- neuralnetwork(
X = xT,
y = yT,
hidden.layers = hidden.layers,
regression = TRUE,
activ.functions = "relu",
n.epochs = 500,
loss.type = "squared",
val.prop = 0
)
pred_T <- predict(NN, xT)$predictions
pred_V <- predict(NN, xV)$predictions
tibble(
Trait = tr,
Rep = r,
Fold = f,
Hidden_Layers_1 = h1,
Hidden_Layers_2 = h2,
R2_Train = round(cor(pred_T, yT), 3),
R2_Valid = round(cor(pred_V, yV), 3),
RMSE_Train = round(sqrt(mean((pred_T - yT)^2)), 3),
RMSE_Valid = round(sqrt(mean((pred_V - yV)^2)), 3),
GEBV_T = list(pred_T),
GEBV_V = list(pred_V)
)
}
stopCluster(cl)
saveRDS(results_cv_NN, file = "output/real_results_cv_NN.rds")Os resultados consolidados da validação cruzada para todos os métodos e traços são apresentados na tabela a seguir. Esta tabela resume as métricas de desempenho (R² e RMSE) para cada método, agrupadas por traço.
# Carregar resultados
results_cv_GBLUP_a <- readRDS("output/real_results_cv_GBLUP_a.rds")
results_cv_GBLUP_ad <- readRDS("output/real_results_cv_GBLUP_ad.rds")
results_cv_GBLUP_ade <- readRDS("output/real_results_cv_GBLUP_ade.rds")
results_cv_MARS_1 <- readRDS("output/real_results_cv_MARS_1.rds")
results_cv_MARS_2 <- readRDS("output/real_results_cv_MARS_2.rds")
results_cv_MARS_3 <- readRDS("output/real_results_cv_MARS_3.rds")
results_cv_DT <- readRDS("output/real_results_cv_DT.rds")
results_cv_Bag <- readRDS("output/real_results_cv_Bag.rds")
results_cv_RF <- readRDS("output/real_results_cv_RF.rds")
results_cv_GBM <- readRDS("output/real_results_cv_GBM.rds")
results_cv_NN <- readRDS("output/real_results_cv_NN.rds")
# Agrupa por Trait e configuração de neurônios
resumo_nn <- results_cv_NN %>%
group_by(Trait, Hidden_Layers_1, Hidden_Layers_2) %>%
summarise(
R2_Valid_Mean = mean(R2_Valid, na.rm = TRUE),
R2_Train_Mean = mean(R2_Train, na.rm = TRUE),
RMSE_Valid_Mean = mean(RMSE_Valid, na.rm = TRUE),
RMSE_Train_Mean = mean(RMSE_Train, na.rm = TRUE),
R2_Valid_SD = sd(R2_Valid, na.rm = TRUE),
R2_Train_SD = sd(R2_Train, na.rm = TRUE),
RMSE_Valid_SD = sd(RMSE_Valid, na.rm = TRUE),
RMSE_Train_SD = sd(RMSE_Train, na.rm = TRUE),
.groups = "drop"
)
# Exibe os melhores por característica com base no maior R²_Valid
melhores_config_por_trait <- resumo_nn %>%
group_by(Trait) %>%
slice_max(order_by = R2_Valid_Mean, n = 1) %>%
arrange(Trait)
print(melhores_config_por_trait)# A tibble: 3 × 11
# Groups: Trait [3]
Trait Hidden_Layers_1 Hidden_Layers_2 R2_Valid_Mean R2_Train_Mean
<chr> <dbl> <dbl> <dbl> <dbl>
1 BM 1 8 0.0749 0.706
2 Cer 8 0 0.0752 0.979
3 Prod.Cor 1 2 0.0739 0.692
# ℹ 6 more variables: RMSE_Valid_Mean <dbl>, RMSE_Train_Mean <dbl>,
# R2_Valid_SD <dbl>, R2_Train_SD <dbl>, RMSE_Valid_SD <dbl>,
# RMSE_Train_SD <dbl># Adicionar coluna Method e combinar todos
combined_results <- bind_rows(
results_cv_GBLUP_a %>% mutate(Method = "GBLUP_A"),
results_cv_GBLUP_ad %>% mutate(Method = "GBLUP_AD"),
results_cv_GBLUP_ade %>% mutate(Method = "GBLUP_ADE"),
results_cv_MARS_1 %>% mutate(Method = "MARS_1"),
results_cv_MARS_2 %>% mutate(Method = "MARS_2"),
results_cv_MARS_3 %>% mutate(Method = "MARS_3"),
results_cv_DT %>% mutate(Method = "DT"),
results_cv_Bag %>% mutate(Method = "BAG"),
results_cv_RF %>% mutate(Method = "RF"),
results_cv_GBM %>% mutate(Method = "GBM")
)
# Resumir por Method e Trait
table_final <- combined_results %>%
group_by(Method, Trait) %>%
summarise(
R2_Valid_Mean = mean(R2_Valid, na.rm = TRUE),
R2_Valid_SD = sd(R2_Valid, na.rm = TRUE),
RMSE_Valid_Mean = mean(RMSE_Valid, na.rm = TRUE),
RMSE_Valid_SD = sd(RMSE_Valid, na.rm = TRUE)
) %>%
ungroup() %>%
bind_rows(
melhores_config_por_trait %>%
mutate(Method = "NN") %>%
dplyr::select(Method, Trait, R2_Valid_Mean, R2_Valid_SD, RMSE_Valid_Mean, RMSE_Valid_SD)
)
# Exibir tabela consolidada
table_final %>%
arrange(Trait, Method) %>%
mutate_if(is.numeric, round, 4) %>%
kable(booktabs = TRUE, caption = "Resultados Consolidados por Método e Traço") %>%
kable_styling(full_width = FALSE)| Method | Trait | R2_Valid_Mean | R2_Valid_SD | RMSE_Valid_Mean | RMSE_Valid_SD |
|---|---|---|---|---|---|
| BAG | BM | -0.0950 | 0.0794 | 0.5299 | 0.0419 |
| DT | BM | -0.0385 | 0.1207 | 0.6761 | 0.0564 |
| GBLUP_A | BM | NaN | NA | 2.1107 | 0.0667 |
| GBLUP_AD | BM | -0.0683 | 0.1041 | 2.1114 | 0.0654 |
| GBLUP_ADE | BM | -0.0805 | NA | 2.1100 | 0.0668 |
| GBM | BM | -0.0886 | 0.1012 | 0.5210 | 0.0435 |
| MARS_1 | BM | -0.0341 | 0.1395 | 0.8444 | 0.1155 |
| MARS_2 | BM | -0.0766 | 0.1189 | 0.8630 | 0.0709 |
| MARS_3 | BM | -0.0431 | 0.1263 | 0.8303 | 0.0608 |
| NN | BM | 0.0749 | 0.1023 | 2.7595 | 1.9837 |
| RF | BM | -0.1135 | 0.0806 | 0.5305 | 0.0420 |
| BAG | Cer | -0.0966 | 0.1355 | 0.5937 | 0.0346 |
| DT | Cer | 0.0207 | 0.1258 | 0.7295 | 0.0513 |
| GBLUP_A | Cer | -0.0287 | 0.0244 | 2.5420 | 0.0399 |
| GBLUP_AD | Cer | 0.0513 | 0.1439 | 2.5484 | 0.0542 |
| GBLUP_ADE | Cer | 0.1348 | 0.2152 | 2.5411 | 0.0411 |
| GBM | Cer | -0.0778 | 0.0997 | 0.5767 | 0.0307 |
| MARS_1 | Cer | -0.0458 | 0.1353 | 0.9167 | 0.1413 |
| MARS_2 | Cer | -0.0255 | 0.0986 | 1.0138 | 0.1573 |
| MARS_3 | Cer | -0.0294 | 0.1014 | 0.9823 | 0.1339 |
| NN | Cer | 0.0752 | 0.1185 | 3.3241 | 2.8533 |
| RF | Cer | -0.0949 | 0.1469 | 0.5939 | 0.0333 |
| BAG | Prod.Cor | -0.0318 | 0.1047 | 3.0902 | 0.3069 |
| DT | Prod.Cor | -0.0826 | 0.1106 | 3.9523 | 0.2953 |
| GBLUP_A | Prod.Cor | -0.0071 | 0.1375 | 5.9770 | 0.4122 |
| GBLUP_AD | Prod.Cor | -0.0196 | 0.1281 | 5.9827 | 0.4052 |
| GBLUP_ADE | Prod.Cor | -0.0202 | 0.1130 | 5.9805 | 0.4064 |
| GBM | Prod.Cor | 0.0155 | 0.0761 | 3.0041 | 0.2998 |
| MARS_1 | Prod.Cor | -0.0229 | 0.1342 | 4.6685 | 0.7683 |
| MARS_2 | Prod.Cor | -0.0283 | 0.1303 | 5.0639 | 0.6063 |
| MARS_3 | Prod.Cor | -0.0146 | 0.1318 | 5.0429 | 0.6617 |
| NN | Prod.Cor | 0.0739 | 0.1045 | 8.4643 | 6.7590 |
| RF | Prod.Cor | -0.0248 | 0.0779 | 3.0785 | 0.3048 |
Nenhum método se sobressaiu em todas as características. Os ensembles de árvores (GBM, BAG, RF) ofereceram consistência e os menores erros médios em BM e Cer, mas não capturaram a variabilidade genética (R² negativas). As redes neurais, por sua vez, foram as únicas a gerar R² positivas em todos os traços, porém sacrificaram precisão, refletida em RMSE elevados. Os modelos GBLUP apresentaram acurácia limitada e alto erro em todas as análises. Esses resultados sugerem que combinações de abordagens — por exemplo, integrando as previsões de árvores e redes neurais ou incorporando covariáveis genômicas adicionais — podem ser necessárias para otimizar a seleção genômica em Coffea arabica.
Neste estudo, avaliamos o desempenho de diferentes métodos de predição genômica em um conjunto de dados real de Cofea arabica. Os resultados indicam que não há um único método que se destaque em todas as características fenotípicas. Os ensembles de árvores, como GBM, BAG e RF, mostraram-se consistentes e com menores erros médios, mas falharam em capturar a variabilidade genética. As redes neurais foram as únicas a gerar R² positivos em todos os traços, embora com RMSE elevados. Os modelos GBLUP apresentaram acurácia limitada e alto erro.
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=Portuguese_Brazil.utf8 LC_CTYPE=Portuguese_Brazil.utf8
[3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C
[5] LC_TIME=Portuguese_Brazil.utf8
time zone: America/Sao_Paulo
tzcode source: internal
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
[4] dplyr_1.1.4 purrr_1.1.0 readr_2.1.5
[7] tidyr_1.3.1 tibble_3.3.0 tidyverse_2.0.0
[10] caret_7.0-1 ggplot2_3.5.2 ANN2_2.3.4
[13] neuralnet_1.44.2 gbm_2.2.2 randomForest_4.7-1.2
[16] rpart_4.1.23 earth_5.3.4 plotmo_3.6.4
[19] plotrix_3.8-4 Formula_1.2-5 sommer_4.4.2
[22] crayon_1.5.3 MASS_7.3-60.2 Matrix_1.7-0
[25] doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2
[28] kableExtra_1.4.0 cvTools_0.3.3 robustbase_0.99-4-1
[31] lattice_0.22-7
loaded via a namespace (and not attached):
[1] pROC_1.18.5 writexl_1.5.4 readxl_1.4.5
[4] rlang_1.1.6 magrittr_2.0.3 git2r_0.36.2
[7] compiler_4.4.1 systemfonts_1.2.3 vctrs_0.6.5
[10] reshape2_1.4.4 pkgconfig_2.0.3 fastmap_1.2.0
[13] utf8_1.2.6 promises_1.3.3 rmarkdown_2.29
[16] prodlim_2025.04.28 tzdb_0.5.0 xfun_0.52
[19] cachem_1.1.0 jsonlite_2.0.0 recipes_1.3.1
[22] later_1.4.2 R6_2.6.1 bslib_0.9.0
[25] stringi_1.8.7 RColorBrewer_1.1-3 parallelly_1.45.0
[28] cellranger_1.1.0 jquerylib_0.1.4 Rcpp_1.1.0
[31] knitr_1.50 future.apply_1.20.0 httpuv_1.6.16
[34] splines_4.4.1 nnet_7.3-19 timechange_0.3.0
[37] tidyselect_1.2.1 rstudioapi_0.17.1 yaml_2.3.10
[40] timeDate_4041.110 codetools_0.2-20 listenv_0.9.1
[43] plyr_1.8.9 withr_3.0.2 evaluate_1.0.4
[46] future_1.58.0 survival_3.6-4 xml2_1.3.8
[49] pillar_1.11.0 whisker_0.4.1 stats4_4.4.1
[52] generics_0.1.4 rprojroot_2.0.4 hms_1.1.3
[55] scales_1.4.0 globals_0.18.0 class_7.3-22
[58] glue_1.8.0 tools_4.4.1 data.table_1.17.8
[61] ModelMetrics_1.2.2.2 gower_1.0.2 fs_1.6.6
[64] grid_4.4.1 ipred_0.9-15 nlme_3.1-168
[67] cli_3.6.5 textshaping_1.0.1 workflowr_1.7.1
[70] viridisLite_0.4.2 svglite_2.2.1 lava_1.8.1
[73] gtable_0.3.6 DEoptimR_1.1-3-1 sass_0.4.10
[76] digest_0.6.37 farver_2.1.2 htmltools_0.5.8.1
[79] lifecycle_1.0.4 hardhat_1.4.1 Weverton Gomes da Costa, Pós-Doutorando, Departamento de Estatística - Universidade Federal de Viçosa, wevertonufv@gmail.com↩︎