Fix SNV density and dN/dS in the sample-set SNV search, and revive its test suite - #20
Open
jbrestel wants to merge 30 commits into
Open
Fix SNV density and dN/dS in the sample-set SNV search, and revive its test suite#20jbrestel wants to merge 30 commits into
jbrestel wants to merge 30 commits into
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 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