ChiLiubio
commented
4 days ago Hi. To compare two types of results and make the steps reproducible, I will use the example data inside the package and shorten the steps as many as possible.
library(microeco)
data(dataset)
data(env_data_16S)
t1 <- trans_env$new(dataset = dataset, add_data = env_data_16S[, 4:11])
t1$cal_ordination(method = "dbRDA", use_measure = "bray")
t1$trans_ordination(adjust_arrow_length = TRUE, max_perc_env = 1.5)
t1$plot_ordination(plot_color = "Group")
# let's check the percentage shown in the axis
raw_dbrda_result <- t1$res_ordination
raw_dbrda_result$CCA$eig/sum(raw_dbrda_result$CCA$eig) * 100
Hi There!
I was performing a dbRDA but notice that the percentage explained by the axis are different doing manually and with the microeco
I will send here the codes and the files
Ordination manual
otu.table.norm <- ps_mesocosmos_norm$otu_table
ps_mesocosmos$sample_table[is.na(ps_mesocosmos$sample_table)] <- 0 amb <- ps_mesocosmos$sample_table[7:21] var_amb <- as.data.frame(scale(amb))
dist_matriz <- vegdist(otu.table.norm, method = "bray")
dist_matriz <- ps_mesocosmos_norm$beta_diversity$bray View(dist_matriz)
dbRDA usando a matriz de distâncias e as variáveis escalonadas
dbrda_result <- dbrda(dist_matriz ~ ., data = var_amb) summary(dbrda_result)
anova(dbrda_result, by = "terms") anova(dbrda_result, by = "axis")
plot(dbrda_result)
site_scores <- scores(dbrda_result, display = "sites") site_scores_df <- as.data.frame(site_scores)
env_scores <- scores(dbrda_result, display = "bp") env_scores_df <- as.data.frame(env_scores)
eigenvalues <- summary(dbrda_result)$cont$importance[2, ] # Linha com % explicada eigenvalues <- eigenvalues * 100 # Converter para porcentagem names(eigenvalues) <- paste0("dbRDA", seq_along(eigenvalues)) # Nomear os eixos
perc_axis <- paste0("dbRDA1 (", round(eigenvalues["dbRDA1"], 2), "%) | dbRDA2 (", round(eigenvalues["dbRDA2"], 2), "%)")
site_scores_df$warming <- ps_mesocosmos$sample_table$warming
test2 <- ggplot(site_scores_df, aes(x = dbRDA1, y = dbRDA2, color = warming)) + geom_point(size = 3) +
geom_text_repel(aes(label = rownames(site_scores_df)), size = 3) +
geom_segment(data = env_scores_df, aes(x = 0, y = 0, xend = dbRDA1, yend = dbRDA2), arrow = arrow(type = "closed", length = unit(0.1, "inches")), color = "black") + geom_text(data = env_scores_df, aes(x = dbRDA1, y = dbRDA2, label = rownames(env_scores_df)), color = "black", size = 3) + labs(x = "db-RDA1", y = "db-RDA2", color = "Waming", subtitle = perc_axis) + theme_minimal() + scale_color_brewer(palette = "Set1")
test2
Teste microeco
ps_mesocosmos_norm$sample_table[is.na(ps_mesocosmos_norm$sample_table)] <- 0 t1 <- trans_env$new(dataset = ps_mesocosmos_norm, env_cols = c(7:21)) t1$cal_ordination(method = "dbRDA", use_measure = "bray") t1$trans_ordination(adjust_arrow_length = TRUE, max_perc_env = 1.5) t1$plot_ordination(plot_color = "warming") t1$cal_ordination_anova() t1$cal_ordination_envfit() t1$res_ordination_terms t1$res_ordination_axis
test1 <- t1$plot_ordination(plot_color = "warming")
Thanks for the attention
tax.txt metadata120.txt otu.txt