BaizeGen 使用手册 手册版本 2026-11-09

示例与技巧#

完整的 DNAseq 流程脚本#

下面的脚本对一个人类全基因组样本运行完整的胚系流程（GPU 版）：比对、排序、标记重复，同时统计 BQSR 校正表和质控指标，然后输出 gVCF 和 VCF。

#!/bin/bash
set -euo pipefail

# ---------- 需要修改的部分 ----------
BAIZEGEN=/opt/baizegen/baizegen-gpu
REF=/data/ref/Homo_sapiens_assembly38.fasta          # 同目录下有 .fai .dict 和 bwa 索引
DBSNP=/data/ref/Homo_sapiens_assembly38.dbsnp138.vcf.gz
MILLS=/data/ref/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
SAMPLE=NA12878
FQ1=/data/fastq/${SAMPLE}_R1.fq.gz
FQ2=/data/fastq/${SAMPLE}_R2.fq.gz
OUT=/data/out/${SAMPLE}
NT=$(nproc)
# ------------------------------------

mkdir -p "$OUT"
cd "$OUT"

# 1. 比对 + 排序 + 标记重复，同时统计 BQSR 校正表和质控指标
/usr/bin/time -v "$BAIZEGEN" fastq2bam --engine bwa-gpu -R "$REF" -t "$NT" -O "$SAMPLE.bam" \
    --recal-table "$SAMPLE.recal.table" --known-sites "$DBSNP" --known-sites "$MILLS" \
    --qc "$SAMPLE" \
    --rg "@RG\tID:$SAMPLE\tSM:$SAMPLE\tPL:ILLUMINA\tLB:lib1" \
    "$FQ1" "$FQ2" -- -K 100000000 \
    2> "$SAMPLE.fastq2bam.log"

# 2. 变异检测：gVCF（用于联合分型）和 VCF，即时应用 BQSR
/usr/bin/time -v "$BAIZEGEN" haplotyper --device gpu -R "$REF" -I "$SAMPLE.bam" -t "$NT" \
    -q "$SAMPLE.recal.table" --emit-mode gvcf -O "$SAMPLE.g.vcf.gz" \
    2> "$SAMPLE.gvcf.log"
/usr/bin/time -v "$BAIZEGEN" haplotyper --device gpu -R "$REF" -I "$SAMPLE.bam" -t "$NT" \
    -q "$SAMPLE.recal.table" -O "$SAMPLE.vcf.gz" \
    2> "$SAMPLE.vcf.log"

在没有 GPU 的机器上，把 BAIZEGEN 改为 baizegen，--engine bwa-gpu 改为 --engine bwa，去掉两处 --device gpu。

从 GATK 流程迁移#

原流程中的命令 BaizeGen 命令
bwa index ref.fa baizegen index --engine bwa ref.fa（已有 bwa 索引可直接使用）
bwa mem -K N -R RG ref.fa r1 r2 | samtools sort + gatk MarkDuplicates baizegen fastq2bam --engine bwa -R ref.fa --rg RG r1 r2 -- -K N（bwa 的 -R RG 改用 --rg RG；BaizeGen 的 -R 是参考基因组）
gatk SortSam baizegen sort（排序顺序与 samtools sort 相同，同一位置上的记录顺序可能与 SortSam 不同）
gatk MarkDuplicates baizegen dedup
gatk BaseRecalibrator baizegen bqsr，或 fastq2bam --recal-table
gatk ApplyBQSR baizegen applybqsr，或 haplotyper -q（不生成新 BAM）
gatk HaplotypeCaller baizegen haplotyper
gatk HaplotypeCaller -ERC GVCF baizegen haplotyper --emit-mode gvcf（也接受 -ERC GVCF）
gatk CombineGVCFs baizegen combinegvcfs
gatk GenotypeGVCFs baizegen gvcftyper（参考基因组用小写 -r，输入用小写 -v），或 combinegvcfs --genotype-out
gatk Mutect2 / FilterMutectCalls baizegen mutect2 / filtermutectcalls
gatk CollectWgsMetrics 等 baizegen qcmetrics，或 fastq2bam --qc
GATK 3 RealignerTargetCreator + IndelRealigner baizegen realign，或 fastq2bam --realign
fastp baizegen fastp，或 fastq2bam --fastp

复现已有的 GATK 结果#

要让 haplotyper 的输出与某一次实际的 GATK HaplotypeCaller 运行逐条相同，需要对齐以下几项：

  1. 相同的 BAM。不同的比对结果当然会得到不同的变异；从 FASTQ 开始复现时，bwa 的 -K 要与原流程实际使用的批大小一致（见 bwa 的批大小）。
  2. 相同的 GATK 版本：--gatk-version，例如 --gatk-version 4.1.8.0。
  3. 相同的 GATK 参数，例如原流程用了 --max-reads-per-alignment-start 0，这里也要加上。
  4. GATK 当时所在的环境：
    • GATK 运行在支持 AVX-512 的 CPU 上：加 --gatk-host-simd avx512（GATK 的原生 PairHMM 在 AVX-512 与 AVX2 上有极少量浮点差异）。
    • GATK 运行在 glibc 2.41 及以后的系统上：加 --gatk-host-glibc 2.43。
    • GATK 使用的是 Java PairHMM（--pair-hmm-implementation LOGLESS_CACHING）：这里也加同样的参数（仅 CPU）。

