Skip to content

Add GATK4 components - #226

Open
lazappi wants to merge 19 commits into
mainfrom
feature/add-gatk4
Open

Add GATK4 components#226
lazappi wants to merge 19 commits into
mainfrom
feature/add-gatk4

Conversation

@lazappi

@lazappi lazappi commented Jul 24, 2026

Copy link
Copy Markdown
Contributor

Description

Add components for functions in the Genome Analysis Toolkit (GATK4):

  • gatk4/gatk4_createsequencedictionary: Create a sequence dictionary (.dict) for a reference FASTA file
  • gatk4/gatk4_markduplicates: Locate and tag duplicate reads in a BAM file
  • gatk4/gatk4_baserecalibrator: Generate a base quality score recalibration (BQSR) table from known variant sites
  • gatk4/gatk4_applybqsr: Apply a base quality score recalibration table to a BAM file
  • gatk4/gatk4_haplotypecaller: Call germline SNVs and indels from a BAM file via local re-assembly of haplotypes
  • gatk4/gatk4_combinegvcfs: Combine per-sample GVCFs into a single multi-sample GVCF
  • gatk4/gatk4_genomicsdbimport: Import per-sample GVCFs into a GenomicsDB workspace for scalable joint genotyping
  • gatk4/gatk4_genotypegvcfs: Perform joint genotyping on GVCFs, combined GVCFs, or a GenomicsDB workspace
  • gatk4/gatk4_variantfiltration: Filter variant calls based on INFO and/or FORMAT annotations
  • gatk4/gatk4_selectvariants: Select a subset of variants from a VCF based on various criteria

Checklist before requesting a review

  • I have performed a self-review of my code

  • Conforms to the Contributing guidelines

  • Proposed changes are described in the CHANGELOG.md

  • I have tested my code with viash ns test --parallel -q <name or namespace>

  • Check the correct box. Does this PR contain:

    • Breaking changes
    • New functionality
    • Major changes
    • Minor changes
    • Documentation
    • Bug fixes

@lazappi
lazappi marked this pull request as ready for review July 24, 2026 12:28
@lazappi
lazappi requested a review from mdgrv July 27, 2026 12:39
type: boolean_true
description: |
Round quantized quality scores down to the nearest quantized value.
Cannot be used with `--quantization_levels`.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

argument changed names:
The real mutually-exclusive partner (per the GATK help dump) is this component's --quantize_quals.

example: ["lowDP"]
- name: "--invalidate_previous_filters"
type: boolean_true
description: |

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rm whitespace after |

log "Validating TEST 2 outputs..."
check_file_exists "$meta_temp_dir/output.g.vcf" "output GVCF file"
check_file_not_empty "$meta_temp_dir/output.g.vcf" "output GVCF file"
check_file_contains "$meta_temp_dir/output.g.vcf" "^##fileformat=VCF" "output GVCF file header"

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Every VCF (VCF, GVCF) matches the "^##fileformat=VCF" header line. Maybe use the "##GVCFBlock header lines to check?


staged_reference=$(stage_reference_trio "$tmp_dir" "$par_reference" "$par_reference_fai" "$par_reference_dict")

# GATK requires a literal `gendb://` prefix to recognize a GenomicsDB

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test cases use GVCF, the gendb:// is not included in the genotypegvfs tests. Did a real run with gendb:// and all runs fine, there is just no test to catch future breaks.

check_file_not_empty "$test2_dir/dedup.bam" "duplicate-removed output BAM"
check_file_exists "$test2_dir/dedup_metrics.txt" "duplication metrics file"
check_file_not_empty "$test2_dir/dedup_metrics.txt" "duplication metrics file"
check_file_contains "$test2_dir/dedup_metrics.txt" "ESTIMATED_LIBRARY_SIZE" "duplication metrics file"

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add a check to see if with --remove_duplicates the BAM output is lower than the output in test 1 (default, w/o removal of duplicates)? We need to have less read than default run in case there are duplicate reads present.

@mdgrv mdgrv left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I reviewed all 10 components statically, ran the full viash ns test -q gatk4 suite, and exercised the entire pipeline end-to-end on real data (SEQC2 HCC1395 chr21, ~30x).

Checked:

  • Static (all 10): correct viash 0.9.x conventions, DRY utils/ layer, meta_* (never par_*), JVM -Xmx from meta memory, multiple:true via split_multiple_to_flags, helper-based staging, genuine 4.6.2.0 help dumps, CHANGELOG cites PR #226, ASCII-clean.
  • Tests: all 10 build + test.sh pass; ASCII build gate green for all 10.
  • Real-data end-to-end: the full chain (markduplicates → baserecalibrator on real bgzipped dbSNP138 → applybqsr → haplotypecaller GVCF → combinegvcfs/genomicsdbimport → genotypegvcfs from both gendb:// and combined GVCF → variantfiltration → selectvariants), every step succeeded. Highlights: gendb vs combined give identical genotypes; paired filter expr/name map correctly; SNP selection works.

Suggestions (up to you, see inline comments):

  • Test-coverage gaps — haplotypecaller GVCF test only checks the generic VCF header (assert <NON_REF>/##GVCFBlock); genotypegvcfs gendb:// branch untested; markduplicates --remove_duplicates doesn't assert dropped reads.
  • Doc edits — applybqsr --round_down_quantized references a nonexistent --quantization_levels (should be --quantize_quals); variantfiltration trailing whitespace after |.

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

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants