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

中秋假期过去了，学习不能停下来。这一期推文继续和大家分享与群体遗传进化相关的知识。这一期主要讲解基因流和Treemix的使用，读完后希望对大家有帮助且有所收获。

什么是基因流(Gene flow)?

基因流（也称基因迁移）是指一个物体中的一些个体从一个群体迁移到另一个群体，这个过程中某些基因或遗传物质会引入到新的群体中，从而产生基因的流动。这个过程会改变群体“基因库”的组成，改变基因的频率。下面的图就是一个非常经典的例子，通过基因流改变了两个鸟群体的基因频率。

通过基因交流向群体中引入新的等位基因，是遗传变异中非常重要的来源，影响群体遗传多样性，可以导致新的性状组合的产生。基因在群体间中流动的水平越大，群体的基因组成相对应的性状就会越均匀或普遍相似，受限制的基因流使群体间发生分化，因为每个群体中都会或多或少的独立发生适应和遗传漂变。群体间不发生基因流可能是因为生殖隔离而没有相互杂交，或因为地理隔离而无法杂交。

Treemix 软件介绍

TreeMix 是由Joseph K. Pickrell和Jonathan K. Pritchard开发，一种推断一组种群历史中种群分化和基因流的工具。在基础模型中，一个物种的现代种群通过祖先种群与共同祖先相关。TreeMix通过从多个种群中获得等位基因频率，返回该种群的最大似然(ML)树，并推断可能发生的杂交事件。

其基本原理可以分为三个要点：

1. 根据基因频率，算出每对群体之间的协方差
2. 根据基因型频率数据，构建最大似然树，利用两个种群在进化树上的关系，计算出协方差的估计值
3. 根据实际值与估计值之间的差的大小，判断两个种群之间是否发生基因流，如果实际值小于估计值，则说明我们构建出来的树夸大了种群之间的差异，则说明种群之间有基因交流，因为基因流会减少种群之间的差异

Treemix安装比较简单，直接下载编辑就能安装：

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

Treemix 软件使用

Treemix的主要输入文件有三个分别是：SNP的VCF文件，群体的分群信息和分组顺序文件。这里我使用我的大豆数据来示范，一共370个个体包括野生种，地方种和栽培种。

群体的分群信息,例子如下：

head -n6 new_pop.cluster

HN001    HN001   Cultivar
HN006    HN006   Landrace
HN007    HN007   Landrace
HN009    HN009   Landrace
HN012    HN012   Landrace
HN_HP025    HN_HP025    Wild-type

分组顺序文件如下：

cat order.list

Cultivar
Landrace
Wild-type

然后由于，Treemix是假设你的SNP是不连锁的，并且它并不喜欢SNP VCF中有缺失的数据，这样我们先对VCF文件进行一些过滤：

###过滤缺失数据
vcftools --vcf SNPs.vcf --max-missing 1 --recode --stdout > SNPs_no_missing.vcf
### 根据LD来过滤连锁的SNP
plink --vcf flooding_no_missing.vcf --indep-pairwise 50 10 0.2 --out tmp.ld --allow-extra-chr --set-missing-var-ids @:# --keep-allele-order
~/biosoft/plink --vcf flooding_no_missing.vcf --extract tmp.ld.prune.in --freq --missing --within new_pop.cluster --out input --allow-extra-chr --set-missing-var-ids @:# --keep-allele-order

压缩频率文件，使用treemix里面的脚本，将freq频率文件转成treemix的输入文件：

gzip input.frq.strat
python2.7 plink2treemix.py input.frq.strat.gz input_treemix.frq.gz

进行treemix基因流分析,假设群体间有0-2次基因流的发生，分别运行：

for i in {0..2}
do
    treemix -i input.frq.treemix.gz -m $i  -root Wild-type  -noss -o out.${i} > treemix_${i}_log &
done

对结果进行可视化：

prefix="out"
library(RColorBrewer)
library(R.utils)
#加载treemix提供的R脚本,https://github.com/lakeseafly/bioconductor_learn/blob/master/plotting_funcs.R
source("plotting_funcs.R")
#绘制残差图
for(edge in 0:2){
 plot_resid(stem=paste0(prefix,".",edge),pop_order="dogs.list")
}

结果简单解读

首先是没有基因流的情况：

由热图可以看出，Wild-type与Cultivar群体之间的Landrace颜色对应的数值在0以上，所以推测有1次基因流动。

现在让我们看看假设有一次基因流的情况：

树上添加了一条箭头，从Wild-type指向Landrace,说明从野生种到地方种存在基因流。

对应的热图标尺中，0SE为最上端的颜色，且热图都没有明显的颜色分化，说明模型对实际群体之间的协方差拟合程度很好。

介绍到这里就结束了，希望大家通过本次推文有所了解什么是基因流，然后怎样通过Treemix来检测基因流，并且运用到自己的数据中来解决对应的生物学问题。

参考资料：

使用Treemix软件进行基因流分析简介（https://www.oebiotech.com/news/news182.html）

基因流与D测验课程纪要 (https://www.bilibili.com/read/cv6041247)