前期准备：

1 beagle填充，还有所谓的phasing（windows 版本）

java -jar beagle.08Feb22.fa4.jar gt=B1B2.vcf out=B1B2-phase.vcf ne=2

2 提取群体样品

vcftools --gzvcf B1B2-phase.vcf.vcf.gz
--keep B1.txt
--recode
--recode-INFO-all
--out B1.beagle

vcftools --gzvcf B1B2-phase.vcf.vcf.gz
--keep B2.txt
--recode
--recode-INFO-all
--out B2.beagle

由于是群体间的，因此需要将两个群体的vcf文件分别提取出来，后面只用其中一个做展示，另一个只是修改输入输出文件名即可
！！！注意：第3-5步只用其中一个群体计算遗传距离即可，是通用的

开始分析：

1 拆分染色体(B1.chr.sh )

for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf B1.beagle.recode.vcf
--recode
--recode-INFO-all
--chr {i};
done

1 拆分染色体(B2.chr.sh )

for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf B2.beagle.recode.vcf
--recode
--recode-INFO-all
--chr {i};
done

2 准备分染色体的map文件(2_map.sh)

for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf B1.chr{k}.MT;
done

3 map计算遗传距离(3_map.distance.sh)

for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk 'BEGIN{OFS=" "} {print 1,".",4}' chr{k}.MT.map.distance;
done

4 计算XPEHH(4_xpehh.sh)

for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do /mnt/d/peng/fst/selscan-master/bin/linux/selscan --xpehh
--vcf B1.chr{k}.recode.vcf
--map chr{k}.B2_B1;
done
得到 chr${k}.EU_JBC文件29个
重测序数据直接计算XPEHH（对每个100K窗口单倍型进行打分，标准化使其正态分布，并统计值>2的SNP所占的比例） 相对芯片来说步骤少点
最终得到的标准化的xpehh值表示的是xpehh>2的SNP在此区间所占的比例
具有方向性，xpehh>2表示目标群体的受选择方向，xpehh<-2表示参考群体的受选择方向

5 第一列改为染色体号(5_add.chr.sh)

并将空格改为tab键分割
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '2,4,6,8}' chr{k}.B2_B1.xpehh.out;
sed -i 's/ /\t/g' Chr{k}',3,5,7,{k}.B2_B1.xpehh.out > Chr{k}.B2_B1.xpehh-1.out > B2_B1.xpehh.out-2;
done

6 加50kb窗口且标准化(6_norm.xpehh.sh)

可按需求修改窗口大小
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do /mnt/d/peng/fst/selscan-master/bin/linux/norm --xpehh
--files Chr${k}.B2_B1.xpehh-1.out
--bp-win --winsize 50000;
done

产生两种文件.100bins.norm文件和.100bins.norm.50kb.windows文件
.100bins.norm.50kb.windows文件的表头为
<win start> <win end> <scores in win> <gt the fraction of XP-EHH scores >2> <lt the fraction of XP-EHH scores < -2> <approx percentile for gt threshold wins> <approx percentile for lt threshold wins> <max score> <min score>

7 提取并合并结果文件(7_extract.merge.all.sh)

for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '1,4,8,{k}.EU_JBC.xpehh.out.norm.50kb.windows > Chr${k}.EU_JBC.chart.xpehh.norm.50kb;
cat ./*.EU_JBC.chart.xpehh.norm.50kb > all.xpehh.norm.50kb;
done