四倍体VCF样本缺失率计算求助：两种缺失类型未被正确统计

一、自定义脚本优化（Python + pysam） #

针对你已能处理位点缺失率的情况，优化脚本即可实现样本缺失率统计：

• 依赖pysam库（安装命令：pip install pysam）
• 核心逻辑：遍历每个变异位点，对每个样本单独检查GT字段是否为四倍体缺失格式（./././.或.:.:.:.），统计每个样本的总位点数和缺失位点数
• 示例代码片段：

import pysam

vcf_path = "your_input.vcf"
sample_miss = {}
sample_total = {}

# 初始化样本统计字典
with pysam.VariantFile(vcf_path) as vcf:
    for sample in vcf.header.samples:
        sample_miss[sample] = 0
        sample_total[sample] = 0

    for record in vcf:
        # 每个位点给所有样本加总计数
        for sample in vcf.header.samples:
            sample_total[sample] += 1
        # 检查每个样本的GT是否缺失
        for sample in vcf.header.samples:
            gt = record.samples[sample]['GT']
            # 四倍体缺失GT的两种情况：全None或全'.'（pysam解析后会转成None）
            if all(allele is None for allele in gt):
                sample_miss[sample] += 1

# 计算缺失率并输出
for sample in sample_miss:
    miss_rate = sample_miss[sample] / sample_total[sample]
    print(f"{sample}\t{miss_rate:.4f}")

• 后续可根据缺失率阈值，用pysam过滤掉不符合要求的样本。

二、利用PLINK 2.0处理多倍体缺失 #

PLINK 2.0原生支持多倍体数据，指定倍性后可准确统计样本缺失率：

1. 将VCF转成PLINK 2.0格式：

plink2 --vcf your_input.vcf --make-pgen --ploidy 4 --out temp_plink

2. 统计样本缺失率：

plink2 --pfile temp_plink --missing --out sample_missing_stats

3. 根据输出的sample_missing_stats.imiss文件中的F_MISS列，过滤缺失率过高的样本：

plink2 --pfile temp_plink --remove sample_missing_stats.imiss --remove-if F_MISS > 0.1 --make-pgen --out filtered_plink

（将0.1替换为你的实际阈值）

三、bcftools自定义统计与格式统一 #

1. 直接用bcftools + awk统计样本缺失率 #

通过bcftools query提取GT字段，再用awk统计两种缺失格式：

bcftools query -f '%CHROM\t%POS[\t%GT]\n' your_input.vcf | \
awk '
BEGIN {
    # 提前用bcftools query -l your_input.vcf获取样本名，保存到samples.txt（每行一个）
    while ((getline line < "samples.txt") > 0) samples[NR] = line
}
{
    for(i=3; i<=NF; i++){
        total[i-2]++
        # 匹配四倍体的两种缺失GT格式
        if ($i ~ /^\.\/\.\/\.\/\./ || $i ~ /^\.:\.:\.:\./) miss[i-2]++
    }
}
END {
    for(j=1; j<=length(samples); j++){
        miss_rate = miss[j] / total[j]
        print samples[j], miss_rate
    }
}' > sample_miss_rates.txt

2. 统一缺失格式后用标准工具统计 #

先将./././.格式的缺失GT替换为.:.:.:.，再用vcftools统计：

# 替换GT格式
bcftools +setGT your_input.vcf -Ov -o unified_miss.vcf -- -t . --mode replace --force -i 'GT="./././."'
# 统计样本缺失率
vcftools --vcf unified_miss.vcf --missing --out sample_missing

输出的sample_missing.imiss文件即为样本缺失统计结果，可据此过滤样本。

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