群体遗传学统计指标——群体间分歧度检验（Fst）

Posted 2023-04-12

参考技术A Fst ：群体间遗传分化指数，是种群分化和遗传距离的一种衡量方法，分化指数越大，差异越大。适用于亚群体间多样性的比较。用于衡量种群分化程度，取值从0到1，为0则认为两个种群间是随机交配的，基因型完全相似；为1则表示是完全隔离的，完全不相似。它往往从基因的多样性来估计，比如SNP。是一种以哈温平衡为前提的种群遗传学统计方法。

在群体遗传学中衡量群体间的遗传分化的程度的指标有许多种，较为常见的就是 遗传分化指数（Fst） ，fst是由F统计量演变而来，F统计量主要有三种（FIS，FIF，FST）。Fst是针对一对等位基因，如果基因座上存在复等位基因，则需要用Gst衡量，基因差异分化系数（gene differentiation coefficient，Gst）。假定有s个地方群体，第k个地方群体相对大小为wk，k地方群体中第i个等位基因频率为qk(i)，杂合体频率观察值为hk，那么，整个群体中观察到的杂合体频率平均值HI，地方群体为理想群体的期望杂合体频率平均值HS，整个群体为理想群体的期望杂合体频率HT，分别为：

FIS，是HI相对于HS减少量的比值，即地方群体的平均近交系数。

FST，是HS相对于HT减少量的比值，即有亲缘关系地方群体间的平均近交系数。

FIT，是HI相对于HT减少量的比值，即整个群体的平均近交系数。

简单来说，两个种群之间遗传差异的基本测量就是统计量FST。在遗传学中，F一词通常代表“近亲繁殖”，它倾向于减少群体中的遗传变异。遗传变异可以用杂合度来衡量，所以F一般表示群体中杂合性的减少。 FST是与它们所属的总群体相比，亚群体中杂合性的减少量。

Fst值的取值范围是【0,1】，最大值为1表明两个群体完全分化，最小值为0表明群体间无分化。

在实际的研究中Fst值为0--0.05时说明群体间遗传分化很小，可以不做考虑；

为0.05--0.15时，表明群体间存在中等程度的遗传分化；

为0.15--0.25时群体间存在较大的遗传分化；

为0.25以上的时候群体间就存在很大的遗传分化了。

1）找出每个亚群的等位基因频率。

2）查找复合群体的平均等位基因频率

3）计算每个亚群的杂合度（2pq）

4）计算这些亚群杂合度的平均值，这是HS。

5）根据总体等位基因频率计算杂合度，这是HT。

6）最后，计算FST =（HT-HS）/ HT

基因SLC24A5是黑色素表达途径的关键部分，其导致皮肤和毛发色素沉着。与欧洲较轻的皮肤色素密切相关的SNP是rs1426654。 SNP有两个等位基因A和G，其中G与轻度皮肤相关，在犹他州的欧裔美国人中，频率为100％。美洲印第安人与美国印第安人混血儿的SNP在频率上有所不同。墨西哥的样本有38％A和62％G;在波多黎各，频率分别为59％A和41％G，查尔斯顿的非裔美国人样本中有19％A和81％G.这个例子中的FST是什么？

看了眼花缭乱的理论知识，感觉云里雾里的应该是正常的感受。总结一下，其实Fst分析就是看两个群体之间分化程度的一种方法，Fst值越大（越接近1）表明两个群体间分化程度越高，亲缘关系越远；Fst值越小（越接近0）表明群体间分化程度越低，亲缘关系越近。

一般我们用于分析的序列都很长，样本量也很大，靠笔算可能比西天取经还难。所以我们用强大的VCF分析工具——vcftools软件来计算Fst值。

SNP单点计算

按照区域（窗口模式）计算

其中--vcf 是输入vcf格式文件，里面保存的是全部样品的位点信息，具体细节大家自己查一查。

--weir-fst-population 这个命令是告诉计算机需要计算的群体ID名，该文件必须是txt格式，每个ID占一行，每次计算两个群之间的关系，因此需要population_1.txt和population_2.txt两个文件。

单点计算就是每个SNP位点单独计算，窗口模式就是按照预先设定的步长，多个SNP位点计算一次，两者的区别就是在后面加上窗口的大小和步长参数：--fst-window-size 500000 --fst-window-step 50000 ，这里窗口设置为500kb，步长设置为50kb，窗口的设置没有一个固定的标准和要求，都是按照自己的需要而定。

数据可视化换句话说就是清（hua）楚（shi）的展示你的结果。回头看看，在Fst分析中，Fst值越大（越接近1）表明两个群体间分化程度越高，亲缘关系越远；Fst值越小（越接近0）表明群体间分化程度越低，亲缘关系越近。那么我们所画的图，应该要清晰的展示两个群体Fst值的大小。因此，我们会选择散点图，它可以展示在整个基因组上，每个区域两个群体Fst值的大小。当然，如果你要比较多个群体之间的关系，也可以把整个基因组看作一个整体，用箱图展示群体之间的Fst值大小关系。

