胚系变异检测（DNAseq）#

本章介绍单个样本从 FASTQ 到 VCF / gVCF 的标准胚系变异检测流程。流程遵循 GATK 的最佳实践（bwa-mem 比对、排序、标记重复、可选的碱基质量校正，然后 HaplotypeCaller），每一步的输出与对应工具一致。

概述#

所需输入#

• 测序数据：双端（两个文件）或单端（一个文件）FASTQ，可以是 gzip 压缩的。
• 参考基因组 FASTA 及其索引：.fai、.dict、bwa 索引，准备方法见参考基因组准备。
• 可选：已知变异位点 VCF（dbSNP、Mills indels 等），用于 BQSR 和 indel 重比对。

流程步骤#

1. 比对、排序和标记重复：fastq2bam
2. 质控指标（可选）：fastq2bam --qc 或 qcmetrics
3. 碱基质量校正 BQSR（可选）：fastq2bam --recal-table 或 bqsr
4. indel 局部重比对（可选，GATK 3 风格流程）：fastq2bam --realign 或 realign
5. 变异检测：haplotyper，输出 VCF 或 gVCF

第 2–4 步都可以在第 1 步的同一个进程中完成，不产生中间 BAM。

分步说明#

1. 比对、排序和标记重复#

fastq2bam 在一个进程中完成比对（bwa-mem）、坐标排序和标记重复，直接写出已排序、已标记重复的 BAM 及其索引。输出与分别运行 align | sort 和 dedup 逐字节相同。

GPU 版：

baizegen-gpu fastq2bam --engine bwa-gpu -R REFERENCE -t NUMBER_THREADS -O DEDUPED_BAM \
    --rg '@RG\tID:GROUP_NAME\tSM:SAMPLE_NAME\tPL:PLATFORM\tLB:LIBRARY' \
    FASTQ_1 [FASTQ_2] -- -K 100000000

CPU 版：

baizegen fastq2bam --engine bwa -R REFERENCE -t NUMBER_THREADS -O DEDUPED_BAM \
    --rg '@RG\tID:GROUP_NAME\tSM:SAMPLE_NAME\tPL:PLATFORM\tLB:LIBRARY' \
    FASTQ_1 [FASTQ_2] -- -K 100000000

需要的输入：

• REFERENCE：参考基因组 FASTA 文件，同目录下要有 bwa 索引。
• NUMBER_THREADS：线程数，一般设为机器的 CPU 线程数。
• FASTQ_1、FASTQ_2：测序数据。给一个文件为单端，两个文件为双端。同一样本测了多次（多对 FASTQ）时，每一端写成逗号分隔的列表（run1_R1.fq.gz,run2_R1.fq.gz run1_R2.fq.gz,run2_R2.fq.gz），结果与先合并文件再比对相同，见多对 FASTQ。
• GROUP_NAME、SAMPLE_NAME、PLATFORM、LIBRARY：读组（read group）信息，用 --rg 给出（等同于 bwa mem -R）。SM 是样本名，haplotyper 用它作为 VCF 中的样本列名，必须提供；LB 用于标记重复时区分文库。
• -- 之后是 bwa mem 自己的参数（如 -k、-M、-Y、-K），写法与 bwa mem 相同（-x、-t、-o 除外），全部参数见 bwa mem 参数。

输出：

• DEDUPED_BAM：按坐标排序、重复 reads 已打上 0x400 标记的 BAM，以及索引 DEDUPED_BAM.bai。染色体长度超过 512 Mb（例如小麦）时自动改写 .csi 索引。
• 输出文件名以 .cram 结尾时写 CRAM 3.1（及 .crai），读取时需要同一个参考基因组。

bwa 的批大小 -K#

bwa-mem 按批读入 reads，并在每一批内估计插入片段长度分布，因此批大小（-K，每批的碱基数）会影响少量 reads 的比对结果。

• 推荐总是明确给出 -K（例如 -- -K 100000000），这样结果与线程数无关。
• 不给 -K 时，与 bwa 一样，批大小为 1000 万 × 线程数，结果会随 -t 变化。
• 要与已有的 bwa 结果逐条对比时，使用对方实际的批大小。例如对方运行的是不带 -K 的 bwa mem -t 48，这里就写 -- -K 480000000。

分步运行#

需要分别得到各步骤的输出时，可以分步运行，最终 BAM 与 fastq2bam 相同：

baizegen align --engine bwa -R REFERENCE -t NUMBER_THREADS \
    --rg '@RG\tID:GROUP_NAME\tSM:SAMPLE_NAME\tPL:PLATFORM\tLB:LIBRARY' \
    FASTQ_1 [FASTQ_2] -- -K 100000000 \
  | baizegen sort -I - -O SORTED_BAM -t NUMBER_THREADS
baizegen dedup -I SORTED_BAM -O DEDUPED_BAM -t NUMBER_THREADS

• align 不给 -O 时把 SAM 写到标准输出，sort -I - 从标准输入读取。
• dedup 同时写出索引（.bai），无需再运行 samtools index；默认只标记重复；加 --remove-duplicates 则直接删除重复记录。

测序数据预处理（fastp）#

加 --fastp "ARGS" 时，fastq2bam 先在同一进程中用 fastp 1.4.0 过滤和修剪 reads，再比对修剪后的 reads。ARGS 是 fastp 自己的参数（"" 表示使用 fastp 的默认设置），结果等于先运行 fastp ARGS -w 1 再运行 fastq2bam，同时写出 fastp 的 JSON / HTML 报告：

