第3章 用Pandas读取VCF：INFO与FORMAT字段解析实战
说实话，VCF文件解析是我做生物信息时最常碰到的活儿。刚开始我也是用循环逐行读文件，后来发现——真没必要。Pandas一出手，这些表格化的数据就能收拾得服服帖帖。
这一章，我就带你把VCF里的INFO和FORMAT字段拆解成DataFrame。你可能会问：直接用文本工具不香吗？当然可以，但一旦数据量上去，批处理和统计就麻烦了。你看完后面的内容就明白了。
3.1 VCF文件结构快速回顾
在动手前，咱们得先搞清楚VCF长啥样。典型的VCF长这样：
##fileformat=VCFv4.2
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE1 SAMPLE2
chr1 100 . A G 50.0 PASS DP=100;AF=0.5 GT:AD:DP 0/1:50,50:100 1/1:0,100:100关键部分我标一下：
- INFO：用分号分隔的键值对，记录变异级别的信息（比如测序深度、等位基因频率）
- FORMAT：冒号分隔的标签列表，定义样品级别的字段
- 样品列：每个样品的数据，用冒号分隔，顺序和FORMAT一一对应
我最初做全外显子组数据时，看到几十兆的VCF文件——嗯，人工解析基本没戏。
3.2 为什么用Pandas而不是自己写循环？
你可能手写过这样的解析代码：
for line in open('input.vcf'):
if line.startswith('#'):
continue
fields = line.strip().split('\t')
info = dict(item.split('=') for item in fields[7].split(';'))
# 然后用info['DP']、info['AF']去做统计……
这代码本身没问题，但可维护性差。数据量大、字段多时，索引容易乱。我习惯用Pandas，主要是因为它：
- 处理大数据集时不卡顿（用底层C数组）
- 能直接做分组、聚合、过滤——这些在生信里太常用
- 导出为Excel、CSV一键搞定
3.3 读取VCF核心列
先跳过注释行。Pandas读不了“#”开头的行，所以得指定跳过注释行：
import pandas as pd
# 定义列名（VCF标准列）
col_names = ['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT']
# 读取VCF，跳过注释行
df = pd.read_csv('input.vcf', sep='\t', comment='#', header=None, names=col_names)
print(df.head())
输出结果里，INFO和FORMAT还是字符串。别着急，下一步就是拆解它们。
3.4 解析INFO字段
INFO字段由分号分隔的键值对构成。比如：
DP=100;AF=0.5;AN=200;AC=100我们把它转为字典：
def parse_info(info_str):
if pd.isna(info_str):
return {}
items = info_str.split(';')
result = {}
for item in items:
if '=' in item:
key, value = item.split('=', 1)
# 尝试转数字
try:
value = float(value) if '.' in value else int(value)
except:
pass
result[key] = value
else:
# 有的是标志位，没有值（比如END=300）
result[item] = True
return result
# 应用函数
df_info = df['INFO'].apply(parse_info).apply(pd.Series)
df = pd.concat([df.drop(columns=['INFO']), df_info], axis=1)
为什么我用split('=', 1)？我之前遇到过一个VCF，INFO里直接写成END=200;SOMETHING_NAME=value，只用split('=')会得到三个元素，然后索引就乱了。按第一个等号分最安全。
解析后的DataFrame变成这样：
| CHROM | POS | DP | AF | AN | AC |
|---|---|---|---|---|---|
| chr1 | 100 | 100 | 0.5 | 200 | 100 |
| chr1 | 200 | 80 | 0.3 | 160 | 48 |
df['DP'].mean()算平均深度了。
3.5 解析FORMAT字段
FORMAT字段稍微复杂点——它描述样品级别的数据。比如：
FORMAT: GT:AD:DP
SAMPLE1: 0/1:50,50:100解析策略：对每个样品列，按冒号分割，再映射到FORMAT里定义的标签。
# 假设样品列从第10列开始（0-based索引）
# 先得到样品名
sample_cols = df.columns[9:].tolist() # 或从原VCF读取的列
def parse_format(format_str, sample_str):
tags = format_str.split(':')
values = sample_str.split(':')
if len(tags) != len(values):
return {}
result = {}
for tag, val in zip(tags, values):
try:
val = float(val) if '.' in val else int(val)
except:
pass
result[tag] = val
return result
for sample in sample_cols:
# 对每个样品，构建其解析后的列
parsed = df.apply(lambda row: parse_format(row['FORMAT'], row[sample]), axis=1)
parsed_df = parsed.apply(pd.Series)
# 添加前缀避免列名冲突
parsed_df.columns = [f"{sample}_{col}" for col in parsed_df.columns]
df = pd.concat([df, parsed_df], axis=1)
解析后，你会看到类似这样的列：
sample1_GT sample1_AD sample1_DP sample2_GT sample2_AD sample2_DP
0/1 50,50 100 1/1 0,100 100
PL标签，它是一个逗号分隔的整数列表——直接用float转换会报错，得单独处理。
3.6 整体流程展示
我把整个流程整合成一张图，你一看就明白：
3.7 实战避坑指南
这里我说几个我亲自踩过的坑：
.（点），Pandas会把它当作字符串。float('.')会报错。解决方案是在try/except里加上if val == '.': val = None。
pd.read_csv(..., chunksize=1000)，然后循环解析每个块，最后合并。我处理过2GB的VCF，不分块直接解析，内存飙到8GB——差点把服务器跑崩。
3.8 小结
通过这一章，你应该能熟练使用Pandas读取VCF，并把INFO和FORMAT字段拆成整齐的列。接下来你就可以对这些数据做各种统计，比如算一下每个位点的等位基因频率分布，或者筛选出深度大于30的位点。
我个人的心得是：解析这一步虽然看起来繁琐，但做好了一劳永逸。把解析函数封装好，以后不管多少VCF文件，一行代码就能搞定。这比每次手动改正则表达式靠谱多了。