PopSizeABC的原理：
1.根据观测数据（真实数据）估算SFS和LD
2.考虑了大量随机种群历史变化情景。对于每个历史变化情景，都会模拟许多序列。这些序列可以反映相应的历史（基于溯祖理论）。
3.然后，估计这些模拟序列的SFS和LD。
4.计算模拟数据和真实数据的汇总统计（SFS和LD）之间的欧氏距离。
5.然后，我们可以将最相似的模拟数据的历史映射到真实数据上，就可以得到种群历史。
参考文献：Boitard, S., Rodriguez, W., Jay, F., Mona, S., & Austerlitz, F. (2016). Inferring population size history from large samples of genome-wide molecular data-an approximate Bayesian computation approach. PLoS genetics, 12(3), e1005877.
软件网址：https://forge-dga.jouy.inra.fr/projects/popsizeabc/
前提是你有一个VCF文件，没有call出来的位点不需要，和reference一样的纯合位点也不需要。
1. split vcf 文件
压缩
bgzip -c 02.Var.F2BV1.vcf > 02.Var.F2BV1.vcf.gz
建立引索
tabix -p vcf 02.Var.F2BV1.vcf.gz
分割
for i in $(<chr_list)
do
{
bcftools view -r ${i} 02.Var.F2BV1.vcf.gz -o ${i}.vcf.gz -Oz
}
done
2. 改脚本（参数和命名）
有很多要改的，这里可能没有列全
for i in *.py;do sed -i 's/Jersey/KN/g' ${i};done
for i in *.py;do sed -i 's/cattle/myanimal/g' ${i};done
for i in *.py;do sed -i 's/n=30/n=32/g' ${i};done
（n：单倍体数量，例如二倍体生物就是个体数量乘以二）
写ped文件
KN A01-1ND 0 0 0 -999
KN A02-3H 0 0 0 -999
KN A04-3H 0 0 0 -999
KN B07-2H 0 0 0 -999
KN B09-H 0 0 0 -999
KN B10-1 0 0 0 -999
KN B11-1H 0 0 0 -999
KN C06-3H 0 0 0 -999
KN C07-3H 0 0 0 -999
替换
simul_data_ex1.py 中的haploid sample size
for i in *.R;do sed -i 's/estim/estim_for_myanimal/g' ${i};done
for i in *.R;do sed -i 's/Jersey/KN/g' ${i};done
for i in *.R;do sed -i 's/cattle/myanimal/g' ${i};done
for i in *.R;do sed -i 's/n=30/n=32/g' ${i};done
for i in *.R;do sed -i 's/n30/n32/g' ${i};done
其中 abc_cv_ex1.R里有一行：
par.estim_for_gibbon=(res$estim)[[1]]
后面这个res$estim这个别改
把abc_cv_ex1.R里的res/ex1_...stat改成ex1_...stat
params同理
记得改generation.R里的generation time
调整prior数值
原来
现在
要注意VCF file里面有多少条染色体就要设置多少个segment，否则差异很大
运行1
mkdir res
python2.7 comp_stat1/stat_from_vcf_ex1.py
运行 2
python2.7 comp_stat1/simul_data_ex1.py
因为我是分了很多个脚本并行模拟的，所以把模拟好的文件直接cat到一起。如果你就模拟了一次就不用cat
cat res/*.params > ex1_simu_all_n32_s1.params
cat res/ex1_simu*.stat > ex1_simu_all_n32_s1_mac6_macld6.stat
运行 3，abc计算并画图
R -f estim/abc_ex1.R
R -f estim/abc_cv_ex1.R
（这个R脚本里的输入文件名也要根据你cat的文件名修改）
结果：
image.png
image.png
图表背后的数据如果要输出的话在R里加一行write.table就行了
网友评论