开始正文之前先了解一下为什么构建群体遗传结构是重要的，并且可以通过哪些方式构建此结构：

对于方法1，我们需要使用Plink构建PCA图

(1) 首先使用plink将vcf格式文件转换成.bed , .bim , .fam文件, 操作如下：

plink --allow-extra-chr -vcf XXX.vcf --make-bed --out filtered_XXX && echo "vcf2bed is done"

(2) 使用plink进行主成分分析，操作如下：

plink -bfile filtered_XXX --allow-extra-chr --pca 10 --out PCA_filtered_snps

产生了两个eigenval eigenvec文件，之后开始绘制PCA散点图，可以使用python脚本进行绘制，代码如下：

#!/usr/bin/python

# -*- coding: UTF-8 -*-

'''绘制PCA散点图'''

import matplotlib.pyplot as plt

import collections

import re

hash = collections.OrderedDict()

eval_file = open("PCA_filtered_snps.eigenval","r")

evec_file = open("PCA_filtered_snps.eigenvec","r")

### 从.eval文件中读取top2主成分

eval_1 = eval_file.readline()

eval_2 = eval_file.readline()

### 从.evec文件中读取在pc上的值,之所以使用re，是因为我的样本名是种+数字编号，例如“PTE01”,但我需要根据种进行聚类，所以使用r'[0-9]+|[A-Za-z]+'进行提取，这一次大家可以自行修改

for i in evec_file:

if re.match(r'#',i.strip()):

continue

tmp = i.strip().split()[1:]

if tmp[-1] in hash.keys():

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]]['1st'].append(eval(tmp[1]))

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]]['2nd'].append(eval(tmp[2]))

else:

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]] = hash.get(re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0],collections.OrderedDict())

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]]['1st'] = hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]].get('1st',[])

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]]['2nd'] = hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]].get('2nd',[])

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]]['1st'].append(eval(tmp[1]))

hash[re.findall(r'[0-9]+|[A-Za-z]+',tmp[0])[0]]['2nd'].append(eval(tmp[2]))

### 绘制散点图

fig = plt.figure(figsize=(12,10))

ax = fig.add_subplot(1, 1, 1) ### 2D figure

#ax = fig.add_subplot(111,projection='3d') ### 3D figure

mark = ['o','v']*10 ### 设置散点性状

col = ['blue','lime','y','red','cyan','fuchsia','orange']*2 ### 设置散点颜色

for m,n in enumerate(hash.keys()): ### 逐点画出

ax.scatter(hash[n]['1st'],hash[n]['2nd'],c=col[m],s=75,marker=mark[m],label=n,alpha=0.7)

ax.legend(loc='best',scatterpoints=1,fontsize='12',framealpha=0)

ax.set_xlabel('PC 1 ({}%)'.format(float('%.2f' % eval(eval_1))),fontsize=13,fontweight='bold')

ax.set_ylabel('PC 2 ({}%)'.format(float('%.2f' % eval(eval_2))),fontsize=13,fontweight='bold')

#### 保存图片

save_FileName = 'PCA_filtered_snps_pc1_pc2.png'

plt.savefig(save_FileName,dpi=400,bbox_inches='tight')

结果如下：

可见pil一类，ppr和peu一类，pte,ptre比较接近，剩下的聚成了一类。

对于方法2，我们使用Admixture构建进化关系

首先交叉验证选择最优的K值：

for K in 10 11 12 13 14 15; do admixture --cv filtered_XXX.bed $K | tee admixture_log${K}.out; done

结果如下：

CV error (K=10): 0.20719

CV error (K=11): 0.17233

CV error (K=12): 0.18846

CV error (K=13): 0.22279

CV error (K=14): 0.26081

CV error (K=15): 0.26839

因此选择K=11 is a sensible modeling choice。选择filtered_snps.11.Q作图，结果如下：

通过Ancestry的结果，可以分成三类：

第一类：ppr

第二类：pang, plas, pnig, plua

第三类：pde, pfre

第四类：ptr

第五类：ptre, pte

第六类：pba

第七类：peu

第八类：pil

作图代码：

1.tbl = read.table("filtered_XXX.11.Q")

2.barplot(t(as.matrix(tbl)),col=rainbow(11),xlab="Individual #",ylab="Ancestry",border=NA)

对于方法3，我们使用phylip和MEGA构建进化关系

首先需要将VCF格式转换成phylip所需的格式

python vcf2phylip.py -i XXX.vcf -o XXX.phy

之后使用MEGA将phy转换成meg格式

之后使用MEGA X建树，参数如下：

(No. of Bootstrap Replications开太多了会严重影响运行速度，开到100-500就可以)

在群体遗传里面，使用邻接法建树最好

等待建树：

结果如下：

根据此结果，结合上述两种方法得到的结果，我们最终可以分成四类：

第一类：pil

第二类：peu, ppr

第三类：pte, ptre

第四类：plas, pnig, pdel, pfre, pang, plau, pba, ptr