【群体遗传】Treemix四倍体分析 & 最优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

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

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

（2）输入文件生成

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

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

image.png

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%的数据时（用该算法计算出来的协方差矩阵作为衡量标准），便可以停止添加migration edge。但是在真实的情况下，99.8%这个数值非常难达到，而且当所使用的数据中非纯合群体数目相较于纯合群体数目更多，Treemix更倾向于缩短纯合群体的branch length，而不是推测真实情况下存在的基因流事件（也就是说，在上述情况下，Treemix更加倾向于将事实简单化，99.8%变成了“永远也达不到的真实”）。

关于这部分，我个人的理解是：当分析的数据中，存在比较古老的基因流事件时，且source为较早分化出来的population，Treemix倾向于缩短该population的branch length（从而导致模型的似然值下降），但于此同时，模型在admixed population中添加migration edges，又增加了模型的似然值，从而达到“拆东墙，补西墙”。。。
【branch length，是covariance matrix中的数值】

OptM则从另一个角度对最优m进行判断（实际上，OptM选择的是给定基因流次数范围内的相对最优m）。作者通过构建 $Δm$ （要求同一m值，进行多次bootstrap运算），利用模型之间似然值的变化速率判断最优m值。
$Δm = \frac{L{''}(m)}{σL(m)}$
当 $Δm$ 达到最大值时，该条件下的m为最优m值。

image.png

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

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()函数得到模型对使用数据的解释。

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

参考文献

[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来研究群体之间的基因流

【群体遗传】Treemix四倍体分析 & 最优m值选择