比较 VCF 时去掉记录命令行等信息的表头行（各类文件的比较方法见结果一致性验证）：

zcat baizegen.vcf.gz | grep -v '^##' | md5sum
zcat gatk.vcf.gz     | grep -v '^##' | md5sum

示例：复现 GATK 4.1.8.0、参数 --max-reads-per-alignment-start 0 的 VCF：

baizegen-gpu haplotyper --device gpu -R ref.fa -I sample.bam -t 32 \
    --gatk-version 4.1.8.0 --max-reads-per-alignment-start 0 -O sample.vcf.gz

外显子与 panel 数据#

  • 质控：在 fastq2bam 中加 --qc PREFIX --qc-bait-intervals baits.bed --qc-target-intervals targets.bed，得到 CollectHsMetrics 的结果。

  • 变异检测：用 -L 限定在捕获区域。GATK 常用 -ip 100 向两侧扩展区间，haplotyper 没有这个参数，可以先扩展 BED：

    bash bedtools slop -i targets.bed -g ref.fa.fai -b 100 | sort -k1,1 -k2,2n | bedtools merge > targets.pad100.bed baizegen-gpu haplotyper --device gpu -R ref.fa -I sample.bam -L targets.pad100.bed -t 32 -O sample.vcf.gz

多 GPU 与选择显卡#

  • 默认使用全部可见、且未被其他进程占用的显卡；结果与卡数无关。
  • 只用某几张卡：CUDA_VISIBLE_DEVICES=0,1 baizegen-gpu …，或用 --gpus N 指定张数。

  • 一台多卡机器同时处理多个样本：每个进程绑定不同的卡，并分配一部分 CPU 线程：

    bash CUDA_VISIBLE_DEVICES=0 baizegen-gpu haplotyper --device gpu -t 24 -R ref.fa -I s1.bam -O s1.vcf.gz & CUDA_VISIBLE_DEVICES=1 baizegen-gpu haplotyper --device gpu -t 24 -R ref.fa -I s2.bam -O s2.vcf.gz & wait

内存受限的机器#

fastq2bam 默认把排序和标记重复的中间数据放在内存中，人类 WGS 的峰值内存约 65–85 GB。内存不够时：

baizegen-gpu fastq2bam --engine bwa-gpu -R ref.fa -t 32 -O sample.bam \
    -m 2g --dedup-max-mem 4g -T /nvme/tmp/sample \
    --run-store disk \
    --rg '@RG\tID:s\tSM:s\tPL:ILLUMINA\tLB:lib1' \
    r1.fq.gz r2.fq.gz -- -K 100000000
  • -m：排序的内存预算，超过后溢出到 -T 指定的位置。
  • --dedup-max-mem：标记重复的内存上限。
  • --run-store disk：GPU 版排序中间数据写到固态盘（内存约减半）。
  • 以上设置都不改变输出。

其他命令：mutect2 减小 --chunk-size 或线程数；combinegvcfs 在主机内存紧张时会等待而不是被系统终止。

CRAM 输入与输出#

  • 所有命令都可以读 CRAM，用 -R 给出的参考基因组解码。
  • fastq2bam -O sample.cram 直接写 CRAM 3.1（人类 WGS 约为 BAM 大小的一半），读取时需要同一个参考基因组。
  • CRAM 中不保存 MD / NM 标签，samtools / htslib 读取时会自动重建；GATK（htsjdk）读取时不重建这两个标签。

非人类物种与多倍体#

  • 倍性：haplotyper --ploidy N；gvcftyper 按每个样本 gVCF 中 GT 的倍性分型，不同倍性的样本可以一起分型。
  • 长染色体（超过 512 Mb，例如小麦）：BAM 和 VCF 的索引自动改用 .csi。
  • 大基因组的 GPU 比对：bwa 索引需要整体放入显存（小麦约 17 GB，需要 24 GB 的卡），否则自动退回 CPU 比对。
  • gVCF 表头中染色体顺序与参考基因组不同不影响 combinegvcfs（与 GATK 相同）；同名染色体长度不一致时才需要加 --disable-sequence-dictionary-validation。

记录运行时间和内存#

用 /usr/bin/time -v 包住命令，日志末尾会给出墙钟时间（Elapsed (wall clock) time）和峰值内存（Maximum resident set size）：

/usr/bin/time -v baizegen-gpu haplotyper --device gpu -R ref.fa -I sample.bam -t 32 -O sample.vcf.gz 2> hc.log

BaizeGen 使用手册 · 手册版本 2026-11-09。命令行参数以所用版本的 baizegen <命令> --help 为准。