bGMYC4 Interactive Workflow

compiler::enableJIT(3) # JIT acceleration / Ускорение JIT
library(ape)
library(bGMYC4)
library(treeio)
library(ggtree)
library(dplyr)
library(plotly)
library(htmlwidgets)
library(future)
library(future.apply)
library(mcmcse)

# EN: Safe numeric scalar input / Безопасный ввод скалярных значений
input_scalar <- function(label, default_val) {
  if (!interactive()) return(default_val)
  val <- readline(sprintf("  %s (default / по умолчанию: %s): ", label, default_val))
  if (nchar(trimws(val)) == 0) return(default_val)
  num <- suppressWarnings(as.numeric(val))
  if (is.na(num)) return(default_val)
  return(num)
}

# EN: Safe numeric vector input / Безопасный ввод векторов
input_vector <- function(label, default_vec) {
  if (!interactive()) return(default_vec)
  val <- readline(sprintf("  %s (comma-separated, default / через запятую, по умолчанию: %s): ", 
                          label, paste(default_vec, collapse = ", ")))
  if (nchar(trimws(val)) == 0) return(default_vec)
  nums <- suppressWarnings(as.numeric(unlist(strsplit(val, "[,;\\s]+"))))
  nums <- nums[!is.na(nums)]
  if (length(nums) != length(default_vec)) return(default_vec)
  return(nums)
}

# EN: Enforce ultrametricity / Обеспечение ультраметричности
fix_ultrametric <- function(tr) {
  if (!is.ultrametric(tr)) {
    tr$edge.length <- round(tr$edge.length, 8)
    if (!is.ultrametric(tr)) stop("❌ Tree remains non-ultrametric / Дерево остаётся неультраметричным.")
  }
  return(tr)
}

consensus_path <- if (interactive()) readline("📥 Path to consensus tree / Путь к консенсусному дереву (.tree): ") else "dummy.tree"
posterior_path <- if (interactive()) readline("📥 Path to posterior trees / Путь к набору деревьев (.trees): ") else "dummy.trees"

tree_beast <- treeio::read.beast(consensus_path)
tree_consensus <- tree_beast@phylo
all_trees <- read.nexus(posterior_path)
class(all_trees) <- "multiPhylo"

# EN: Analysis Mode Selection / Выбор режима анализа
mode_str <- if (interactive()) readline("Mode 1 (Single) or 2 (Multi)? / Режим 1 или 2? [2]: ") else "2"
analysis_mode <- ifelse(nchar(trimws(mode_str)) == 0, 2, as.integer(mode_str))

# EN: Tree sampling with burnin / Выборка деревьев с учетом burnin
if (analysis_mode == 2) {
  n_total <- length(all_trees)
  burnin_idx <- floor(n_total * 0.10) # 10% burnin
  n_sample <- input_scalar("Number of trees to sample / Кол-во деревьев", 10)
  set.seed(42)
  trees_sample <- all_trees[sample((burnin_idx + 1):n_total, n_sample)]
  class(trees_sample) <- "multiPhylo"
} else {
  trees_sample <- NULL
}

# EN: Outgroup removal via regex patterns / Удаление аутгрупп по паттернам
if (interactive() && tolower(readline("Drop outgroups? / Удалять аутгруппы? (y/n): ")) == "y") {
  og_str <- readline("Enter patterns (comma-separated) / Введите паттерны: ")
  og_patterns <- trimws(unlist(strsplit(og_str, "[,;]+")))
  matching_tips <- unique(unlist(lapply(og_patterns, function(p) {
    pattern <- paste0("(^|[^A-Za-z0-9])", p, "($|[^A-Za-z0-9])")
    tree_consensus$tip.label[grepl(pattern, tree_consensus$tip.label, ignore.case = TRUE, perl = TRUE)]
  })))
  
  if (length(matching_tips) > 0) {
    tree_consensus <- drop.tip(tree_consensus, matching_tips)
    if (analysis_mode == 2) {
      trees_sample <- lapply(trees_sample, drop.tip, tip = matching_tips)
      class(trees_sample) <- "multiPhylo"
    }
    tree_beast <- treeio::drop.tip(tree_beast, matching_tips)
  }
}

