我终于会算Fst了

估计是过节结束了，最近公众号的粉丝增加终于回到了正常水平了，尽管依旧很慢。不过不着急，反正很多人加了也是处于不会打开的状态，我还是安心更新吧。

上游分析流程得到的变异信息总需要进一步的数据挖掘才能找到有趣的故事，但是首先我们得知道可以用什么工具进行什么类型的分析。这篇主要介绍stacks/populations

对了，由于我也是新手，希望大家能给我推荐一些经典的群体遗传学教程。

什么是Fst

2年前（也就是2016年,时间过的真快。。），我在现在的实验室轮转的时候，师姐想让我帮她算Fst。但是那个时候啥都不懂，未知的东西实在是太多了，所以到最后我都没算出来，那个时候我最大的贡献就是帮她装了几盘土种苗。

首先我得知道什么叫做Fst。Fst是群体遗传学的一个统计值，是F-统计值家族的一员，和F检验没啥关系。根据维基百科

Fst(Fixation index,固定指数)是衡量因遗传结构(genetics structure)而引起群体差异(population differentiation), 是Sewall Wright提出的F-统计值的特例，可以说是目前群体遗传学最常用的统计值。

Sewall Wright是现代群体遗传学的奠基者之一，其余两人是Ronald Fisher（费希尔）和 J.B.S. Haldane。群体遗传学好多概念就是他提出的，比如说近交系数，遗传漂变等。

之前只用了stacks的populations进行snp过滤和数据导出为常用格式，但其实它还可以计算群体遗传统计值如期望/观测杂合率，π(核苷酸多样性)， Fis（个体,individual相对亚群subpopulation的近交系数），还可以逐对比较所有群体计算Fst，用来判断两个群体间的差异，还可以计算基于单倍体的群体遗传统计值，如单倍型多样性， ΦST, and FST’.

第一步，获取VCF文件。其实stacks是可以根据sstacks处理后的数据进行分析，但是VCF可以说是目前许多文件格式的中转站，用的也比较广泛，为了方便后续分析和协作分析，最好还是先得到一个VCF吧。

min_samples=0.80min_maf=0.05max_obs_het=0.80populations -P 03-stacks-analysis/ref-based/ -r $min_samples --min_maf $min_maf \--max_obs_het $max_obs_het --vcf &> populations_ref_based.oepopulations -P 03-stacks-analysis/de-novo/ -r $min_samples --min_maf $min_maf \--max_obs_het $max_obs_het --vcf &> populations_de_novo.oe

最后会在03-stacks-analysis/ref-based/和03-stacks-analysis/de-novo/下都有一个batch_1.vcf，这就是后续用来分析的基础。分别改名成de_novo.vcf和ref_based.vcf便于区分。

populations支持导出的格式非常的多，如fasta, fasta_strict, vcf, vcf_haplotypr, structure,phase,fastphase,beagle,beagle_phased,plink,hzar,phylip,phylip_var,phylip_var_all,treemix

第二步，声明用于配对比较的群体. 这一步生成-M,--popmap输入文件。比如说比较淡水鱼(pcr,wj)和海水鱼（cs,sj)，

sed -r 's/\t(cs|sj)/\tocceanic/; s/\t(pcr|wc)/\tfreshwater/;' info/popmap.tsv \> info/popmap.oceanic_freshwater.tsv

第三步，计算群体遗传统计值

popmap=info/popmap.oceanic_freshwater.tsvpopulations -V ref_based.vcf -M $popmap -O ./ -p 2 \--fstats -k --sigma 100000

对于有参考基因，可以使用-k计算kernel-smoothed π, FIS, FST, FST’, and ΦST，至于有啥意义，我目前还不懂。此外还可计算重抽样的几个统计值，对应--bootstrap_xx。启用后计算时间会明显增加。最后会得到如下文件

• ref_based.fst_occeanic-freshwater.tsv

• ref_based.fst_summary.tsv

• ref_based.haplotypes.tsv

• ref_based.hapstats.tsv

• ref_based.phistats.tsv

• ref_based.phistats_occeanic-freshwater.tsv

• ref_based.sumstats.tsv

• ref_based.sumstats_summary.tsv

第四步，R可视化展示. 结果必须可视化才能便于我们找到数据里有意思的地方。这一步要用到上一步得到的ref\_based.phistats\_occeanic-freshwater.tsv

x = read.delim('ref_based.fst_occeanic-freshwater.tsv')x.vii = subset(x, Chr=='groupVII')plot(x.vii$BP, x.vii$AMOVA.Fst,pch=3, cex = 0.5,xaxt='n',xaxs= 'i',yaxs = 'i',xlab='groupVII Position (Mb)',ylab='Fst')lines(x.vii$BP, x.vii$Smoothed.AMOVA.Fst,col='blue')axis(1,at=seq(0,ceiling(max(x.vii$BP)/1000000)*1000000,1000000),labels = seq(0,ceiling(max(x.vii$BP)/1000000),1))

由于VCF文件还是比较小的，所以有兴趣的小伙伴可以到我分享的百度云盘下载。链接: https://pan.baidu.com/s/1gho2l7H 密码: q8gi

题外话，昨天Jimmy和我说有人在他出租的服务器里跑流程，感觉自己的教程没白写呀。当然他的服务器是64G内存，跑起来也是非常OK。