VCF变异注释图表构建与生物学结果分析 - 第1节

第1章 VCF格式基础：先搞懂文件里到底写了什么

各位做生信的朋友，咱们在做变异检测的时候，VCF文件几乎是天天打交道的。但说实话，我见过不少人直接用工具处理VCF，却连里面字段是什么含义都说不清。这样搞可视化，很容易画出一堆没有意义的图表。

VCF，全称是Variant Call Format，说白了就是记录基因组上变异信息的一种标准化格式。它由VCF规范定义，现在主流版本是VCF 4.x。我个人习惯把它看作一张表格，只不过是带了很多“元信息”的特殊表格。

这一章，咱们就把VCF的“底裤”扒干净，搞清楚每一行、每一列到底在记录什么。

VCF文件的组成结构

一个标准的VCF文件，说白了就是三部分：卷首的元数据、一条表头、以及成千上万条数据行。我刚开始接触时总觉得元数据行是废话，后来才发现，很多注释工具的关键信息就藏在这里面。

1. 元数据行（以 ## 开头）

这部分是对整个文件的描述。包括文件格式版本、参考基因组版本、以及后续数据行会用到的所有字段定义。你想想看，没有这些定义，光看INFO那列的一大串字母数字，谁知道是什么意思？

我在项目中遇到过最坑的事：不同工具生成的VCF，INFO字段的标签可能同名但含义完全不同。所以，做可视化之前一定要先读一下元数据行，确认每个标签的真实含义。

常见的元数据行示例：

##fileformat=VCFv4.2
##fileDate=20231015
##source=Platypus
##reference=file:///ref/hg38.fa
##contig=<ID=chr1,length=248956422>
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic Depth">

2. 表头行（以 #CHROM 开头）

这行是数据列的“索引”，固定有 8 个必填列，后面可以跟任意数量的样本列。嗯，这里要注意，不同的VCF版本可能会略有区别，但万变不离其宗。

表头的列名如下：

3. 数据行的核心逻辑

搞懂数据行之前，我们先看一条真实数据长什么样。我记得自己刚入行那会儿，看到这一长串用制表符分割的数字，完全懵了。

chr1    12345   .       A       G       58.2    PASS   DP=234;AF=0.33;AN=2;AC=1   GT:AD:DP:GQ   0/1:120,114:234:99

注意，VCF每一列之间是制表符分隔的，我见过有人用逗号或空格处理，结果解析全错了。你写脚本处理VCF时，一定要按制表符分割。

这里我们拆解一下这条数据：

• CHROM (chr1)：变异发生在1号染色体。
• POS (12345)：位置是第12345号碱基。
• ID (.)：点号表示这个变异在dbSNP里没有对应的rsID。
• REF (A)：参考基因组在这个位置上是A。
• ALT (G)：检测到的变异是G碱基，这是一个SNP。
• QUAL (58.2)：质量值，数值越高越可靠。说实话，58.2属于中等偏上水平。
• FILTER (PASS)：通过了所有的过滤条件。
• INFO：
  • DP=234：这个位点的总覆盖深度是234条reads。
  • AF=0.33：变异等位基因频率是33%。
  • AN=2：这个位点的有效等位基因总数是2（通常是二倍体样本的两个拷贝）。
  • AC=1：携带变异等位基因的等位基因数是1。
• FORMAT (GT:AD:DP:GQ)：定义了样本列中各字段的顺序。
• 样本列 (0/1:120,114:234:99)：
  • GT=0/1：杂合变异（一个参考等位基因，一个变异等位基因）。
  • AD=120,114：分别支持参考和变异等位基因的reads数。
  • DP=234：样本在该位点的深度。
  • GQ=99：基因型质量值，接近满分。

关键字段解析：INFO 的“秘密”

INFO 列是整个VCF注释的核心战场。我过去做可视化时，80%的图表数据都来自这里。不同的注释工具（VEP、SnpEff、ANNOVAR等）会在INFO里塞进各种标签。

常见的 INFO 标签如下：

常见注释信息的读取方式

在实际项目中，VCF中的注释信息远不止上面那些。像VEP产生的CSQ、SnpEff产生的ANN，都是非常冗长的字符串。我个人的做法是：先把这些字段拆解成表格，再做分析。

以 VEP 的 CSQ 为例，它的格式是：

##INFO=<ID=CSQ,Number=.,Type=String,Description=
"Consequence annotations from VEP. Format: 
Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|
BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|
Protein_position|Amino_acids|Codons|Existing_variation|
DISTANCE|STRAND|FLAGS|VARIANT_CLASS|SYMBOL_SOURCE|
HGNC_ID|...">

你看，这个字段里用竖线分隔了超多信息。如果想提取“这个变异影响了哪个基因”，你就需要按竖线分割，然后找到 SYMBOL 对应的位置（VEP里通常是第4个字段）。

我在处理这类数据时，会用 Python 写一个简单的解析函数：

def parse_csq(csq_string):
    fields = [
        "Allele", "Consequence", "IMPACT", "SYMBOL",
        "Gene", "Feature_type", "Feature", "BIOTYPE",
        "EXON", "INTRON", "HGVSc", "HGVSp"
    ]
    values = csq_string.split("|")
    # 只取需要的字段
    result = {}
    for i, key in enumerate(fields):
        if i < len(values):
            result[key] = values[i]
    return result

这个函数可以把一大串CSQ解析成字典，然后用 pd.DataFrame 收集起来，方便后续画图。嗯，在实际应用中，我还会加上缺失值的判断，有些VCF的CSQ字段可能只填写了前几个值。

grep "^##INFO=