(1) F S T F_{ST} FST是什么?含义是什么?

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

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

(2) F S T F_{ST} FST如何计算?

示例: F S T F_{ST} FST计算原理
AAAaaa
Pop1125250125
Pop2503020
Pop3100500400
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 ∗ 1 1000 \frac{125*2 + 250*1}{1000} 10001252+2501,即0.5

  • a allele实际的频率为 1- 0.5 = 0.5

Pop2中,

  • A allele实际的频率为 2 ∗ 50 + 30 ∗ 1 200 \frac{2*50 + 30*1}{200} 200250+301,即0.65

  • a allele实际的频率为 1- 0.65 = 0.35

Pop3中,

  • A allele的频率为 100 ∗ 2 + 500 2000 \frac{100*2 + 500}{2000} 20001002+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

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

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

Pop1为0.5,Pop2为0.455,Pop3为0.455

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

6、计算A allele的频率均值

p ‾ = 2 ∗ 125 + 250 + 2 ∗ 50 + 30 + 2 ∗ 100 + 500 1000 + 200 + 2000 \overline{p} = \frac{2*125 + 250 + 2*50 + 30 + 2*100 + 500}{1000 + 200 + 2000} p=1000+200+20002125+250+250+30+2100+500,即0.4156

7、计算a allele的频率均值

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

8、计算the global heterozygosity indices

1.首先使用 H o b s H_{obs} Hobs计算 H I H_{I} HI

H I = H o b s 1 ∗ N 1 + H o b s 2 ∗ N 2 + H o b s 3 ∗ N 3 N t o t a l H_{I}=\frac{H_{obs1}*N_{1} + H_{obs2}*N2 + H_{obs3}*N3}{N_{total}} HI=NtotalHobs1N1+Hobs2N2+Hobs3N3,带入数值,即0.4875

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

H S = H e x p 1 ∗ N 1 + H e x p 2 ∗ N 2 + H e x p 3 ∗ N 3 N t o t a l H_{S} = \frac{H_{exp1}*N_{1} + H_{exp2}*N2 + H_{exp3}*N3}{N_{total}} HS=NtotalHexp1N1+Hexp2N2+Hexp3N3,带入数值,即0.4691

3.计算global heterozygosity indicex的期望值

H T = 1 − ∑ ( p ‾ 2 + q ‾ 2 ) = 1 − ( 0.414 6 2 + 0.584 4 2 ) 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.计算 F I S = H S − H I H S F_{IS} = \frac{H_{S} - H_{I}}{H_{S}} FIS=HSHSHI,即-0.0393
2.计算 F S T = H T − H S H T F_{ST} = \frac{H_{T} - H_{S}}{H_{T}} FST=HTHTHS,即0.0344
3.计算 F I T = H T − H I H T F_{IT} = \frac{H_{T} - H_{I}}{H_{T}} FIT=HTHTHI,即-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) F S T F_{ST} FST计算完了之后该干啥?

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

当然,对 F S T F_{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
image.png

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