GATK4_GATHERPILEUPSUMMARIES concatenates scattered pileup tables out of genome order

[Hope2925]

Description of the bug

Bug report: GATK4_GATHERPILEUPSUMMARIES concatenates scattered pileup tables out of genome order; CalculateContamination crashes / silently produces wrong contamination

Environment

• nf-core/sarek 3.10.0 (commit 8ccac7a)
• Nextflow 25.10.4
• GATK 4.6.2.0, container community.wave.seqera.io/library/gatk4_gcnvkernel:edb12e4f0bf02cd3
• Executor: Google Batch (googlebatch profile), but the mechanism is executor-independent (see root cause)
• Reference: GRCh38 analysis-set (Broad Homo_sapiens_assembly38.fasta/.dict)

What happened:

A tumor/normal WES pair failed at BAM_VARIANT_CALLING_SOMATIC_MUTECT2:CALCULATECONTAMINATION:

Process NFCORE_SAREK:SAREK:BAM_VARIANT_CALLING_SOMATIC_ALL:BAM_VARIANT_CALLING_SOMATIC_MUTECT2:CALCULATECONTAMINATION (sample1_T_vs_sample1_N) terminated with an error exit status (3)
...
java.lang.IllegalArgumentException: Invalid interval. Contig:chr22 start:37986297 end:37650688
at org.broadinstitute.hellbender.utils.Utils.validateArg(Utils.java:804)
at org.broadinstitute.hellbender.utils.SimpleInterval.validatePositions(SimpleIl.java:60)
at org.broadinstitute.hellbender.utils.SimpleInterval.(SimpleInterval.jav
at org.broadinstitute.hellbender.tools.walkers.contamination.ContaminationSegmeambda$findContigSegments$4(ContaminationSegmenter.java:71)
...
at org.broadinstitute.hellbender.tools.walkers.contamination.ContaminationModel>(ContaminationModel.java:61)
at org.broadinstitute.hellbender.tools.walkers.contamination.CalculateContaminaoWork(CalculateContamination.java:124)

CalculateContamination's ContaminationSegmenter seems to assume its input pileup table is sorted in genome order and builds intervals from consecutive rows; a start>end interval means two adjacent the table go backward in genomic position.

Direct evidence: the gathered table is out of order

Pulled the actual sampleID_T.mutect2.pileups.table fed to CalculateContamination (output of GATHERPILEUPSUMMARIES_TUMOR) and looked at contig transitions in file order:

row 15474: chr3 -> chr22
row 16168: chr22 -> chr3
row 63753: chr21 -> chr22
row 64775: chr22 -> chrX

chr22 appears in two disjoint blocks. The reference .dict order is chr3 (record 4), chr21 (22),chr22 (23),chrX (24), so the second chr22 block (between chr21 and chrX) is correctly placed, but the first one (spliced between two chr3 blocks) is not. This produced a negative-length interval when ContaminationSegmenter walks the file assuming monotonic genome order.

Root cause: two independent sarek behaviors that combine badly

1. subworkflows/local/prepare_intervals/main.nf deliberately reorders interval chunks by estimated runtime before handing them to every scatter/gather consumer:

intervals_bed = intervals_bed.flatten()
.map{ intervalFile -> [ duration, intervalFile ] }
.toSortedList({ a, b -> b[0] <=> a[0] }) // longest job first, NOT genome order
.flatten().collate(2).map{ _duration, intervalFile -> intervalFile }.collect()

This is fine for consumers whose merge step re-sorts (e.g. MergeVcfs), but it means every process fed from this channel receives chunks in "biggest first" order, not chromosome order.

2. modules/nf-core/gatk4/gatherpileupsummaries/main.nf concatenates whatever order it's given, with no re-sort and no validation:

def input_list = pileup.collect { pileup_ -> "--I ${pileup_}" }.join(' ')

GATK GatherPileupSummaries is documented/behaves as a straight concatenation tool,it trusts the caller to supply inputs in coordinate order.

Contrast with subworkflows/local/bam_variant_calling_somatic_mutect2/main.nf's VCF path, which merges through GATK4_MERGEVCFS Picard's MergeVcfs does a real coordinate-sorting merge, so it tolerates arbitrary input order. GatherPileupSummaries does not have that tolerance, which is the asymmetry that makes this bug pileup-summary-specific.

Reproduction recipe

• --tools mutect2 with --germline_resource/--germline_resource_tbi set (gnomAD AF-only VCF)
• --nucleotides_per_second set low enough that CREATE_INTERVALS_BED (modules/local/create_intervals_bed/main.nf) splits the WES target intervals into enough chunks (24 in our run) that the runtime-descending resort in prepare_intervals/main.nf interleaves chunks belonging to the same chromosome.
• Any tumor/normal WES pair

Possible fixes?

Either:

• Sort chunks fed into the pileup-summary scatter/gather path back into genome order before GATHERPILEUPSUMMARIES (independent of the runtime-based scatter order used elsewhere), or
• Make GATK4_GATHERPILEUPSUMMARIES's wrapper sort its --I inputs by each file's first data row's (contig rank via .dict, position) before invoking GATK

Command used and terminal output

Relevant files

No response

System information

No response