斯坦的系统发育推断
phylostan的Python项目详细描述
利用stan进行系统发育推断
简介
phylostan是一个用python编写的工具，用于从核苷酸数据集推断系统发生树。 它使用stan语言生成各种系统发育模型。 通过pystan库，phylostan可以访问stan的变分推理和采样（nuts和hmc）引擎。
功能
系统发育模型组成：
- 核苷酸替代模型：JC69，hky，gtr
- 利率异质性：离散威布尔分布和一般离散分布
- 拓扑上具有一致先验的无时钟约束树
- 时间树：
- 同步序列：相同的采样日期
- 异时序列：在不同时间点采样的序列
- 分子钟：
- 严格
- Autocorrelated
- Uncorrelated：记录正常的分层优先级
- 合并模型：
斯坦提供的算法：
- 变分推理：
- 平均场分布
- 全秩分布
- 无U形回转取样器（NUTS）
- 哈密顿蒙特卡罗（hmc）
先决条件
| Program/Library | Version | Description |
|---|---|---|
| python | Tested on python 2.7, 3.5, 3.6, 3.7 | |
| pystan | >=2.19 | API for Stan |
| dendropy | Library for manipulating trees and alignments | |
| numpy | >=1.7 |
您可以使用PIP安装Phylostan
pip install phylostan
您也可以在本地运行它
python -m phylostan.phylostan <COMMAND>
其中<COMMAND>是build或run命令。
命令行用法
phylostan被分解为两个子命令：
- build：创建一个stan文件：包含模型的文本文件。
- run：使用数据运行stan文件。
这两个步骤是分开的，因此用户可以编辑stan模型。主要原因是要修改prior。
要获得有关build或run命令的帮助：
phylostan build --help phylostan run --help
快速启动
我们将使用fluA.fa对齐和fluA.tree树文件。这个数据集包含了从1981年到1998年间分离的69个甲型流感病毒血凝素核苷酸序列。
首先，需要使用build命令生成stan脚本：
cd examples/fluA phylostan build -s fluA-GTR-W4.stan -m HKY -C 4\ --heterochronous --estimate_rate --clock strict --coalescent constant
此命令将使用以下模型创建stan文件fluA-GTR-W4.stan：
- 长谷川、kishino和yano（hky）核苷酸替代模型
- 使用weibull分布的4种利率类别的利率异质性
- 假设采集的序列是不同的时间点（非同步）
- 恒定有效种群大小
- 估计替代率
在第二步中，我们使用数据编译并运行脚本
phylostan run -s fluA-GTR-W4.stan -m HKY -C 4\ --heterochronous --estimate_rate --clock strict --coalescent constant \ -i fluA.fa -t fluA.tree -o fluA -q meanfield
run命令需要数据（树和对齐）和输出参数。 它还需要为build命令提供的参数。 输出将包含4个文件：
fluA：这个文件是stan的输出文件。它包含从变分分布中提取的样本（或mcmc样本）。fluA.diag：这个文件也是由stan生成的，它包含一些信息，比如每次迭代时的elbo。fluA.trees：这个文件是一个包含树的nexus文件。它可以用FigTree等程序打开，也可以使用BEAST或BEAST2中的treeannotator进行摘要。fluA-GTR-W4.pkl：stan脚本被编译成这个二进制文件。此文件可以由phylostan自动重用，除非必须重新编译，否则可以使用--compile选项。
在运行结束时，phylostan将在屏幕上打印感兴趣参数的平均可信区间和95%：
Weibull (shape) mean: 0.488 95% CI: (0.383,0.616)
Strict clock (rate) mean: 0.00499 95% CI: (0.00432,0.00577)
Constant population size (theta) mean: 4.03 95% CI: (3.14,5.05)
HKY (kappa) mean: 5.58 95% CI: (4.37 7.039)
Root height mean: 18.96 95% CI: (18.36 19.74)
在这个例子中，我们使用平均场分布（^ {CD10>}）来使用变量来逼近后部。推理。 stan模型已经编译好了，因此我们可以运行nuts算法，而无需重新生成脚本文件，只需发出以下命令：
phylostan run -s fluA-GTR-W4.stan -m HKY -C 4\ --heterochronous --estimate_rate --clock strict --coalescent constant \ -i fluA.fa -t fluA.tree -o fluA -a nuts
nuts算法比变分推理慢得多（而且更精确），因此应该在小数据集上使用。
参考
马修·福门特和亚伦·E·达林。系统发育学中概率规划和快速变分贝叶斯推理的评价。{EM1}$BioXiV
欢迎加入QQ群-->： 979659372
