#-----------------------------------------------------------------------# Pacotes.library(reactable)library(segmented)library(tidyverse)library(broom)library(janitor)#-----------------------------------------------------------------------# Tema.custom_theme <-function() {theme_linedraw() +theme(panel.grid =element_blank(),strip.background =element_rect(fill ="gray90",color ="black"),strip.text =element_text(color ="black") )}# Configure globalmente o tema personalizadotheme_set(custom_theme())
2 Importação e preparo dos dados
2.1 Importação
Código
#-----------------------------------------------------------------------# Importação e preparo dos dados.# Importa todas as abas.url <-"Cmel_nym_chr.xlsx"tbs <-map(1:3, ~readxl::read_excel(url, sheet = .))# Converte variáveis para fator.tbs <-map(tbs, ~mutate_at(., c("isolate", "species"), factor))# str(tbs)# NOTE: Cada experimento avaliou apenas uma espécie. Juntar em uma única# tabela.tb <-bind_rows(tbs)# Doses usadas.tb |>distinct(dose) |>mutate(dose_log10 =log10(dose))
# Ordena isolados dentro de espécie.tb <- tb |>arrange(species, isolate) |>mutate(isolate =factor(isolate, levels =unique(isolate)))
2.2 Transformação da escala da dose para ajuste dos modelos
Código
#-----------------------------------------------------------------------# Transformação na escala da dose para ajuste das curvas.# Experimenta transformações para ficar próximo de um equidistante pois# isso é interessante para o ajuste de modelos polinomiais.shift <-0.01tb |>distinct(dose) |>mutate(log_dose =log10(dose + shift) +2)
# IMPORTANT: Praticamente equidistante. O valor `0` é a testemunha. Os# valores de dose ficam positivos ao somar 2. Assim o intercepto será o# crescimento na dose 0. Sabendo isso, é resolver em x para saber qual# valor em x que retorna 0.5 * y(x = 0).# Cria uma nova coluna com a dose na escala logarítmica.tb <- tb |>mutate(log_dose =log10(dose + shift) +2)
Figura 1. Análise exploratória para o crescimento de isolados de Colletotrichum em função da dose de fungicida. Curvas ajustadas considerando um polinômio de terceiro grau.
#-----------------------------------------------------------------------# Ajusta os modelos e calcula EC50.# Calcula ec50 para cada isolado usando o modelo linear.k <-1tb_ec50_1 <- tb |>group_by(species, isolate) |>transmute(x = log_dose, y = growth) |>summarise(result =list(find_ec50_polyroot(tibble(x = x, y = y), k = k)),.groups ="drop") |>unnest(result)tb_ec50_1
# A tibble: 71 × 6
species isolate x y y0 model
<fct> <fct> <dbl> <dbl> <dbl> <list>
1 C. melonis 107_J 4.74 0.552 0.552 <lm>
2 C. melonis 123_J 3.34 0.856 0.856 <lm>
3 C. melonis 128_J 5.47 1.02 1.02 <lm>
4 C. melonis 129_J 2.62 0.870 0.870 <lm>
5 C. melonis 133_J 7.03 0.933 0.933 <lm>
6 C. melonis 137_J 5.36 2.26 2.26 <lm>
7 C. melonis 21_J 2.92 0.555 0.555 <lm>
8 C. melonis 28_J 1.77 0.216 0.216 <lm>
9 C. melonis 29_SP 5.42 0.785 0.785 <lm>
10 C. melonis 32_J 2.39 0.317 0.317 <lm>
# ℹ 61 more rows
Código
# Calcula ec50 para cada isolado usando o modelo cúbico.k <-3tb_ec50_3 <- tb |>group_by(species, isolate) |>transmute(x = log_dose, y = growth) |>summarise(result =list(find_ec50_polyroot(tibble(x = x, y = y), k = k)),.groups ="drop") |>unnest(result)tb_ec50_3
# A tibble: 147 × 6
species isolate x y y0 model
<fct> <fct> <dbl> <dbl> <dbl> <list>
1 C. melonis 107_J 4.03 0.572 0.572 <lm>
2 C. melonis 123_J 3.38 0.901 0.901 <lm>
3 C. melonis 128_J -3.66 1.01 1.01 <lm>
4 C. melonis 129_J -1.26 0.867 0.867 <lm>
5 C. melonis 129_J 2.25 0.867 0.867 <lm>
6 C. melonis 129_J 4.48 0.867 0.867 <lm>
7 C. melonis 133_J 4.82 0.964 0.964 <lm>
8 C. melonis 137_J -6.21 2.35 2.35 <lm>
9 C. melonis 21_J -1.41 0.543 0.543 <lm>
10 C. melonis 21_J 2.62 0.543 0.543 <lm>
# ℹ 137 more rows
Código
# Juntar os dois ajustes para fazer o restante dos processos.tb_ec50 <-bind_rows(list("1"= tb_ec50_1,"3"= tb_ec50_3),.id ="order")# str(tb_ec50, max.level = 2, give.attr = FALSE)#-----------------------------------------------------------------------# Elimina soluções inválidas para EC50.# Valida os valores de EC50 na escala em que foi estimada, ou seja, `x`.tb_ec50 <- tb_ec50 |>group_by(order, species, isolate) |># y e y0 devem ser próximos e x deve estar entre 0 e 8.mutate(x =if_else(between(x, 0, 8), x, NA_real_)) |># Pega o menor valor quando houver vários válidos.mutate(x =if_else(x ==min(x, na.rm =TRUE), x, NA_real_)) |>ungroup()# Aplica a função inversa para obter a EC50.tb_ec50 <- tb_ec50 |>mutate(ec50 = x -2,ec50 =10^ec50 - shift)# Quais isolados não tiveram EC50 estimada?tb_ec50 |>group_by(order, species, isolate) |>summarise(n =n(),n_ec50 =sum(!is.na(ec50)),.groups ="drop") |>filter(n_ec50 ==0) |># droplevels() |># pull(isolate) |># as.character() |># dput() |>identity()
# A tibble: 4 × 5
order species isolate n n_ec50
<chr> <fct> <fct> <int> <int>
1 1 C. crisophylum 55_PR 1 0
2 3 C. melonis 128_J 1 0
3 3 C. melonis 137_J 1 0
4 3 C. nymphae 20_RS 1 0
Figura 2. EC50 para os isolados de Colletotrichum a partir dos ajuste dos modelos polinomiais de order 1 (linear) e 3 (cúbico) aos dados de crescimento em função da dose de fungicida. Linhas verticais indicam o valor da EC50 e linhas horizontais indicam o crescimento estimado quando a dose é a EC50.
3.4 Ajustes finais
Código
#-----------------------------------------------------------------------# Pegar o ajuste mais apropriado.# Escolhe o ajuste.tb_ec50_final <- tb_ec50 |>drop_na(ec50) |>group_by(species, isolate) |>arrange(order) |>do({if (nrow(.) <=1) {# Se tem uma linha ou menos, retornar isso. . } else {# Usa o LRT para decidir qual modelo usar. p_value <-anova(.$model[[1]], .$model[[2]])[2, "Pr(>F)"]if (p_value <0.05) { .[2, ] } else { .[1, ] } } }) |>ungroup()tb_ec50_final
# A tibble: 71 × 8
order species isolate x y y0 model ec50
<chr> <fct> <fct> <dbl> <dbl> <dbl> <list> <dbl>
1 1 C. melonis 107_J 4.74 0.552 0.552 <lm> 545.
2 1 C. melonis 123_J 3.34 0.856 0.856 <lm> 21.6
3 1 C. melonis 128_J 5.47 1.02 1.02 <lm> 2954.
4 1 C. melonis 129_J 2.62 0.870 0.870 <lm> 4.15
5 1 C. melonis 133_J 7.03 0.933 0.933 <lm> 108166.
6 1 C. melonis 137_J 5.36 2.26 2.26 <lm> 2309.
7 1 C. melonis 21_J 2.92 0.555 0.555 <lm> 8.30
8 3 C. melonis 28_J 1.82 0.208 0.208 <lm> 0.654
9 3 C. melonis 29_SP 5.37 0.830 0.830 <lm> 2369.
10 3 C. melonis 32_J 3.19 0.308 0.308 <lm> 15.6
# ℹ 61 more rows
Código
# Filtra os pontos.tb_ec50_pred <-left_join(tb_ec50_final[, c("order", "species", "isolate")], tb_ec50_pred)# Exibir ordenando pela EC50 dentro de cada espécie.tb_ec50_final <- tb_ec50_final |>arrange(species, ec50) |>mutate(isolate =as.character(isolate),isolate =factor(isolate,levels = isolate))tb <- tb |>mutate(isolate =factor(isolate,levels = tb_ec50_final$isolate))# Cria os rótulos para exibir no gráfico.tb_ec50_final <- tb_ec50_final |>mutate(ec50_label =sprintf(fmt =case_when(ec50 <0.1~"EC[50] == %0.4f", ec50 <1~"EC[50] == %0.2f", ec50 <100~"EC[50] == %0.1f",TRUE~"EC[50] == %0.0f"), ec50))
Figura 3. EC50 para os isolados de Colletotrichum de acordo com o modelo polinomial mais apropriado quando ajustado aos dados de crescimento em função da dose de fungicida. Linhas verticais indicam o valor da EC50 e linhas horizontais indicam o crescimento estimado quando a dose é a EC50.
Figura 4. EC50 para os isolados de Colletotrichum de acordo com o modelo polinomial mais apropriado quando ajustado aos dados de crescimento em função da dose de fungicida. Linhas verticais indicam o valor da EC50 e linhas horizontais indicam o crescimento estimado quando a dose é a EC50.
Tabela 2. Estimativas de EC50 para os isolados de Colletotrichum acompanhados de medidas de ajuste para os modelos ajustados aos dados de crescimento em função da dose de fungicida. ec50 é a estimativa da concentração para obter a metade do crescimento estimado na ausência de fungicida (EC50), growth50 é o crescimento estimado para a EC50, order corresponde a ordem do polinômio empregado (1 = linear, 3 = cúbico), r_squared é o coeficiente de determinação ajustado, p_value é o valor-p do teste F para o hipótese nula de que o modelo não difere do modelo trivial e stars são sinais gráficos para indicar a classe de significância.
