一 首先用BayesScan删除受选择的位点
1.输入文件:PGDSpider将VCF 转为GESTE/BayScan格式,即下面的bates.xml
2.运行命令
cd /home/yl/biosoft/BayeScan-master/my_file/
/home/lsq/BayeScan2.1/asia-usnp-R3/BayeScan2.1_linux64bits /home/yl/biosoft/BayeScan-master/input_examples/**bates.xml** -n 100 -thin 50 -burn 200 -pr_odds 10000 -threads 8
3.输出文件在XXX_fst.txt中,q值小于0.01定义为受选择位点
二 利用easySFS生成SFS文件
1.位点谱频(SFS),又称等位基因频谱,用来显示特定基因座上各等位基因所占的频率。
2.unfold 还是 fold: 如果祖先状态已知为unfold(我理解是有参考基因组或者有外类群);反之为fold(比如简化基因组测序得到的数据)
3.命令:
安装详见 >https://github.com/xiaoming-liu/stairway-plot-v2
# 运行命令
conda activate easySFS ### 好像是运行前先改变环境吧?要在终端前面是(easySFS),如下一行所示才能运行
### (easySFS) yl@zhab001:~/biosoft/easySFS-master/BSZ-ZY$ 终端前面是(easySFS)
python /home/yl/biosoft/easySFS-master/easySFS.py -i zy_bsz.vcf -p zy_pops.txt --preview -a
python /home/yl/biosoft/easySFS-master/easySFS.py -i zy_bsz.vcf -p zy_pops.txt -a --proj=38 # --proj=38是根据上一步preview的结果得到的,比如我的个体是20个,一般就可以默认 --proj=40,但是上一步preview结果显示序号后面的数字在后面急剧变小的话,那就设定为38。简单来说:通常可设定个体n的2倍,2n;若后面存在剧降趋势则选择开始下降时的值。建议若区别不大,还是用自己个体的2倍吧,不然之后stairwploti一步,sfs的数字个数多次没弄对,想死,如果这里减了,之后nseq也得减
三 stairWay-plot画图
输入文件详见 https://github.com/xiaoming-liu/stairway-plot-v2/blob/master/READMEv2.1.pdf
使用easySFS中生成fastsimcoal2中的结果,对于unfoled和folded,easySFS都会生成2n+1个值,unfolded同样需要去掉第一个和最后一个,foled因为祖先状态未知,n+1之后的数值为0,用于stairway plot 2需要去掉第一个和所有的0。
cd /home/yl/biosoft/stairway_plot_v2.1.1/
java -cp stairway_plot_es Stairbuilder bsz-zy.blueprint
#生成.sh文件用于下一步
bash bsz-zy.blueprint.sh
参考
https://mp.weixin.qq.com/s/GS_g1aPYOSVIcVOljrOzJA
https://www.jianshu.com/p/73202f89f487
https://www.jianshu.com/p/08efb118cd00
吐槽一句,折腾一个周差点要了我老命,还好有师姐指点,qu愿我的下午能画出图来