bcftools-variant-manipulation

bcftools — VCF/BCF 变异操作工具包

概述

bcftools 是 HTSlib 生态系统中用于处理 VCF（Variant Call Format，变异调用格式）和 BCF（Binary Call Format，二进制调用格式）文件的标准命令行工具包。它涵盖变异调用后的完整工作流：格式转换、质量筛选、变异规范化、多样本合并、使用外部数据库进行注释、基因型提取和 QC 统计。bcftools 的设计采用流式处理——大多数命令从标准输入读取并写入标准输出，因此非常适合在大型队列上构建内存高效的流水线。

适用场景

• 在变异调用后，按质量（QUAL、DP、AF）筛选变异
• 将多个样本的 VCF 文件合并为联合调用集
• 为变异调用添加 rsID 或基因注释
• 将特定字段（基因型、等位基因深度）提取为表格输出
• 规范化插入缺失变异的表示形式并拆分多等位基因记录
• 根据 pileup 输出调用变异（mpileup + call）
• 计算每个样本及整体的 VCF QC 统计数据
• 在人类样本中使用局部重比对调用变异时，改用 GATK HaplotypeCaller
• 进行群体遗传学统计（Fst、LD、Hardy-Weinberg）时，改用 VCFtools
• 在 HTSlib 流水线中使用 bcftools；使用 picard 标记重复项并计算文库指标

先决条件

• 安装：bcftools 1.17+（属于包含 samtools 的 HTSlib 套件）
• 输入要求：用于区域查询的 VCF 或经 BGzip 压缩并建立 tabix 索引的 VCF（.vcf.gz + .vcf.gz.tbi）
• 配套工具：用于 BAM 处理的 samtools；用于 VCF 索引的 tabix

安装前请检查：当前环境中可能已提供该工具（例如，在 pixi / conda 环境中）。先运行 command -v bcftools，如果返回路径，则跳过下面的安装命令。在 pixi 项目中运行时，请使用 pixi run bcftools 调用该工具，而不是直接运行 bcftools。

# Bioconda (recommended — installs HTSlib suite)
conda install -c bioconda bcftools

# Homebrew (macOS)
brew install bcftools

# Verify
bcftools --version | head -1
# bcftools 1.20

# Index a VCF for region queries
bcftools index -t variants.vcf.gz   # creates .tbi
bcftools index -c variants.vcf.gz   # creates .csi (for chromosomes > 512 Mb)

快速入门

# Typical post-calling workflow: normalize → filter → annotate → extract
bcftools norm -d any -f reference.fa variants.vcf.gz \
  | bcftools filter -i 'QUAL>20 && DP>10' \
  | bcftools annotate -a dbSNP.vcf.gz -c ID \
  | bcftools view -O z -o final.vcf.gz

# Index the output
bcftools index -t final.vcf.gz

# Count variants at each stage
bcftools stats final.vcf.gz | grep "^SN"

核心 API

模块 1：VCF/BCF 输入输出与格式转换

在文本 VCF 和二进制 BCF 之间转换；压缩并建立索引，以支持随机访问。

# VCF → compressed BCF (fastest format for piping)
bcftools view -O b -o variants.bcf variants.vcf

# BCF → VCF (for human-readable output)
bcftools view -O v -o variants.vcf variants.bcf

# VCF → bgzipped + indexed (standard archive format)
bcftools view -O z -W -o variants.vcf.gz variants.vcf
# -W automatically creates .tbi index after writing

# Extract specific samples
bcftools view -s sample1,sample2 -O z -o subset.vcf.gz variants.vcf.gz

# Exclude samples (prefix with ^)
bcftools view -s ^outlier_sample -O z -o cleaned.vcf.gz variants.vcf.gz

# Extract by region (fast; requires index)
bcftools view -r chr1:1000000-2000000 variants.vcf.gz -O v -o chr1_region.vcf

# Streaming pipeline: no intermediate files
samtools mpileup -Ou input.bam | bcftools call -m -Oz -o calls.vcf.gz

模块 2：变异筛选

应用质量阈值和基于 FLAG 的筛选条件，以保留高置信度的调用结果。

# Expression-based filter (include)
bcftools filter -i 'QUAL>20 && DP>10' variants.vcf.gz -O z -o filtered.vcf.gz

# Expression-based filter (exclude)
bcftools filter -e 'QUAL<10 || DP<5' variants.vcf.gz -O v -o filtered.vcf

# Soft filter: mark but keep (sets FILTER field to label)
bcftools filter -s LowQual -e 'QUAL<20' variants.vcf.gz -O z -o soft_filtered.vcf.gz
# Variants with QUAL<20 get FILTER="LowQual"; others get FILTER=PASS

