如何将位点间成对系统发育距离输出为dist类矩阵？

数据与需求 #

先加载系统发育树和位点类群有无矩阵：

set.seed(020123)
tree_ <- phytools::pbtree(b = 1, n = 7)

matrix_ <- structure(
  c(1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0,
    0, 0, 0, 1, 0, 0, 1, 0, 0, 1),
  dim = c(3L, 7L),
  dimnames = list(c("site_1", "site_2", "site_3"),
  c("t1", "t2", "t3", "t4", "t5", "t6", "t7"))
)

需求是计算所有位点间的成对系统发育距离，并将结果存储为规范的dist类对象。

已实现的距离计算函数 #

已完成两个位点间系统发育距离的计算函数：

get_community_taxa <- function(x) {
  list(
    colnames(x)[x[1, ] == 1L],
    colnames(x)[x[2, ] == 1L]
  )
}

phylo_dist <- function(x, tree) {
  distance_mtx <- ape::cophenetic.phylo(tree)
  communities <- get_community_taxa(x)
  phylo_dist <- distance_mtx[
    communities[[1]],
    communities[[2]]
  ]
  if (length(phylo_dist) == 1L) {
    return(phylo_dist)
  }
  mean(c(
    apply(phylo_dist, 1, function(x) min(x, na.rm = TRUE)),
    apply(phylo_dist, 2, function(x) min(x, na.rm = TRUE))
  ))
}

已完成的尝试 #

通过combn结合apply，已得到所有位点对的距离向量：

dist_vec <- apply(
  combn(nrow(matrix_), 2),
  MARGIN = 2,
  function(x) phylo_dist(matrix_[x, ], tree_)
)
#> [1] 1.276455 2.191359 3.045831

转换为dist类对象 #

有两种简便方式将向量转换为符合要求的dist对象：

方法一：直接构造dist对象 #

利用structure给向量添加dist类属性：

site_names <- rownames(matrix_)
site_dist <- structure(
  dist_vec,
  class = "dist",
  Size = length(site_names),
  Labels = site_names,
  Diag = FALSE,
  Upper = FALSE
)

# 查看结果
site_dist
#>         site_1   site_2
#> site_2 1.276455         
#> site_3 2.191359 3.045831

方法二：先构建对称矩阵再转换 #

先创建空对称矩阵，填充下三角后转为dist对象：

n_sites <- length(site_names)
dist_mtx <- matrix(0, nrow = n_sites, ncol = n_sites, dimnames = list(site_names, site_names))
dist_mtx[lower.tri(dist_mtx)] <- dist_vec

site_dist <- as.dist(dist_mtx, upper = FALSE, diag = FALSE)

两种方法均可得到目标格式的dist类对象。

内容的提问来源于stack exchange，提问作者Junitar