bwa-mem2-dna-aligner

BWA-MEM2 — DNA ショートリードアライナー

概要

BWA-MEM2 は、BWT-FM インデックスを使用して、短い DNA リード（Illumina、50～250 bp）を参照ゲノムにアラインメントします。全ゲノムシーケンス（WGS）、全エクソームシーケンス（WES）、ChIP-seq、ATAC-seq の DNA アラインメントで標準的に使用されるアライナーです。BWA-MEM2 は元の BWA-MEM より 2 倍高速でありながら、同一の結果を生成します。GATK HaplotypeCaller と Picard ツールで必要となる適切なリードグループ（@RG）ヘッダーを含む SAM 形式で出力します。ペアエンドリードの場合は、適切なペアを識別し、キメラ／分割リードを補足アラインメントに解決します。

使用する場面

• バリアントコール（SNP、indel、SV）のために、WGS または WES の Illumina リードを参照ゲノムにアラインメントする場合
• ChIP-seq または ATAC-seq の DNA をアラインメントし、MACS3 でピークコールするための BAM ファイルを生成する場合
• @RG リードグループタグを含む GATK 互換の BAM ファイルを生成する場合
• 50 bp 以上のリードをアラインメントする場合。短いリード（50 bp 未満）には、BWA-backtrack の方が適していることがあります
• 更新された参照ゲノムアセンブリに、旧形式の FASTQ ファイルを再アラインメントする場合
• スプライスジャンクションをまたぐ RNA-seq リードには、代わりに STAR を使用します
• ローカルアラインメントの代替として、またはインデックスのサイズを最小限に抑える必要がある場合は Bowtie2 を使用します

前提条件

• ソフトウェア：bwa-mem2（conda またはコンパイル済みバイナリ）、samtools
• 参照データ：ゲノム FASTA（例：GRCh38、hg19、mm10）
• RAM：ヒトゲノムのインデックスには約 28 GB、マウスには 6～8 GB

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

# Install with conda (recommended)
conda install -c bioconda bwa-mem2 samtools

# Or download pre-compiled binary
wget https://github.com/bwa-mem2/bwa-mem2/releases/download/v2.2.1/bwa-mem2-2.2.1_x64-linux.tar.bz2
tar -jxf bwa-mem2-2.2.1_x64-linux.tar.bz2
export PATH="$PWD/bwa-mem2-2.2.1_x64-linux:$PATH"

# Verify
bwa-mem2 version
# 2.2.1

クイックスタート

# 1. Build genome index (~30 min, run once)
bwa-mem2 index GRCh38.fa

# 2. Align paired-end reads and sort
bwa-mem2 mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA" \
    GRCh38.fa sample1_R1.fastq.gz sample1_R2.fastq.gz \
    | samtools sort -@ 8 -o sample1.sorted.bam

# 3. Index the BAM
samtools index sample1.sorted.bam
echo "Aligned reads: $(samtools view -c -F 4 sample1.sorted.bam)"

ワークフロー

ステップ 1：参照ゲノムをダウンロードする

対象のアセンブリに対応する参照ゲノム FASTA ファイルを取得します。

# Download GRCh38 primary assembly (human)
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
mv GCA_000001405.15_GRCh38_no_alt_analysis_set.fna GRCh38.fa

# Or use ENSEMBL/GENCODE
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/GRCh38.primary_assembly.genome.fa.gz
gunzip GRCh38.primary_assembly.genome.fa.gz

echo "Reference size: $(du -sh GRCh38.fa)"

ステップ 2：BWA-MEM2 インデックスを作成する

参照ゲノムのインデックスを作成します。ゲノムごとに 1 回必要で、ヒトの場合は約 25～35 分かかります。

# Build index (~28 GB RAM required for human genome)
bwa-mem2 index GRCh38.fa

# This creates: GRCh38.fa.0123, GRCh38.fa.amb, GRCh38.fa.ann,
#               GRCh38.fa.bwt.2bit.64, GRCh38.fa.pac
echo "Index files: $(ls GRCh38.fa.* | wc -l) created"
ls -lh GRCh38.fa.*

ステップ 3：ペアエンドリードをアラインメントする

GATK との互換性に必要なリードグループヘッダーを指定して、FASTQ リードをアラインメントします。

# Align with read group (required for GATK)
# @RG fields: ID (run ID), SM (sample name), PL (platform), LB (library), PU (flowcell)
bwa-mem2 mem \
    -t 16 \
    -R "@RG\tID:sample1_run1\tSM:sample1\tPL:ILLUMINA\tLB:lib1\tPU:flowcell1" \
    GRCh38.fa \
    sample1_R1.fastq.gz \
    sample1_R2.fastq.gz \
    | samtools sort -@ 8 -m 2G -o sample1.sorted.bam

