gatk-variant-calling

GATK — 生殖細胞系列バリアントコールパイプライン

概要

GATK（Genome Analysis Toolkit）は、Illumina WGS および WES データから SNP とインデルをコールするための GATK Best Practices ワークフローを実装しています。このパイプラインでは、サンプルごとに HaplotypeCaller を実行して GVCF ファイルを生成し、GenomicsDBImport で GVCF を統合し、GenotypeGVCFs で統合ジェノタイピングを行った後、VQSR（Variant Quality Score Recalibration）またはハードフィルターでバリアントをフィルタリングします。GATK の入力には、BWA-MEM2 でアラインメントし、重複マーキングと塩基品質スコア再キャリブレーション（BQSR）を行った BAM ファイルが必要です。前処理と後処理には Picard tools、samtools、bcftools を連携して使用します。GATK4 ワークフローは、研究および臨床ゲノミクスにおける生殖細胞系列バリアントコールの NIH/ENCODE 標準です。

使用する場面

• 集団遺伝学や臨床バリアント解析のため、WGS または WES サンプルから生殖細胞系列 SNP とインデルをコールする場合
• コホート規模の研究（家族、症例対照）で複数サンプルの統合ジェノタイピングを行う場合
• HaplotypeCaller の実行前に塩基品質スコア再キャリブレーション（BQSR）を適用して、バリアントコールの精度を向上させる場合
• スケーラブルなコホート拡張に向けて GVCF ファイルを生成する場合：既存サンプルを再処理せずに新しいサンプルを追加できます
• Ensembl VEP、ANNOVAR、SnpEff による下流のアノテーション用に、バリアントコールセットを作成する場合
• 比較可能な精度を保ちながら、より高速なディープラーニング手法を使うには、代わりに DeepVariant（Google）を使用してください
• アセンブリベースの局所再アラインメントを行わずに迅速にバリアントコールするには、代わりに bcftools call を使用してください

前提条件

• ソフトウェア：GATK4、Java 17+、samtools、BWA-MEM2
• 参照ファイル：ゲノム FASTA ＋既知バリアント VCF（dbSNP、1000G、Mills indels）
• 入力：@RG リードグループヘッダーを含む、重複マーキング済みのソート済み BAM（BWA-MEM2 由来）

インストール前に確認してください：現在の環境（例：pixi / conda 環境内）で、ツールをすでに利用できる場合があります。まず command -v gatk を実行し、パスが返された場合は以下のインストールコマンドをスキップしてください。pixi プロジェクト内で実行する場合は、gatk をそのまま実行するのではなく、pixi run gatk を使ってください。

# Install GATK4
wget https://github.com/broadinstitute/gatk/releases/download/4.6.0.0/gatk-4.6.0.0.zip
unzip gatk-4.6.0.0.zip
export GATK="$PWD/gatk-4.6.0.0/gatk"

# Or with conda
conda install -c bioconda gatk4

# Verify
gatk --version
# GATK v4.6.0.0

# Download GATK resource bundle files (GRCh38)
# From gs://gcp-public-data--broad-references/hg38/v0/ (requires gsutil or Broad FTP)

事前確認

解析コードを書く前に、ユーザーと以下の点を確認してください。

decisions:
  - id: D1
    param: referenceBundle
    kind: derived
    source: upstream
    ask: "Which reference FASTA were the BAMs aligned to?"
    default: "carried from the alignment stage; must match exactly"

  - id: D2
    param: callingIntervals
    kind: required
    source: user
    ask: "Call across the whole genome, or restrict to capture targets or regions of interest?"
    default: null

  - id: D3
    param: cohortMode
    kind: required
    source: user
    ask: "Call each sample on its own, or emit GVCFs now and genotype the cohort jointly afterwards?"
    default: "joint genotyping when more than one sample is being analysed together"

  - id: D4
    param: filteringStrategy
    kind: required
    source: user
    depends_on: [D3]
    ask: "Filter calls with fixed hard thresholds, or train a variant-recalibration model on the cohort?"
    default: "hard filters; recalibration needs a cohort large enough to train on"

  - id: D5
    param: knownSites
    kind: optional
    source: user
    ask: "Annotate calls against a known-variant resource such as dbSNP?"
    default: "none"

  - id: D6
    param: callConfidence
    kind: optional
    source: user
    ask: "How confident must a genotype be before the variant is emitted at all?"
    default: 30

  - id: D7
    param: threadsAndHeap
    kind: never_ask
    source: data
    reason: "Pair-HMM threads and JVM heap size affect runtime and memory, not the calls"
    default: "4 threads, heap sized to the genome"

