PopSizeABC：一百万年至今的种群历史模拟 demogra

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就行了