# Keep only PASS variants
bcftools view -f PASS variants.vcf.gz -O z -o pass_only.vcf.gz

# SNP-only output
bcftools view --type snps variants.vcf.gz -O z -o snps.vcf.gz

# Indel-only output
bcftools view --type indels variants.vcf.gz -O z -o indels.vcf.gz

# Filter by allele frequency and depth
bcftools filter -i 'AF>0.1 && DP>20 && MQ>40' variants.vcf.gz -O z -o confident.vcf.gz

# Remove SNPs within 3 bp of indels
bcftools filter --SnpGap 3 variants.vcf.gz -O z -o gapfiltered.vcf.gz

模块 3：VCF 查询与提取

将 VCF 内容转换为表格文本，以便进行下游分析。

# Extract chrom, position, ref, alt, quality
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\n' variants.vcf.gz > variants.txt

# With header row (-H adds #-prefixed column names)
bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\n' variants.vcf.gz > variants.tsv

# Per-sample genotypes and allele depths
bcftools query -f '[%SAMPLE\t%GT\t%AD\n]' variants.vcf.gz > genotypes.txt
# Output: sample1  0/1  25,18  (ref_depth,alt_depth)

# Rare variants (AF < 1%)
bcftools query -i 'AF<0.01' -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\n' \
    variants.vcf.gz > rare_variants.txt

# Count variants per chromosome
bcftools query -f '%CHROM\n' variants.vcf.gz | sort | uniq -c | sort -rn

# Extract genotype matrix across all samples
bcftools query -f '%CHROM:%POS\t[%GT\t]\n' -H variants.vcf.gz > genotype_matrix.tsv

模块 4：多文件操作

合并来自多个样本的 VCF 文件（merge），或拼接来自多个染色体的 VCF 文件（concat）。

# Merge: join VCFs from DIFFERENT sample sets (same variants)
bcftools merge sample1.vcf.gz sample2.vcf.gz sample3.vcf.gz \
    -O z -o cohort.vcf.gz

# Merge with auto-indexing and threading
bcftools merge -O b -W --threads 4 sample*.vcf.gz > cohort.bcf

# Concat: join VCFs from SAME sample set (different chromosomes or batches)
bcftools concat chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz -O z -o full.vcf.gz

# Concat with overlap handling (from batched calling)
bcftools concat -a --threads 4 batch*.vcf.gz -O z -o concat.vcf.gz

# Pipeline: merge → filter → normalize
bcftools merge sample1.vcf.gz sample2.vcf.gz \
    | bcftools filter -i 'QUAL>20' \
    | bcftools norm -d any -f genome.fa \
    | bcftools view -O z -o merged_clean.vcf.gz
bcftools index -t merged_clean.vcf.gz

# Extract genotype matrix from merged cohort
bcftools merge cohort*.vcf.gz | bcftools query -f '[%GT\t]\n' > gt_matrix.tsv

模块 5：变异注释

为 VCF 记录添加标识符、基因注释或外部数据。

# Add rsIDs from dbSNP
bcftools annotate -a dbSNP.vcf.gz -c ID variants.vcf.gz -O z -o rsid_annotated.vcf.gz

# Annotate with BED file (adds gene names)
bcftools annotate -a genes.bed.gz \
    -h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description="Gene name">') \
    -c CHROM,FROM,TO,GENE \
    variants.vcf.gz -O z -o gene_annotated.vcf.gz

# Remove unwanted INFO fields
bcftools annotate -x INFO/AC,INFO/AN,INFO/MQ variants.vcf.gz -O v -o stripped.vcf

# Normalize: left-align indels, split multi-allelic records
bcftools norm -f reference.fa variants.vcf.gz -O v -o normalized.vcf

# Split multi-allelic sites into separate records
bcftools norm -m -any variants.vcf.gz -O v -o split.vcf

# Deduplicate overlapping records
bcftools norm -d any variants.vcf.gz -O z -o deduped.vcf.gz

# Full normalize pipeline
bcftools norm -m -any variants.vcf.gz | \
    bcftools norm -d any -f reference.fa | \
    bcftools view -O z -o normalized_split.vcf.gz

模块 6：统计与 QC

生成汇总指标和每个样本的变异计数。

# Full VCF statistics report
bcftools stats variants.vcf.gz > qc.stats.txt

# Extract Summary Numbers section only
grep "^SN" qc.stats.txt | cut -f3,4
# number of records:    45231
# number of SNPs:       38941
# number of indels:     6290
# ...

# Per-sample stats (PSC = per-sample counts)
bcftools stats -s - variants.vcf.gz | grep "^PSC" > per_sample.txt
# cols: id  sample  hom_RR  het  hom_AA  ts  tv  indel  missing  singleton

