种群遗传学树状混交分析和最佳 m 值确定

之前已经有生信菜鸟团前辈写过关于Treemix的分析，但值得关注的是，现在的软件大部分都只适用于二倍体，那今天介绍一下四倍体怎么进行Treemix分析以及如何选择最优m值。

首先，还是对Treemix做一个简要的介绍：Joseph Pickrell和Jonathan Pritchard于2012年，发表在Nature Precedings上的文章《Inference of population splits and mixtures from genome-wide allele frequency data》中的一个算法 & 软件，利用allele frequency构建最大似然树。

（1）软件安装

# 简单粗暴一点吧
conda install treemix

当然，也可以选择自行下载安装

wget https://bitbucket.org/nygcresearch/treemix/downloads/treemix-1.13.tar.gz
./configure
make
make install

对一些参数做个简要介绍：

-k         # 将长度为n的窗口内的所有SNP，当作一个整体用于分析
-root      # 设定最大似然树的根（一般是外类群）
-g         # 提供vertices.gz和edges.gz文件，作为目前分析的基础，即在之前分析的基础上，继续提升基因流的次数
-bootstrap # 取子集，进行最大似然树的构建

（2）输入文件生成

这一步：如何成为一个及格过线的调包侠。

调用ScanTools工具包中的conversionTreemixMajda.py（这个脚本适用于二倍体 & 四倍体混合群体），就可以得到用于Treemix分析的输入数据。由于ScanTools工具包在分析过程中，会对VCF文件进行格式转换，而且无法直接使用ScanTools，因此就需要手动进行一些操作，大致流程为：

python conversionTreemixMajda.py -i "allele_table/" -o "allele_table/" 
# 最终得到treemix_input.table.gz

（3）Treemix运行

for i in {0..10}
    do
    for j in {1..10}
        do
        treemix -i treemix_input.table.gz -root Outgroup -k 500 -m ${i} -bootstrap -o migration_${i}_bootstrap_${j}
    done
done

多个m值以及需要多次bootstrap的情况下，还是推荐并行运行，使用for循环所需的运行时间太久了。

（4）使用OptM判断最优m值

在进行Treemix分析时，Treemix作者推荐当构建的model能够解释99.8%的数据时（用该算法计算出来的协方差矩阵作为衡量标准），就已经说明模型已经足够符合真实情况。但是在真实的情况下，99.8%这个数值非常难达到，而且当所使用的数据中非纯合群体（admixed population）数目相较于纯合群体（unadmixed population）数目更多，Treemix更倾向于将纯合群体的枝条长度缩短，而不是推测真实情况下存在的基因流事件（也就是说，在上述情况下，Treemix更加倾向于将事实简单化，忽略了unadmixed population向admixed population的基因流，99.8%变成了“永远也达不到的真实”）。

OptM则从另一个角度对最优m进行判断（实际上，OptM选择的是给定基因流次数范围内的相对最优m）。作者通过构建Δm（要求同一m值情况下，进行多次bootstrap运算）来判断最优m值（L''(m)和L(m)具体含义参考OptM文章附件）

当Δm达到最大值时，该条件下的m为最优m值。

在R中，基本使用方法如下：

# install.packages("OptM")
library(OptM)

# 读入数据
dat <- optM("treemix_results")   # treemix_results为Treemix分析结果文件夹，包含每一次运行的llik、modelcov.gz等文件

# 判断哪一个m是最优的
plot_optM(dat, method = "Evanno", plot = TRUE, pdf = "OptM.pdf")

（5）Treemix结果图绘制 & 解读

Treemix结果可视化，使用Treemix的自带R脚本：plotting_funcs.R

library(RColorBrewer)
library(R.utils)
source("plotting_funcs.R")

# 绘制Treemix graph
for (i in 1:10){
  pdf(paste('migration_m_bt_', i, '.pdf', sep = ""), width = 14, height = 7)
  plot_tree(paste('treemix_results/migration_m_bt_', i, sep=""))
  plot_resid(paste('treemix_results/migration_m_bt_', i, sep=""), 'poplist.txt')
  dev.off()
}

此处对残差值进行一个简要说明：

•当残差为正值时（observed value＞estimated value），表示对应population之间的关系被模型所低估（两个population之间有introgression可能符合真实情况）•当残差为负值（observed value＜estimated value），表示Treemix构建的最大似然树高估了对应population之间的关系（两个population没有introgression更符合真实情况）

同时，也可以使用get_f()函数得到模型对使用数据的解释。

对于残差热图绘制的解释，已经有前辈做了非常好的介绍，我这边就不再赘述。

群体遗传系列之：一文了解和使用Treemix来研究群体之间的基因流

参考文献

[1] Inference of population splits and mixtures from genome-wide allele frequency data

[2] OptM : estimating the optimal number of migration edges on population trees using Treemix

[3] Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study

[4] 群体遗传系列之：一文了解和使用Treemix来研究群体之间的基因流