Как добавить загрузочные значения узлов кластера (дерева) в формате NEWICK в R

Вот одно из решений: объекты класса phylo имеют доступный слот с именем node.label, который, соответственно, дает вам метку узла. Вы можете использовать его для хранения значений начальной загрузки. В вашем файле Newick File в соответствующем месте будет записано, как вы можете видеть в коде .write.tree2:

> .write.tree2
function (phy, digits = 10, tree.prefix = "") 
{
    brl <- !is.null(phy$edge.length)
    nodelab <- !is.null(phy$node.label)

...

    if (is.null(phy$root.edge)) {
        cp(")")
        if (nodelab) 
            cp(phy$node.label[1])
        cp(";")
    }
    else {
        cp(")")
        if (nodelab) 
            cp(phy$node.label[1])
        cp(":")
        cp(sprintf(f.d, phy$root.edge))
        cp(";")
    }

...

Настоящая трудность состоит в том, чтобы найти правильный порядок узлов. Я искал и искал, но не смог найти способ найти правильный порядок апостериори.... так что это означает, что мы должны будем получить эту информацию во время преобразования из объект класса hclust в объект класса phylo.

И, к счастью, если вы посмотрите на функцию as.phylo.hclust, там есть вектор, содержащий индексы узлов в их правильном порядке по отношению к предыдущему объекту hclust:

> as.phylo.hclust
function (x, ...) 
{
    N <- dim(x$merge)[1]
    edge <- matrix(0L, 2 * N, 2)
    edge.length <- numeric(2 * N)
    node <- integer(N)              #<-This one
...

Это означает, что мы можем сделать наш собственный as.phylo.hclust с параметром nodenames, если он находится в том же порядке, что и узлы в объекте hclust (что имеет место в вашем примере, поскольку pvclust поддерживает согласованный внутренний порядок, т.е. порядок узлов в hclust такое же, как и в таблице, в которой вы ковыряли бутстрапы):

# NB: in the following function definition I only modified the commented lines
as.phylo.hclust.with.nodenames <- function (x, nodenames, ...) #We add a nodenames argument
{
    N <- dim(x$merge)[1]
    edge <- matrix(0L, 2 * N, 2)
    edge.length <- numeric(2 * N)
    node <- integer(N)
    node[N] <- N + 2L
    cur.nod <- N + 3L
    j <- 1L
    for (i in N:1) {
        edge[j:(j + 1), 1] <- node[i]
        for (l in 1:2) {
            k <- j + l - 1L
            y <- x$merge[i, l]
            if (y > 0) {
                edge[k, 2] <- node[y] <- cur.nod
                cur.nod <- cur.nod + 1L
                edge.length[k] <- x$height[i] - x$height[y]
            }
            else {
                edge[k, 2] <- -y
                edge.length[k] <- x$height[i]
            }
        }
        j <- j + 2L
    }
    if (is.null(x$labels)) 
        x$labels <- as.character(1:(N + 1))
    node.lab <- nodenames[order(node)] #Here we define our node labels
    obj <- list(edge = edge, edge.length = edge.length/2, tip.label = x$labels, 
        Nnode = N, node.label = node.lab) #And you put them in the final object
    class(obj) <- "phylo"
    reorder(obj)
}

В конце концов, вот как вы могли бы использовать эту новую функцию в своем тематическом исследовании:

bootstraps <- (round(y$edges,2)*100)[,1:2]
yy<-as.phylo.hclust.with.nodenames(y$hclust, nodenames=bootstraps[,2])
write.tree(yy,tree.names=TRUE,digits=2)
[1] "((x5:0.27,x8:0.27)100:0.24,((x7:0.25,(x4:0.14,(x1:0.13,x3:0.13)61:0.014)99:0.11)100:0.23,(x2:0.46,x6:0.46)56:0.022)61:0.027)100;"
#See the bootstraps    ^^^ here for instance
plot(yy,show.node.label=TRUE) #To show that the order is correct
plot(y) #To compare with (here I used the yellow value)