如何从VCF文件中按读取顺序统计各群体的个体数量？

嘿，我来帮你搞定这个需求！要按群体在VCF样本中首次出现的顺序统计每个群体的个体数量，核心是先把VCF里的样本顺序和它们的群体对应起来，再按顺序计数。下面是两种实用的方法：

前提准备：群体-个体映射文件 #

VCF的头部（你提供的内容）里没有存储群体信息，所以你需要一个群体映射文件（比如命名为sample_pop_map.txt），每行格式是[个体ID] [群体名称]，而且个体的顺序必须和VCF文件里样本列的顺序完全一致。

如果还没有这个文件，先从VCF里提取样本列表：

# 用bcftools提取VCF的样本ID列表，顺序和VCF里一致
bcftools query -l your_data.vcf > samples.txt

然后给每个样本标注所属群体，保存成上面说的映射文件即可。

方法1：用awk快速统计（命令行工具） #

这个方法适合快速处理，而且能严格保持群体的首次出现顺序：

1. 创建一个awk脚本文件count_pops.awk，内容如下：

BEGIN {
    # 用来记录群体首次出现的顺序
    pop_index = 0
}
{
    # 如果这个群体还没被记录过，加入顺序列表
    if (!($2 in pop_order)) {
        pop_order[pop_index++] = $2
    }
    # 统计每个群体的个体数
    pop_count[$2]++
}
END {
    # 按群体首次出现的顺序输出结果
    for (i = 0; i < pop_index; i++) {
        print pop_order[i], pop_count[pop_order[i]]
    }
}

2. 运行脚本：

awk -f count_pops.awk sample_pop_map.txt

输出结果的群体顺序就是它们在VCF样本中首次出现的顺序，完全符合你的要求。

方法2：用Python脚本（灵活定制） #

如果需要更复杂的处理（比如直接从VCF读取样本并关联群体），Python脚本会更灵活：

def get_vcf_samples(vcf_path):
    """从VCF文件中提取样本ID列表，顺序和VCF一致"""
    samples = []
    with open(vcf_path, 'r') as vcf_file:
        for line in vcf_file:
            if line.startswith('#CHROM'):
                # VCF从第10列开始是样本ID
                samples = line.strip().split('\t')[9:]
                break
    return samples

def get_population_list(map_path):
    """读取群体映射文件，返回按样本顺序排列的群体列表"""
    pop_list = []
    with open(map_path, 'r') as map_file:
        for line in map_file:
            if line.strip():
                _, pop = line.strip().split()
                pop_list.append(pop)
    return pop_list

def count_populations_by_order(pop_list):
    """按群体首次出现的顺序统计个体数量"""
    pop_counts = {}
    pop_order = []
    for pop in pop_list:
        if pop not in pop_order:
            pop_order.append(pop)
        pop_counts[pop] = pop_counts.get(pop, 0) + 1
    
    # 按顺序输出结果
    print("群体\t个体数量")
    for pop in pop_order:
        print(f"{pop}\t{pop_counts[pop]}")

if __name__ == '__main__':
    # 替换成你的文件路径
    VCF_FILE = "your_data.vcf"
    POP_MAP_FILE = "sample_pop_map.txt"

    # 如果已经确认映射文件顺序正确，可以跳过样本提取步骤
    # samples = get_vcf_samples(VCF_FILE)
    pop_list = get_population_list(POP_MAP_FILE)
    count_populations_by_order(pop_list)

运行脚本：

python count_vcf_pops.py

这个脚本会清晰输出每个群体的个体数量，顺序严格遵循群体在VCF样本中首次出现的先后。

关键注意点 #

• 群体映射文件的个体顺序必须和VCF样本列顺序完全一致，否则统计的顺序会出错。
• 如果你的VCF是压缩格式（.vcf.gz），可以用bcftools query -l your_data.vcf.gz提取样本列表，无需解压。

内容的提问来源于stack exchange，提问作者Ella Bowles