# Transition/transversion ratio (genome-wide QC)
bcftools stats variants.vcf.gz | grep "Ts/Tv"
# Ts/Tv ratio: 2.06 (healthy WGS; <1.8 or >2.2 suggests quality issues)

# Check for sample contamination (F-statistic per sample)
bcftools stats -s - variants.vcf.gz | grep "^PSC" | awk '{print $2, $9}'
# sample  F_missing (high = poor sample quality)

# Variant calling (mpileup → call pipeline)
samtools mpileup -Ou -f genome.fa *.bam | bcftools call -m -v -Oz -o calls.vcf.gz
bcftools stats calls.vcf.gz | grep "^SN"

关键概念

输出格式标记

-O v  → VCF text (uncompressed)       default for human inspection
-O z  → bgzipped VCF (.vcf.gz)        standard for archiving
-O b  → binary BCF (compressed)        fastest for piping
-O u  → binary BCF (uncompressed)      fastest output (no compression)

规则：中间流水线步骤使用 -O b 或 -O u（不产生 I/O 开销）。对于需要存储或使用 tabix 建立索引的文件，使用 -O z。

筛选表达式语法

表达式使用 INFO 和 FORMAT 字段以及比较运算符：

# INFO fields (one value per variant)
QUAL>20          # quality score
DP>10            # total depth
AF<0.05          # allele frequency
MQ>40            # mapping quality

# FORMAT fields (per-sample; use [] to iterate)
GT=="1/1"        # homozygous alternate
AD[1]>5          # alt allele depth > 5
GQ>20            # genotype quality

# Combined
QUAL>20 && DP>10 && AF>0.01
(GT=="0/0" || GT=="1/1") && GQ>30

常见工作流

工作流 1：变异 QC 与筛选流水线

目标：对原始变异调用集进行规范化、筛选和注释，以供下游分析使用。

#!/bin/bash
VCF="raw_calls.vcf.gz"
REF="reference.fa"
DBSNP="dbSNP_hg38.vcf.gz"
FINAL="variants_filtered_annotated.vcf.gz"

# 1. Normalize: left-align indels, split multi-allelic, deduplicate
bcftools norm -m -any $VCF \
    | bcftools norm -d any -f $REF \
    | bcftools view -O z -o normalized.vcf.gz
bcftools index -t normalized.vcf.gz

# 2. Filter by quality and depth
bcftools filter -i 'QUAL>20 && DP>10' normalized.vcf.gz \
    | bcftools filter --SnpGap 3 \
    | bcftools view -f PASS -O z -o filtered.vcf.gz
bcftools index -t filtered.vcf.gz

# 3. Annotate with rsIDs
bcftools annotate -a $DBSNP -c ID filtered.vcf.gz -O z -o $FINAL
bcftools index -t $FINAL

# Report variant counts
echo "Final variant count:"
bcftools stats $FINAL | grep "number of records"

工作流 2：多样本队列合并与基因型提取

目标：将每个样本的 VCF 合并为队列 VCF；为 GWAS 提取基因型矩阵。

#!/bin/bash
SAMPLES=(sample1 sample2 sample3 sample4 sample5)

# 1. Ensure all VCFs are indexed
for s in "${SAMPLES[@]}"; do
    bcftools index -t ${s}.vcf.gz
done

# 2. Merge into cohort VCF (only sites present in ALL samples: -m none)
bcftools merge -m none "${SAMPLES[@]/%/.vcf.gz}" \
    -O z --threads 8 -o cohort.vcf.gz
bcftools index -t cohort.vcf.gz

# 3. Filter: PASS, SNPs only, MAF > 1%
bcftools view -f PASS --type snps cohort.vcf.gz \
    | bcftools filter -i 'AF>0.01 && AF<0.99' \
    | bcftools view -O z -o cohort_snps_filtered.vcf.gz

# 4. Extract numeric genotype matrix (for plink/R/Python)
bcftools query -H \
    -f '%CHROM\t%POS\t%REF\t%ALT\t[%GT\t]\n' \
    cohort_snps_filtered.vcf.gz > genotype_matrix.tsv
echo "Genotype matrix: $(wc -l < genotype_matrix.tsv) variants x ${#SAMPLES[@]} samples"

关键参数

参数	命令	默认值	范围/选项	作用
-O	大多数	v	v,z,b,u	输出格式：VCF、bgzip-VCF、BCF、未压缩 BCF
-r	大多数	—	chr:pos-end	区域筛选（需要 tabix 索引）
-s	大多数	全部	样本名称	包含指定样本（添加前缀 ^ 可排除）
-i	filter	—	表达式	包含符合表达式的变异
-e	filter	—	表达式	排除符合表达式的变异
-f	query	—	格式字符串	自定义输出格式字符串
--threads	大多数	1	1–N	压缩/解压线程数
-a	annotate	—	文件路径	注释来源（BED、VCF、TSV）
-m	norm	none	-any、+any	拆分（−）或合并（+）多等位记录
-d	norm	—	all、any、snps	去重策略
-W	view	—	标志	写入后自动创建索引
-v	call	—	标志	仅输出变异位点（跳过参考位点）

