VCF文件处理：bcftools过滤/注释/统计

1
# 提取 CHROM, POS, REF, ALT, QUAL
2
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\n' variants.vcf.gz | head
3

4
# 提取样本的 GT 和 DP
5
bcftools query -f '%CHROM\t%POS\t[%GT\t%DP]\n' variants.vcf.gz | head
6
# [] = 对每个样本重复，每个样本输出 GT 和 DP
7

8
# 单个样本的 GT
9
bcftools query -f '%CHROM\t%POS\t[%SAMPLE %GT\n]' variants.vcf.gz | head
10

11
# 输出到 TSV 做下游分析
12
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%INFO/DP\t%INFO/AF\n' \
13
    variants.vcf.gz > variants.tsv

1
# 生成统计报告
2
bcftools stats variants.vcf.gz > stats.txt
3

4
# 看几个关键指标
5
cat stats.txt | grep -E "^SN"

1
plot-vcfstats -p vcf_stats_plot/ stats.txt
2
# 生成一套 PNG 图：QUAL 分布、DP 分布、Indel 长度分布等

1
# 两个 VCF 的交集（两个都有的变异）
2
bcftools isec -n=2 -p isec_out sample1.vcf.gz sample2.vcf.gz
3

4
# 差集（只在 sample1 中的变异）
5
bcftools isec -C sample1.vcf.gz sample2.vcf.gz -O z -o sample1_only.vcf.gz
6

7
# 多文件
8
bcftools isec -n+3 -p isec_3way s1.vcf.gz s2.vcf.gz s3.vcf.gz
9
# -n+3 = 至少在3个文件中出现的变异

1
# 纵向合并（同样的位点，不同样本）—— concat
2
bcftools concat -O z -o merged.vcf.gz chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz
3

4
# 横向合并（不同位点或不同样本）—— merge
5
bcftools merge -O z -o combined.vcf.gz sample1.vcf.gz sample2.vcf.gz
6

7
# 注意：concat 要求 VCF header 一致（样本顺序不同可通过重排序解决）

1
# 先下载 dbSNP VCF（或用已有的）
2
# 用 bcftools annotate 加注释
3
bcftools annotate \
4
    -a dbsnp.vcf.gz \
5
    -c ID \
6
    -O z -o annotated.vcf.gz \
7
    variants.vcf.gz
8

9
# -a = annotation source
10
# -c ID = 把 annotation source 的 ID 列复制过来

1
# 检查是不是 bgzip
2
file file.vcf.gz
3
# gzip compressed data  ← 这是 gzip
4
# BGZF compressed data   ← 这才是 bgzip
5

6
# 修复
7
zcat file.vcf.gz | bgzip -c > file_bgzip.vcf.gz
8
tabix -p vcf file_bgzip.vcf.gz

1
# 先把多等位基因拆分为双等位基因
2
bcftools norm -m -both variants.vcf.gz -O z -o split.vcf.gz
3
# -m -both = 把多等位拆成多行
4

5
# 或者只保留双等位
6
bcftools view -m2 -M2 variants.vcf.gz -O z -o biallelic.vcf.gz
7
# -m2 -M2 = 恰好2个等位基因

1
# -i 是 include（保留），-e 是 exclude（排除）
2
# 保留 QUAL ≥ 30：
3
bcftools filter -i 'QUAL>=30' input.vcf.gz   # ✅ 正确
4

5
# 排除 QUAL < 30：
6
bcftools filter -e 'QUAL<30' input.vcf.gz    # ✅ 也正确，效果一样
7

8
# 常见错误：用 -e 但写了保留条件
9
bcftools filter -e 'QUAL>=30' input.vcf.gz   # ❌ 把高质量的全排除了

1
# 查看 VCF 的染色体名
2
bcftools query -f '%CHROM\n' variants.vcf.gz | sort -u
3

4
# 如果需要统一，用 sed 改 header 和数据行
5
# （注意：这很危险，最好从源头统一参考基因组命名）

1
# 先统一样本列表
2
bcftools view -S sample_list.txt chr1.vcf.gz -O z -o chr1_reorder.vcf.gz
3
# -S = 只保留这些样本，按文件中的顺序重新排列
4

5
# 如果样本不同，用 merge 而非 concat
6
bcftools merge chr1.vcf.gz chr2.vcf.gz -O z -o merged.vcf.gz
7
# merge 自动处理样本不同和缺失

1
# VCF 规定 INFO 是分号分隔的键值对：
2
# DP=100;AF=0.5     ← 正确
3
# DP=100 AF=0.5     ← 错误
4

5
# 检查
6
bcftools view -H variants.vcf.gz | head -3 | cut -f8
7
# 看输出里分隔符是 ; 还是空格
8

9
# 修复（慎用）
10
bcftools view -H variants.vcf.gz | awk 'BEGIN{OFS="\t"}{gsub(/ /, ";", $8); print}' \
11
    | cat <(bcftools view -h variants.vcf.gz) - | bcftools view -O z -o fixed.vcf.gz

1
# 检查索引是否存在
2
ls -lh variants.vcf.gz.tbi
3

4
# 确保所有 .vcf.gz 都有索引
5
for vcf in *.vcf.gz; do
6
    if [ ! -f "${vcf}.tbi" ]; then
7
        echo "Missing index for $vcf, building..."
8
        tabix -p vcf "$vcf"
9
    fi
10
done

1
# 通常的过滤策略：
2
# 位点级别：QUAL ≥ 30
3
# 样本级别：GQ ≥ 20（对于群体分析可以适当放宽）
4
bcftools filter -i 'QUAL>=30 && FMT/GQ>=20' variants.vcf.gz