samtools index sample1.sorted.bam
echo "Alignment complete."
echo "Total reads:   $(samtools view -c sample1.sorted.bam)"
echo "Mapped reads:  $(samtools view -c -F 4 sample1.sorted.bam)"

ステップ 4：PCR 重複を識別する

バリアントコールの前に、光学重複と PCR 重複を削除するか、識別します。

# Option A: samtools markdup (fast)
samtools fixmate -m sample1.sorted.bam sample1.fixmate.bam
samtools sort -@ 8 -o sample1.fixmate.sorted.bam sample1.fixmate.bam
samtools markdup -@ 8 sample1.fixmate.sorted.bam sample1.markdup.bam
samtools index sample1.markdup.bam

echo "Duplication rate:"
samtools flagstat sample1.markdup.bam | grep "duplicate"

# Option B: Picard MarkDuplicates (GATK best practices)
picard MarkDuplicates \
    INPUT=sample1.sorted.bam \
    OUTPUT=sample1.markdup.bam \
    METRICS_FILE=sample1.dupmetrics.txt \
    REMOVE_DUPLICATES=false \
    CREATE_INDEX=true
cat sample1.dupmetrics.txt | grep -A2 "ESTIMATED"

ステップ 5：アラインメント品質を評価する

アラインメント統計を生成し、主要な品質指標を確認します。

# Full alignment statistics
samtools flagstat sample1.markdup.bam > sample1.flagstat.txt
cat sample1.flagstat.txt

# Coverage statistics
samtools coverage sample1.markdup.bam | head -30

# Parse key metrics with Python
python3 - << 'EOF'
from pathlib import Path

flagstat = Path("sample1.flagstat.txt").read_text()
for line in flagstat.splitlines():
    if any(kw in line for kw in ["total", "mapped", "properly paired", "duplicate"]):
        print(line)
EOF

ステップ 6：WGS/WES パイプラインを完了する → バリアントコール

BWA-MEM2 の出力を GATK HaplotypeCaller に直接パイプします。

#!/bin/bash
# Complete WGS alignment → variant calling pipeline
GENOME="GRCh38.fa"
SAMPLE="sample1"
R1="data/${SAMPLE}_R1.fastq.gz"
R2="data/${SAMPLE}_R2.fastq.gz"
THREADS=16
OUTDIR="results/${SAMPLE}"
mkdir -p "$OUTDIR"

# Step 1: Align + sort
bwa-mem2 mem -t $THREADS \
    -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA\tLB:lib1\tPU:run1" \
    $GENOME $R1 $R2 \
    | samtools sort -@ 8 -o $OUTDIR/${SAMPLE}.sorted.bam
samtools index $OUTDIR/${SAMPLE}.sorted.bam

# Step 2: Mark duplicates
samtools fixmate -m $OUTDIR/${SAMPLE}.sorted.bam - \
    | samtools sort -@ 8 \
    | samtools markdup -@ 8 - $OUTDIR/${SAMPLE}.markdup.bam
samtools index $OUTDIR/${SAMPLE}.markdup.bam

# Step 3: GATK variant calling
gatk HaplotypeCaller \
    -R $GENOME \
    -I $OUTDIR/${SAMPLE}.markdup.bam \
    -O $OUTDIR/${SAMPLE}.g.vcf.gz \
    -ERC GVCF \
    --native-pair-hmm-threads 4

echo "Pipeline complete: $OUTDIR/${SAMPLE}.g.vcf.gz"

主なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
-t	1	1–64	CPU スレッド数。本番実行では 8～16 を使用します
-R	—	@RG\tID:...\tSM:...	リードグループ文字列。GATK との互換性に必須です
-k	19	10–28	最小シード長。値を小さくすると、短いリードに対する感度が高くなります
-w	100	50–500	Smith-Waterman アラインメントのバンド幅
-M	オフ	フラグ	分割／補足リードをセカンダリアラインメントとしてマークします（BWA-MEM 方式）。Picard との互換性に必要です
-a	オフ	フラグ	シングルエンドリードのすべてのアラインメントを出力します（シーディング用。ファイルサイズが増加します）
-p	オフ	フラグ	入力をインターリーブされたペアエンド FASTQ として扱います
-c	500	100–10000	MEM カウントがしきい値を超える場合、アラインメントをスキップします（マルチマッパーのノイズを低減します）
-T	30	20–60	アラインメントスコアの最小しきい値。値を小さくすると、低品質のアラインメントがより多く報告されます
-Y	オフ	フラグ	補足アラインメントにソフトクリッピングを使用します（GATK で推奨）

よく使うレシピ

レシピ 1：複数のサンプルを一括でアラインメントする

#!/bin/bash
# Align all samples in parallel using GNU parallel or sequential loop
GENOME="GRCh38.fa"
SAMPLES=(ctrl_1 ctrl_2 treat_1 treat_2)
THREADS=12

