Treemix预测基因流原理和方法

什么是基因流

在群体遗传学上，基因流（也称基因迁移）是指从一个物种的一个种群向另一个种群引入新的遗传物质，从而改变群体“基因库”的组成。通过基因交流向群体中引入新的等位基因，是遗传变异一个非常重要的来源，影响群体遗传多样性，产生新的性状组合。基因流会减少种群之间的差异。
下图就形象地展示了基因交流的过程。在某个时期某处发生地质事件，形成一座很高的山峰将一群鸟类群体分隔开来，且鸟类无法自由飞越，一段时间后，由于鸟群生活环境发生巨变，造成山脉分割的两个群体羽毛颜色发生改变（一个群体是显性蓝色HH，一个群体是阴性红色hh），其它羽毛颜色的鸟类不适应环境逐渐灭绝。假设在这两个群体鸟类还未产生生殖隔离前，由于某种原因，如人类活动、鸟类迁徙等，红蓝两个群体的鸟都引入了对方颜色的个体，从基因层面看就是两个亚群都引入的新的等位基因，所以两个群体的基因型频率便会随之改变 ，宏观层面看就是红色鸟类中出现了蓝色的鸟，同样蓝色鸟中也出现了红色的鸟。

自然界中，基因流案例很多，即使在人类社会中也是很常见的：
例如跨国婚姻中不同肤色、眼睛颜色、头发颜色等差异的种族间的结合等。

Treemix软件介绍

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

• 使用基因频率数据可以计算出每对群体之间的协方差，这是实际的协方差（Real value）；
• 使用基因型频率数据可以构建最大似然树，利用两个种群在树上的关系，可以计算出协方差的估计值（Estimated value）；
• 通过实际值与估计值之间的差的大小，判断两个种群之间是否发生基因流，即如果实际值小于估计值很多时，则说明我们构建出来的树夸大了种群之间的差异，则说明种群之间有基因交流，因为基因流会减少种群之间的差异。

1. 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")
}

Treemix结果展示

假设种群间没有基因流：

假设种群间有一次基因流：

当热图中所有种群之间值为0时，表明预测的最大似然树模型与实际相吻合，因此可判断种群间有一次基因流。