【群体遗传】Fst（群体间分化指数）

（1）FSTF_{ST}FST是什么？含义是什么？

FSTF_{ST}FST，全称为fixation index，是一种用于衡量群体间分化程度的统计检验量（由Wright’s F-statistics衍生而来）。一般从SNP或microsatellites数据计算得到，且一般用在群体遗传学分析中。

microsatellites，即微卫星序列，是在一种串列重复序列 —— https://en.wikipedia.org/wiki/Microsatellite
但是现在WGS和RAD-Seq都已经非常普遍了，使用的是否还多我也不了解，就略了~

（2）FSTF_{ST}FST如何计算？

示例：FSTF_{ST}FST计算原理

	AA	Aa	aa
Pop1	125	250	125
Pop2	50	30	20
Pop3	100	500	400

1、统计每一个群体的等位基因数量

每一个Pop对应的基因型（genotype | genotyped individuals）数量为：

500
100
1000

每一个Pop的等位基因数量（the number of allele）为：

1000
200
2000

这边是biallelic类型（A or a），因此等位基因数量为基因型数量的2倍。

2、计算每一个群体实际的等位基因频率

Pop1中，

A allele实际的频率为125∗2+250∗11000\frac{125*2 + 250*1}{1000}1000125∗2+250∗1，即0.5
a allele实际的频率为 1- 0.5 = 0.5

Pop2中，

A allele实际的频率为2∗50+30∗1200\frac{2*50 + 30*1}{200}2002∗50+30∗1，即0.65
a allele实际的频率为 1- 0.65 = 0.35

Pop3中，

A allele的频率为100∗2+5002000\frac{100*2 + 500}{2000}2000100∗2+500，即0.35
a allele的频率为 1- 0.35 = 0.65

3、计算每一个群体期望的基因型数量 & 差值

【标注】期望，即服从HD平衡理论，可以看看北京大学生物演化课程

Pop1中，

AA genotype期望的频率为125
Aa genotype期望的频率为250
aa genotype期望的频率为125

因此，Pop1中对应的基因型数量均无偏差。

Pop2中，

AA genotype期望的频率为42.25
Aa genotype期望的频率为45.5
aa genotype期望的频率为12.25

对应基因型数量的差值为+7.25, -15.5, +7.75。

Pop3中，

AA genotype期望的频率为122.5
Aa genotype期望的频率为455
aa genotype期望的频率为422.5

对应基因型数量的差值为-22.5, +45, -22.5。

对计算结果的理解，Pop1与计算得到的期望数值一样，服从HD平衡；Pop2实际纯合基因型数目与期望纯合基因型数目差值为正，表明存在inbreeding（近亲繁殖）事件；Pop3实际纯合基因型数目与期望纯合基因型数目差值为负，表明存在outbred事件，即亚群之间的isolation（生殖隔离）被打破，导致亚群之间能够产生后代。

4、统计每一个群体实际的杂合基因型占比

Pop1为0.5，Pop2为0.3，Pop3为0.5

【公式标注】Hobs=杂合基因型数目总个体数H_{obs} = \frac{杂合基因型数目}{总个体数}Hobs=总个体数杂合基因型数目

5、计算每一个群体期望的杂合基因型占比

Pop1为0.5，Pop2为0.455，Pop3为0.455

【公式标注】Hexp=1−∑(p2+q2)H_{exp} = 1-\sum(p^2 + q^2)Hexp=1−∑(p2+q2)

6、计算A allele的频率均值

p‾=2∗125+250+2∗50+30+2∗100+5001000+200+2000\overline{p} = \frac{2*125 + 250 + 2*50 + 30 + 2*100 + 500}{1000 + 200 + 2000}p=1000+200+20002∗125+250+2∗50+30+2∗100+500，即0.4156

7、计算a allele的频率均值

1−p‾=q‾1 - \overline{p} = \overline{q}1−p=q，即0.5844

8、计算the global heterozygosity indices

1.首先使用HobsH_{obs}Hobs计算HIH_{I}HI

HI=Hobs1∗N1+Hobs2∗N2+Hobs3∗N3NtotalH_{I}=\frac{H_{obs1}*N_{1} + H_{obs2}*N2 + H_{obs3}*N3}{N_{total}}HI=NtotalHobs1∗N1+Hobs2∗N2+Hobs3∗N3，带入数值，即0.4875

2.使用H_{exp}计算H_{S}

HS=Hexp1∗N1+Hexp2∗N2+Hexp3∗N3NtotalH_{S} = \frac{H_{exp1}*N_{1} + H_{exp2}*N2 + H_{exp3}*N3}{N_{total}}HS=NtotalHexp1∗N1+Hexp2∗N2+Hexp3∗N3，带入数值，即0.4691

3.计算global heterozygosity indicex的期望值

HT=1−∑(p‾2+q‾2)=1−(0.41462+0.58442)H_{T} = 1 - \sum(\overline{p}^2 + \overline{q}^2) = 1 - (0.4146^2 + 0.5844^2)HT=1−∑(p2+q2)=1−(0.41462+0.58442)，即0.4845

9、计算the global F-statistics

1.计算FIS=HS−HIHSF_{IS} = \frac{H_{S} - H_{I}}{H_{S}}FIS=HSHS−HI，即-0.0393
2.计算FST=HT−HSHTF_{ST} = \frac{H_{T} - H_{S}}{H_{T}}FST=HTHT−HS，即0.0344
3.计算FIT=HT−HIHTF_{IT} = \frac{H_{T} - H_{I}}{H_{T}}FIT=HTHT−HI，即-0.0036

10、计算结果说明了什么？

群体间分化的程度达到了3.4%

示例：vcftools计算F_{ST}

【标注】只适用于二倍体。

vcftools --gzvcf input.vcf.gz --weir-fst-pop pop1_sample_id.txt --weir-fst-pop pop2_sample_id.txt --fst-window-size 10000 --fst-window-step 10000 --out pop1_pop2# 参数说明
--gzvcf            # 要求输入为.gz格式的vcf文件
--weir-fst-pop     # 输入VCF文件中的sample，为一个文本文件，每一行一个sample
--fst-window-size  # 设置计算Fst的窗口大小，根据自己的数据进行设置，看看别人文章里怎么用的
--fst-window-step  # 设置计算Fst的步长长度，根据自己的数据进行设置

（4）FSTF_{ST}FST计算完了之后该干啥？

在对两个群体之间进行不同区段的FSTF_{ST}FST计算之后，需要判断哪一些区段，是“真正”受到了选择压力，根据近期看的文章，得到可以选择前5%的作FSTF_{ST}FST为一个阈值，对区域进行划分，高于该阈值的被认为受到了选择压力的影响，进一步就可以得到是受到影响的是哪些SNP，最终即可得到受到影响的是哪些gene。

当然，对FSTF_{ST}FST的计算结果可视化，当然也是非常重要的一部分，但是这篇文章主要想写的是计算原理以及如何使用vcftools进行计算。

参考资料

[1] https://en.wikipedia.org/wiki/Fixation_index
[2] http://www.uwyo.edu/dbmcd/popecol/maylects/fst.html
[3] The genome of oil-Camellia and population genomics analysis provide insights into seed oil domestication