ln-πratio对数频率比率

Posted on Edited on Views: Waline: Word count in article: 3.9k Reading time ≈ 7 mins.

群体遗传选择信号分析之ln-πratio

在群体基因组学研究中，ln-piratio 是一个常用的统计量，用于评估群体间的遗传分化程度。ln-piratio 统计量结合了两个关键的概念：π（pi）和θ（theta）。其中 π 表示群体内的平均遗传多样性，而 θ 则是从群体内随机抽取两个个体并计算其遗传差异的一种估计。ln-piratio 通常用于比较不同群体或不同地理区域之间的遗传分化。

ln-piratio 的计算

ln-piratio 通常定义为群体间 π 与群体内 π 的比值的自然对数。具体来说，对于两个群体A和B，ln-piratio 可以这样计算：

1
\text{ln-piratio} = \ln\left(\frac{\pi_{AB}}{\pi_A}\right)
  • 正值：如果 ln-piratio 的值为正，表示群体间遗传多样性高于群体内遗传多样性，表明群体间存在遗传分化。
  • 负值：如果 ln-piratio 的值为负，表示群体间遗传多样性低于群体内遗传多样性，表明群体间遗传分化较小。

start.sh用于分析两群体间单点ln-ratio

1
bash start.sh --vcf /home/liuxiao/1254/filter/1254_chrall.vcf.recode.vcf.gz --popdir /home/liuxiao/1254/group1/group_id/qun --hu Hu.txt --outdir /home/liuxiao/1254/group1/ln-ratio

start.sh对应

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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#!/bin/bash
function usage() {
    echo "Usage: bash $0 --vcf <vcf> --popdir <popdir> --hu <hu> --outdir <outdir>"
    echo "required options"
    echo "-v|--vcf       vcf file"
    echo "-d|--popdir    群体txt文件目录"
    echo "-h|--hu        Hu群体txt文件"
    echo "-o|--outdir    输出目录"
    exit 1;
}

vcf=""
popdir=""
hu=""
outdir=""

while [[ $# -gt 0 ]]
do
  case "$1" in
    -v|--vcf )
        vcf=$2 ; shift 2 ;;
    -d|--popdir )
        popdir=$2 ; shift 2 ;;
    -h|--hu )
        hu=$2 ; shift 2 ;;
    -o|--outdir )
        outdir=$2 ; shift 2 ;;
    *) echo "Option error!";
       usage    
       ;;
  esac
done

if [ -z $vcf ] || [ -z $popdir ] || [ -z $hu ] || [ -z $outdir ]; then
    echo "Option --vcf and --popdir and --hu and --outdir not specified" >&2
    usage
fi

function process_population() {
    local hu_basename=$1
    local pop_file=$2
    local outdir=$3
    local vcf=$4

    local pop_basename=$(basename $pop_file .txt)

    if [ "$pop_basename" != "$hu_basename" ]; then

        vcftools --gzvcf $vcf --freq --keep $pop_file --out ${outdir}/${pop_basename}

        python ln_ratio_per_site.py --group1 ${outdir}/${hu_basename}.frq --group2 ${outdir}/${pop_basename}.frq --outprefix ${outdir}/${hu_basename}_vs_${pop_basename}
    fi
}

function main() {
    hu_basename=$(basename $hu .txt)
    mkdir -p $outdir

    vcftools --gzvcf $vcf --freq --keep ${popdir}/${hu} --out ${outdir}/${hu_basename}

    for pop_file in ${popdir}/*.txt; do
        process_population $hu_basename $pop_file $outdir $vcf &
    done

    wait
}

main

对应的ln_ratio_per_site.py

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
30
31
32
33
34
35
36
37
38
39
40
import pandas as pd
import numpy as np
import argparse

def extract_maf(row):
    # 提取等位基因频率并找到次常见等位基因频率
    freqs = [float(af.split(':')[1]) for af in row['{ALLELE:FREQ}'].split()]
    maf = min(freqs) 
    return maf

def calculate_ln_ratio(group1_file, group2_file, outprefix):
    # 读取数据文件
    group1 = pd.read_csv(group1_file, delim_whitespace=True, comment='#')
    group2 = pd.read_csv(group2_file, delim_whitespace=True, comment='#')
    
    # 提取 MAF 列
    group1['MAF_1'] = group1.apply(extract_maf, axis=1)
    group2['MAF_2'] = group2.apply(extract_maf, axis=1)
    
    # 合并数据文件
    merged = pd.merge(group1[['CHROM', 'POS', 'MAF_1']], group2[['CHROM', 'POS', 'MAF_2']], on=['CHROM', 'POS'])
    
    # 计算 ln_ratio
    merged['ln_ratio'] = np.log(merged['MAF_1'] / merged['MAF_2'])
    
    # 保存结果
    merged.to_csv(f'{outprefix}_ln_ratio.txt', sep='\t', index=False)

def main():
    parser = argparse.ArgumentParser(description='Calculate ln ratio of Pi values per site.')
    parser.add_argument('--group1', type=str, required=True, help='Group 1 Pi values file')
    parser.add_argument('--group2', type=str, required=True, help='Group 2 Pi values file')
    parser.add_argument('--outprefix', type=str, required=True, help='Output prefix for the results')

    args = parser.parse_args()
    calculate_ln_ratio(args.group1, args.group2, args.outprefix)

if __name__ == "__main__":
    main()