1. 用pixy计算群体的pi，dxy，fst

pixy可以基于多组学的vcf格式文件来计算群体的pi，dxy，fst值。

1.1. 安装pixy

1
2
3

conda install --yes -c conda-forge pixy
conda install --yes -c bioconda htslib
pixy --help

1.2. 准备vcf文件

pixy需要的vcf文件与只包含snp的vcf文件不同，是同时包含变异位点（variant sites）和非变异位点（invariant sites）的全位点（all sites）文件。

获取all-sites的vcf文件

gatk --java-options "-Xmx32g" HaplotypeCaller -R ref.fa -I input.bam -O output.g.vcf -ERC GVCF # 生成单个样本gvcf文件
gatk --java-options "-Xmx32g" GenomicsDBImport -V output1.g.vcf -V output2.g.vcf ... --genomicsdb-workspace-path <allsamples_genomicsdb>
for chr in $(cat chr.list);do gatk --java-options "-Xmx32g" GenotypeGVCFs -R ref.fa -V samples.g.vcf -L ${chr} -all-sites -O ${chr}_allsites.vcf.gz;done  # joint-calling这里一定要加参数-all-sites(把non-variant位点也保存到vcf文件中，后面才能进行pixy的分析）。必须使用-L指定染色体，否则有报错风险。耗时步骤。
ls *_allsites.vcf.gz >vcf.list # 把生成的染色体的vcf文件名保存到vcf.list
gatk MergeVcfs -I vcf.list -O samples_allsites.vcf.gz # 合并所有染色体vcf文件到一个vcf文件，耗时步骤

分别获取invariant和variant位点，过滤后，再合并

vcftools --gzvcf samples_allsites.vcf.gz --max-maf 0 --minQ 30 --max-missing 0.9 --min-meanDP 10 --max-meanDP 43 --remove filter.indv --recode --stdout | bgzip -c > samples_invariant_filtered.vcf.gz & # 提取invariant位点并过滤，同时生成out.log名称的log文件
vcftools --gzvcf samples_allsites.vcf.gz --mac 1 --min-alleles 2 --max-alleles 2 --minQ 30 --max-missing 0.9 --min-meanDP 10 --max-meanDP 43 --remove filter.indv --recode --stdout | bgzip -c > samples_variant_filtered.vcf.gz & # 提取variant位点并过滤，同时也生成out.log名称的log文件，最好与上一步骤分开在两个文件夹运行
tabix samples_invariant_filtered.vcf.gz
tabix samples_variant_filtered.vcf.gz
bcftools concat --allow-overlaps samples_variant.vcf.gz samples_invariant.vcf.gz -O z -o samples_allsites_filtered.vcf.gz # 合并invariant和variant位点到一个文件

1.3. pixy计算dxy，pi和fst

popfile.txt群体文件

popfile.txt文件包含样本个体分组信息，两列内容，tab分隔。第一列样本ID，第二列分组名称（常用物种作为分组依据）。

ERS223827   BFS
ERS223759   BFS
ERS223750   BFS
ERS223967   AFS
ERS223970   AFS
ERS223924   AFS
ERS224300   AFS
ERS224168   KES
ERS224314   KES

运行

1
2
3

bgzip samples_allsites_filtered.vcf
tabix samples_allsites_filtered.vcf.gz
nohup pixy --stats pi fst dxy --vcf sample_allsites_filtered.vcf.gz --populations popfile.txt --window_size 10000 --n_cores 16 &

1.4. pixy的结果

根据–stats的设定，生成三个统计量的结果文件。

pixy_pi.txt：居群内的核苷酸多样性（nucleotide diversity (pi)）

pop：popfile.txt群体文件中群体ID
chromosome：染色体或contig的ID
window_pos_1：基因组window的起始位置
window_pos_2：基因组window的终止位置
avg_pi：window的位点的平均核苷酸多样性。更具体地说，pixy 计算window中所有位点的每个位点的加权平均核苷酸多样性，其中权重由每个位点的基因分型样本数决定。
no_sites：window中至少具有一种有效基因型的位点总数。此统计信息被提供给用户，不会直接用于任何计算。
count_diffs：window中所有基因型之间成对差异的原始数量。这是计算 avg_pi 的分子。
count_comparisons：window中所有基因型之间非缺失成对比较的原始数量（即比较两种基因型且均有效的情况）。这是计算 avg_pi 的分母。
count_missing：window中所有基因型之间缺失配对比较的原始数量（即比较两种基因型且至少缺失一种的情况）。

