使用 treemix 进行基因渗入分析 & 使用 hierfstat 计算成对 Fst

发表于 | 更新于 | 分类于

treemix 是利用等位基因频率来推断群体间分化和杂合（基因流动或者基因渗入）的软件。

格式转换

通常，我们的变异数据储存格式为vcf,首先将其转换为tped格式。

1
vcftools --vcf test.vcf --plink-tped --out test

转换后的tfam文件如下：第一列为FID,第二列为IID，我们需要修改第一列数据为群体信息。

1
2
3
4
5
wHAIPI052555-68 wHAIPI052555-68 0 0 0 -9
wHAIPI052556-77 wHAIPI052556-77 0 0 0 -9
wHAIPI052557-76 wHAIPI052557-76 0 0 0 -9
wHAIPI052558-63 wHAIPI052558-63 0 0 0 -9
wHAIPI052559-71 wHAIPI052559-71 0 0 0 -9

修改好的tfam文件如下：XM和GD是2个群体

1
2
3
4
5
XM wHAIPI052555-68 0 0 0 -9
XM wHAIPI052556-77 0 0 0 -9
XM wHAIPI052557-76 0 0 0 -9
GD wHAIPI052558-63 0 0 0 -9
GD wHAIPI052559-71 0 0 0 -9

然后再编辑一个群体信息文件(pop.cov)，格式如下：

1
2
3
4
5
XM wHAIPI052555-68	XM
XM wHAIPI052556-77	XM
XM wHAIPI052557-76	XM
GD wHAIPI052558-63	GD
GD wHAIPI052559-71	GD

编辑好之后，计算等位基因组频率：

1
plink --tfile test --freq  --within pop.cov

然后使用treemix 自带的脚本 plink2treemix.py进行格式转换，转换前需要对等位基因频率文件压缩。

1
python2 plink2treemix.py plink.frq.strat.gz treemix_in.gz

基因渗入分析

1
treemix -se -bootstrap -k 1000 -m 1 -i treemix_in.gz

或者进行多个m和重复：

1
2
3
4
5
6
7
8
9
10
11
12
for m in {1..5}
   do
   for i in {1..20}
      do
      treemix -se -bootstrap \
         -i treemix_in.gz \
         -o TreeMix.${i}.${m} \
         -global \
         -m ${m} \
         -k 1000
      done 
done

其中m参数为 the number of migration edges，中文翻译过来其实就是基因渗入方向的个数，默认为0，可以尝试1-10，并在每个m中重复10次。

最佳迁移边缘个数选择

将生成的llik、cov、modelcov 文件放置同一文件夹，使用 R 包OptM进行分析：

1
2
3
library(OptM)
linear = optM("./result")
plot_optM(linear)

生成图中，当Δm值最小时的 migration edges 为最佳迁移边缘个数。

基因渗入作图

确定最佳迁移边缘个数之后，我们使用 m=最佳迁移边缘个数的结果文件作图。

1
2
source("plotting_funcs.R") #treemix scr文件夹中R脚本
plot_tree("TreeMix") #TreeMix为结果文件前缀

pairwise Fst

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(vcfR)
library(adegenet)
library(hierfstat)
kiwipang<-read.vcfR("my.vcf") #读取vcf
df<-vcfR2genind(kiwipang) #vcf2genind
ploidy(df)<-2 #倍型为2
pop.Kiwipang<-read.table("pop.txt",sep="\t",header=F) #读取群体结构文件，第一列为样本名，第二列为群体名，顺序与vcf文件一致
all(colnames(kiwipang@gt)[-1] == pop.Kiwipang$V1) #检查文件名
pop(df)<-pop.Kiwipang$V2 #指定群体
df$pop #查看群体
mydf<-genind2hierfstat(df) #genind2hierfstat
popFst1<-genet.dist(mydf,diploid=F) #计算pairwise Fst
popFst<-as.matrix(popFst1) #转为矩阵
pheatmap::pheatmap(popFst,cluster_rows = F,cluster_cols = F) #作图