结果一致性验证#

BaizeGen 的每个工具都以一个标准工具为准绳。本章说明怎样自己验证这一点：每个命令对照哪个工具的哪个版本，比较之前要对齐哪些条件，各类输出文件怎样比较，以及结果不同时从哪里查起。

GPU 版与 CPU 版之间、不同线程数和不同卡数之间的比较，方法相同（结果应逐字节相同），另见 GPU 自检。

对照的标准工具与版本#

比较之前#

逐条相同的前提是两边做的是同一件事。比较之前确认以下几点，大部分“不一致”都出在这里：

1. 同一份输入文件：同一个 FASTQ / BAM 文件，而不是同一样本的另一次比对或另一份拷贝；同一个参考基因组 FASTA 和索引（bwa 的 .alt 文件在不在也会改变比对结果）。从不同的 BAM 出发，变异结果当然不同。
2. 同一个引擎：align 和 fastq2bam 的默认引擎是 minibwa，与 bwa 比较时必须加 --engine bwa 或 --engine bwa-gpu。
3. 同样的参数：标准工具命令中的每个参数都原样带到 BaizeGen（BaizeGen 接受 GATK 的写法，例如 -ERC GVCF）；-L 区间也要相同。不同 GATK 版本的默认值不同，用 --gatk-version 选择对应版本。
4. bwa 的批大小：两边给出相同的 -K；都不给 -K 时，批大小随线程数变化，两边的 -t 也要相同，见 bwa 的批大小。
5. 排序工具：BaizeGen 的坐标排序与 samtools sort 相同（同一位置上正链在前，其余保持输入顺序）。Picard SortSam 在同一位置、同一链上再按 read 名排序，顺序可能不同；MarkDuplicates 按文件顺序挑选保留哪条记录，所以下游的重复标记和变异结果也可能随之不同。与用 SortSam 排序的流程比较时，对它的 SortSam 输出运行 baizegen dedup，再与它的 MarkDuplicates 输出比较。
6. GATK 的运行方式：HaplotypeCaller 的 VCF 模式在整个运行中共用一个随机数生成器（下采样和 QD 抖动），按区间拆成多个进程运行再合并的结果，可能与单个进程的结果不同。与这样的结果比较时，BaizeGen 也按同样的区间分别运行。
7. GATK 当时所在的环境：GATK 运行在支持 AVX-512 的 CPU 上、glibc 2.41 及以后的系统上，或使用 Java PairHMM 时，加上对应的参数，见复现已有的 GATK 结果。

比较方法#

总原则：比较解压后的内容，不比较文件本身。BAM、.vcf.gz、.fq.gz 的压缩字节随压缩库和压缩级别变化，内容相同时文件 md5 也可能不同；索引文件（.bai、.tbi、.csi）不必比较。另外，各工具都会在文件中写入自己的程序名、命令行、版本或运行时间，比较时去掉这些行：

SAM / BAM#

适用于比对、排序、标记重复、fastq2bam、applybqsr、realign 和 postalt 的输出：

# 记录：两个 md5 应相同
samtools view -@ 8 BAIZEGEN.bam | md5sum
samtools view -@ 8 STANDARD.bam | md5sum

# 表头：没有输出即相同
diff <(samtools view -H BAIZEGEN.bam | grep -E '^@(SQ|RG)') \
     <(samtools view -H STANDARD.bam | grep -E '^@(SQ|RG)')

• align 输出 SAM 时直接比较文本：grep -v '^@' BAIZEGEN.sam | md5sum。
• CRAM 加 -T REFERENCE 解码。CRAM 与 BAM 互相比较时，两边的 MD / NM 标签可能不同（CRAM 不保存这两个标签，读取时重建），建议同一种格式之间比较。
• realign 另外比较区间文件：cmp BAIZEGEN.intervals GATK3.intervals。
• dedup 的输出只比较重复标记时：samtools view -c -f 0x400 FILE 给出两边的重复记录数，samtools view FILE | awk '{print int($2/1024)%2}' | md5sum 比较每条记录的标记。

fastp#

# 输出的 FASTQ：比较解压后的内容
zcat BAIZEGEN_R1.fq.gz | md5sum
zcat FASTP_R1.fq.gz    | md5sum

# JSON 报告：没有输出即相同
diff <(grep -v '"command"' BAIZEGEN.json) <(grep -v '"command"' FASTP.json)

fastp 多线程运行时，重复率估计等部分结果依赖线程之间的先后，每次运行可能不同；BaizeGen 在任何线程数下都等于 fastp -w 1 的结果，所以对照的 fastp 要以 -w 1 运行。fastq2bam --fastp "ARGS" 等于先运行 fastp ARGS -w 1，再对它的输出运行 fastq2bam。

BQSR 校正表#

GATK 的校正表中没有日期、路径和版本，直接比较整个文件：

cmp BAIZEGEN.recal.table GATK.recal.table && echo 相同

两边的 --known-sites 和 -L 要相同。

质控指标#

Picard 的指标文件从 ## METRICS CLASS 一行开始比较；DepthOfCoverage 和 flagstat 的输出比较整个文件：

for m in alignment_summary_metrics insert_size_metrics wgs_metrics quality_yield_metrics duplicate_metrics; do
    if diff -q <(sed -n '/^## METRICS CLASS/,$p' QC_PREFIX.$m) <(sed -n '/^## METRICS CLASS/,$p' PICARD.$m) > /dev/null
    then echo "$m 相同"; else echo "$m 不同"; fi
