vcftools计算fst选择信号

Posted on Edited on In , Views: Waline: Word count in article: 3.1k Reading time ≈ 6 mins.

fst（F-statistics）是一种用于选择信号检测的统计方法，用于测量和比较不同群体之间的遗传差异。它可以帮助确定是否存在选择压力导致的遗传变异。通过计算不同群体之间的遗传差异（fst 值）并进行统计检验，fst 方法可以帮助检测到与自然选择相关的遗传变异。这种方法在种群遗传学、进化生物学和生态学研究中得到了广泛应用。

fst 的选择信号原理如下：

  1. 遗传差异的衡量：fst 通过计算群体间和群体内的遗传差异来衡量基因组水平上的遗传变异。它使用基因频率或基因型频率作为指标来评估群体间的遗传差异。
  2. 群体间遗传差异：fst 使用群体间的遗传差异来测量选择信号。群体间遗传差异是指在不同群体中，相同基因座上基因型或等位基因频率的差异。如果在某个基因座上，不同群体之间的遗传差异显著高于群体内的遗传差异，可能表明该基因座受到选择的影响。
  3. 统计计算：fst 的计算通常基于基因频率或基因型频率的方差分析。常用的 fst 计算方法包括 Weir and Cockerham 方法、Hudson’s unbiased estimator 和 AMOVA (Analysis of Molecular Variance) 方法。这些方法使用群体间和群体内的遗传差异来计算 fst 值。
  4. 统计显著性检验：计算得到的 fst 值可以与理论空模型进行比较来进行统计显著性检验。常用的检验方法包括置换检验、Bootstrap 方法和模拟方法。

在群体遗传学中衡量群体间的遗传分化的程度的指标有许多种，较为常见的就是遗传分化指数（Fst），fst是由F统计量演变而来，F统计量主要有三种（FIS，FIF，FST）。Fst是针对一对等位基因，如果基因座上存在复等位基因，则需要用Gst衡量，基因差异分化系数（gene differentiation coefficient，Gst）。假定有s个地方群体，第k个地方群体相对大小为wk，第k个地方群体中第i个等位基因频率为qk(i)，杂合体频率观察值为hk，那么整个群体中观察到的杂合体频率平均值HI，地方群体为理想群体的期望杂合体频率平均值HS，整个群体为理想群体的期望杂合体频率HT，分别为：

FIS，是HI相对于HS减少量的比值，即地方群体的平均近交系数。
FST，是HS相对于HT减少量的比值，即有亲缘关系地方群体间的平均近交系数。
FIT，是HI相对于HT减少量的比值，即整个群体的平均近交系数。
Fst值的取值范围是【0,1】，最大值为1表明两个群体完全分化，最小值为0表明群体间无分化。

在实际的研究中Fst值为0–0.05时说明群体间遗传分化很小，可以不做考虑；
0.05–0.15时，表明群体间存在中等程度的遗传分化；
0.15–0.25时群体间存在较大的遗传分化；
0.25以上的时候群体间就存在很大的遗传分化了。

目前主要使用的是vcftools来计算：

1
2
3
4
5
6
7
8
conda install -c bioconda vcftools
#单点计算
vcftools --vcf test.vcf --weir-fst-pop popa.txt --weir-fst-pop popb.txt --out a_b.fst
#窗口计算
vcftools --vcf test.vcf --weir-fst-pop popa.txt --weir-fst-pop popa.txt --out a_b.fst --fst-window-size 150000
#滑动窗口计算
vcftools --vcf test.vcf --weir-fst-pop popa.txt --weir-fst-pop popa.txt --out a_b.fst --fst-window-size 150000 --fst-window-step 75000

计算

1群体较少

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#! /bin/bash
#crate by lx
if [ $# -ne 6 ]; then
 echo "error.. need args"
 echo "command:$0 <VCF> <Pop1> <Pop2> <Win> <Step> <Out>"
 exit 1
fi
VCF=$1
Pop1=$2
Pop2=$3
Win=$4
Step=$5
Out=$6

if [[ "${file##*.}" = "vcf.gz" ]]; then
  vcftools --gzvcf ${VCF} --weir-fst-pop ${Pop1} --weir-fst-pop ${Pop2} --fst-window-size ${Win} --fst-window-step ${Step} --out ${Out} --max-missing 0.9 --maf 0.05
else
  vcftools --vcf ${VCF} --weir-fst-pop ${Pop1} --weir-fst-pop ${Pop2} --fst-window-size ${Win} --fst-window-step ${Step} --out ${Out} --max-missing 0.9 --maf 0.05
fi

2流水线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#!/bin/bash
#crate by lx
vcf_file="/home/liuxiao/1254/miss0.8.recode.vcf"
hu_file="/home/liuxiao/1254/group/group/txt/Hu.txt"
group_folder="/home/liuxiao/1254/group/group/txt/"
out_folder='/home/liuxiao/1254/result/fst_snp/'

txt_files=$(ls "$group_folder"/*.txt | grep -v "Hu.txt")
pids=()
for hu_group_file in $hu_file; do
  hu_group_name=$(basename "$hu_group_file" .txt)
  for txt_file in $txt_files; do
    group_name=$(basename "$txt_file" .txt)
    output_file="${hu_group_name}_${group_name}_8miss.fst"
    vcftools --vcf "$vcf_file" --weir-fst-pop "$hu_group_file" --weir-fst-pop "$txt_file" --out "$out_folder/$output_file"  &
    pids+=($!)
  done
done

for pid in "${pids[@]}"; do
  wait "$pid"
done

echo "所有计算完成"