Консенсусные сети

Рисуем деревья, используем fancy палитры и радуемся жизни

Библиотеки

library(readr)
library(stylo)
library(tidyverse)
library(purrr)
library(ape) 
library(phangorn)
library(TreeTools)
library(paletteer)

Код

Для начала необходимо представить полученный текстовый файл в табличном формате:

# читаем текстовый файл. Преобразуем в табличный формат.
data <- read.table("sholohov.txt", header = TRUE, row.names = 1)
# транспонируем матрицу, чтобы ряды и столбцы поменялись местами.
sholohov <- t(data)
# превращаем в ДФ
sholohov <- as.data.frame.matrix(sholohov)

Строим консенсусное дерево через stylo и метод BCT:

bct_result <- stylo(gui = FALSE, 
                    frequencies = sholohov,
                    analysis.type = "BCT",
                    mfw.min = 100,
                    mfw.max = 500,
                    mfw.incr = 100,
                    distance.measure = "wurzburg",
                    colors.on.graphs = "greyscale",
                    write.png.file = FALSE,
                    consensus.strength = 0.5,
                    plot.custom.width = 8, 
                    plot.custom.height = 8)

Теперь пишем функцию, которая получает деревья:

# пишем функцию, которая получает деревья
get_tree <- function(df) {
  X <- df[ , sample(3000, replace = F, size = 100)]
  distmx <- dist(scale(X))
  tr <- as.phylo(hclust(distmx))
  tr
}

Создаем зерно, проходимся функцией по данным и красим будущую сеть:

#задаем зерно
set.seed(20022025)
#проходим функцией по ДФ
trees_result <- map(1:15, ~get_tree(sholohov))

#загружаем палитру для будущего графика
pal <- paletteer_d("ggprism::floral")
 
# вычисляем консенсус
mph <- as.multiPhylo(trees_result)
cons.nw <- consensusNet(mph, prob = 0.3, rooted = FALSE)

# раскрашиваем сеть
cons.nw$col <- str_remove_all(cons.nw$tip.label, "_.+")

col_tbl <- tibble(label = unique(cons.nw$col),
       
                             col = pal)
color_group <- tibble(label = cons.nw$col) |> 
  left_join(col_tbl)

cons.nw$col <- color_group$col

Рисуем сеть, применяя загруженную ранее палитру:

# рисуем
set.seed(2002)
par(mar = c(1, 1, 1, 1)) # задаем поля
plot(cons.nw, type = "2D", 
     direction = "horizontal", # названия будут писаться горизонтально
     tip.color = cons.nw$col,#используем палитру
     edge.color = "grey50",
     edge.width = 1, 
     cex = 0.6)

Получившаяся сеть (сделала текст поменьше, иначе получалась совсем уж каша)