done
cmp QC_PREFIX.flagstat SAMTOOLS.flagstat
cmp QC_PREFIX_depth.sample_summary GATK_DEPTH.sample_summary

• 对照的 Picard / GATK 用 --gatk-version 所选版本的 jar（默认 4.6.2.0），以默认参数或与 qcmetrics 相同的参数运行；DepthOfCoverage 的统计区域（-L）、覆盖度阈值和分箱与 qcmetrics 的设置相同。
• flagstat 的格式随 samtools 版本不同（1.13 之前没有 primary 相关的行），用 --samtools-version 指定对照的 samtools 版本。
• 各文件对应的标准工具见 qcmetrics。

VCF / gVCF#

适用于 haplotyper、combinegvcfs、gvcftyper、variantfiltration 和 selectvariants：

# 记录（含 #CHROM 行）：两个 md5 应相同
zcat -f BAIZEGEN.vcf.gz | grep -v '^##' | md5sum
zcat -f GATK.vcf.gz     | grep -v '^##' | md5sum

# 表头：没有输出即相同
diff <(zcat -f BAIZEGEN.vcf.gz | grep '^##' | grep -v '^##GATKCommandLine') \
     <(zcat -f GATK.vcf.gz     | grep '^##' | grep -v '^##GATKCommandLine')

combinegvcfs 加 --merge-ties-by-input 时，同一位置记录的合并顺序与 GATK 不同（见加速选项），不与 GATK 比较；这时可以比较不同机器、不同卡数的运行之间是否相同。

甲基化#

bsalign --aligner bowtie2 与 Bismark 比较时，两边使用同一个 bowtie2（--path_to_bowtie2）和同样的选项，输出的文件名也相同：

• BAM：按 SAM / BAM 的方法比较记录；表头除 @PG 行外相同。
• 比对报告：diff <(sed '1d;$d' BAIZEGEN_PE_report.txt) <(sed '1d;$d' BISMARK_PE_report.txt)（去掉第一行和最后一行）。
• 甲基化提取（bsextract）：各 context 的调用文件去掉第一行比较（zcat -f FILE | sed 1d | md5sum）；splitting report 去掉 Bismark Extractor Version 行；其余文件比较解压后的全部内容。
• 报告中写有输入文件路径的（例如去重报告），两边的输入路径相同时才会逐字节相同。

统计一致的输出#

mutect2 等体细胞工具、bsalign --aligner baize 是按算法实现的，与标准工具的结果高度一致，但不保证逐条相同。这类输出比较一致的程度，而不是是否相同。

体细胞 VCF：用 GATK Concordance 以 GATK 的结果为“真值”，统计 SNP 和 indel 的灵敏度与精确度：

gatk Concordance --truth GATK.vcf.gz --evaluation BAIZEGEN.vcf.gz --summary concordance.tsv

或者先拆分多等位位点、左对齐，再看两边各自独有的位点：

for f in GATK BAIZEGEN; do
    bcftools norm -m -any -f REFERENCE $f.vcf.gz -Oz -o $f.norm.vcf.gz && bcftools index -t $f.norm.vcf.gz
done
bcftools isec -p isec_out GATK.norm.vcf.gz BAIZEGEN.norm.vcf.gz
# isec_out/0000.vcf：仅 GATK 有；0001.vcf：仅 BaizeGen 有；0002.vcf / 0003.vcf：两边共有

比对结果：两个比对器对同一批 reads 的输出（按输入顺序，未排序），可以用 samcompare 逐条比较，按 MAPQ 分层统计不同的记录。

甲基化 --aligner baize：比较相同比对的 read 对比例和共有 CpG 位点甲基化水平的相关性，并与 Bismark 用不同随机数种子运行两次之间的差异对照，见甲基化分析。

同一个统计一致的工具，-t N 与 -t 1、GPU 与 CPU 之间仍然逐字节相同，可以用上面逐条相同的方法比较。

结果不同时#

先定位第一处不同。大文件不要直接 diff（全基因组的 SAM 文本有上百 GB），用 cmp 找到第一个不同的行号，再取出那一行：

cmp <(samtools view BAIZEGEN.bam) <(samtools view STANDARD.bam)
# 输出例如 "... differ: byte 123456789, line 4567"
samtools view BAIZEGEN.bam | sed -n '4567p;4567q'
samtools view STANDARD.bam | sed -n '4567p;4567q'

VCF 同理，把 samtools view FILE 换成 zcat -f FILE | grep -v '^##'。有索引时，也可以按染色体或区间（samtools view FILE chr20、tabix FILE chr20）缩小范围。然后依次检查：

1. 比较之前的各项条件，尤其是输入文件、比对引擎、bwa 的 -K、排序工具，以及 GATK 是否拆分运行。
2. 对照工具的版本：bwa 0.7.17 与 0.7.19、不同版本的 GATK、不同版本的 samtools（flagstat）的输出本身就可能不同。
3. 已知的细微差别：GATK 在 AVX-512 与 AVX2 CPU 上、不同 glibc 版本上的少量 PL / QUAL 差异，以及 GATK 4.0.12.0 自身的随机差异，见常见问题。

仍然不同时请反馈给我们，附上两边的完整命令、标准工具的版本、第一处不同的两行，以及反馈时请附上中列出的信息。