废话太多，先上个某文章的经典图片：

文章对它的描述如下：A major genomic signature of modern processing tomatoes and three causative variants for pink fruit. FST values for all SNP sites between tomatoes for fresh consumption and modern processing tomatoes. Blue dots above the horizontal dashed line indicate highly divergent SNPs (top 1%; FST = 0.4464).

我是用R画的，当然你也可以直接把数据导入EXCEL，用EXCEL画也是非常漂亮的。

我的数据整理后长这样，用到的是窗口模式返回结果的WEIGHTED_FST列数据（这里我重新命名为Sample）：

分享一下我写得一个R：（大家需要自己根据自己的数据就行调整，但是万变不离其中，你们可以的！）

Vcftools Manual

Genomic analyses provide insights into the history of tomato breeding

Fst的计算原理与实战

重测序（RADseq）做群体遗传分析套路

参考技术A

构建的群体，或自然群体，如各地方品种。

提取DNA后，构建文库，简要步骤如下：
① 限制性内切酶TaqI酶切；
② 连接P1接头；
③ DNA随机打断片断化；
④ 目的片段回收与末端修复；
⑤ 连接P2接头；
⑥ RAD片段富集；
⑦ 上机测序。
参考： Rapid and cost-effective polymorphism identification and genotyping using restriction site associated DNA (RAD) markers

根据识别标签序列得到每个个体的测序reads，使用trimmomatic进行过滤（其他质控软件，如fastqc，multiQC等）
设置过滤参数为：SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50。 过滤标准：两端质量低于5的碱基进行切除，并以5bp为窗口进行滑动过滤，对平均质量低于20的窗口进行切除。

BWA （其他比对软件如bowtie2/soap2/MAQ等）将过滤后的个体clean reads比对到参考基因组序列上。样本比对率反映的是样本测序数据与参考基因组的相似性，覆盖深度和覆盖度能够直接反映测序数据的均一性与参考序列的同源性。

使用GATK（或samtools+bcftools)Haplotype Caller模块进行变异检测，获得群体变异集文件（VCF 格式）。对变异进行过滤：过滤参数为缺失率小于或等于0.2、杂合率小于或等于0.2、最小等位基因频率（MAF） 大于或等于0.05，最终得到高质量的基因型数据。

在获得高质量的标记数据以后，利用vcftools将vcf文件处理得到plink.ped和plink.map文件（整理为plink软件所需格式）。

使用plink 软件随机选择连锁不平衡（LD）小于0.1，且相邻间隔在300kb以上的SNP位点，最后得到一个包含3420 个SNP位点的标记集，一般是生成.bed文件。

1.祖先成分堆叠图
使用ADMIXTURE对此 SNP位点集（bed文件）进行群体结构分析（Structure），利用交叉验证过程确定确定合适的祖先数或亚群（K值）。若不知道理想的K值，可用ADMIXTURE计算，一般当cross-validation error值最低时所对应的K值为最合适的K值。

考虑到样本所归属的分类单元，即看看哪几个物种聚在一起，对合适的K值利用Structure软件（速度慢，其他软件如frappe，ADMIXTURE也可做群体结构图，并且很快）聚类图，一些R包如hapmap也是可以做群体结构图的。

2.PCA
利用GCTA对SNP数据集进行样本的PCA分析（其他软件如EIGENSOFT中的smartpca）。GCTA可以直接读取.bed , .bim , .fam文件，利用–make-grm 生成个体对之间的遗传关系矩阵，并将GRM的下三角元素保存为二进制文件.grm.id , .grm.bin ， .grm.N.bin。使用 –pca 设置要生成主成分的数目，一般来说就可以刻画出群体结构。这一步会生成 .eigenval 和 .eigenvec 两个文件。.eigenval文件为各主成分可解释遗传信息的比例，.eigenvec文件为每个样本在top4主成分上的分解值。

3.系统发育树
构树的方法有非加权分组平均法（UPGMA，已经很少用）、最小进化法（ME）、邻接法（NJ）、最大简约法（MP）、最大似然法（ML）等。

构树软件如FastTree/MEGA/cluster X/phylip，美化可以用FigTree/ggtree/treeview/GraPhIAn。

NJ法是基于最小进化原理经常被使用的一种算法，它不检验所有可能的拓扑结构，能同时给出拓扑结构和分支长度。

GWAS的群体遗传分析也是包含这三个图，RADseq毕竟是简化基因组，得到的SNP有限，做这种群体分析效果肯定没有GWAS好。