tree_consensus <- fix_ultrametric(tree_consensus)
ntips <- length(tree_consensus$tip.label)

params <- list(
  mcmc = 10000, burnin = 1000, thinning = 10,
  py1 = 0, py2 = 1.5, pc1 = 0, pc2 = 2,
  t1 = 2, t2 = min(35, ntips - 1),
  scale = c(20, 10, 5), start = c(1, 1, floor(ntips/3))
)

# EN: Diagnostic tuning loop / Цикл диагностики и настройки
repeat {
  res_single <- bgmyc.singlephy(
    phylo = tree_consensus, mcmc = params$mcmc, burnin = params$burnin, 
    thinning = params$thinning, py1 = params$py1, py2 = params$py2, 
    pc1 = params$pc1, pc2 = params$pc2, t1 = params$t1, t2 = params$t2, 
    scale = params$scale, start = params$start
  )
  
  # EN: Convergence checks / Проверка сходимости
  ar <- res_single$accept
  cat(sprintf("Acceptance rates: py=%.3f | pc=%.3f | th=%.3f\n", ar[1], ar[2], ar[3]))
  if (requireNamespace("mcmcse", quietly = TRUE)) {
    ess_vals <- sapply(1:4, function(col) round(mcmcse::ess(res_single$par[, col])))
    cat(sprintf("ESS: py=%d | pc=%d | th=%d | logL=%d [Target > 200]\n", 
                ess_vals[1], ess_vals[2], ess_vals[3], ess_vals[4]))
  }
  
  plot(res_single)
  
  if (!interactive() || tolower(readline("Accept parameters? / Принять параметры? (y/n): ")) != "n") break
  
  # EN: Update parameters / Обновление параметров
  params$mcmc <- input_scalar("mcmc", params$mcmc)
  params$burnin <- input_scalar("burnin", params$burnin)
  params$thinning <- input_scalar("thinning", params$thinning)
  params$scale <- input_vector("scale", params$scale)
  params$start <- input_vector("start", params$start)
}

if (analysis_mode == 2) {
  # EN: Parallel execution / Параллельное выполнение
  n_workers <- min(parallel::detectCores(logical = FALSE) - 1, length(trees_sample))
  plan(multisession, workers = max(1, n_workers))
  
  final_res <- future_lapply(seq_along(trees_sample), function(i) {
    bgmyc.singlephy(
      phylo = trees_sample[[i]], mcmc = params$mcmc, burnin = params$burnin, 
      thinning = params$thinning, py1 = params$py1, py2 = params$py2, 
      pc1 = params$pc1, pc2 = params$pc2, t1 = params$t1, t2 = params$t2, 
      scale = params$scale, start = params$start
    )
  }, future.seed = TRUE)
  class(final_res) <- "multibgmyc"
  
  # EN: Gelman-Rubin R-hat / Статистика Гелмана-Рубина
  if (requireNamespace("mcmcse", quietly = TRUE) && length(final_res) > 1) {
    chains_list <- lapply(final_res, function(res) res$par[, 3])
    names(chains_list) <- paste0("Tree_", seq_along(final_res))
    rhat <- round(mcmcse::gelman(chains_list)$Rhat, 3)
    cat(sprintf("Gelman-Rubin R-hat: %.3f [Target < 1.05]\n", rhat))
  }
} else {
  final_res <- list(res_single)
  class(final_res) <- "multibgmyc"
}

probmat <- spec.probmat(final_res)

# EN: Synchronize tip order between tree and matrix / Синхронизация порядка таксонов
p_tree <- suppressWarnings(ggtree(tree_beast, layout = "rectangular"))
tips_data <- p_tree$data %>% filter(isTip) %>% arrange(y)
tip_order <- tips_data$label
probmat <- probmat[tip_order, tip_order]

