Skip to content

hap.py benchmarking silently produces empty results for GRCh37/b37 (non-chr-prefixed) references #324

Description

@mbiob

Description of the bug

conf/modules.config hardcodes --fixchr for the HAPPY_PREPY process:

withName: "HAPPY_PREPY" {
    ext.args   = {"--fixchr --filter-nonref --bcftools-norm"}
}

hap.py 0.3.15's pre.py --fixchr unconditionally force-prepends chr to
chromosome names before running bcftools norm against --reference. For a
GRCh37/b37-style reference (contigs named 1, 2, ... with no chr prefix
— e.g. human_g1k_v37.fasta), this makes bcftools norm look for chr1 in
the fasta index, fail, and get silently swallowed. The net effect: pre.py
emits an essentially empty VCF (header only), and every downstream hap.py
summary metric reports TRUTH.TP=0 / QUERY.TOTAL=0 across the board — with no
error surfaced to the user. rtg-tools (run in the same invocation) is
unaffected since it doesn't go through this code path, so the failure is easy
to miss if you're only glancing at one benchmarking tool's output.

I confirmed this by running the generated pre.py command manually in the
quay.io/biocontainers/hap.py:0.3.15--py27hcb73b3d_0 container:

  • With --fixchr (as hardcoded): 0 variants in output, warning
    Empty input or not tabix indexed
  • With --fixchr removed, everything else identical: 30,140 variants in
    output (correct, non-empty result)

I see #313 / #316 already removed HAPPY_PREPY entirely on dev (routing
regions into hap.py's own -f/-T instead, which doesn't have this issue) —
but that fix isn't in any tagged release yet. master and the latest release
(1.5.0) both still have the hardcoded --fixchr, so this silently breaks
hap.py benchmarking today for anyone using a chr-less reference (b37/hg19
without chr, which is still common, e.g. GATK's legacy b37 bundle). Since
the failure is silent rather than a hard error, it seems worth a callout
(and/or a backport/patch release) rather than only relying on the dev-branch
rewrite, since anyone on a released version could be misreading empty results
as "this caller has 0% recall" rather than "the tool didn't run".

Command used and terminal output

NXF_VER=25.04.0 nextflow run nf-core/variantbenchmarking -r 1.5.0 -profile docker \
  -params-file params_HG001_small.yaml

(params: analysis=germline, variant_type=small, method=happy,rtgtools,
preprocess=normalize,deduplicate,prepy, fasta/fai=local b37 human_g1k_v37,
truth_vcf=GIAB NISTv4.2.1 HG001 GRCh37, regions_bed=GIAB high-confidence
BED, targets_bed=exome capture BED)

hap.py summary.csv output (both test callers identical pattern):

Type,Filter,TRUTH.TOTAL,TRUTH.TP,TRUTH.FN,QUERY.TOTAL,QUERY.FP,QUERY.UNK,...
INDEL,ALL,1753,0,1753,0,0,0,...
SNP,ALL,28474,0,28474,0,0,0,...

.command.log from the HAPPY_PREPY task:

[I] Total VCF records:         35413
[I] Non-reference VCF records: 35413
WARNING  Empty input or not tabix indexed

Relevant files

conf/modules.config, HAPPY_PREPY withName block (line ~438-444 in 1.5.0 / master)

System information

  • nf-core/variantbenchmarking 1.5.0 (also present on master @ HEAD)
  • Nextflow 25.04.0
  • Docker, macOS (Apple Silicon, linux/amd64 emulation)
  • Reference: GRCh37/b37 (human_g1k_v37.fasta), no chr prefix

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions