1

是否有可能0/1从像左边那样的分叉系统发育树生成像下面右图所示的字符矩阵。矩阵中的1表示存在联合进化枝的共享特征。

这段代码生成了很好的随机树，但我不知道从哪里开始将结果转换为字符矩阵。

library(ape) # Other package solutions are acceptable

forest <- rmtree(N = 2, n = 10, br = NULL)
plot(forest)

为了清楚起见，我可以使用以下代码生成随机矩阵，然后绘制树。

library(ape)
library(phangorn)

ntaxa <- 10
nchar <- ntaxa - 1

char_mat <- array(0, dim = c(ntaxa, ntaxa - 1))

for (i in 1:nchar) {
  char_mat[,i] <- replace(char_mat[,i], seq(1, (ntaxa+1)-i), 1)
}

char_mat <- char_mat[sample.int(nrow(char_mat)), # Shuffle rows 
                     sample.int(ncol(char_mat))] # and cols

# Ensure all branch lengths > 0
dist_mat <- dist.gene(char_mat) + 0.5
upgma_tree <- upgma(dist_mat)
plot.phylo(upgma_tree, "phylo")

我想要的是生成随机树，然后从这些树中制作矩阵。此解决方案无法生成正确类型的矩阵。

为清楚起见进行编辑：我正在生成二进制字符矩阵，学生可以使用它来使用简单的简约法绘制系统发育树。该1字符代表将分类群联合成进化枝的同源性。因此，所有行必须共享一个字符（一1列中的所有行），并且某些字符必须仅由两个分类单元共享。（我在打折自体变形。）

例子：

在此处输入图像描述

4

2 回答 2

2

您可以看一下非常简单的rTraitDisc功能：ape

library(ape)
## You'll need to simulate branch length!
forest <- rmtree(N = 2, n = 10)

## Generate on equal rate model character
(one_character <- rTraitDisc(forest[[1]], type = "ER", states = c(0,1)))
# t10  t7  t5  t9  t1  t4  t2  t8  t3  t6 
#   0   0   0   1   0   0   0   0   0   0 
# Levels: 0 1

## Generate a matrix of ten characters
(replicate(10, rTraitDisc(forest[[1]], type = "ER", states = c(0,1))))

#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# t10 "0"  "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  
# t7  "0"  "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  
# t5  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t9  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t1  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t4  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t2  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t8  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t3  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t6  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"

要将其应用于多个树，最好的方法是创建一个 lapply 函数，如下所示：

## Lapply wrapper function
generate.characters <- function(tree) {
    return(replicate(10, rTraitDisc(tree, type = "ER", states = c(0,1))))
}

## Generate 10 character matrices for each tree
lapply(forest, generate.characters)

# [[1]]
#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# t10 "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  
# t7  "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  
# t5  "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  
# t9  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t1  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t4  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t2  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t8  "0"  "0"  "0"  "1"  "0"  "1"  "0"  "0"  "0"  "1"  
# t3  "0"  "0"  "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  
# t6  "0"  "0"  "0"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  

# [[2]]
#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# t7  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t9  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t5  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t2  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t4  "0"  "1"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  
# t6  "0"  "1"  "0"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  
# t10 "0"  "1"  "1"  "0"  "1"  "1"  "0"  "0"  "0"  "1"  
# t8  "0"  "1"  "1"  "0"  "1"  "0"  "0"  "0"  "0"  "0"  
# t3  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  
# t1  "0"  "1"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"

另一种选择是使用sim.morpho包中的dispRity。该函数重用了该rTraitDisc函数，但实现了更多模型，并允许将速率作为分布提供，从中进行采样。它还允许角色在没有太多不变数据的情况下看起来更“逼真”，并确保生成的角色“看起来”像一个真实的形态学角色（比如具有适量的同质性等......）。

library(dispRity)
## You're first tree
tree <- forest[[1]]
## Setting up the parameters
my_rates = c(rgamma, rate = 10, shape = 5)
my_substitutions = c(runif, 2, 2)

## HKY binary (15*50)
matrixHKY <- sim.morpho(tree, characters = 50, model = "HKY",
     rates = my_rates, substitution = my_substitutions)

## Mk matrix (15*50) (for Mkv models)
matrixMk <- sim.morpho(tree, characters = 50, model = "ER", rates = my_rates) 

## Mk invariant matrix (15*50) (for Mk models)
matrixMk <- sim.morpho(tree, characters = 50, model = "ER", rates = my_rates,
     invariant = FALSE)

## MIXED model invariant matrix (15*50)
matrixMixed <- sim.morpho(tree, characters = 50, model = "MIXED",
     rates = my_rates, substitution = my_substitutions,  invariant = FALSE,
     verbose = TRUE)

我建议您阅读该功能以获取有关模型如何工作的正确参考资料，或阅读 dispRity 包手册sim.morpho中的相关部分。

于 2019-01-28T03:19:00.663 回答
0

我想出了如何使用Descendantsphangorn包制作矩阵。我仍然需要使用合适的节点标签对其进行调整，以匹配原始问题中的示例矩阵，但框架就在那里。

library(ape)
library(phangorn)

ntaxa <- 8
nchar <- ntaxa - 1

tree <- rtree(ntaxa, br = NULL)

# Gets descendants, but removes the first ntaxa elements,
# which are the individual tips
desc <- phangorn::Descendants(tree)[-seq(1, ntaxa)]

char_mat <- array(0, dim = c(ntaxa, nchar))

for (i in 1:nchar) {
  char_mat[,i] <- replace(char_mat[,i], y <- desc[[i]], 1)
}

rownames(char_mat) <- tree$tip.label
char_mat
#>    [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> t6    1    1    0    0    0    0    0
#> t3    1    1    1    0    0    0    0
#> t7    1    1    1    1    0    0    0
#> t2    1    1    1    1    1    0    0
#> t5    1    1    1    1    1    0    0
#> t1    1    0    0    0    0    1    1
#> t8    1    0    0    0    0    1    1
#> t4    1    0    0    0    0    1    0

plot(tree)

由reprex 包（v0.2.1）于 2019 年 1 月 28 日创建

于 2019-01-28T15:52:30.580 回答