第4章 VCF注释信息提取与清洗：处理缺失值、过滤低质量变异

终于到了真正动手处理数据的环节了。

VCF文件经过注释后，信息量暴增。但信息多不代表能用，里面掺杂着大量缺失值和低质量位点。如果不做清洗，后续的可视化就是“垃圾进，垃圾出”。

这一章，我会带你把这堆原始数据变成干净、可用、有意义的表格。

4.1 为什么必须做数据清洗？

VCF注释后的数据，说白了就是一张巨大的稀疏矩阵。

你想想看，一个变异位点在某个样本里可能没有覆盖到，或者覆盖深度不够；另一个位点可能在群体里就是纯合参考——这些都会被记录为缺失或低质量。

我刚开始做分析时犯过一个错：直接把所有位点丢进可视化工具。结果出来的图全是离群信号，查了半天才发现是未过滤的假阳性位点捣乱。

4.2 质量过滤：砍掉不可靠的变异

过滤的第一步，是干掉那些质量不过关的位点。

VCF文件里有两个核心质量指标：QUAL（位点质量值）和FILTER（过滤标记）。

我个人习惯先看FILTER列。如果标注了PASS，说明通过了GATK等软件的默认过滤；如果标注了LowQual或SbCluster，直接丢弃。

但仅靠软件的默认过滤不够。我们需要自己设定阈值：

看一个实际的过滤代码：

# 使用bcftools进行质量过滤
bcftools filter -i 'QUAL > 30 && FILTER="PASS"' \
    -O z -o filtered_high_qual.vcf.gz \
    input_annotated.vcf.gz

# 过滤后生成索引
tabix -p vcf filtered_high_qual.vcf.gz

4.3 缺失值处理：别让NaN毁了你的图

质量过滤完了，但数据里还有一堆缺失值。缺失值不等于没有变异，它意味着不可知——我们不能直接扔掉它。

常见的缺失值类型：

• ./.：基因型缺失，最常见
• 0/0：纯合参考（技术上不是缺失，但对可视化来说也是一种“无信息”状态）
• AD=0：等位基因深度为0，表示没有read支持该等位基因

我的处理原则很简单：把缺失值区分为“随机缺失”和“系统性缺失”。

# Python脚本：统计每个位点的缺失率
import pysam

vcf = pysam.VariantFile("filtered_high_qual.vcf.gz")
for record in vcf:
    dp = record.samples["SAMPLE01"].get("DP")
    gt = record.samples["SAMPLE01"].get("GT")
    
    if gt is None or dp is None:
        # 记录缺失信息
        print(f"{record.chrom}:{record.pos}\tGT缺失\tDP缺失")
    elif gt == (0, 0) and dp < 5:
        # 深度低的参考位点，标记为低置信
        print(f"{record.chrom}:{record.pos}\t低深度参考位点")
    else:
        # 有效位点
        pass

对于可视化来说，我建议根据不同的绘图任务做不同处理：

• 曼哈顿图：直接剔除缺失位点
• 热图：用均值/中位数填充，或标记为特殊颜色
• PCA/聚类：必须用完整矩阵，建议做LD pruning后再填充

4.4 信息提取：从VCF到可读表格

过滤和补缺之后，我们要把VCF变成能直接画图的格式。最常见的格式就是TSV或CSV。

需要的字段通常包括：

1. 坐标信息：染色体、位置
2. 等位基因：REF、ALT
3. 注释信息：基因名、功能影响、保守性分数
4. 群体频率：gnomAD、ExAC等
5. 功能预测：CADD、SIFT、PolyPhen

写个Python函数把关键字段提取出来：

def extract_vcf_to_df(vcf_path):
    import pandas as pd
    import pysam
    
    vcf = pysam.VariantFile(vcf_path)
    records = []
    
    for record in vcf:
        # 基本信息
        chrom = record.chrom
        pos = record.pos
        ref = record.ref
        alt = ','.join([str(a) for a in record.alts])
        qual = record.qual
        
        # 注释信息（假设已用SnpEff注释）
        annotation = record.info.get('ANN', [''])[0]
        gene = annotation.split('|')[3] if annotation else 'NA'
        impact = annotation.split('|')[2] if annotation else 'NA'
        
        # 频率信息
        af_gnomad = record.info.get('AF_gnomAD', 'NA')
        if af_gnomad != 'NA':
            af_gnomad = float(af_gnomad)
        
        records.append({
            'CHROM': chrom,
            'POS': pos,
            'REF': ref,
            'ALT': alt,
            'QUAL': qual,
            'GENE': gene,
            'IMPACT': impact,
            'AF_gnomAD': af_gnomad
        })
    
    return pd.DataFrame(records)

# 实际使用
df = extract_vcf_to_df("filtered_high_qual.vcf.gz")
df.to_csv("cleaned_variants.tsv", sep='\t', index=False)

4.5 核心逻辑：清洗工作流

前面的步骤不算少。我把它们串成一张图，方便你记住整个流程。

这张图概括了从原始VCF到可绘图数据的完整流程。个人觉得最容易被忽略的是中间那个分支——缺失值的分类处理。绝大部分教程只教你怎么过滤，不教你怎么区分缺失类型。这里你得多留个心眼。

4.6 实战环节：跑一个完整的清洗管线

纸上谈兵够了，咱们实际跑一遍。

#!/bin/bash
# 完整清洗管线

# 步骤1：基础质量过滤
bcftools filter \
    -i 'QUAL > 30 && FILTER="PASS"' \
    -O z -o step1_qual.vcf.gz \
    input_annotated.vcf.gz
tabix -p vcf step1_qual.vcf.gz

# 步骤2：按样本深度过滤（只保留深度≥10的基因型）
bcftools view \
    -i 'FMT/DP >= 10' \
    -O z -o step2_dp.vcf.gz \
    step1_qual.vcf.gz
tabix -p vcf step2_dp.vcf.gz

# 步骤3：过滤缺失率超过20%的位点
bcftools filter \
    -i 'F_MISSING < 0.2' \
    -O z -o step3_missing.vcf.gz \
    step2_dp.vcf.gz
tabix -p vcf step3_missing.vcf.gz

# 步骤4：提取关键信息到表格
python3 extract_vcf.py step3_missing.vcf.gz cleaned_variants.tsv

echo "清洗完成！输出文件：cleaned_variants.tsv"

嗯，到这里，你已经拿到了干净、结构化的数据。可以安心进入下一阶段的可视化了。