在Rstudio中使用ape包的ace函数时出错：object "phy" is not of class "phylo"

Question

我正在尝试运行 ASR 分析，并尝试将二进制数据组合到先前生成的树的叶子中。从 ape 包运行 ace 函数时，我不断收到相同的错误：对象“phy”不属于“phylo”类。以下是我正在运行的完整代码。

install.packages("ape")
install.packages("phytools")
install.packages("tidyverse")
library(tidyverse)
library(ape)
library(phytools)

# Read the tree
tree_path <- "/Users/USERNAME/Downloads/aeruginosa_galaxy_tree"
tree <- read.tree(tree_path)

# Read the CSV file
csv_path <- "/Users/USERNAME/Desktop/Pf4_gene_results_modified.csv"
gene_data <- read.csv(csv_path, row.names = 1)  # Assuming the first column contains row names

# Use only common identifiers
common_identifiers <- intersect(rownames(gene_data), tree$tip.label)

# Filter the gene data and tree
gene_data_filtered <- gene_data[common_identifiers, , drop = FALSE]
tree_filtered <- drop.tip(tree, setdiff(tree$tip.label, common_identifiers))


# List identifiers in gene_data that are not in the tree's tip labels
missing_in_tree <- setdiff(rownames(gene_data_filtered), tree$tip.label)

# List identifiers in the tree's tip labels that are not in gene_data
missing_in_gene_data <- setdiff(tree$tip.label, rownames(gene_data_filtered))

# Output the number of missing identifiers
length(missing_in_tree)
length(missing_in_gene_data)

class(tree)

# Example: Using a Parsimony model for binary data
fit <- ace(tree, gene_data_filtered, type="discrete", model="parsimony")

# Plot the tree with ancestral states
plot(fit, tree)

我认为这是 ape 包及其与我的 Rstudio 环境交互的问题。我已经重启过环境，在新的环境下尝试过，并重新安装了多次ape包。当我使用 ape 包中的内置数据集时，也会出现此问题。

> data(bird.orders)
> fit_example <- ace(bird.orders$tree, bird.orders$trait, model="ER")
    Error in ace(bird.orders$tree, bird.orders$trait, model = "ER") : 
      object "phy" is not of class "phylo"

Answer 1

ace()

不断抛出错误，因为它期望将

phylo

类的对象传递给

phy

参数。如果您查看

ace()

的帮助文件，您会发现前两个参数（

x

和

phy

）应该分别是字符状态向量和树。

例如：

set.seed(1)
x <- round(rnorm(23, 30, 5))
ans <- ace(x, bird.orders)

par(mar=c(0, 0, 0, 0))
plot(bird.orders, type="tidy", FALSE, label.offset=2)
tiplabels(text=x, cex=0.7, adj=-0.1)
nodelabels(text=round(ans$ace, 1), cex=0.7)