2. 多批次VCF数据准备：收集、版本检查、压缩与索引
好，我们正式开始讲故事。
做多批次VCF一致性分析，最怕什么？最怕数据还没看，工具就开始报错。我见过太多同学，兴致勃勃跑流程，结果第一步就被卡住——VCF格式不对，版本不兼容，甚至文件都没压缩。嗯，这些坑我都踩过，而且不止一次。
这一章，我们就来聊聊怎么把原始数据整理成“可以上战场”的状态。说白了，就是数据准备四步走：收集、版本检查、压缩、索引。每一步我都给你掰开揉碎讲清楚。
2.1 收集数据，别做无头苍蝇
很多人觉得收集数据嘛，扔到一个文件夹就行了。我建议你不要这样。我个人的习惯是，在项目启动时就定好一套命名规则。比如：
# 我常用的命名格式
project_batch1_20250301.vcf
project_batch2_20250315.vcf
project_batch3_20250328.vcf为什么要这么干？你想想看，三个月后你可能手头有30个批次，没有命名规则，光靠文件名猜内容，能把你搞疯。我在一个项目中就看到过同事用“final_v1.vcf”、“final_v2.vcf”、“real_final.vcf”这种命名，结果根本分不清哪个是哪个版本。
| 批次编号 | 文件名 | 生成日期 | 样本数 | 备注 |
|---|---|---|---|---|
| Batch 1 | proj_b1_20250301.vcf | 2025-03-01 | 48 | 基准批次 |
| Batch 2 | proj_b2_20250315.vcf | 2025-03-15 | 52 | 新增样本 |
| Batch 3 | proj_b3_20250328.vcf | 2025-03-28 | 50 | 重测数据 |
2.2 版本检查，别掉进格式陷阱
VCF格式有v4.1、v4.2、v4.3三个主要版本。为什么非要检查这个？因为不同版本的关键字段定义不一样啊！我曾经遇到过，两个VCF文件看起来数据都对，结果一合并，QUAL字段全乱了——因为一个是v4.1，一个是v4.2，QUAL的精度定义不同。
怎么检查版本？用grep看文件头就行：
# 查看VCF版本信息
head -5 sample.vcf | grep "##fileformat"
# 输出示例
##fileformat=VCFv4.2
##fileDate=20250301
##source=GATK4.3.0.0
##reference=GRCh38
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE1- fileformat — 格式版本（必须一致）
- reference — 参考基因组（hg19 vs hg38混用会出大事）
- 样本ID — 列名重复？不要有同名样本
- FORMAT字段 — 定义是否一致（GT、DP、GQ等）
我个人建议写一个脚本自动检查所有待分析VCF的文件头。拿Python写也不复杂：
import gzip
def check_vcf_header(filepath):
with gzip.open(filepath, 'rt') if filepath.endswith('.gz') else open(filepath, 'r') as f:
for line in f:
if line.startswith('##fileformat'):
print(f"版本: {line.strip()}")
if line.startswith('##reference'):
print(f"参考基因组: {line.strip()}")
if line.startswith('#CHROM'):
samples = line.strip().split('\t')[9:]
print(f"样本ID: {samples}")
break2.3 压缩必须是bgzip，别用错
好，版本确认没问题了，接下来就是压缩。很多同学觉得，压缩嘛，gzip不就行了？嗯，还真不行。
为什么？因为我们要做索引。bgzip和gzip的区别在于，bgzip可以按块压缩，支持随机访问。你在tabix索引的时候，需要能在不拆开整个文件的情况下快速定位到某个区域。gzip做不到这一点，所以VCF工具链（bcftools、GATK）都要求用bgzip。
# 正确的做法 —— 用bgzip
bgzip -c input.vcf > output.vcf.gz
# 或者直接原地压缩（推荐）
bgzip input.vcf
# 会生成 input.vcf.gzgunzip -t output.vcf.gz 就能测完整性。
还有一个坑：有些同学用Python的 gzip.open() 来写压缩文件，这个方法默认会生成gzip格式，不是bgzip格式。虽然文件名都是 .gz，但是BCFtools是不认的！
所以，老老实实用 bgzip 命令行工具来完成压缩。如果你非要程序里处理，可以用 pysam.bgzip：
import pysam
# pysam自带的bgzip压缩
pysam.tabix_compress("input.vcf", "output.vcf.gz")2.4 索引，让工具能“找到”数据
压缩之后必须做索引。索引的作用是什么？说白了，就是让工具能快速定位到基因组上的某个位置，而不需要把整个VCF都读到内存里。你需要同时查第1号染色体和第22号染色体？没有索引的话，工具就得从头扫两遍。有了索引，几毫秒就定位到了。
索引使用 tabix 工具：
# 对bgzip压缩后的VCF做索引
tabix -p vcf input.vcf.gz
# 完成后会生成 input.vcf.gz.tbi 文件
# 查看索引内容
tabix -l input.vcf.gz-p vcf— 指定输入格式是VCF（也可以是bcf、gff、bed等）- 输出文件
.tbi必须和.vcf.gz放在同一目录 - 如果VCF文件非常大（比如超过500M），可以考虑用
--csi参数生成CSI索引 - 批量索引时，使用
for f in *.vcf.gz; do tabix -p vcf $f; done
我个人的经验：索引之后，一定要验证能不能正常读取。随便挑一个位置试一下：
# 验证索引有效
tabix input.vcf.gz chr1:1000000-2000000 | head -5如果输出正常，说明索引生效了。如果报错说“index not found”或者“file is not bgzip compressed”，那就是前两步错了，回头检查。
2.5 统一数据管理的一站式方案
到现在为止，每个VCF都已经经过了收集、版本检查、压缩、索引四道工序。项目进行到这里，我建议你做一个数据管理脚本，把上述步骤全串起来。这是我当年做大项目时写的一个简化版脚本：
#!/bin/bash
# 多批次VCF数据准备一键脚本（简化版）
INPUT_DIR="./raw_vcf"
OUTPUT_DIR="./prepared_vcf"
mkdir -p $OUTPUT_DIR
for vcf in $INPUT_DIR/*.vcf; do
echo "=== 处理文件: $vcf ==="
# 1. 版本检查
head -1 $vcf | grep -q "^##fileformat=VCFv4.2"
if [ $? -ne 0 ]; then
echo "警告: $vcf 不是VCFv4.2格式"
continue
fi
# 2. 检查参照基因组
head -5 $vcf | grep -q "GRCh38"
if [ $? -ne 0 ]; then
echo "警告: $vcf 不是GRCh38参考基因组"
continue
fi
# 3. 检查样本ID
SAMPLE_COUNT=$(head -1 $vcf | awk -F'\t' '{print NF-9}')
echo "样本数: $SAMPLE_COUNT"
# 4. bgzip压缩
bgzip -c $vcf > $OUTPUT_DIR/$(basename $vcf).gz
# 5. 创建索引
tabix -p vcf $OUTPUT_DIR/$(basename $vcf).gz
# 6. 验证
tabix -l $OUTPUT_DIR/$(basename $vcf).gz | head -3
echo "=== 完成: $vcf ==="
doneconda install -c bioconda bcftools tabix。这样既能隔离环境，又能保证版本一致性。
2.6 一个真实的踩坑记录
最后分享一个我经历过的真实案例。有一次，我从合作方拿到了5个批次的VCF文件，对方保证说“格式都处理好啦”。我没多想，直接做了一致性分析。结果跑了两天才发现，第二个批次的样本ID里有空格，第四个批次的QUAL精度不对，第五个批次的索引文件是旧的。
从那次以后，我发誓：数据到手后，第一步必做版本检查，第二步必重新压缩和索引，绝不做“拿来主义”。
你想想看，一个一致性分析流程，如果数据本身都有问题，结果怎么可能可靠？所以，数据准备这一关，值得你花200%的精力去过。