for sample in "${SAMPLES[@]}"; do
    echo "=== Aligning $sample ==="
    bwa-mem2 mem -t $THREADS \
        -R "@RG\tID:${sample}\tSM:${sample}\tPL:ILLUMINA\tLB:lib1\tPU:run1" \
        $GENOME \
        data/${sample}_R1.fastq.gz \
        data/${sample}_R2.fastq.gz \
        | samtools sort -@ 4 -m 2G -o results/${sample}.sorted.bam
    samtools index results/${sample}.sorted.bam
    
    MAPPED=$(samtools view -c -F 4 results/${sample}.sorted.bam)
    TOTAL=$(samtools view -c results/${sample}.sorted.bam)
    echo "$sample: $MAPPED / $TOTAL reads mapped"
done

レシピ 2：Python でアラインメント指標を解析する

import subprocess
import pandas as pd
from pathlib import Path

samples = ["ctrl_1", "ctrl_2", "treat_1", "treat_2"]
metrics = []

for sample in samples:
    bam = f"results/{sample}.sorted.bam"
    result = subprocess.run(
        ["samtools", "flagstat", bam], capture_output=True, text=True
    )
    stats = {}
    for line in result.stdout.splitlines():
        if "total" in line:
            stats["total"] = int(line.split()[0])
        elif "mapped" in line and "%" in line:
            stats["mapped"] = int(line.split()[0])
            stats["pct_mapped"] = float(line.split("(")[1].split("%")[0])
        elif "properly paired" in line:
            stats["properly_paired"] = int(line.split()[0])
    stats["sample"] = sample
    metrics.append(stats)

df = pd.DataFrame(metrics).set_index("sample")
print(df[["total", "mapped", "pct_mapped", "properly_paired"]])
df.to_csv("alignment_metrics.tsv", sep="\t")

期待される出力

出力	形式	説明
*.sorted.bam	BAM	座標順にソートされたアラインメント済みリードです。samtools index でインデックスを作成します
*.sorted.bam.bai	BAI	BAM インデックスです。GATK と IGV でランダムアクセスするために必要です
*.flagstat.txt	テキスト	アラインメントの概要です。合計数、マッピング数、ペア数、重複数とその割合を示します
*.markdup.bam	BAM	重複マーク済み BAM です。GATK HaplotypeCaller への入力として使用します
*.dupmetrics.txt	テキスト	推定ライブラリサイズを含む Picard の重複メトリクスです

トラブルシューティング

問題	原因	解決策
マッピング率が低い（< 85%）	ゲノムの不一致、コンタミネーション、または品質の低いリード	ゲノムアセンブリがサンプルと一致することを確認します。FastQC を実行し、Trim Galore でアダプターをトリミングします
GATK で @RG ヘッダーが見つからないエラー	アラインメント時に -R フラグを指定していない	-R "@RG\tID:...\tSM:...\tPL:ILLUMINA" を指定して再アラインメントします
インデックス作成中にメモリ不足になる	ゲノムインデックスに対して RAM が不足している	BWA-MEM2 はヒトゲノムで約 28 GB を必要とします。必要に応じて、より少ない RAM で使用できる従来の bwa を使います
ペアエンドのリード数が一致しない	インターリーブされた FASTQ、またはファイルの不一致	インターリーブの場合は -p を追加します。zcat r1.fq.gz | wc -l で R1/R2 のリード数を確認します
[E::bwa_idx_load_from_disk] エラー	インデックスファイルがない、またはプレフィックスが間違っている	bwa-mem2 index genome.fa を再実行します。すべての .0123 および .bwt.2bit.64 ファイルが存在することを確認します
アラインメント速度が遅い	スレッド数が少ない、または I/O が遅い	-t 16 以上を使用します。データを SSD に保存し、samtools sort に直接パイプします
補助アラインメントが問題を引き起こす	下流ツールでリードが分割される	分割リードをセカンダリとしてマークするため、-M フラグを追加します（Picard との互換性）
GATK の塩基品質スコア再較正に失敗する	既知のバリアント VCF がない	NCBI または GATK resource bundle から、使用するゲノムアセンブリ用の dbSNP VCF をダウンロードします

参考資料

• BWA-MEM2 GitHub: bwa-mem2/bwa-mem2 — ソースコード、コンパイル済みバイナリ、ベンチマーク
• Vasimuddin M et al. (2019) "Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems" — IPDPS 2019. DOI:10.1109/IPDPS.2019.00041
• Li H & Durbin R (2009) "Fast and accurate short read alignment with Burrows-Wheeler Aligner" — Bioinformatics 25(14):1754-1760. DOI:10.1093/bioinformatics/btp324
• GATK Best Practices for germline short variant discovery — BWA-MEM2 によるアラインメントを推奨する GATK ワークフロー