Fix SNV density and dN/dS in the sample-set SNV search, and revive its test suite - #20
Merged
Conversation
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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.
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.
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>
…s 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>
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>
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>
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>
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>
…-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>
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.
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.
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>
…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>
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>
…ents Fixes the two miscomputed statistics in the HSSS gene characteristics search (SNV density was neither CDS nor coding; the nonsyn/syn ratio had no site normalization), repairs the test suite that had been silently dead since 2014, and accepts 'noncoding' as an SNP class. Requires ApiCommonModel c3193504 - the plugin now advertises span_snp_density and validateColumns rejects it unless the model declares that wsColumn.
This was referenced Aug 10, 2026
* 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>
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.
Fixes two statistics in the
GenesByNgsSnpssearch (SNV Characteristics Within a Group of Samples) that had never computed what their labels claimed, and repairs the test suite that should have caught it.Requires
ApiCommonModelPR (dnaseq-merge-experiments) to merge together. The plugin now advertises aspan_snp_densitycolumn, andPluginExecutor.validateColumnsrejects it unless the model declares thatwsColumn. A half-deploy fails at run time on every run of this search.What was wrong
Density.
hsssGeneCharacteristicsFiltercomputed$snpsCount / ((end - start) / 1000)— all variants over genomic span, including introns and UTRs — while the column, both param prompts and the help text all said "SNPs per KB (CDS)". Wrong in both numerator and denominator.The nonsyn/syn ratio was a raw count ratio with no site normalization, so it carried the genome's codon bias. The search's own PlasmoDB description apologised for this in prose and promised normalized ratios "in subsequent releases".
snp_class"Non-Coding" was unreachable —legalParamsomitted"noncoding"(and listed"coding"twice), so selecting it threwPluginUserExceptionbefore the script ran, even though both the perl filter and the generator supported it.Approach
Numerators stay sample-set-dependent and are still computed by HSSS. The denominators — CDS length and Nei-Gojobori synonymous/nonsynonymous site counts — are gene properties already derived once from the genetic code in
apidbtuning.GeneVariationSummary, so the plugin reads them there viaLEFT JOINand passes them throughgeneLocations.txt. Both this search and the precomputedGenesByVariantCharacteristicsnow rest on one definition of a synonymous site.The old gene-span density is retained as a new
span_snp_densitycolumn rather than discarded.Test suite repair
The suite had been silently dead since 2014. Four generator call sites were three arguments short, so every argument shifted left and the geneChars stage received the literal string
5where its gene locations file should have been — which is why itsdiffwas commented out and its expected file never existed. Three other expected files had been stale sinced3771afandf1ac0d9(2014) changed output formatting and the product-class encoding without re-baselining.Fixed the arguments, extended the fixture to exercise all four product classes (it previously produced only
synandnon-syn, sononCodingCountandnonsenseCountwere structurally zero), and re-baselined.Still red at the
majorAllelesstage — a separate pre-existing failure, deliberately untouched and recorded under "Known remaining breakage" in the plan.Verification
Whole-cohort run on a dev instance against
unidb_shu_a, 5,590 pfal3D7 transcripts:(nonsyn/nonsyn_sites)/(syn/syn_sites)exactly — 31/31 sampled genes, max absolute error0.000000.PF3D7_0100100, 1,644 coding variants / 6,492 bp = 253.23, reported 253.23. Its span density is 228.86 — the figure the old code reported under the CDS label.GeneVariationSummary's independently computed piN/piS median of 0.512. The un-normalized ratio sat at 2.0, implying genome-wide positive selection.Full evidence in
docs/superpowers/specs/2026-08-08-hsss-gene-stats-validation.md; design and plan alongside it.Deployment
The plugin jar is not rebuilt by
wb model. It needsbld ApiCommonWebService/WSFPluginplus a webapp reload, and the perl filter must be installed to$GUS_HOME/bin— the test suite resolves tools offPATH, not the checkout, so an un-installed edit yields a false pass.Note on scope
This branch carries 16 prior commits of dnaseq HSSS plumbing (variation ID conventions, param renames, data paths) authored before this work. They are not part of this change but are on the branch.
🤖 Generated with Claude Code
Related PRs
The variants work spans four repos. These four merge together —
ApiCommonModel's record rename is a contract the other three depend on by name:ApiCommonModelRename the Variation record to Variant; add a precomputed SNV characteristics search ApiCommonModel#212 — record rename + precomputed SNV characteristics searchApiCommonWebServiceFix SNV density and dN/dS in the sample-set SNV search, and revive its test suite #20 — density and dN/dS fixes, test-suite repairApiCommonWebsitePrefix HSSS result IDs with Variant_ instead of NGS_SNP. ApiCommonWebsite#310 — HSSS result-ID prefix (NGS_SNP.→Variant_)web-monorepoShow the Variant record id without the record-type prefix web-monorepo#1841 — record-heading override, filename tracking the renameCoupling worth knowing when sequencing the merges: #1841's customization is resolved by
recordClass.fullName, so without #212 it matches nothing and the override silently stops applying — no error, just the default heading back. #310's prefix and #20's separator together produce the ids #212's record class resolves; any two without the third leaves HSSS emitting ids the record class rejects.Related but independent, cherry-picked onto
master— merge before the four above:EbrcModelCommonOptional dataset gate for presenter template injection (partial-appDb builds) EbrcModelCommon#131 — optional dataset gate for presenter template injection. Not part of the variants work; needed to build a site against a partial appDb.The release note originally in VEuPathDB/ApiCommonWebsite#310 has been deferred (reverted in
965009fdd) because itsdateattribute was an unresolvable placeholder. Recover the drafted text withgit show a039aeee9in ApiCommonWebsite and re-PR it once the release date is known.