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 倍，且结果完全相同。它输出 SAM 格式，并包含 GATK HaplotypeCaller 和 Picard 工具所需的读组（@RG）标头。对于双端读段，它会标记正确配对的读段，并将嵌合／拆分读段解析为补充比对。

适用场景

• 将 WGS 或 WES Illumina 读段比对到参考基因组，以进行变异检测（SNP、indel、SV）
• 比对 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 项目中运行时，请使用 pixi run bwa-mem2 调用该工具，而不要直接使用 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 索引

为参考基因组建立索引——每个基因组只需执行一次，人类基因组约需 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：比对双端读段

比对 FASTQ 读段，并添加 GATK 兼容性所需的读组标头。

# 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	Text	比对摘要：总数、已比对、成对和重复读段的数量及百分比
*.markdup.bam	BAM	标记重复读段的 BAM；用作 GATK HaplotypeCaller 的输入
*.dupmetrics.txt	Text	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 资源包下载适用于你的基因组组装的 dbSNP VCF

参考资料

• BWA-MEM2 GitHub：bwa-mem2/bwa-mem2 — 源代码、预编译二进制文件和基准测试
• Vasimuddin M 等（2019）《面向多核系统的 BWA-MEM 高效架构感知加速》— IPDPS 2019。DOI:10.1109/IPDPS.2019.00041
• Li H 和 Durbin R（2009）《使用 Burrows-Wheeler Aligner 快速准确地比对短读段》— Bioinformatics 25(14):1754-1760。DOI:10.1093/bioinformatics/btp324
• GATK 胚系短变异发现最佳实践 — 推荐使用 BWA-MEM2 进行比对的 GATK 工作流