2019-12-20Fst的计算原理与实战

来源https://www.jianshu.com/p/b73a8d6233be

概念回顾

Fst:群体间遗传分化指数,是种群分化和遗传距离的一种衡量方法,分化指数越大,差异越大。适用于亚群体间多样性的比较。

用于衡量种群分化程度,取值从0到1,为0则认为两个种群间是随机交配的,基因型完全相似;为1则表示是完全隔离的,完全不相似。它往往从基因的多样性来估计,比如SNP或者microsatellites(串联重复序列一种,长度小于等于10bp)。是一种以哈温平衡为前提的种群遗传学统计方法。

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

具体可以下面的公式表示:

Fst= (Ht-Hs)/ Ht

Hs:亚群体中的平均杂合度

Ht:复合群体中的平均杂合度

理论上计算Fst的步骤

理论上要估算FST,需要以下步骤:

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

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

计算每个亚群的杂合度(2pq)

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

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

最后,计算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是什么?

手工计算的步骤如下:

在得到每个群体中等位基因的频率后:

首先第一步是计算每个亚群的杂合度H,这一步简单H=2pq (p和q是对应群体中等位基因的频率)

第二步计算这些亚群杂合度的平均值HS,取三个不同种群中的H值,除去平均数3,得到结果

第三步,计算群体等位基因频率,和上述方法类似,取每个等位基因在对应种群中的频率叠加除去平均数。

第四步,计算所有群体中的杂合度HT,使用2pq,p和q是群体中的等位基因频率,在第三步骤已经算好了。

最后,根据计算公式算出Fst。

在实际中,当然不可能用手工对每个位点进行这样的计算,怎样在电脑中分析FST等值,会在下次内容中介绍。

实战

当然了,实际中我们并不需要像上述例子那样逐个计算。等位基因频率的信息都藏在snp calling后的vcf中,再使用恰当的工具就可以快速计算出Fst。下面给大家介绍一种比较popular的计算Fst的方法(还有其他方法不限于一种),使用vcftools:

##对每一个SNP变异位点进行计算vcftools--vcf test.vcf--weir-fst-pop1_population.txt--weir-fst-pop2_population.txt--outp_1_2—single##按照区域来计算vcftools--vcf test.vcf--weir-fst-pop1_population.txt--weir-fst-pop2_population.txt--outp_1_2_bin--fst-window-size500000--fst-window-step50000# test.vcf是SNP calling 过滤后生成的vcf 文件;# p_1_2_3 生成结果的prefix# 1_population.txt是一个文件包含同一个群体中所有个体,一般每行一个个体。个体名字要和vcf的名字对应。# 2_population.txt 包含了群体二中所有个体。#计算的窗口是500kb,而步长是50kb (根据你的需其可以作出调整)。我们也可以只计算每个点的Fst,去掉参数(--fst-window-size 500000 --fst-window-step 50000)即可。

分别对不同结果进行图形绘制

##图1library(ggplot2)data<-read.table("test1.out.windowed.weir.fst",header=T)sc3=subset(data,CHROM=="Gm01")p<-ggplot(sc3,aes(x=BIN_END/1000000,y=WEIGHTED_FST))+geom_point(size=0.5,colour="blue")+xlab("Physical distance (Mb)")+ylab("Fst")+ylim(-1,1)p+theme_bw()##图2library(ggplot2)data<-read.table("test.out.weir.fst",header=T)sc3=subset(data,CHROM=="Gm01")p<-ggplot(sc3,aes(x=POS,y=WEIR_AND_COCKERHAM_FST))+geom_point(size=0.5,colour="blue")+xlab("Physical distance (Mb)")+ylab("Fst")+ylim(-1,1)p+theme_bw()

可以看到两幅图还是基本对应的,然后对一些接近1,与大部分点偏离比较高的点可以将其与功能注释相结合,还有一些选择压力分析的工具相结合,寻找出其对应的基因,观察该基因是否是被选择。Pi也是选择分析中一个很要的参数,这一部分内容有时间这一点以后再补回来。

参考文章:

https://blog.csdn.net/u014182497/article/details/52672308

http://www.omicshare.com/forum/thread-3688-1-1.html

作者:lakeseafly

链接:https://www.jianshu.com/p/b73a8d6233be

来源:简书

著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,222评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,455评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,720评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,568评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,696评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,879评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,028评论 3 409
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,773评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,220评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,550评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,697评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,360评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,002评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,782评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,010评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,433评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,587评论 2 350

推荐阅读更多精彩内容