最佳实践

1. 区域查询前先建立索引：bcftools view -r chr1:... 需要 .tbi 或 .csi 索引。创建任何 bgzip 压缩的 VCF 后，都应运行 bcftools index -t output.vcf.gz。

2. 合并或注释前先进行标准化：不同的变异检测工具对同一插入缺失变异的表示方式可能不同。合并前运行 bcftools norm -m -any | bcftools norm -d any -f ref.fa，以避免重复记录。

3. 将 -O u 用于管道中间步骤：未压缩的 BCF 输出（-O u）可消除多步管道中的压缩/解压开销，通常比 -O z 快 2–3 倍。

4. 变异检测后检查 Ts/Tv 比率：bcftools stats variants.vcf.gz | grep Ts/Tv。人类全基因组测序（WGS）的预期值为 2.0–2.1；外显子组测序为 2.5–3.0。超出这些范围可能表示存在质量问题。

5. 大型队列合并前先筛选：合并前筛选每个样本的 VCF，可降低内存和 I/O 开销。运行 bcftools merge 前，对每个样本应用位点级质量控制（QUAL>20 && DP>5）。

6. 检查染色体命名是否一致：合并使用 chr1 格式和 1 格式的 VCF 时，bcftools 可能静默失败。使用 bcftools view -h file.vcf.gz | grep "^##contig" 检查。

常用示例

示例：统计筛选前后的变异数量

echo "Before:" $(bcftools view -c 1 raw.vcf.gz | wc -l)
echo "PASS only:" $(bcftools view -f PASS raw.vcf.gz | bcftools view -c 1 | wc -l)
echo "SNPs PASS:" $(bcftools view -f PASS --type snps raw.vcf.gz | bcftools view -c 1 | wc -l)

示例：提取一个样本的杂合位点

bcftools view -s SAMPLE_A variants.vcf.gz \
    | bcftools filter -i 'GT="0/1"' \
    | bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' > het_sites.txt
echo "$(wc -l < het_sites.txt) heterozygous sites"

示例：比较两个 VCF 文件的一致性

# Find variants unique to each file and shared
bcftools isec -p isec_dir file1.vcf.gz file2.vcf.gz
ls isec_dir/
# 0000.vcf: private to file1
# 0001.vcf: private to file2
# 0002.vcf: shared (from file1 perspective)
# 0003.vcf: shared (from file2 perspective)
echo "Shared: $(wc -l < isec_dir/0002.vcf) variants"

故障排除

问题	原因	解决方案
Missing index 错误	VCF 未建立索引（-r 区域查询需要索引）	运行 bcftools index -t file.vcf.gz
[E::vcf_parse_format] 解析错误	VCF FORMAT 或 INFO 字段格式错误	验证：bcftools view file.vcf.gz 2>&1 | head -5；检查来源工具版本
筛选结果为空	表达式过于严格或字段缺失	测试表达式：bcftools view -h file.vcf.gz | grep "##INFO=<ID=QUAL"
合并：样本重复	输入 VCF 中存在重复样本名称	合并前使用 bcftools reheader -s new_names.txt sample.vcf.gz 重命名
Ts/Tv 比率不正确（<1.8）	变异检测质量较低或覆盖度不足	应用更严格的质量筛选（QUAL>30 && DP>15）；检查比对质量
concat 因重叠失败	输入 VCF 中存在重叠区域	使用 -a 标志：bcftools concat -a file1.vcf.gz file2.vcf.gz
注释不匹配	chr 命名冲突（chr1 与 1）	检查：bcftools view -h file.vcf.gz | grep contig；使用 bcftools annotate --rename-chrs 重命名
query 返回空字段	样本的 FORMAT 字段未填充	检查 VCF 头部：bcftools view -h file.vcf.gz | grep FORMAT

相关技能

• samtools-bam-processing — 为 bcftools 变异检测工作流提供输入的 BAM 处理
• bedtools-genomic-intervals — 将 VCF 变异与基因组特征（基因、区域）求交集
• gget-genomic-databases — 查询 Ensembl/NCBI，为变异的基因背景添加注释

参考资料

• bcftools 文档 — 完整的命令参考
• GitHub：samtools/bcftools — 源代码、发行版本、问题跟踪器
• Danecek 等人（2021）“SAMtools 和 BCFtools 的十二年历程” — GigaScience 10(2)
• VCF 文件格式规范 — 字段定义和编码规则