pop	chromosome	window_pos_1	window_pos_2	avg_pi	no_sites	count_diffs	count_comparisons	count_missing
BFS	HiC_scaffold_1	1	10000	0.0	6	0	6	0
AFS	HiC_scaffold_1	1	10000	0.0	6	0	6	0
... ...

pixy_dxy.txt：居群间的核苷酸差异（nucleotide divergence (dxy)）

pop1：popfile.txt群体文件中群体ID，比较的第一个居群
pop2：popfile.txt群体文件中群体ID，比较的第二个居群
chromosome：染色体或contig的ID
window_pos_1：基因组window的起始位置
window_pos_2：基因组window的终止位置
avg_dxy：window的平均每个位点的核苷酸差异
no_sites：window中两个群体中至少具有一种有效基因型的位点总数。此统计信息提供给用户，不会直接用于任何计算。
count_diffs：所有基因型之间的成对交叉种群差异的原始数量。这是计算 avg_dxy 的分子。
count_comparisons：window中所有基因型之间非缺失成对交叉种群比较的原始数量（即比较两种基因型且均有效的情况）。这是计算 avg_dxy 的分母。
count_missing：window中所有基因型之间缺失的成对跨种群比较的原始数量（即两个基因型被比较而至少有一个缺失的情况）。这个统计数字提供给用户，不直接用于任何计算。

1
2
3

pop1	pop2	chromosome	window_pos_1	window_pos_2	avg_dxy	no_sites    count_diffs	count_comparisons	count_missing
BFS	AFS	HiC_scaffold_1	1	10000	0.0	6	24	0
CFS	DES	HiC_scaffold_1	1	10000	0.0	6	24	0

pixy_fst.txt：Weir and Cockerham’s FST

pop1：popfile.txt群体文件中群体ID，比较的第一个居群
pop2：popfile.txt群体文件中群体ID，比较的第二个居群
chromosome：染色体或contig的ID
window_pos_1：基因组window的起始位置
window_pos_2：基因组window的终止位置
avg_wc_fst：window的平均每个SNP（不是每个位点）的 Weir and Cockerham’s FST
no_snps：window内变异位点（SNPs）的总数

1
2
3

pop1	pop2	chromosome	window_pos_1	window_pos_2	avg_wc_fst	no_snps
BFS	AFS	HiC_scaffold_1	30001	40000	NA	67
CFS	DES	HiC_scaffold_1	30001	40000	0.0	67

结果的后续处理

计算平均值的注意事项

值得注意的是，如果想用pixy的结果继续计算更大的window或者整个基因组的平均值不能直接用window的pi/dxy值来计算。

正确的计算方法是使用原始计数重新计算differences/comparisons ratios，公式如下：

$$(window 1 count_diffs + window 2 count_diffs + … + window n count_diffs) / (window 1 comparisons + window 2 comparisons + … + window n comparisons)$$

计算整个基因组的pi值

cat pixy_pi.txt |cut -f1,2 |sort|uniq >pops.list
cat pops.list | while read line
do
    a=$(echo $line | awk '{print $1}');
    b=$(echo $line | awk '{print $2}');
    cat pixy_dxy.txt |awk '$1 == "'$a'" && $2 == "'$b'" {suma+=$8;sumb+=$9} END {print "'$a'", "'$b'", suma, sumb, suma/sumb}' >> dxy_average.txt
done

计算整个基因组的dxy平均值

cat pixy_dxy.txt |cut -f1 |sort|uniq >pop.list
cat pop.list | while read line
do
    a=$(echo $line | awk '{print $1}');
    cat pixy_pi.txt |awk '$1 == "'$a'" {suma+=$7;sumb+=$8} END {print "'$a'", suma, sumb, suma/sumb}' >> pi_average.txt
done