问题描述
我有 96 个氨基酸序列，我用 MAFFT 比对并手动修剪（FASTA 格式），用 Prottest 选择氨基酸替换模型（LG+I+G 模型），用 MEGAX（ML 方法， bootstrap 测试 1000 次重复，Newick 格式的树）和使用 PAML 的祖先重建，总共 664 个最终氨基酸位置。但是，我的对齐方式有插入缺失。我用字母（A 到 T）和相应的酰胺酸位置范围命名每个 indel：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出现在祖先序列中的概率是多少。有人告诉我，“猿 - 系统发育和进化分析”包上的 R-studio“ace”功能可以执行此任务。我已经安装了“ape”和“ggtree”。我检查了这个网页 https://www.rdocumentation.org/packages/ape/versions/3.0-1/topics/ace，但是，我不知道如何构建脚本。我是一名生物学家，也是 R 的新手。
有人可以帮忙吗？将不胜感激，谢谢。
解决方法
很难从示例中确切地确定您需要什么，但以下内容可能符合总体思路：
1 - 在 R
中加载您的树
对于这一步，您可以使用 read.tree 或 read.nexus 函数，具体取决于您的树格式：即您的系统发育软件是否输出 NEXUS 文件（通常这些文件中的第一行是 {{1} } 并且最后一行是 #NEXUS 或 end;) 或一个新的输出（通常，第一行直接以 END; 这样的系统发育开始并以 ((my_species... 结束）。你可以找到这个文件，然后在 R 中使用：
;2 - 在 ## Loading the package
library(ape)
## Reading the tree
my_tree <- read.tree("<the_path_to_your_file>")
中加载特征数据
然后您需要以 R 或 matrix 的形式加载特征数据（例如您在上面列出的 indels 位置）。最简单的方法是将它们设为 data.frame 格式（“逗号分隔值”），然后您可以使用 .csv 函数在 R 中读取：
read.csv3 - 运行祖先性格估计
最后，您可以使用 ## Reading the variables as a matrix
my_variables <- read.csv("<the_path_to_your_file>")
包中的 ace 函数对每个变量进行祖先特征估计：
ape当然还有更多内容，但希望这可以帮助您入门。