bcftools-variant-manipulation

bcftools — VCF/BCF バリアント操作ツールキット

概要

bcftools は、HTSlib エコシステムにおける VCF（Variant Call Format）および BCF（Binary Call Format）ファイル処理の標準的なコマンドラインツールキットです。バリアントコール後のワークフロー全体に対応しており、フォーマット変換、品質フィルタリング、バリアントの正規化、複数サンプルのマージ、外部データベースによるアノテーション、遺伝子型の抽出、QC 統計の算出が行えます。bcftools はストリーミングを前提として設計されており、ほとんどのコマンドは標準入力から読み込み、標準出力に書き込むため、大規模なコホートを扱うメモリ効率の高いパイプラインに適しています。

使用する場面

• バリアントコール後に、品質（QUAL、DP、AF）に基づいてバリアントをフィルタリングする
• 複数サンプルの VCF ファイルをマージして、統合コールセットを作成する
• バリアントコールに rsID や遺伝子アノテーションを追加する
• 特定のフィールド（遺伝子型、アレル深度）を表形式で抽出する
• indel 表現を正規化し、複数アレルのレコードを分割する
• 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 プロジェクト内で実行する場合は、bcftools 単体ではなく pixi run 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 ファイルをマージするか、複数染色体の VCF ファイルを連結します。

# 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 を使用します（入出力のオーバーヘッドがありません）。保存するファイルや 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	Most	v	v,z,b,u	出力形式：VCF、bgzip-VCF、BCF、非圧縮BCF
-r	Most	—	chr:pos-end	領域フィルター（tabix インデックスが必要）
-s	Most	All	サンプル名	特定のサンプルを含めます（除外するには接頭辞 ^ を付けます）
-i	filter	—	式	式に一致するバリアントを含めます
-e	filter	—	式	式に一致するバリアントを除外します
-f	query	—	書式文字列	カスタム出力書式文字列
--threads	Most	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 インデックスが必要です。bgzipped 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 の前に、各サンプルにサイトレベルの QC（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)

レシピ：1 サンプルのヘテロ接合サイトを抽出する

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"

レシピ：2 つの 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）「Twelve years of SAMtools and BCFtools」— GigaScience 10(2)
• VCF ファイル形式仕様 — フィールド定義とエンコード規則