Comparison of Bayesian and maximum-likelihood inference of population genetic parameters
Migrate-n 方法原始文献，2006 年发表在 Bioinformatics 上。作者为 Peter Beerli。
使用 coalescence theory，结合了最大似然法（ML）和贝叶斯推断，推断群体大小、migration rates。
Parameters are scaled effective population sizes Θi (4× effective population size × mutation rate per site per generation), and scaled immigration rates Mji (immigration rate divided by mutation rate).
θ = 4Neμ，M = m/μ。m 是每一代迁徙的个体的比例，μ 是每个位点的突变率，Ne 是有效种群大小。
记一下我总结的 migrate-n 的用法：
migrate-n 是为桑格测序数据设计的一款进行种群大小、迁移率推断的软件，可以首先设定几个可能的模型（例如两个种群，可以设置 full migration model 即两个种群间有双向基因流，也可以设置单向基因流的模型），运行 migrate-n 得到每个模型的 marginal likelihood。
对于二代测序数据，migrate-n 有将 vcf 转换为 migration 格式的脚本，见 tutorial。但是这个转换格式的脚本不支持过多的染色体数目，例如我自己的数据是 scaffold level 的参考基因组，就无法完整地进行转换。而如果是染色体水平的参考基因组，虽然这一步可以成功转换，但到了运行的步骤也很可能失败；例如我选择了最大的 20 条 scaffolds，无法成功运行，日志文件里并不明说理由，程序中止，很可能是因为数据量太大。在我自己的例子里，减小数据量可以解决问题。可以随机地选取一些不太大的 scaffolds。
设置模型：
我自己的二群体模型，简单的三种情况。
软件的 tutorial 里详细地介绍了模型的设置方法，可以非常复杂。除了migration model 以外还有 divergence model。我这里是 migration model。
运行 migrate-n，在 menu 里设置好输入文件、parameter 等，运行。
我遇到了一个很奇怪的 bug，有几个个体的序列就是读不进去。后来改了个体名称就好了。
获得的 outfile。可以将不同模型的 outfile 分别命名为 model1.out 等。作者在 tutorial 里写到用以下命令获取最佳模型：
grep "All" model*.out | sort -n -k 4,4 | python bf.py
注意这个 python 脚本不在软件包里，需要用 tutorial 里的链接下载。此外在我自己的情况里，上述脚本不能直接用于计算，因为包含 All 的行不仅有模型总结参数，还有 θ 和 M 的估计，所以我自己又在 sort 之后加上了 head -n N，这个 N 是比较的模型的数目。获得的结果如下表：
软件 tutorial 里给的示例
LBF 是贝叶斯因子的 log 值。等于 0 时意味着贝叶斯因子为 1，在其他地方读到说这意味着完全否定原假设 。这里我想需要进一步了解 migrate-n 的内容才能读明白了。总之这里给出了最佳模型。可以拿最佳模型的 outfile，选取出估计得到的 θ 和 M，进一步得到有效种群大小和迁移率。
可以读一读其他用到这一软件的文章。
Genetic diversity and structure of the Chinese lake gudgeon ( Sarcocheilichthys sinensis)
一种鱼的线粒体 DNA 分析。用了许多模拟的方法，包括 ABC、S-DIVA、migrate-n、TMRCA，研究种群的地理历史。证明了这个物种在长江、岷江间的扩散历史。
写到用 migrate-n 估计了 mutation-scaled population size θ = 2Neμ，以及 mutation-scaled immigration rate M = m/μ。（对于线粒体的数据）
mitochondrial DNA θ = 2Neμ (where Ne is the effective population size of females, and μ is the mutation rate).
网友评论