baizegen-gpu fastq2bam --engine bwa-gpu -R REFERENCE -t NUMBER_THREADS -O DEDUPED_BAM \
    --fastp "--detect_adapter_for_pe" --fastp-json SAMPLE.fastp.json --fastp-html SAMPLE.fastp.html \
    --rg '@RG\tID:GROUP_NAME\tSM:SAMPLE_NAME\tPL:PLATFORM\tLB:LIBRARY' \
    FASTQ_1 FASTQ_2 -- -K 100000000

2. 质控指标（可选）#

在 fastq2bam 中加 --qc PREFIX，在写 BAM 的同时统计质控指标，几乎不增加耗时：

baizegen-gpu fastq2bam ... --qc QC_PREFIX

对已有的 BAM 单独统计：

baizegen qcmetrics -I DEDUPED_BAM -R REFERENCE -O QC_PREFIX -t NUMBER_THREADS

输出文件（每个文件的内容与所选 GATK 版本自带的 Picard / GATK，以及所选版本的 samtools 逐字节相同）：

外显子或 panel 数据再加 --qc-bait-intervals BAITS --qc-target-intervals TARGETS（interval_list 或 BED），额外输出 CollectHsMetrics 的结果。参数详见 qcmetrics。

3. 碱基质量校正 BQSR（可选）#

BQSR 分两步：统计校正表（BaseRecalibrator），再按表校正碱基质量（ApplyBQSR）。

统计校正表，二选一：

# 在 fastq2bam 中同时统计（GPU 版在 GPU 上统计）
baizegen-gpu fastq2bam ... --recal-table RECAL_TABLE \
    --known-sites DBSNP --known-sites KNOWN_INDELS

# 或对已有的 BAM 单独统计
baizegen bqsr -I DEDUPED_BAM -R REFERENCE -k DBSNP -k KNOWN_INDELS -O RECAL_TABLE -t NUMBER_THREADS

• DBSNP、KNOWN_INDELS：已知变异位点 VCF，可重复给出多个，至少一个。
• RECAL_TABLE：GATKReport 格式的校正表，与 GATK BaseRecalibrator 的输出逐字节相同。

应用校正表，二选一：

# 推荐：变异检测时即时校正，不生成新的 BAM
baizegen-gpu haplotyper --device gpu -R REFERENCE -I DEDUPED_BAM -q RECAL_TABLE -t NUMBER_THREADS -O VARIANT_VCF

# 或写出校正后的 BAM
baizegen applybqsr -I DEDUPED_BAM -q RECAL_TABLE -O RECALED_BAM -t NUMBER_THREADS

4. indel 局部重比对（可选）#

GATK 4 的 HaplotypeCaller 自带局部组装，一般不需要这一步。需要复现 GATK 3 风格流程时，可以运行 GATK 3.8-1 的 RealignerTargetCreator + IndelRealigner，目标区间和输出 BAM 的每条记录都与 GATK 3.8-1 相同：

# 在 fastq2bam 中，写 BAM 之前直接重比对（需要 -t 2 以上）
baizegen-gpu fastq2bam ... --realign --known-indels KNOWN_INDELS

# 或对已有的 BAM 单独运行
baizegen realign -I DEDUPED_BAM -R REFERENCE -O REALIGNED_BAM --known-indels KNOWN_INDELS -t NUMBER_THREADS

在 fastq2bam 中同时使用 --realign 和 --recal-table / --qc 时，BQSR 和质控统计的是重比对之后的记录。

5. 变异检测#

haplotyper 是 GATK HaplotypeCaller 的实现，默认复现 GATK 4.6.2.0 的结果。

输出 VCF：

baizegen-gpu haplotyper --device gpu -R REFERENCE -I DEDUPED_BAM -t NUMBER_THREADS \
    [-q RECAL_TABLE] -O VARIANT_VCF

输出 gVCF（用于之后的多样本联合分型）：

baizegen-gpu haplotyper --device gpu -R REFERENCE -I DEDUPED_BAM -t NUMBER_THREADS \
    [-q RECAL_TABLE] --emit-mode gvcf -O GVCF

CPU 版把 baizegen-gpu 换成 baizegen 并去掉 --device gpu，结果逐字节相同。

需要的输入：

• REFERENCE：参考基因组 FASTA。
• DEDUPED_BAM：已排序、已标记重复并带索引的 BAM 或 CRAM。
• VARIANT_VCF / GVCF：输出文件。以 .gz 结尾时写 bgzip 压缩的文件和 tabix 索引（.tbi；染色体超过 512 Mb 时为 .csi）。

常用的可选参数：

• -L INTERVAL：只分析指定区间，可以是 chr20:1000000-2000000 这样的区间，也可以是 .bed、.interval_list、.list 文件；可以重复给出。
• --ploidy N：样本倍性，默认 2。
• --gatk-version V：要复现的 GATK 版本，见下文。
• -A StrandBiasBySample：VCF 模式下也输出 FORMAT 中的 SB 字段。
• --gpus N|all：GPU 版使用的显卡数，默认使用全部空闲的卡；多卡时按染色体分配，结果与单卡相同。

一个 BAM 中有多个样本（读组的 SM 不同）时，haplotyper 输出多样本 VCF，与 GATK HaplotypeCaller 的行为一致。

复现指定的 GATK 版本#

--gatk-version 支持 GATK 4 各个主版本线的最后一个版本，输出与该版本的 HaplotypeCaller 逐条记录相同：

要逐字节复现某一次实际的 GATK 运行，还需要使用与那次运行相同的参数，并告诉程序 GATK 当时所在的环境，详见复现已有的 GATK 结果。

参考耗时#

人类 30× WGS（约 6.6–8.1 亿条 reads，2×150 bp），1 × RTX 5090 + Threadripper 3960X（24 核），-t 48，FASTQ / BAM 在本地 NVMe：

耗时随数据、硬件和存储变化很大，仅供参考；从 NFS 直接读写时整体会被网速拖慢。与 bwa、GATK 的对比见性能对比。