# EN: Extract posterior probabilities for branches / Извлечение posterior для ветвей
post_col <- intersect(c("posterior", "prob", "Posterior"), colnames(p_tree$data))[1]
node_posterior <- p_tree$data[[post_col]]
names(node_posterior) <- p_tree$data$node

# EN: Custom color gradient function / Функция градиента цвета
get_pp_color <- function(pp) {
  if(is.na(pp)) return("#CCCCCC")
  pp <- max(0, min(1, pp))
  colors <- list(c(0.0, 1.0, 0.0, 0.0), c(0.25, 1.0, 0.5, 0.0), 
                 c(0.5, 1.0, 1.0, 0.0), c(0.75, 0.5, 1.0, 0.0), c(1.0, 0.0, 0.7, 0.0))
  for(i in 1:(length(colors)-1)) {
    if(pp >= colors[[i]][1] && pp <= colors[[i+1]][1]) {
      t <- (pp - colors[[i]][1]) / (colors[[i+1]][1] - colors[[i]][1])
      r <- colors[[i]][2] + t * (colors[[i+1]][2] - colors[[i]][2])
      g <- colors[[i]][3] + t * (colors[[i+1]][3] - colors[[i]][3])
      b <- colors[[i]][4] + t * (colors[[i+1]][4] - colors[[i]][4])
      return(sprintf("#%02X%02X%02X", round(r*255), round(g*255), round(b*255)))
    }
  }
  return("#CCCCCC")
}

# EN: Build Plotly Tree / Построение дерева в Plotly
edges <- p_tree$data %>% filter(!is.na(parent))
fig_tree <- plot_ly()
for(i in 1:nrow(edges)) {
  child <- edges[i, ]
  parent <- p_tree$data %>% filter(node == child$parent)
  branch_color <- get_pp_color(node_posterior[as.character(child$node)])
  fig_tree <- fig_tree %>% add_segments(
    x = parent$x, xend = child$x, y = child$y, yend = child$y,
    line = list(color = branch_color, width = 5), showlegend = FALSE)
  fig_tree <- fig_tree %>% add_segments(
    x = parent$x, xend = parent$x, y = parent$y, yend = child$y,
    line = list(color = branch_color, width = 5), showlegend = FALSE)
}

# EN: Build Plotly Heatmap / Построение тепловой карты
fig_heat <- plot_ly(
  z = probmat, x = 1:nrow(probmat), y = 1:ncol(probmat), type = "heatmap",
  colorscale = list(list(0.0, "#F7FCF5"), list(0.5, "#41AB5D"), list(1.0, "#00441B")),
  zmin = 0, zmax = 1, showscale = TRUE
)

# EN: Combine 1:1 / Объединение 1:1
fig_combined <- subplot(fig_tree, fig_heat, nrows = 1, widths = c(0.5, 0.5), shareY = TRUE) %>%
  layout(yaxis = list(autorange = "reversed", showticklabels = FALSE),
         xaxis = list(showticklabels = FALSE),
         xaxis2 = list(showticklabels = FALSE),
         yaxis2 = list(autorange = "reversed", showticklabels = FALSE))

htmlwidgets::saveWidget(fig_combined, "bGMYC_interactive_heatmap.html", selfcontained = TRUE)

# EN: Export clusters at different thresholds / Экспорт кластеров на разных порогах
for (p in c(0.05, 0.01)) {
  out <- bgmyc.point(probmat, ppcutoff = p)
  df <- data.frame(
    Sequence = unlist(out),
    MOTU_bGMYC = rep(seq_along(out), lengths(out)),
    stringsAsFactors = FALSE
  )
  write.table(df, file = sprintf("Delimitation_bGMYC_%.2f.csv", p),
              row.names = FALSE, sep = ";", dec = ".", quote = FALSE, fileEncoding = "UTF-8")
}

# EN: Export full probability matrix / Экспорт полной матрицы вероятностей
spec_out <- bgmyc.spec(final_res)
write.csv(spec_out$specprobs, "bGMYC_full_probs.csv", row.names = FALSE)