如何将我的数据和问题放入 Rstudio 中 ape pakage 的 ace 函数的脚本中？

Question

我有96个氨基酸序列，我用MAFFT比对并手动修剪（FASTA格式），用ProtTest选择氨基酸替换模型（LG+I+G模型），用MEGAX做了系统发育重建（ML方法， bootstrap 测试 1000 次重复，Newick 格式的树）和使用 PAML 的祖先重建，总共 664 个最终氨基酸位置。但是，我的对齐有插入缺失。我用字母 (A 到 T) 命名每个插入缺失，并且各自的酰胺酸位置范围：A:89-92, B:66-67, C:181-186, D:208-208, E:214-219 , F:244-250, G:237-296, H:278-280, I:295-295, J:329-334, K:345-349, L:371-375, M:390-425, N ：432-433，O：440-443，P：480-480，Q：500-500，R：541-544，S：600-600。序列的首尾部分变化很大，因此从位置 0 到 34（首字母）和 600 到 664（最终），每个氨基酸位置可能代表一个插入缺失。

我想知道，在每个祖先节点上，每个 indel 出现在祖先序列中的概率是多少。有人告诉我，包“ape - 系统发育和进化分析”上的 R-studio“ace”功能可以执行此任务。我已经安装了“ape”和“ggtree”。我检查了这个网页https://www.rdocumentation.org/packages/ape/versions/3.0-1/topics/ace，但是，我不知道如何构建脚本。我是一名生物学家，也是 R. 的新手。

有人可以帮忙吗？将不胜感激，谢谢。

Answer 1

很难从您的示例中准确找出您需要什么，但以下内容可能符合总体思路：

1 - 将你的树加载到R

对于这一步，您可以使用函数read.tree或read.nexus，具体取决于您的树格式：即您的系统发育软件是否输出NEXUS文件（通常这些文件中的第一行是#NEXUS，最后一行是end; 或END;）或newick 输出（通常，第一行直接以((my_species... 之类的系统发育开始并以; 结束）。您可以找到该文件，然后在 R 中使用：

## Loading the package
library(ape)
## Reading the tree
my_tree <- read.tree("<the_path_to_your_file>")

2 - 在R 中加载您的特征数据

然后您需要将特征数据（例如您上面列出的插入缺失位置）加载为matrix 或data.frame。最简单的方法是将它们设为.csv 格式（“逗号分隔值”），然后您可以使用read.csv 函数读取R：

## Reading the variables as a matrix
my_variables <- read.csv("<the_path_to_your_file>")

3 - 运行祖先字符估计

最后，您可以使用 ape 包中的 ace 函数为每个变量运行祖先字符估计：

## Selecting the variable of interest (e.g. the first column of the dataset)
one_variable <- my_variables[, 1]
## Running the ancestral character estimation for this variable
my_ace <- ace(x = one_variable, phy = my_tree, type = "discrete")
## Looking at the results
my_ace

当然还有更多内容，但希望这可以让你开始。