基因流分析——Treemix

2021-11-27 本文已影响0人 EwanH

基本概念：

基因流（也称基因迁移）

是指从一个物种的一个种群向另一个种群引入新的遗传物质，从而改变群体“基因库”的组成。通过基因交流向群体中引入新的等位基因，是遗传变异一个非常重要的来源，影响群体遗传多样性，产生新的性状组合。

Treemix介绍

treemix由Joseph K. Pickrell和Jonathan K. Pritchard开发，文章Inference of population splits and mixtures from genome-wide allele frequency data. 通过从多个种群中获得等位基因频率，返回该种群的最大似然树，并推断可能发生的杂交事件。
Treemix软件使用全基因组的等位基因频率数据，推断多个群体的分化和混合的模式。该软件输入数据为多个群体的等位基因频率数据，可以生成这些群体的最大似然树，并且可以推测群体混合事件。软件的示例数据是53个人类基因组数据，得到如下结果：

FigureA. Maximum likelihood human tree

人种最大似然树。群体的颜色代表其地理位置。下方的比例尺展示了样品协方差矩阵中元素的10倍平均标准差。

FigureB. Residual fit from tree

残差拟合热图。通过图A的最大似然树得到的残差拟合值。将每对群体（群体i和群体j）之间的残差协方差值除以所有样品对之间的平均标准差，使用这个标准化后的残差绘制该图。右侧为颜色标尺。白色（0点）以上的残差表示对应群体之间的关系比最大似然树上的关系更紧密，暗示这些群体之间有基因渗入事件。

FigureC. ML tree of 53 human populations with inferred migration edges
在最大释然树上展示渗入事件（箭头）。

Treemix基本原理

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

软件下载：

依赖软件：

最新版Boost
最新版gsl
下载

wget -c http://ftp.club.cc.cmu.edu/pub/gnu/gsl/gsl-latest.tar.gz

安装

tar -zxf gsl-latest.tar.gz
cd gsl-2.x
./configure
make
make install
#添加环境变量
export PATH=$PATH:/usr/local/bin
export C_INCLUDE_PATH=$C_INCLUDE_PATH:/usr/local/include
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib
export GSL_LD=/usr/local/lib

下载Treemix

下载地址：Treemix
安装

tar -zxvf treemix-1.13.tar.gz 
cd treemix-1.13/
#默认安装需要root权限
./configure 
make
sudo make install

下载需要用到格式转换脚本：

conversion script
plink2treemix.py
把这些脚本都放到环境变量即可。

Treemix的使用

Treemix有官方的说明书，写得非常详细，下面我只列举一些我用到的参数。

输入文件

1、存储每个个体基因型的vcf文件（YourFileName.vcf）；
2、样品分组文件（第一列和第二列都为样本名称，第三列为分类单元名称，population.txt）；
3、分组的排序文件（poporder.txt）。
前两个文件是构建进化树所必须的。第三个文件规定了热图中的分组顺序。

运算

导入vcf文件
File=YourFileName.vcf
#导入clust文件
clust=population.txt
#LD过滤
plink --vcf $File --indep-pairwise 100 20 0.2 -out ${File%%.*}_LD100_20_02 --make-bed
plink --bfile ${File%%.*}_LD100_20_02 --extract ${File%%.*}_LD100_20_02.prune.in --out ${File%%.*}_LD100_20_02_vcf --recode vcf-iid
#使用文件格式转换脚本生成treemix输入文件
gzip ${File%%.*}_LD100_20_02_vcf.vcf
vcf2treemix.sh ${File%%.*}_LD100_20_02_vcf.vcf.gz $clust
#运行treemix
for i in {0..15}
do
        treemix -i ${File%%.*}_LD100_20_02_vcf.treemix.frq.gz -m $i -o ${File%%.*}_LD100_20_02_vcf.$i -root 1 -bootstrap -k 200 > treemix_${i}_log &
done

参数解释：

-root 指定外群（outgroup）
-k 以滑动窗口的方式选择SNP位点构建树
-m 指定预估可能有几次基因流事件。比如根据经验推测可能有两次基因流事件，则-m参数设置为2（{0..15}执行-m等于0-15）
-g 指定先验进化树
-bootstrap 进行bootstrap replicate
-noss 关闭样本数量矫正

数据可视化

#R环境
source("plotting_funcs.R")
poporder="poporder.txt"
#绘制-m等于11的结果
outstem="YourFileName_100_20_02_vcf.11"
#绘制FigureC
plot_tree(outstem)
#绘制FigureB（残差图）
plot_resid(outstem,poporder)