D4 は D3 を前提としています。バリアントの再キャリブレーションはコホートを対象とする手法であり、 単一サンプルでは学習に使えるデータがありません。そのため、どのフィルタリング方法が適切かは、 コホートに関する質問への回答が得られて初めて判断できます。

クイックスタート

GENOME="GRCh38.fa"
DBSNP="dbsnp_146.hg38.vcf.gz"

# Run HaplotypeCaller in GVCF mode
gatk HaplotypeCaller \
    -R $GENOME \
    -I sample1.markdup.bam \
    -O sample1.g.vcf.gz \
    -ERC GVCF \
    --dbsnp $DBSNP \
    --native-pair-hmm-threads 4

echo "GVCF: sample1.g.vcf.gz"

ワークフロー

ステップ 1：塩基品質スコア再キャリブレーション（BQSR）

バリアントコールの前に、塩基品質スコアの系統誤差を補正します。

GENOME="GRCh38.fa"
KNOWN_SITES="dbsnp_146.hg38.vcf.gz Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
KNOWN_FLAGS=$(printf -- '--known-sites %s ' $KNOWN_SITES)

# Step 1a: Build recalibration table
gatk BaseRecalibrator \
    -R $GENOME \
    -I sample1.markdup.bam \
    $KNOWN_FLAGS \
    -O sample1.recal.table

# Step 1b: Apply recalibration
gatk ApplyBQSR \
    -R $GENOME \
    -I sample1.markdup.bam \
    --bqsr-recal-file sample1.recal.table \
    -O sample1.bqsr.bam

echo "BQSR BAM: sample1.bqsr.bam"
samtools flagstat sample1.bqsr.bam | head -3

ステップ 2：HaplotypeCaller でバリアントをコール（GVCF モード）

サンプルごとにバリアントコールを実行し、統合ジェノタイピング用の中間 GVCF を生成します。

# HaplotypeCaller in GVCF mode (recommended for cohort analysis)
gatk HaplotypeCaller \
    -R GRCh38.fa \
    -I sample1.bqsr.bam \
    -O gvcfs/sample1.g.vcf.gz \
    -ERC GVCF \
    --dbsnp dbsnp_146.hg38.vcf.gz \
    --native-pair-hmm-threads 4

# For WES: specify target intervals
# gatk HaplotypeCaller ... -L exome_targets.interval_list --interval-padding 100

echo "GVCF: gvcfs/sample1.g.vcf.gz"
zcat gvcfs/sample1.g.vcf.gz | grep -v "^#" | wc -l

ステップ 3：GenomicsDBImport で GVCF を統合

効率的な統合ジェノタイピングのため、サンプルごとの GVCF をマージします。

# Create sample map file: sample_name\tpath_to_gvcf
printf "ctrl_1\tgvcfs/ctrl_1.g.vcf.gz\n" > sample_map.txt
printf "ctrl_2\tgvcfs/ctrl_2.g.vcf.gz\n" >> sample_map.txt
printf "treat_1\tgvcfs/treat_1.g.vcf.gz\n" >> sample_map.txt
printf "treat_2\tgvcfs/treat_2.g.vcf.gz\n" >> sample_map.txt

# Import GVCFs into GenomicsDB for each chromosome
for CHR in chr1 chr2 chr3; do
    gatk GenomicsDBImport \
        --sample-name-map sample_map.txt \
        --genomicsdb-workspace-path genomicsdb/${CHR} \
        -L $CHR \
        --reader-threads 4
done

echo "GenomicsDB created for $(ls genomicsdb/ | wc -l) chromosomes"

ステップ 4：GenotypeGVCFs で統合ジェノタイピング

GenomicsDB 全体で、すべてのサンプルを同時にジェノタイピングします。

# Joint genotype all samples
mkdir -p vcfs
for CHR in chr1 chr2 chr3; do
    gatk GenotypeGVCFs \
        -R GRCh38.fa \
        -V gendb://genomicsdb/${CHR} \
        --dbsnp dbsnp_146.hg38.vcf.gz \
        -O vcfs/cohort_${CHR}.vcf.gz
done

