Read the merged VCF directly to build Variant record tables - #21
Merged
Conversation
The published merged.ann.vcf.gz declares only CANN, CA and DFS - no GT, DP, AD/RO/QR/AO/QA or contig lines - and every INFO field opens with a bare '.'. This spike settles what htsjdk does with that before any code depends on it. htsjdk already resolves to 5.0.0 via the parent gus-project-pom:1.0.0 (dependency:tree on WSFPlugin shows com.github.samtools:htsjdk:jar:5.0.0:compile with no local pin), so no pom.xml change was needed. The spike passed on cedar against the real file - no lenient-codec workaround required. DP and CA read correctly despite being absent from the header, and CANN reads correctly despite the header's dangling leading '.' in INFO. GT is asserted against the specific expected value (GenotypeType.HOM_REF, matching the independently-verified 0:178:... call at this locus) rather than merely non-null: Genotype.getType() never returns null - a degraded parse falls back to NO_CALL (alleles present but all '.') or UNAVAILABLE (allele list empty), both non-null, so only an exact-value check catches either failure mode. That said, GT was never actually at risk from the header gap this spike targets: htsjdk's AbstractVCFCodec.createGenotypeMap() special-cases the GT key by literal match against the record's own FORMAT column, independent of any header FORMAT declaration for GT. So this test's genuine header-tolerance finding is about DP, CANN and CA, not GT - the GT assertion is a correctness guard on the sample call, included for rigor, not evidence that the GT-specific hazard was ever reachable. Skipped unless -Dvariants.vcf names the file, which lives only on the webserver.
CANN is an INFO field; CA is the per-sample FORMAT field naming which CANN entries describe that sample's genotype. Slots split on / or |, transcript keys within a slot on ';', bare 'r' means reference-with-no-CDS-annotation. Exceeds the plan's 8-test floor by two, closing coverage gaps found in review: the truncated-entry defensive branch in parse() and the '.' amino-acid filter in aminoAcidsFor() were otherwise never exercised. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
One row per sample regardless of ploidy - a diploid het collapses to an IUPAC ambiguity code, matching processSequenceVariations.jl's gt_to_base. Read frequency reports support for the CALLED allele, so a reference call shows RO/(RO+AO_total) rather than the pipeline's alt-centric 0.0, and a coverage-filled call (GT and DP only, every other FORMAT field '.') shows 100. The alt branch indexes the AO entry for the allele actually carried, by VCF-wide allele index, rather than summing all AO - summing would dilute the frequency with reads for an alt this sample doesn't carry on a multi-allelic record. The published file is biallelic by construction (write_vcf_entry splits multi-allelic sites into one record per alt: measured 40,000 loci / 548,040 alt calls, zero multi-alt, zero samples with >1 nonzero AO), so this is defensive rather than hot; synthetic multi-alt fixture records and tests exercise it directly, including both the out-of-range fallback and the in-range-but-sentinel fallback (AO correctly sized for the alt count, but '.' sits exactly at the called slot). intsOf preserves per-slot position for a multi-value FORMAT field instead of dropping '.' entries: dropping would shift later slots left and misalign an index-based lookup (a sample called for alt2 could silently read alt3's count). A '.'/unparseable slot becomes a MISSING (-1) sentinel; sums skip it, and the indexed lookup treats a sentinel at the called index as missing rather than returning it as a negative read count. A field that is entirely missing (bare '.', or every comma slot blank) still collapses to an empty array, so the coverage-filled 'no AO/RO data at all' check - which tests array length, not sentinel content - is unaffected. The fixture reproduces the real file's under-declared header on purpose. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The locus query returns sequence, position, name_for_filenames and the EDA sample table suffix in one round trip; verified against genomicsdb_071n. Metadata is keyed on provider_label ['country'], not the VAR_<hash> stable_id, which is a digest of the label and therefore site-specific. SampleMetadataLookup distinguishes two edaSuffix states rather than one: a null/blank suffix (no dnaseq EDA study for this organism) yields an empty map, while a non-null suffix that fails the identifier guard is a hard failure. The SQLException catch is narrowed to SQLState 42P01 (undefined_table, walking the next-exception/cause chain for wrapped drivers) so a missing attributevalue_* table still renders an empty table, but a genuine connection/pool/timeout fault is no longer swallowed into the same silent empty-map outcome. That walk carries an identity-based visited set, since a mutual getCause() cycle or a cyclic getNextException() chain (real pooled drivers produce the latter) would otherwise recurse forever. VariantLocusResolver.vcfPath now validates WEBSERVICEMIRROR and the build number before building the path, rather than silently concatenating a null into "build-null" and failing later as an inscrutable file-not-found.
No-call samples get a row rather than being dropped - the record already advertises no_call_strain_count, so omitting them would contradict the panel above the table. A sample resolving to several transcripts' amino acids stays on one row with the values comma-joined.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every row of the new Strains / Samples table rendered an empty Amino Acid, at loci where CANN was populated. Found by verifying the live page against the real VCF, not by any test. CANN is declared Number=. in the header, so htsjdk parses it into a List<String> - one element per comma-separated entry. Calling getAttributeAsString on a List-valued attribute returns the LIST'S toString: bracketed and ", "-separated, "[r0|TGT|C|..., k0|AGT|S|...]". CannIndex then split that on ',', producing a first key of "[r0" and later keys with a leading space, so no key ever matched a CA value and aminoAcidsFor returned empty for every sample. getAttributeAsStringList is NOT the fix on its own: for a genuine List it returns the elements unchanged, but for a String-valued attribute CommonInfo .getAttributeAsList wraps the whole string as a singleton without splitting on ','. That would hand CannIndex several entries glued together and corrupt the '|'-split fields at the join. So parseCann branches on the raw attribute's runtime type instead of assuming either shape. Why no test caught it: the Task 1 spike asserted only that CANN was non-null, which held either way; no MergedVcfReaderTest case asserted on aminoAcids; and the composer test builds SampleCall by hand. Nothing covered the reader -> CannIndex seam. Tests added at that seam, plus one pinning the bracketed-list form as garbage so nobody later "fixes" it by teaching the parser to tolerate brackets rather than fixing the caller. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The dnaseq pipeline's write_vcf_entry emits one VCF record per unique ALT, so a locus routinely has several records at the same CHROM/POS. MergedVcfReader.read() took only the first record and stopped, silently rendering every sample whose alt lived on a later record as reference. Measured against the real merged.ann.vcf.gz: - 9,384 of 47,189 positions (19.9%) have more than one record; one had seven. - 13.9% of sample-locus alt calls involve a sample carrying alts on more than one record. - Concrete example, Pf3D7_01_v3:538376 (7 records, same REF, different indel ALTs): 43 samples split 22/9/3/2/3/3/1 across the records. Reading only record 1 rendered 21 of those 43 as reference. read() now consumes every record at the position and merges them per sample in MergedVcfReader.mergeSample(): - noCall only if the genotype is no-call on every record. - a "contributing" record is one where the genotype is present and NOT hom-ref - a hom-ref call on a record whose alt the sample doesn't carry is the pipeline's per-ALT splitting artifact, not a reference call. - chromosomeAlleles (the aggregation weight): no contributing record but a real call -> the first real record's alleles; exactly one contributing record -> that record's alleles as-is; more than one -> the non-ref allele from each contributing record. - aminoAcids: union across ALL records where the sample's CA resolves, each against its OWN record's CannIndex - CANN keys (r0, k0, ...) are scoped per record and must never be merged into one map, or a sample's CA would resolve against the wrong record's key. - depth/readFrequency: from a contributing record when one exists (the first, if several), else the first record with a real call. - coverageFilled: only when there is no contributing record and the reference call is coverage-filled by the existing rule. - genotype: raw GT when exactly one record contributed (or plain reference/no-call), else the slash-joined raw GTs of the contributing records - no longer shown on the page, but still used by downloads. - ploidy: max across records where the sample has a real call. LocusCalls.cann is removed: a single merged CannIndex cannot exist now that keys are per-record, and the field was never read by anything. altAlleles is now the union of ALTs across records, in record order. Allele display is unified (confirmed with the user) across the within-record het case and the new across-record case: 1 distinct chromosome allele -> that allele; 2 distinct single-character alleles -> IUPAC ambiguity code (unchanged); otherwise -> slash-joined distinct alleles (e.g. "ATT/ATATT"). This changes MergedVcfReaderTest#complexHetSlashJoinsDistinctAllelesNotIupac (previously ...ReturnsFirstNonRefAlleleNotIupac): a single-record het of A/AG now displays "A/AG" instead of dropping the reference and showing "AG" alone. Every other single-record-locus test is unchanged. Extended the fixture with a two-record locus at chr1:700 (REF=A, ALT=AT and ALT=ATT, each with its own CANN whose k0 maps to a different amino acid) and added tests for: an alt on the second record rendering instead of reference; alts split across two records rendering slash-joined with both in chromosomeAlleles; amino acids unioned and resolved per-record (proving no key collision); hom-ref on every record staying a reference call; no-call on every record staying a no-call. Also, since GT is no longer displayed (hidden via variantRecords.xml in the ApiCommonModel commit), a no-call sample must show "No call" in the Allele column instead of blank - VariantLocusComposer.strainRows and its test are updated accordingly. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Ploidy-weighted per processSequenceVariations.jl: one unit per chromosome slot, so a diploid het contributes to both alleles and the denominator is the sum of ploidies rather than a sample count. Haploid P. falciparum makes the two coincide; TriTryp and Fungi aneuploids do not. The denominator is each country's own chromosomes, so a row reads on its own. Samples with no country are excluded entirely. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Equal strain counts fell back to the VCF's sample order - deterministic per file but arbitrary to a reader. Sort is now (-strainCount, country). Two tests passed vacuously. The sort-order test grouped Mali first, so LinkedHashMap insertion order alone satisfied it; the calls are reordered so Gambia is encountered first and both positions are asserted. Verified by deleting the sort: exactly those two ordering tests go red. Also covers a real rank-3 allele, which only had its empty-string form pinned. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The two row builders each had their own idea of 'has a country' - getOrDefault
in one, a null/empty check in the other. They now share countryOf(), which also
trims: SampleMetadataLookup's SQL rejects NULL but not blank, and nothing
normalises EDA free text, so ' ' was becoming its own country group and
'Mali ' was splitting Mali across two rows.
CountryRow's allele fields hold rendered cells ('C (0.3333)'), not the bare
allele that SampleCall.allele and StrainRow.allele hold two files away, so they
are renamed *WithFrequency and the record now documents the mixed units:
strainCount counts samples, the frequencies count chromosomes.
Ranks beyond 3 are dropped for a three-column table; that is now commented and
covered by a test rather than left for a reader to discover.
Verified by mutation: removing the trim reddens exactly the three
normalisation tests.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Also normalise country values at the SampleMetadataLookup boundary so a blank or untrimmed EDA value cannot split one country into two rows for any future consumer of the map. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The two VCF-backed Variant record-page tables duplicated ~45 lines of control flow carrying real invariants: which conditions are a benign empty table and at what log level, the pk-column prefix, and the exception wrapping. A third consumer is coming, so the flow moves into AbstractVariantLocusTablePlugin and each plugin becomes columns + label + cell composition. The country plugin's early return on empty sample metadata - which skips the VCF open entirely - is preserved deliberately, expressed as requiresSampleMetadata(). Also: - Catch TribbleException/UncheckedIOException in the base class. Files.exists() checks the VCF but not the .tbi, and MergedVcfReader opens with requireIndex = true, so a missing/stale index or truncated bgzf block threw unchecked past the WdkModelException catch with neither the variant nor the path attached. Now wrapped with both, and still an error rather than a silently empty table. - Memoize SampleMetadataLookup per (study, attribute). The ~537-row join is keyed on the study, not the variant, yet ran on every uncached record-page view, twice per page. Cached values are immutable copies; the 42P01 path is deliberately not cached, since the eda tables can appear from a data load without a webapp restart. - Reword the empty-metadata log, which conflated three causes and printed "Study null"; move a misplaced comment and correct a javadoc in VariantLocusComposer; document why countryOf() keeps its own trim(). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This was referenced Aug 12, 2026
jbrestel
added a commit
that referenced
this pull request
Aug 13, 2026
…s test suite (#20) * Add the HSSS variation plumbing design spec Makes the HighSpeedSnpSearch plugins emit variation record IDs and read the variation directory layout, so the ported HSSS searches can work at all. Two mismatches, both invisible until a search runs: the reconstruct script hardcodes a "." between sequence and location where the variation source_id needs "_", and getSearchDir returns /highSpeedSnpSearch where the real files live under /dnaseq. The ID one is the dangerous half -- it would return rows matching no variation record, so the search yields zero results and no error. Editing the plugins in place rather than subclassing is justified in the spec: the assembled model references no HSSS plugin and defines no SnpQuestions or SnpChipQuestions, because the imports are commented out in the shared apiCommonModel.xml -- so this holds for every project. Spans two repos; the idPrefix is a Conifer cohort variable in ApiCommonWebsite, not a template in this one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add the HSSS variation plumbing plan; correct the spec's test claims The spec claimed ApiCommonWebService/Test was a working JUnit harness whose tests would need updating. Both HSSS harnesses are in fact already broken, independently of this change: the JUnit module references FindPolymorphismsPlugin.PARAM_STRAIN_LIST, which is defined nowhere, so it does not compile (and is absent from build.xml's default depends); hsssTestSuite passes 8 positional args where the generator now consumes 5 standard plus 6 extras, leaving reconstructCmdName and idPrefix undefined. So no green test run gates this change, and the spec and plan now say so. Also adds a fourth production edit the spec had missed: getStrainFilterParamName returns "ngsSnp_strain_meta", which is a contract with the consuming model XML. Renaming it to variation_sample_meta now, while the plugin has no consumer, keeps snp vocabulary out of new variation XML. The plan's one real test runs the reconstruct script straight from the source tree against a checked-in fixture -- no build, database, or webserver -- and the sed that rewrites the 11 baked-in fixture IDs was verified on copies. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Build variation IDs with an underscore separator hsssReconstructSnpId hardcoded a '.' between sequence and location, so no combination of the configurable idPrefix/idSuffix could produce a VariationRecordClass source_id (Variant_Pf3D7_01_v3_29514). Both the STDOUT and STDERR joins now use '_'. Verified against the checked-in fixture: contig 80/location 896 with prefix Variant_ now yields Variant_a80_896, and the resulting shape matches a real row in apidbtuning.VariationAttributes. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add plan Task 1b: the second ID-composition site Task 1's Step 5 grep was written expecting no other script to compose source_ids the dotted way. hsssGenomicLocationsFilter does, at two branches, and it is a live alternative pipeline tail rather than dead chip code -- FindSnpsByGeneIdsPlugin overrides getGenerateScriptName to route through it, so VariationsByGeneIds would have shipped with the exact silent zero-results failure this change prevents. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Build variation IDs with an underscore in the locations filter too hsssGenomicLocationsFilter composes source_ids the same dotted way hsssReconstructSnpId did, at both of its output branches. It is not dead chip code: HsssGenomicLocationFilterScriptGenerator returns it as the final command, so it substitutes for the reconstruct script rather than following it, and FindSnpsByGeneIdsPlugin overrides getGenerateScriptName to route through it. Without this, VariationsByGeneIds would still emit dotted IDs matching no variation record -- zero results, no error -- while the isolate-group and location searches worked, since those inherit the reconstruct path. Found by Task 1's Step 5 grep, which was written expecting no second site. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Fix a plan grep that could only fail, and record two Task 1b findings Task 1b's Step 4 used grep without -F on a pattern containing ${...}, where BRE treats the braces as interval syntax and matches nothing. It reported 0 even for correct code, so it was a check an implementer could only satisfy by breaking something. Now uses -F. Also records two findings from Task 1b as explicit non-goals: the separator is now hardcoded in three places and wants a shared helper eventually, and hsssReconstructSnpId echoes every composed ID to stderr bypassing the sequence/location filter. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Read HSSS variation data from the dnaseq directory The variation HSSS files are laid out as <organism>/dnaseq/readFreq<N>/, not <organism>/highSpeedSnpSearch/readFreq<N>/, so findOrganismDir would throw 'Organism dir does not exist'. The chip plugins keep their own override. Not independently verifiable -- findOrganismDir is only reached when a search invokes the plugin. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Record hsssCopyFilesToWebSvcDir as a deliberate non-goal Flagged during Task 2: it writes into <organism>/highSpeedSnpSearch, which now disagrees with what the plugin reads. Left alone because it is a run-once snp-era prototype copier rather than a deployment path -- hardcoded prototype source dir, hardcoded organism table, dies if the target exists -- and is not what produced the dnaseq directories. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Name the strain filter param for variations, not snps getStrainFilterParamName is a contract with the model XML: the consuming filterParam must carry this exact name or the plugin rejects the request as missing a required parameter. Renaming it now, while the plugin has no consumer, keeps snp vocabulary out of new variation model XML. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Defer FindMajorAllelesPlugin's param rename to its own search's spec Flagged during Task 3: FindMajorAllelesPlugin hardcodes ngsSnp_strain_meta_a and _m (plus _wiz variants in sharedParams.xml) as its own required-param contract, with the same failure mode Task 3 fixed. Deferred rather than done now: the two-isolate-groups search it serves needs per-strain data that does not exist yet, so naming its param family before designing it would be a guess. Recorded with a note that the rename must not be forgotten, and that the _a/_m asymmetry looks like a typo for _a/_b. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Update HSSS fixtures to the variation ID convention Renames the fixture search dir to dnaseq and rewrites the 11 baked-in IDs in the three expected files from NGS_SNP.<contig>.<loc> to Variant_<contig>_<loc>. Fixes no test: both HSSS harnesses are already broken independently of this change -- the JUnit module references a constant that exists nowhere, and hsssTestSuite passes the wrong argument count to hsssGeneratePolymorphismScript. This only means a future revival starts from fixtures matching the current convention. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Correct Task 6's conifer invocation and config path Two errors found while executing it. The generated highSpeedSnpSearch-config.xml lives directly in gus_home/config/, not the per-project config/PlasmoDB/ subdirectory, so the verification grep would have reported a missing file and read as 'Task 5 failed'. And conifer configure needs --cohort, --project, --webapp-ctx and --tomcat-webapp-ctx. Without --cohort it refuses to run; without --tomcat-webapp-ctx it gets most of the way, regenerates the file being checked, then fails on log4j2.json -- a half-success where the check passes while one config silently goes unregenerated. The step now says to read the PLAY RECAP. Also drops the separate conifer install step: bld ApiCommonWebService already installs ApiCommonWebsite's conifer vars. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Mark the HSSS variation plumbing spec implemented Six commits across two repos; 14 files, 17 insertions, 17 deletions of production code. The ID fix is verified twice -- source tree and installed copy both emit Variant_a80_896 -- and the format is confirmed against a real apidbtuning.VariationAttributes row. bld succeeded, conifer configure landed idPrefix: Variant_ in the generated config, the model still loads, and the error logs stayed silent. The search directory and filter param name are correct by inspection but were never exercised: nothing invokes the plugin until a variation search exists. This is deliberately not described as end-to-end verified. Records the four plan/spec errors execution exposed -- most importantly a Step 5 grep whose false assumption uncovered a second, live ID-composition site that would have made VariationsByGeneIds return zero results silently. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Rename the chromosome param contract for variation searches FindPolymorphismsWithSeqFilterPlugin serves VariationsByLocation now. Its sibling contract, the strain filter, was already renamed to variation_sample_meta; leaving this one as chromosomeOptionalForNgsSnps would give a single processQuery two differently-named eras of the same plugin and invite the next reader to re-derive that the snp name is meaningless. The chip plugin's own PARAM_CHROMOSOME (chromosomeOptional) is untouched: it serves live chip-snp searches. * Rename the two sample-group param contracts for variation searches FindMajorAllelesPlugin serves VariationsByTwoIsolateGroups now. It extends HighSpeedSnpSearchAbstractPlugin directly rather than FindPolymorphismsPlugin, so these two names are its own constants and were not covered by the earlier strain-filter rename. The odd _m becomes _b: nothing in the plugin distinguishes it beyond being the second group, and its prompts already read Set B. With this, no snp-era param name survives in any variation search. * Accept noncoding as an SNP class in FindGenesWithSnpCharsPlugin The enumParam geneParams.snp_class offers "Non-Coding" and hsssGeneCharacteristicsFilter branches on it, but the plugin's legalParams whitelist omitted it (and listed "coding" twice), so choosing it threw PluginUserException before the script ever ran. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Design: fix the two broken statistics in the HSSS gene characteristics search Specs the real fix for SNV density (neither CDS nor coding today) and the nonsyn/syn ratio (no site normalization, worth 1.43x in pfal), plus honest labels for both and for the mislabelled "non-coding" class. Key decision recorded: numerators stay sample-set-dependent and are computed by HSSS; denominators are gene properties already derived from the genetic code in GeneVariationSummary and are read from there, so both searches rest on one definition of a synonymous site. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Plan: implement the HSSS gene statistics fix Ten tasks across ApiCommonWebService, ApiCommonModel and ApiCommonWebsite, starting by giving hsssGeneCharacteristicsFilter the regression test it never had - the suite's diff was commented out and its expected file absent. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * test: pass the missing reconstruct/prefix/suffix args in hsssTestSuite HsssScriptGenerator::extractArgs consumes the five standard args and hands the rest to getFinalCommandString, which unpacks the polymorphism args -- including reconstructCmdName, idPrefix and idSuffix -- BEFORE any stage-specific ones. Every call site in the suite supplied those three too few, so every argument after strains_list_file shifted three positions left. Effects, by stage: - polymorphism (both variants): the reconstruct command name was empty, so the generated pipeline read "... | <dir>/contigIdToSourceId.dat 1", i.e. it tried to execute the data file. Permission denied, suite aborts under set -e. - genomicLocations: same shift; the filter file landed on idPrefix. - geneChars: worst case -- geneFilters.txt was swallowed as reconstructCmdName and the gene locations file became the literal string "5". The filter has therefore never been given its gene locations at all, which is why its diff was commented out and test/expected/geneCharsFilter.txt was never created. Variant_/NULL are the values the committed expected files were built with: with this fix the ID column of polymorphismSearchWithSourceIds.txt matches exactly. The remaining columns of the three filter-stage expected files are still stale (they predate f1ac0d9); re-baselining them is a separate, reviewed change. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Plan: add Task 0 and record the test suite breakage left unfixed The original Task 1 assumed hsssTestSuite exercised the geneChars filter and discarded the result. It never reached it: four generator call sites were three args short, so the gene locations file was replaced by the literal "5". Also records three defects found and deliberately not fixed - the dead majorAlleles stage, an orphaned self-contradictory fixture, and a usage text that documents the wrong product-class encoding. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * test: cover all four product classes and re-baseline the HSSS fixtures The fixture only ever produced syn and non-syn SNPs, so the suite could not have caught a regression in the non-coding or nonsense paths -- the two this branch is about to change. Extends the fixture first, then baselines once. Fixture (6 added rows, no existing row touched): - contig 99 @ 1500, product byte 0 in strain3/strain4 and referenceGenome. productClass stays 0 because nonSyn is never set (updateCounts requires product > 0) and refProduct is not > 0 -> non-coding. Note 'X' does NOT work here: hsssFindPolymorphic maps 'X' to -1 only on the strain side (line 64), while refProduct is read raw, so 'X' (88) on the reference scores as a real product and yields syn. - contig 103 @ 30500, product byte 42 ('*') in strain3/strain4 and referenceGenome -> nonsense, and since refProduct > 0 the class is 1 negated to -1 rather than 0 negated to 0 (0 * -1 is still non-coding, so a nonsense SNP must also be coding to be reported as one). Both are polymorphic (2 of 4 strains differ from the reference) and both fall inside a gene span, so geneChars now reports non-zero nonCodingCount (g1) and nonsenseCount (g4) instead of columns that were structurally always 0. strain1/strain2 are deliberately untouched so expected/mergeStrains.txt, which is built from those two alone, keeps its meaning. Baselines regenerated. Two blessed pre-existing behaviour changes are absorbed here, neither of which had ever been re-baselined: - d3771af (2014-07-26) "%d" -> "%.1f" on the percentage columns: 100 -> 100.0. Every value is numerically unchanged; this is formatting only. - f1ac0d9 (2014-08-19) replaced the boolean y/blank product column with the four-value label: blank -> syn, y -> non-syn, plus non-coding and "has stop codon" which the old encoding could not express. expected/mergeStrainsConsensus.txt is also updated: it is derived from referenceGenome.dat, so the two new reference rows appear in it. Verified by re-running with the fixture reverted, where that stage still matched. geneCharsFilter.txt is new -- the geneChars diff in hsssTestSuite had been commented out and this file never existed, so the filter had no coverage at all. Its assertion is now enabled. geneChars needs the widened arguments "all 0 -1 0 -1 0 -1"; the previous "coding 2 5 .1 .9 3 1000" selects no genes and produces an empty file, which would assert nothing. Not addressed here, deliberately: expected/majorAlleles.txt is stale for an unrelated reason (its product letters are each one higher than the code now emits -- J vs I, U vs T, O vs N) and still fails, so the suite is not yet green end to end. hsssTestSuite:135 also calls hsssGenerateMajorAllelesScript with 12 of the 14 arguments it requires. Both are pre-existing and confirmed unaffected by this change. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Plan: correct the exit-code assertion and record the majorAlleles off-by-one hsssTestSuite runs under set -e, so it dies at the diff and every exit -1 after one is dead code; the plan's expected exit=255 was wrong. A pre-existing majorAlleles failure also means the suite exits non-zero even when geneChars passes, so checks must grep for 'matched' rather than test $?. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * test: add normalizer columns to the geneChars fixture cdsLen, synSites, nonsynSites per gene, covering the normal case, an all-zero non-coding gene, and a gene with zero synonymous sites. The filter ignores them until the next commit, so output is unchanged here. * Compute CDS density and site-normalized dN/dS in the geneChars filter Density was total variants over GENOMIC span while claiming coding variants over CDS length; it is now the latter, with the old value retained as a separate span-density column. The nonsyn/syn ratio had no site normalization and so carried the genome's codon bias - worth 1.43x in pfal3D7, where the pooled synonymous-site fraction is 17.49% rather than the textbook ~25%. Normalizers arrive per gene in geneLocations.txt; the filter treats an empty one as "statistic not defined" and excludes the gene only when the matching filter has actually been narrowed. Output grows from 8 fields to 9. The Java that parses it changes in the next commit; the two must deploy together. * Plan: warn that hsssTestSuite tests the installed copy, not the checkout Every tool the suite calls resolves off PATH into GUS_HOME/bin, so editing a filter in the checkout and re-running the suite tests the old code and can print an undeserved 'matched'. Cost most of a task to discover. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Feed per-gene normalizers to the geneChars filter and read its ninth column The gene locations query gains cds_length, syn_sites and nonsyn_sites from apidbtuning.GeneVariationSummary via LEFT JOIN, so both this search and GenesByVariantCharacteristics rest on one definition of a synonymous site rather than two that can drift. makeResultRow now expects 9 fields and maps the new span_snp_density column. This commit and the previous one must deploy together. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Plan: correct Task 5 scope and record the dormant chip-plugin coupling * Record the whole-cohort validation of the corrected geneChars statistics dN/dS reproduces (nonsyn/nonsyn_sites)/(syn/syn_sites) exactly on 31 sampled genes, max error 0.000000 - the design doc expected only a correlation, but both sides draw site counts from the same tuning table so the normalization is exactly reproducible. Median dN/dS is 0.4718 with 85% below 1, against GeneVariationSummary's independently computed piN/piS median of 0.512, where the un-normalized ratio sat at 2.0 implying genome-wide positive selection. All 254 no-CDS transcripts return with a blank CDS density and a populated span density, so non-coding genes survive. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Read the merged VCF directly to build Variant record tables (#21) * Prove htsjdk reads the merged VCF's under-declared header The published merged.ann.vcf.gz declares only CANN, CA and DFS - no GT, DP, AD/RO/QR/AO/QA or contig lines - and every INFO field opens with a bare '.'. This spike settles what htsjdk does with that before any code depends on it. htsjdk already resolves to 5.0.0 via the parent gus-project-pom:1.0.0 (dependency:tree on WSFPlugin shows com.github.samtools:htsjdk:jar:5.0.0:compile with no local pin), so no pom.xml change was needed. The spike passed on cedar against the real file - no lenient-codec workaround required. DP and CA read correctly despite being absent from the header, and CANN reads correctly despite the header's dangling leading '.' in INFO. GT is asserted against the specific expected value (GenotypeType.HOM_REF, matching the independently-verified 0:178:... call at this locus) rather than merely non-null: Genotype.getType() never returns null - a degraded parse falls back to NO_CALL (alleles present but all '.') or UNAVAILABLE (allele list empty), both non-null, so only an exact-value check catches either failure mode. That said, GT was never actually at risk from the header gap this spike targets: htsjdk's AbstractVCFCodec.createGenotypeMap() special-cases the GT key by literal match against the record's own FORMAT column, independent of any header FORMAT declaration for GT. So this test's genuine header-tolerance finding is about DP, CANN and CA, not GT - the GT assertion is a correctness guard on the sample call, included for rigor, not evidence that the GT-specific hazard was ever reachable. Skipped unless -Dvariants.vcf names the file, which lives only on the webserver. * Add CANN parsing and CA -> amino acid resolution CANN is an INFO field; CA is the per-sample FORMAT field naming which CANN entries describe that sample's genotype. Slots split on / or |, transcript keys within a slot on ';', bare 'r' means reference-with-no-CDS-annotation. Exceeds the plan's 8-test floor by two, closing coverage gaps found in review: the truncated-entry defensive branch in parse() and the '.' amino-acid filter in aminoAcidsFor() were otherwise never exercised. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Read per-sample calls from the merged VCF by tabix seek One row per sample regardless of ploidy - a diploid het collapses to an IUPAC ambiguity code, matching processSequenceVariations.jl's gt_to_base. Read frequency reports support for the CALLED allele, so a reference call shows RO/(RO+AO_total) rather than the pipeline's alt-centric 0.0, and a coverage-filled call (GT and DP only, every other FORMAT field '.') shows 100. The alt branch indexes the AO entry for the allele actually carried, by VCF-wide allele index, rather than summing all AO - summing would dilute the frequency with reads for an alt this sample doesn't carry on a multi-allelic record. The published file is biallelic by construction (write_vcf_entry splits multi-allelic sites into one record per alt: measured 40,000 loci / 548,040 alt calls, zero multi-alt, zero samples with >1 nonzero AO), so this is defensive rather than hot; synthetic multi-alt fixture records and tests exercise it directly, including both the out-of-range fallback and the in-range-but-sentinel fallback (AO correctly sized for the alt count, but '.' sits exactly at the called slot). intsOf preserves per-slot position for a multi-value FORMAT field instead of dropping '.' entries: dropping would shift later slots left and misalign an index-based lookup (a sample called for alt2 could silently read alt3's count). A '.'/unparseable slot becomes a MISSING (-1) sentinel; sums skip it, and the indexed lookup treats a sentinel at the called index as missing rather than returning it as a negative read count. A field that is entirely missing (bare '.', or every comma slot blank) still collapses to an empty array, so the coverage-filled 'no AO/RO data at all' check - which tests array length, not sentinel content - is unaffected. The fixture reproduces the real file's under-declared header on purpose. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Resolve variant locus, VCF path and EDA sample metadata The locus query returns sequence, position, name_for_filenames and the EDA sample table suffix in one round trip; verified against genomicsdb_071n. Metadata is keyed on provider_label ['country'], not the VAR_<hash> stable_id, which is a digest of the label and therefore site-specific. SampleMetadataLookup distinguishes two edaSuffix states rather than one: a null/blank suffix (no dnaseq EDA study for this organism) yields an empty map, while a non-null suffix that fails the identifier guard is a hard failure. The SQLException catch is narrowed to SQLState 42P01 (undefined_table, walking the next-exception/cause chain for wrapped drivers) so a missing attributevalue_* table still renders an empty table, but a genuine connection/pool/timeout fault is no longer swallowed into the same silent empty-map outcome. That walk carries an identity-based visited set, since a mutual getCause() cycle or a cyclic getNextException() chain (real pooled drivers produce the latter) would otherwise recurse forever. VariantLocusResolver.vcfPath now validates WEBSERVICEMIRROR and the build number before building the path, rather than silently concatenating a null into "build-null" and failing later as an inscrutable file-not-found. * Compose per-strain rows from VCF calls and EDA metadata No-call samples get a row rather than being dropped - the record already advertises no_call_strain_count, so omitting them would contradict the panel above the table. A sample resolving to several transcripts' amino acids stays on one row with the values comma-joined. * Add the VariantStrains WSF plugin Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Fix CANN parsing: htsjdk hands back a List, not a String Every row of the new Strains / Samples table rendered an empty Amino Acid, at loci where CANN was populated. Found by verifying the live page against the real VCF, not by any test. CANN is declared Number=. in the header, so htsjdk parses it into a List<String> - one element per comma-separated entry. Calling getAttributeAsString on a List-valued attribute returns the LIST'S toString: bracketed and ", "-separated, "[r0|TGT|C|..., k0|AGT|S|...]". CannIndex then split that on ',', producing a first key of "[r0" and later keys with a leading space, so no key ever matched a CA value and aminoAcidsFor returned empty for every sample. getAttributeAsStringList is NOT the fix on its own: for a genuine List it returns the elements unchanged, but for a String-valued attribute CommonInfo .getAttributeAsList wraps the whole string as a singleton without splitting on ','. That would hand CannIndex several entries glued together and corrupt the '|'-split fields at the join. So parseCann branches on the raw attribute's runtime type instead of assuming either shape. Why no test caught it: the Task 1 spike asserted only that CANN was non-null, which held either way; no MergedVcfReaderTest case asserted on aminoAcids; and the composer test builds SampleCall by hand. Nothing covered the reader -> CannIndex seam. Tests added at that seam, plus one pinning the bracketed-list form as garbage so nobody later "fixes" it by teaching the parser to tolerate brackets rather than fixing the caller. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Merge all VCF records at a locus, not just the first The dnaseq pipeline's write_vcf_entry emits one VCF record per unique ALT, so a locus routinely has several records at the same CHROM/POS. MergedVcfReader.read() took only the first record and stopped, silently rendering every sample whose alt lived on a later record as reference. Measured against the real merged.ann.vcf.gz: - 9,384 of 47,189 positions (19.9%) have more than one record; one had seven. - 13.9% of sample-locus alt calls involve a sample carrying alts on more than one record. - Concrete example, Pf3D7_01_v3:538376 (7 records, same REF, different indel ALTs): 43 samples split 22/9/3/2/3/3/1 across the records. Reading only record 1 rendered 21 of those 43 as reference. read() now consumes every record at the position and merges them per sample in MergedVcfReader.mergeSample(): - noCall only if the genotype is no-call on every record. - a "contributing" record is one where the genotype is present and NOT hom-ref - a hom-ref call on a record whose alt the sample doesn't carry is the pipeline's per-ALT splitting artifact, not a reference call. - chromosomeAlleles (the aggregation weight): no contributing record but a real call -> the first real record's alleles; exactly one contributing record -> that record's alleles as-is; more than one -> the non-ref allele from each contributing record. - aminoAcids: union across ALL records where the sample's CA resolves, each against its OWN record's CannIndex - CANN keys (r0, k0, ...) are scoped per record and must never be merged into one map, or a sample's CA would resolve against the wrong record's key. - depth/readFrequency: from a contributing record when one exists (the first, if several), else the first record with a real call. - coverageFilled: only when there is no contributing record and the reference call is coverage-filled by the existing rule. - genotype: raw GT when exactly one record contributed (or plain reference/no-call), else the slash-joined raw GTs of the contributing records - no longer shown on the page, but still used by downloads. - ploidy: max across records where the sample has a real call. LocusCalls.cann is removed: a single merged CannIndex cannot exist now that keys are per-record, and the field was never read by anything. altAlleles is now the union of ALTs across records, in record order. Allele display is unified (confirmed with the user) across the within-record het case and the new across-record case: 1 distinct chromosome allele -> that allele; 2 distinct single-character alleles -> IUPAC ambiguity code (unchanged); otherwise -> slash-joined distinct alleles (e.g. "ATT/ATATT"). This changes MergedVcfReaderTest#complexHetSlashJoinsDistinctAllelesNotIupac (previously ...ReturnsFirstNonRefAlleleNotIupac): a single-record het of A/AG now displays "A/AG" instead of dropping the reference and showing "AG" alone. Every other single-record-locus test is unchanged. Extended the fixture with a two-record locus at chr1:700 (REF=A, ALT=AT and ALT=ATT, each with its own CANN whose k0 maps to a different amino acid) and added tests for: an alt on the second record rendering instead of reference; alts split across two records rendering slash-joined with both in chromosomeAlleles; amino acids unioned and resolved per-record (proving no key collision); hom-ref on every record staying a reference call; no-call on every record staying a no-call. Also, since GT is no longer displayed (hidden via variantRecords.xml in the ApiCommonModel commit), a no-call sample must show "No call" in the Allele column instead of blank - VariantLocusComposer.strainRows and its test are updated accordingly. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Aggregate per-country allele frequencies from the VCF Ploidy-weighted per processSequenceVariations.jl: one unit per chromosome slot, so a diploid het contributes to both alleles and the denominator is the sum of ploidies rather than a sample count. Haploid P. falciparum makes the two coincide; TriTryp and Fungi aneuploids do not. The denominator is each country's own chromosomes, so a row reads on its own. Samples with no country are excluded entirely. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Tie-break country rows by name, and make the ordering tests bite Equal strain counts fell back to the VCF's sample order - deterministic per file but arbitrary to a reader. Sort is now (-strainCount, country). Two tests passed vacuously. The sort-order test grouped Mali first, so LinkedHashMap insertion order alone satisfied it; the calls are reordered so Gambia is encountered first and both positions are asserted. Verified by deleting the sort: exactly those two ordering tests go red. Also covers a real rank-3 allele, which only had its empty-string form pinned. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Normalise country lookup once, and say what the country row holds The two row builders each had their own idea of 'has a country' - getOrDefault in one, a null/empty check in the other. They now share countryOf(), which also trims: SampleMetadataLookup's SQL rejects NULL but not blank, and nothing normalises EDA free text, so ' ' was becoming its own country group and 'Mali ' was splitting Mali across two rows. CountryRow's allele fields hold rendered cells ('C (0.3333)'), not the bare allele that SampleCall.allele and StrainRow.allele hold two files away, so they are renamed *WithFrequency and the record now documents the mixed units: strainCount counts samples, the frequencies count chromosomes. Ranks beyond 3 are dropped for a three-column table; that is now commented and covered by a test rather than left for a reader to discover. Verified by mutation: removing the trim reddens exactly the three normalisation tests. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add the VariantCountrySummary WSF plugin Also normalise country values at the SampleMetadataLookup boundary so a blank or untrimmed EDA value cannot split one country into two rows for any future consumer of the map. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Extract a shared base plugin for the variant locus tables The two VCF-backed Variant record-page tables duplicated ~45 lines of control flow carrying real invariants: which conditions are a benign empty table and at what log level, the pk-column prefix, and the exception wrapping. A third consumer is coming, so the flow moves into AbstractVariantLocusTablePlugin and each plugin becomes columns + label + cell composition. The country plugin's early return on empty sample metadata - which skips the VCF open entirely - is preserved deliberately, expressed as requiresSampleMetadata(). Also: - Catch TribbleException/UncheckedIOException in the base class. Files.exists() checks the VCF but not the .tbi, and MergedVcfReader opens with requireIndex = true, so a missing/stale index or truncated bgzf block threw unchecked past the WdkModelException catch with neither the variant nor the path attached. Now wrapped with both, and still an error rather than a silently empty table. - Memoize SampleMetadataLookup per (study, attribute). The ~537-row join is keyed on the study, not the variant, yet ran on every uncached record-page view, twice per page. Cached values are immutable copies; the 42P01 path is deliberately not cached, since the eda tables can appear from a data load without a webapp restart. - Reword the empty-metadata log, which conflated three causes and printed "Study null"; move a misplaced comment and correct a javadoc in VariantLocusComposer; document why countryOf() keeps its own trim(). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds the server side of two VCF-backed tables on the PlasmoDB Variant record page. Rows are read live from the merged VCF with htsjdk rather than loaded into the database — record-page tables are uncached in WDK (
SingleRecordAnswerValueroutes non-SQL table fields throughgetUncachedResults()), so a plugin runs per page view and no cache table exists to join against. Aggregation is therefore done in Java.What is here
MergedVcfReader— tabix-indexed random access, merging all VCF records at a locus. Positions commonly carry more than one record (measured: 19.9% of 47,189), and reading only the first silently rendered alt-carrying samples as reference.CannIndex/CannEntry— theCANNINFO field, resolved per sample via theCAFORMAT field.CANNisNumber=., so htsjdk parses it to aList;getAttributeAsStringon it yields[r0|...and every product came out empty.VariantLocusComposer— per-strain rows, and per-country allele frequencies that are ploidy-weighted (one unit per chromosome slot, denominator = that country's own chromosomes) perprocessSequenceVariations.jl:1820. They deliberately differ from the locus-widesnp_major_allele_frequency.AbstractVariantLocusTablePlugin— the shared resolve → lookup → read → compose skeleton, extracted once both plugins existed; it also catches theTribbleExceptionthat a missing or stale.tbithrows past theWdkModelExceptioncatch.SampleMetadataLookup— EDA metadata by study, memoized per(study, attribute)since the values are static between builds, with the42P01path deliberately not cached so an eda table appearing after a data load is picked up without a restart.Verification
42 unit tests. Verified live against
Variant_Pf3D7_01_v3_100057(537 rows in 344ms, strain/country/depth/read-frequency cross-checked against the raw VCF) and against a 7-record locus atPf3D7_01_v3:538376, where all 43 alt-carrying samples now render.🤖 Generated with Claude Code