# Merge per-chromosome VCFs
gatk MergeVcfs \
    $(ls vcfs/cohort_chr*.vcf.gz | sed 's/^/-I /') \
    -O vcfs/cohort_all.vcf.gz

echo "Joint genotyping complete: vcfs/cohort_all.vcf.gz"
gatk CountVariants -V vcfs/cohort_all.vcf.gz

ステップ 5：バリアントのフィルタリング（ハードフィルター）

VQSR の検出力が不足する小規模コホートには、ハードフィルターを適用します。

# Separate SNPs and indels
gatk SelectVariants -V vcfs/cohort_all.vcf.gz --select-type-to-include SNP -O vcfs/snps.vcf.gz
gatk SelectVariants -V vcfs/cohort_all.vcf.gz --select-type-to-include INDEL -O vcfs/indels.vcf.gz

# Apply hard filters: SNPs
gatk VariantFiltration \
    -V vcfs/snps.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    --filter-expression "MQ < 40.0" --filter-name "MQ40" \
    --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
    -O vcfs/snps_filtered.vcf.gz

# Apply hard filters: Indels
gatk VariantFiltration \
    -V vcfs/indels.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 200.0" --filter-name "FS200" \
    --filter-expression "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20" \
    -O vcfs/indels_filtered.vcf.gz

# Merge filtered SNPs + indels
gatk MergeVcfs -I vcfs/snps_filtered.vcf.gz -I vcfs/indels_filtered.vcf.gz \
    -O vcfs/cohort_filtered.vcf.gz

echo "PASS variants: $(bcftools view -f PASS vcfs/cohort_filtered.vcf.gz | grep -v '^#' | wc -l)"

ステップ 6：Python で VCF の結果を解析

バリアントを抽出し、遺伝子情報でアノテーションを付けて、DataFrame を作成します。

import subprocess
import pandas as pd
import io

# Use bcftools query to extract fields from filtered VCF
result = subprocess.run(
    ["bcftools", "query",
     "-f", "%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t%INFO/QD\t%INFO/FS\n",
     "-i", "FILTER='PASS'",
     "vcfs/cohort_filtered.vcf.gz"],
    capture_output=True, text=True
)

cols = ["CHR", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "QD", "FS"]
df = pd.read_csv(io.StringIO(result.stdout), sep="\t", names=cols)
df["QUAL"] = pd.to_numeric(df["QUAL"], errors="coerce")
df["QD"] = pd.to_numeric(df["QD"], errors="coerce")

print(f"PASS variants: {len(df)}")
print(f"SNPs:  {(df['REF'].str.len() == 1) & (df['ALT'].str.len() == 1)).sum()}")
print(f"Indels: {((df['REF'].str.len() > 1) | (df['ALT'].str.len() > 1)).sum()}")
print(df.head())
df.to_csv("pass_variants.tsv", sep="\t", index=False)

主要なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
-ERC	NONE	GVCF, BP_RESOLUTION	リファレンス信頼度モードを出力します。コホートのワークフローでは GVCF を指定します
--native-pair-hmm-threads	4	1–32	HaplotypeCaller の pair-HMM のスレッド数（最も CPU 負荷の高いステップ）
-L	ゲノム全体	インターバルファイルまたは chr	呼び出し対象をインターバルに限定します（エクソームのターゲット、BED 領域）
--dbsnp	—	VCF パス	出力に rsID アノテーションを付与するための dbSNP VCF
--stand-call-conf	30	0–100	バリアントコールを出力するための最小ジェノタイプ品質スコア
-G	StandardAnnotation	アノテーショングループ	適用するアノテーションモジュール（HaplotypeCaller では StandardHCAnnotation）
--sample-name-map	—	TSV ファイル	GenomicsDBImport 用のサンプルと GVCF の対応付け
--reader-threads	1	1–16	GenomicsDBImport の読み込みに使うスレッド数
--filter-expression	—	JEXL 式	ハードフィルター式（例："QD < 2.0"）
--java-options	—	-Xmx4g	Java ヒープサイズ。大きなゲノムでは -Xmx16g を使用します

よく使うレシピ

レシピ 1：単一サンプルのバリアントコール（コホートなし）

# For a single sample, skip GenomicsDBImport and call directly
gatk HaplotypeCaller \
    -R GRCh38.fa \
    -I sample.bqsr.bam \
    -O sample.vcf.gz \
    --dbsnp dbsnp_146.hg38.vcf.gz \
    --native-pair-hmm-threads 8

# Hard filter the single-sample VCF
gatk VariantFiltration \
    -V sample.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    -O sample_filtered.vcf.gz

echo "PASS variants: $(bcftools view -f PASS sample_filtered.vcf.gz | grep -v '^#' | wc -l)"

レシピ 2：Snakemake でパイプライン全体を実行する

# Snakefile — GATK BQSR + HaplotypeCaller
configfile: "config.yaml"
SAMPLES = config["samples"]
GENOME  = config["genome"]
DBSNP   = config["dbsnp"]
KNOWN   = config["known_sites"]

rule all:
    input: expand("gvcfs/{sample}.g.vcf.gz", sample=SAMPLES)

rule bqsr:
    input: bam="markdup/{sample}.markdup.bam"
    output: bam="bqsr/{sample}.bqsr.bam"
    shell:
        """
        gatk BaseRecalibrator -R {GENOME} -I {input.bam} \
             --known-sites {KNOWN} -O bqsr/{wildcards.sample}.recal.table &&
        gatk ApplyBQSR -R {GENOME} -I {input.bam} \
             --bqsr-recal-file bqsr/{wildcards.sample}.recal.table \
             -O {output.bam}
        """

rule haplotypecaller:
    input: bam="bqsr/{sample}.bqsr.bam"
    output: gvcf="gvcfs/{sample}.g.vcf.gz"
    threads: 4
    shell:
        "gatk HaplotypeCaller -R {GENOME} -I {input.bam} -O {output.gvcf} "
        "-ERC GVCF --dbsnp {DBSNP} --native-pair-hmm-threads {threads}"

想定される出力

出力	形式	説明
*.g.vcf.gz	GVCF	リファレンス信頼度ブロックを含むサンプルごとの GVCF。GenomicsDBImport への入力です
cohort_all.vcf.gz	VCF	複数サンプルを統合ジェノタイプした VCF。フィルタリング前です
cohort_filtered.vcf.gz	VCF	フィルタリング済み VCF。合格したバリアントは FILTER=PASS です
*.recal.table	BQSR	BaseRecalibrator による塩基品質再較正テーブル
*.bqsr.bam	BAM	再較正済み BAM。HaplotypeCaller の入力として使用します
genomicsdb/	ディレクトリ	統合ジェノタイピング用の染色体ごとの GenomicsDB ワークスペース

トラブルシューティング

問題	原因	解決策
SAM/BAM file has no @RG header	BWA-MEM2 のリードグループがありません	-R "@RG\tID:...\tSM:...\tPL:ILLUMINA" を指定して再アラインメントします
Java OutOfMemoryError	ヒープサイズが不足しています	--java-options "-Xmx16g" 以上を追加します
HaplotypeCaller が非常に遅い	HMM がシングルスレッドで実行されています	--native-pair-hmm-threads 8 を追加します。インターバルリストを使って chr ごとに並列化します
GVCF の出力が空です	インターバルが誤っているか、領域内にリードがありません	samtools view -c sample.bam chrN でリード数を確認します
VQSR が失敗する（バリアントが 10k 未満）	トレーニングに使えるバリアントが少なすぎます	小規模なコホートやエクソームでは VQSR ではなくハードフィルターを使用します
既存のパスに対する GenomicsDB のインポートが失敗します	GenomicsDB ワークスペースがすでに存在します	再実行する前に、既存のワークスペースを削除します：rm -rf genomicsdb/chr1
IndexOutOfBoundsException	BAM とリファレンス間で染色体名が一致していません	ゲノム FASTA と BAM で同じ染色体名（chr1 と 1 のどちらか）が使われていることを確認します
BCFtools/tabix のインデックスがありません	Tabix インデックス（.tbi）が作成されていません	gatk IndexFeatureFile -I file.vcf.gz または tabix -p vcf file.vcf.gz を実行します

参考資料

• GATK ドキュメント — 公式のパイプラインガイドとツールドキュメント
• 生殖細胞系列の短いバリアントに関する GATK ベストプラクティス — 手順を追ったプロトコル
• GATK GitHub：broadinstitute/gatk — ソースコード、リリース、リソースバンドル
• Van der Auwera GA & O'Connor BD（2020）Genomics in the Cloud — O'Reilly Media。GATK4 の包括的なガイド