diff --git a/HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter b/HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter index db36e31e..649bb358 100644 --- a/HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter +++ b/HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter @@ -20,7 +20,8 @@ close(C); open(F, $geneLocationsFile) || die "Can't open gene locations file '$geneLocationsFile'\n"; my $geneLocationLine = ; chomp $geneLocationLine; -my ($filterContigId, $filterStart, $filterEnd, $filterGeneId) = split(/\t/, $geneLocationLine); +my ($filterContigId, $filterStart, $filterEnd, $filterGeneId, + $filterCdsLen, $filterSynSites, $filterNonsynSites) = split(/\t/, $geneLocationLine); my $withinFilter = 0; my $snpsCount = 0; @@ -65,7 +66,8 @@ while() { $geneLocationLine = ; chomp $geneLocationLine; last unless $geneLocationLine; - ($filterContigId, $filterStart, $filterEnd, $filterGeneId) = split(/\t/, $geneLocationLine); + ($filterContigId, $filterStart, $filterEnd, $filterGeneId, + $filterCdsLen, $filterSynSites, $filterNonsynSites) = split(/\t/, $geneLocationLine); } # if this SNP is inside the next gene, update counts @@ -86,8 +88,25 @@ sub processGene { return unless $snpsCount; my $nonCodingCount = $snpsCount - $codingCount; - my $dnds = $synCount? $nonSynCount / $synCount : undef; - my $density = $snpsCount / (($filterEnd - $filterStart) / 1000); + + # Densities. cdsDensity is coding variants over coding length, which is what this + # search has always CLAIMED to report; spanDensity is what it actually reported, kept + # under a name that admits it. cdsLen is empty for a gene with no coding sequence. + my $cdsDensity = $filterCdsLen ? 1000 * $codingCount / $filterCdsLen : undef; + my $spanDensity = 1000 * $snpsCount / ($filterEnd - $filterStart); + + # dN/dS, each count normalized by the number of sites of its class (Nei-Gojobori, + # computed from the genetic code in apidbtuning.GeneVariationSummary). Without this + # normalization the ratio carries the genome's codon bias: the pooled synonymous-site + # fraction in pfal3D7 is 17.49%, not the textbook ~25%, worth 1.43x on every gene. + # + # defined($dn), NOT $dn: a gene with zero nonsynonymous variants has dN = 0 and a real + # ratio of 0, which is a strong purifying-selection signal and exactly what someone + # filtering a low range wants. Truth-testing $dn would silently drop those genes. + # $ds IS truth-tested, because zero there is a division by zero, not a result. + my $dn = $filterNonsynSites ? $nonSynCount / $filterNonsynSites : undef; + my $ds = $filterSynSites ? $synCount / $filterSynSites : undef; + my $dnds = (defined($dn) && $ds) ? $dn / $ds : undef; if ($snpClass && ($snpsMin || $snpsMax != -1)) { if ($snpClass eq 'coding') { @@ -106,18 +125,29 @@ sub processGene { } if ($dndsMin || $dndsMax != -1) { - if ($synCount == 0 && $nonSynCount != 0) { - return 0 unless $dndsMax == -1; + # An undefined ratio cannot be shown to be in range, so the gene is excluded - but + # only because the user narrowed this filter. Leaving it alone (min 0, max -1) skips + # this block entirely, which is what keeps non-coding genes in the result. + if (!defined($dnds)) { + return 0; } else { return 0 if $dnds < $dndsMin || ($dndsMax != -1 && $dnds > $dndsMax); } } + # Filters CDS density only. Span density is reported but not filterable, to avoid a + # fifteenth and sixteenth param on a form that already carries fourteen. if ($densityMin || $densityMax != -1) { - return 0 if ($density < $densityMin || ($densityMax != -1 && $density > $densityMax)); + return 0 if (!defined($cdsDensity)); + return 0 if ($cdsDensity < $densityMin || ($densityMax != -1 && $cdsDensity > $densityMax)); } - print STDOUT join("\t", $filterGeneId, sprintf("%.2f",$density), $synCount ? sprintf("%.2f",$dnds) : undef, $synCount, $nonSynCount, $nonCodingCount, $nonsenseCount, $snpsCount) . "\n"; + print STDOUT join("\t", + $filterGeneId, + defined($cdsDensity) ? sprintf("%.2f", $cdsDensity) : '', + sprintf("%.2f", $spanDensity), + defined($dnds) ? sprintf("%.4f", $dnds) : '', + $synCount, $nonSynCount, $nonCodingCount, $nonsenseCount, $snpsCount) . "\n"; } sub usage { @@ -128,19 +158,19 @@ Usage: hsssGeneCharacteristicsFilter contig_id_file gene_locations_filter_file s Where: - contig_id_file: tab delimited, two columns, first column contig index (1,2,...); second column contig source_id - - gene_locations_filter_file: tab delimited: contig_source_id, start, end, gene_source_id. Must be sorted by location. + - gene_locations_filter_file: tab delimited: contig_source_id, start, end, gene_source_id, cds_length, syn_sites, nonsyn_sites. Must be sorted by location. The last three may be empty for a gene with no coding sequence; the statistics that need them are then reported empty. - snp_class: all, coding, noncoding, synonymous, nonsynonymous, nonsense - - snps_min: min percent of SNPs in the gene that belong to the specified class - - snps_min: max percent of SNPs in the gene that belong to the specified class - - dnds_min: min dn/ds ratio - - dnds_min: max dn/ds ratio - - density_min: min SNPs density - - density_max: max SNPs density + - snps_min: min NUMBER of SNPs in the gene that belong to the specified class + - snps_max: max NUMBER of SNPs in the gene that belong to the specified class + - dnds_min: min site-normalized dN/dS ratio + - dnds_max: max site-normalized dN/dS ratio + - density_min: min coding SNPs per kb of CDS + - density_max: max coding SNPs per kb of CDS - snp_search_result: tab_delimited where first column is contig index and second is gene location. Replaces the first two columns of snp_search_result with a single column that is the concatenation of the contig_source_id-location, ie, a snp source id. -Outputs these columns (tab delim): geneId density dndsRatio synCount nonSynCount nonCodingCount nonsenseCount snpsCount +Outputs these columns (tab delim): geneId cdsDensity spanDensity dndsRatio synCount nonSynCount nonCodingCount nonsenseCount snpsCount "; } diff --git a/HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter b/HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter index a4cbed7c..7920f182 100644 --- a/HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter +++ b/HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter @@ -48,7 +48,7 @@ while() { elsif ($contigSourceId eq $filterContigId && $location >= $filterStart && $location <= $filterEnd) { $idPrefix = $idPrefix=~/^NULL$/ ? '' : $idPrefix; $idSuffix = $idSuffix =~/^NULL$/ ? '' : $idSuffix; - print STDOUT join("\t", $idPrefix."$contigSourceId.$location".$idSuffix, @fields) . "\n"; + print STDOUT join("\t", $idPrefix."${contigSourceId}_${location}".$idSuffix, @fields) . "\n"; } # read next filter if beyond current filter, and print if within that next filter @@ -64,7 +64,7 @@ while() { if ($contigSourceId eq $filterContigId && $location >= $filterStart && $location <= $filterEnd) { $idPrefix = $idPrefix=~/^NULL$/ ? '' : $idPrefix; $idSuffix = $idSuffix =~/^NULL$/ ? '' : $idSuffix; - print STDOUT join("\t", $idPrefix."$contigSourceId.$location".$idSuffix, @fields) . "\n"; + print STDOUT join("\t", $idPrefix."${contigSourceId}_${location}".$idSuffix, @fields) . "\n"; } } } diff --git a/HighSpeedSnpSearch/bin/hsssReconstructSnpId b/HighSpeedSnpSearch/bin/hsssReconstructSnpId index dcf32c8a..9042a680 100644 --- a/HighSpeedSnpSearch/bin/hsssReconstructSnpId +++ b/HighSpeedSnpSearch/bin/hsssReconstructSnpId @@ -39,8 +39,8 @@ while() { die "Can't map contigIndex '$contigIndex' in stdin" unless $contigSourceId; $prefix = $prefix=~/^NULL$/ ? '' : $prefix; $suffix = $suffix =~/^NULL$/ ? '' : $suffix; - print STDERR join("\t", $prefix."$contigSourceId.$location".$suffix, @fields) . "\n" ; - print STDOUT join("\t", $prefix."$contigSourceId.$location".$suffix, @fields) . "\n" unless ($seqFilter && ($contigSourceId ne $seqFilter || $location < $minLoc || $location > $maxLoc)); + print STDERR join("\t", $prefix."${contigSourceId}_${location}".$suffix, @fields) . "\n" ; + print STDOUT join("\t", $prefix."${contigSourceId}_${location}".$suffix, @fields) . "\n" unless ($seqFilter && ($contigSourceId ne $seqFilter || $location < $minLoc || $location > $maxLoc)); } sub usage { diff --git a/HighSpeedSnpSearch/bin/hsssTestSuite b/HighSpeedSnpSearch/bin/hsssTestSuite index 8597fce5..98d323a8 100644 --- a/HighSpeedSnpSearch/bin/hsssTestSuite +++ b/HighSpeedSnpSearch/bin/hsssTestSuite @@ -38,7 +38,7 @@ echo "matched" # generate findPolymorphism script echo -e "1\n2\n3\n4" > strainsList.txt -hsssGeneratePolymorphismScript $testDir $testDir 1 runPolymorphismSearch polymorphismSearch_result.txt 20 1 strainsList.txt +hsssGeneratePolymorphismScript $testDir $testDir 1 runPolymorphismSearch polymorphismSearch_result.txt 20 1 strainsList.txt hsssReconstructSnpId Variant_ NULL chmod +x runPolymorphismSearch # run that script and compare output with expected. @@ -56,7 +56,7 @@ echo "" # generate findPolymorphism script with single genomic location filter echo -e "1\n2\n3\n4" > strainsList.txt -hsssGeneratePolymorphismScript $testDir $testDir 1 runPolymorphismSearchWithFilter polymorphismSearchWithFilter_result.txt 20 1 strainsList.txt f100 21 25 +hsssGeneratePolymorphismScript $testDir $testDir 1 runPolymorphismSearchWithFilter polymorphismSearchWithFilter_result.txt 20 1 strainsList.txt hsssReconstructSnpId Variant_ NULL f100 21 25 chmod +x runPolymorphismSearchWithFilter # run that script and compare output with expected. @@ -73,7 +73,7 @@ echo "" # generate findPolymorphism script with genomic locations filter echo -e "1\n2\n3\n4" > strainsList.txt -hsssGenerateGenomicLocationsScript $testDir $testDir 1 runGenomicLocations genomicLocations_result.txt 20 1 strainsList.txt $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/textData/genomicLocationFilters.txt +hsssGenerateGenomicLocationsScript $testDir $testDir 1 runGenomicLocations genomicLocations_result.txt 20 1 strainsList.txt hsssReconstructSnpId Variant_ NULL $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/textData/genomicLocationFilters.txt chmod +x runGenomicLocations # run that script and compare output with expected. @@ -90,19 +90,19 @@ echo "" # generate findPolymorphism script with genes filter echo -e "1\n2\n3\n4" > strainsList.txt -hsssGenerateGeneCharsScript $testDir $testDir 1 runGeneChars geneChars_result.txt 20 1 strainsList.txt $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/textData/geneFilters.txt coding 2 5 .1 .9 3 1000 +hsssGenerateGeneCharsScript $testDir $testDir 1 runGeneChars geneChars_result.txt 20 1 strainsList.txt hsssReconstructSnpId Variant_ NULL $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/textData/geneFilters.txt all 0 -1 0 -1 0 -1 chmod +x runGeneChars # run that script and compare output with expected. ./runGeneChars -#echo "Comparing expected runGeneChars output with result..." -#diff $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt geneChars_result.txt -#diffStat=$? -#if [ $diffStat != 0 ]; then -# exit -1 -#fi -#echo "matched" +echo "Comparing expected runGeneChars output with result..." +diff $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt geneChars_result.txt +diffStat=$? +if [ $diffStat != 0 ]; then + exit -1 +fi +echo "matched" echo "" # test making a consensus from a merged output diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/1 b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/1 similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/1 rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/1 diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/2 b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/2 similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/2 rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/2 diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/3 b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/3 similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/3 rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/3 diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/4 b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/4 similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/4 rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/4 diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/contigIdToSourceId.dat b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/contigIdToSourceId.dat similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/contigIdToSourceId.dat rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/contigIdToSourceId.dat diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/referenceGenome.dat b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/referenceGenome.dat similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/referenceGenome.dat rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/referenceGenome.dat diff --git a/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/strainIdToName.dat b/HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/strainIdToName.dat similarity index 100% rename from HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/readFreq80/strainIdToName.dat rename to HighSpeedSnpSearch/test/TestDB/Hsapiens123/dnaseq/readFreq80/strainIdToName.dat diff --git a/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt b/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt new file mode 100644 index 00000000..41d5a555 --- /dev/null +++ b/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt @@ -0,0 +1,4 @@ +g1 0.83 1.00 0.0000 1 0 1 0 2 +g2 1.67 1.44 0.0000 1 0 0 0 1 +g3 0.03 0 1 0 0 1 +g4 2.22 0.20 0 1 0 1 2 diff --git a/HighSpeedSnpSearch/test/expected/genomicLocationFilter.txt b/HighSpeedSnpSearch/test/expected/genomicLocationFilter.txt index 3bd51a43..0713f7bc 100644 --- a/HighSpeedSnpSearch/test/expected/genomicLocationFilter.txt +++ b/HighSpeedSnpSearch/test/expected/genomicLocationFilter.txt @@ -1,2 +1,4 @@ -NGS_SNP.e99.2011 100 50 -NGS_SNP.h103.30021 100 25 y +Variant_e99_1500 100.0 50.0 non-coding +Variant_e99_2011 100.0 50.0 syn +Variant_h103_30021 100.0 25.0 non-syn +Variant_h103_30500 100.0 50.0 has stop codon diff --git a/HighSpeedSnpSearch/test/expected/mergeStrainsConsensus.txt b/HighSpeedSnpSearch/test/expected/mergeStrainsConsensus.txt index c6f3962c..f14f63f7 100644 --- a/HighSpeedSnpSearch/test/expected/mergeStrainsConsensus.txt +++ b/HighSpeedSnpSearch/test/expected/mergeStrainsConsensus.txt @@ -2,10 +2,12 @@ 86 13441 1 80 0 2 80 0 7500 2500 0 88 150 1 76 0 0 0 0 10000 0 0 90 876 2 67 0 0 0 0 10000 0 0 +99 1500 2 0 0 0 0 0 10000 0 0 99 2011 2 73 0 0 0 0 10000 0 0 100 23 3 76 0 4 77 0 4000 4000 1 102 4334 3 84 0 1 84 0 7500 2500 0 103 30021 2 69 0 1 67 0 7500 2500 0 +103 30500 2 42 0 0 0 0 10000 0 0 104 3002 3 78 0 2 78 0 7500 2500 0 201 54 3 73 0 0 0 0 10000 0 0 302 91 4 81 0 0 0 0 10000 0 0 diff --git a/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIds.txt b/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIds.txt index a8ba36f7..8d716e72 100644 --- a/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIds.txt +++ b/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIds.txt @@ -1,8 +1,10 @@ -NGS_SNP.a80.896 100 25 y -NGS_SNP.b86.13441 100 50 -NGS_SNP.e99.2011 100 50 -NGS_SNP.f100.23 100 20 -NGS_SNP.g102.4334 100 50 y -NGS_SNP.h103.30021 100 25 y -NGS_SNP.i104.3002 100 50 -NGS_SNP.j201.54 100 20 +Variant_a80_896 100.0 25.0 non-syn +Variant_b86_13441 100.0 50.0 syn +Variant_e99_1500 100.0 50.0 non-coding +Variant_e99_2011 100.0 50.0 syn +Variant_f100_23 100.0 20.0 syn +Variant_g102_4334 100.0 50.0 non-syn +Variant_h103_30021 100.0 25.0 non-syn +Variant_h103_30500 100.0 50.0 has stop codon +Variant_i104_3002 100.0 50.0 syn +Variant_j201_54 100.0 20.0 syn diff --git a/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIdsAndSeqFilter.txt b/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIdsAndSeqFilter.txt index 169a8806..410d5b8d 100644 --- a/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIdsAndSeqFilter.txt +++ b/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIdsAndSeqFilter.txt @@ -1 +1 @@ -NGS_SNP.f100.23 100 20 +Variant_f100_23 100.0 20.0 syn diff --git a/HighSpeedSnpSearch/test/textData/geneFilters.txt b/HighSpeedSnpSearch/test/textData/geneFilters.txt index b17a8c76..826e0920 100644 --- a/HighSpeedSnpSearch/test/textData/geneFilters.txt +++ b/HighSpeedSnpSearch/test/textData/geneFilters.txt @@ -1,5 +1,5 @@ -e99 1000 3000 g1 -f100 5 700 g2 -g102 3001 40000 g3 -h103 30021 40000 g4 -j201 20 50 g5 +e99 1000 3000 g1 1200 300 900 +f100 5 700 g2 600 150 450 +g102 3001 40000 g3 0 0 0 +h103 30021 40000 g4 900 0 675 +j201 20 50 g5 300 75 225 diff --git a/HighSpeedSnpSearch/test/textData/referenceGenome.txt b/HighSpeedSnpSearch/test/textData/referenceGenome.txt index 10a2c801..e0d4ad89 100644 --- a/HighSpeedSnpSearch/test/textData/referenceGenome.txt +++ b/HighSpeedSnpSearch/test/textData/referenceGenome.txt @@ -2,11 +2,13 @@ 86 13441 1 80 88 150 1 76 90 876 2 67 +99 1500 2 0 99 2011 2 73 100 23 4 77 102 4334 3 84 103 7 2 69 103 30021 2 69 +103 30500 2 42 104 3002 3 78 201 54 3 73 302 91 4 81 diff --git a/HighSpeedSnpSearch/test/textData/strain3.txt b/HighSpeedSnpSearch/test/textData/strain3.txt index 1dc24a22..5aaab9d2 100644 --- a/HighSpeedSnpSearch/test/textData/strain3.txt +++ b/HighSpeedSnpSearch/test/textData/strain3.txt @@ -1,8 +1,10 @@ 86 13441 2 80 +99 1500 1 0 99 2011 1 73 100 23 3 76 102 4334 2 71 103 7 0 0 +103 30500 1 42 104 3002 2 78 201 54 1 73 201 54 3 73 diff --git a/HighSpeedSnpSearch/test/textData/strain4.txt b/HighSpeedSnpSearch/test/textData/strain4.txt index 90dc5802..74e20751 100644 --- a/HighSpeedSnpSearch/test/textData/strain4.txt +++ b/HighSpeedSnpSearch/test/textData/strain4.txt @@ -1,5 +1,7 @@ 80 896 2 80 +99 1500 1 0 99 2011 1 73 100 23 3 76 102 4334 2 71 103 7 0 0 +103 30500 1 42 diff --git a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java index 46248c9b..7ab61c24 100644 --- a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java +++ b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java @@ -25,8 +25,12 @@ */ public class FindGenesWithSnpCharsPlugin extends FindPolymorphismsPlugin { + // Must stay in step with the enumList of geneParams.snp_class in ApiCommonModel and + // with the classes hsssGeneCharacteristicsFilter branches on. An entry missing here + // fails at run time, not at build time - which is how "noncoding" was unreachable + // despite both the other two layers supporting it. private static final Set legalParams = new HashSet(Arrays.asList(new String[] { "coding", - "nonsynonymous", "synonymous", "nonsense", "all", "coding" })); + "noncoding", "nonsynonymous", "synonymous", "nonsense", "all" })); private static final String geneLocationsFileName = "geneLocations.txt"; @@ -45,6 +49,7 @@ public class FindGenesWithSnpCharsPlugin extends FindPolymorphismsPlugin { public static final String COLUMN_PROJECT_ID = "project_id"; public static final String COLUMN_MATCHED_RESULT = "matched_result"; public static final String COLUMN_DENSITY = "cds_snp_density"; + public static final String COLUMN_SPAN_DENSITY = "span_snp_density"; public static final String COLUMN_DNDS = "ngs_dn_ds_ratio"; public static final String COLUMN_SYN = "ngs_num_synonymous"; public static final String COLUMN_NONSYN = "ngs_num_non_synonymous"; @@ -82,8 +87,11 @@ protected void initForBashScript(File jobDir, Map params, File o FileWriter w = new FileWriter(filtersFile); bw = new BufferedWriter(w); if (snpClass.equals("unit test")) { - String[] testFilters = new String[] { "e99\t1000\t3000\tg1", "f100\t500\t700\tg2", - "h103\t30021\t40000\tg3", "j201\t20\t50\tg4" }; + String[] testFilters = new String[] { + "e99\t1000\t3000\tg1\t1200\t300\t900", + "f100\t500\t700\tg2\t600\t150\t450", + "h103\t30021\t40000\tg3\t\t\t", // no coding sequence: normalizers empty + "j201\t20\t50\tg4\t900\t0\t675" }; // zero synonymous sites: ratio undefined for (String filter : testFilters) { bw.write(filter); bw.newLine(); @@ -96,8 +104,21 @@ protected void initForBashScript(File jobDir, Map params, File o String organism = removeSingleQuotes(params.get(PARAM_ORGANISM)); // can interpolate organism into sql w/o fear of injection because it came from a vocabulary param - String sql = "select g.sequence_id, g.start_min, g.end_max, g.source_id" + newline + - "from webready.GeneAttributes_p g " + newline + "where g.source_id is not null" + newline + + // LEFT JOIN, never inner: GeneVariationSummary holds one row per gene that has + // COHORT variants (5,579 of 5,720 annotated pfal3D7 genes), and a gene missing + // from it must still get a locations line and still report its counts. Only its + // normalized statistics come back empty. + // + // These three columns are gene properties derived from the genetic code, not + // from any sample set, which is why reading them here is sound: the numerators + // stay sample-set-dependent and HSSS still computes them. + String sql = "select g.sequence_id, g.start_min, g.end_max, g.source_id," + newline + + " gvs.cds_length, gvs.syn_sites, gvs.nonsyn_sites" + newline + + "from webready.GeneAttributes_p g " + newline + + "left join apidbtuning.GeneVariationSummary gvs" + newline + + " on gvs.gene_source_id = g.source_id" + newline + + " and gvs.project_id = g.project_id" + newline + + "where g.source_id is not null" + newline + " and g.organism = '" + organism + "'"; ResultSet rs = null; @@ -110,7 +131,13 @@ protected void initForBashScript(File jobDir, Map params, File o String start = rs.getString(2); String end = rs.getString(3); String geneId = rs.getString(4); - bw.write(seqId + "\t" + start + "\t" + end + "\t" + geneId); + // getString returns null for a SQL NULL; the filter tests these for truth, + // so an empty string reads as "no normalizer" exactly like a zero would. + String cdsLen = rs.getString(5) == null ? "" : rs.getString(5); + String synSites = rs.getString(6) == null ? "" : rs.getString(6); + String nonsynSites = rs.getString(7) == null ? "" : rs.getString(7); + bw.write(seqId + "\t" + start + "\t" + end + "\t" + geneId + "\t" + + cdsLen + "\t" + synSites + "\t" + nonsynSites); bw.newLine(); } @@ -150,8 +177,8 @@ protected void initForBashScript(File jobDir, Map params, File o @Override public String[] getColumns(PluginRequest request) { - return new String[] { COLUMN_GENE_SOURCE_ID, COLUMN_PROJECT_ID, COLUMN_DENSITY, COLUMN_DNDS, - COLUMN_SYN, COLUMN_NONSYN, COLUMN_NONCODING, COLUMN_NONSENSE, COLUMN_TOTAL }; + return new String[] { COLUMN_GENE_SOURCE_ID, COLUMN_PROJECT_ID, COLUMN_DENSITY, COLUMN_SPAN_DENSITY, + COLUMN_DNDS, COLUMN_SYN, COLUMN_NONSYN, COLUMN_NONCODING, COLUMN_NONSENSE, COLUMN_TOTAL }; } @Override @@ -192,22 +219,23 @@ protected String getGenerateScriptName() { @Override protected String[] makeResultRow(String[] parts, Map columns, String projectId) throws PluginModelException { - if (parts.length != 8) - throw new PluginModelException("Wrong number of columns in results file. Expected 8, found " + + if (parts.length != 9) + throw new PluginModelException("Wrong number of columns in results file. Expected 9, found " + parts.length); - String[] row = new String[11]; + String[] row = new String[12]; row[columns.get(COLUMN_GENE_SOURCE_ID)] = parts[0]; row[columns.get(COLUMN_SOURCE_ID)] = null; row[columns.get(COLUMN_PROJECT_ID)] = projectId; row[columns.get(COLUMN_MATCHED_RESULT)] = "Y"; row[columns.get(COLUMN_DENSITY)] = parts[1]; - row[columns.get(COLUMN_DNDS)] = parts[2]; - row[columns.get(COLUMN_SYN)] = parts[3]; - row[columns.get(COLUMN_NONSYN)] = parts[4]; - row[columns.get(COLUMN_NONCODING)] = parts[5]; - row[columns.get(COLUMN_NONSENSE)] = parts[6]; - row[columns.get(COLUMN_TOTAL)] = parts[7]; + row[columns.get(COLUMN_SPAN_DENSITY)] = parts[2]; + row[columns.get(COLUMN_DNDS)] = parts[3]; + row[columns.get(COLUMN_SYN)] = parts[4]; + row[columns.get(COLUMN_NONSYN)] = parts[5]; + row[columns.get(COLUMN_NONCODING)] = parts[6]; + row[columns.get(COLUMN_NONSENSE)] = parts[7]; + row[columns.get(COLUMN_TOTAL)] = parts[8]; return row; } } diff --git a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindMajorAllelesPlugin.java b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindMajorAllelesPlugin.java index a72dbf30..47fe7a86 100644 --- a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindMajorAllelesPlugin.java +++ b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindMajorAllelesPlugin.java @@ -17,11 +17,11 @@ public class FindMajorAllelesPlugin extends HighSpeedSnpSearchAbstractPlugin { // required parameter definition - public static final String PARAM_STRAIN_FILTER_A = "ngsSnp_strain_meta_a"; + public static final String PARAM_STRAIN_FILTER_A = "variation_sample_meta_a"; public static final String PARAM_MIN_PERCENT_KNOWNS_A = "MinPercentIsolateCalls"; public static final String PARAM_MIN_PERCENT_MAJOR_ALLELES_A = "MinPercentMajorAlleles"; public static final String PARAM_READ_FREQ_PERCENT_A = "ReadFrequencyPercent"; - public static final String PARAM_STRAIN_FILTER_B = "ngsSnp_strain_meta_m"; + public static final String PARAM_STRAIN_FILTER_B = "variation_sample_meta_b"; public static final String PARAM_MIN_PERCENT_KNOWNS_B = "MinPercentIsolateCallsTwo"; public static final String PARAM_MIN_PERCENT_MAJOR_ALLELES_B = "MinPercentMajorAllelesTwo"; public static final String PARAM_READ_FREQ_PERCENT_B = "ReadFrequencyPercentTwo"; diff --git a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java index 166cf318..0ed8ddc8 100644 --- a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java +++ b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java @@ -40,7 +40,7 @@ public FindPolymorphismsPlugin() { @Override protected String getStrainFilterParamName() { - return "ngsSnp_strain_meta"; + return "variation_sample_meta"; } diff --git a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsWithSeqFilterPlugin.java b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsWithSeqFilterPlugin.java index cda8be13..033fd5a4 100644 --- a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsWithSeqFilterPlugin.java +++ b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsWithSeqFilterPlugin.java @@ -14,7 +14,7 @@ public class FindPolymorphismsWithSeqFilterPlugin extends FindPolymorphismsPlugin { // required parameter definition - public static final String PARAM_CHROMOSOME = "chromosomeOptionalForNgsSnps"; + public static final String PARAM_CHROMOSOME = "chromosomeOptionalForVariations"; public static final String PARAM_SEQUENCE = "sequenceId"; public static final String PARAM_START_POINT = "start_point"; public static final String PARAM_END_POINT = "end_point"; diff --git a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java index 7ad4df65..de9d747b 100644 --- a/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java +++ b/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java @@ -194,7 +194,7 @@ String removeSingleQuotes(String inputText) { } protected String getSearchDir() { - return "/highSpeedSnpSearch"; + return "/dnaseq"; } File findOrganismDir(Map params, String projectId) throws PluginModelException, PluginUserException { diff --git a/docs/superpowers/plans/2026-08-05-hsss-variation-plumbing.md b/docs/superpowers/plans/2026-08-05-hsss-variation-plumbing.md new file mode 100644 index 00000000..3c969715 --- /dev/null +++ b/docs/superpowers/plans/2026-08-05-hsss-variation-plumbing.md @@ -0,0 +1,922 @@ +# HSSS Variation Plumbing Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Make the HighSpeedSnpSearch (HSSS) plugins emit `variation` record IDs and read +the variation HSSS directory layout, so the ported HSSS searches can work at all. + +**Architecture:** Four small edits across two repos — one Perl line (the ID separator), two +Java one-liners (the search directory and the strain-filter param name), and one Conifer +variable (the ID prefix). Plus fixture hygiene. No new files, no new classes: the plugins +are variation-only from here, so they are edited in place. + +**Tech Stack:** Java (WSF plugins, built with `bld`), Perl (the `hsss*` scripts installed +into `$GUS_HOME/bin`), Conifer/Ansible (site config generation), `bld` and +`bin/veup-build.sh` from the `agentic-veupath-dev` control plane. + +**Spec:** `docs/superpowers/specs/2026-08-05-hsss-variation-plumbing-design.md`. Section +references below prefixed `spec §` point there. + +--- + +## Read this before starting + +**There is no test suite to lean on.** Both HSSS harnesses are already broken, +independently of this change (spec §4): the JUnit module references a constant +(`FindPolymorphismsPlugin.PARAM_STRAIN_LIST`) that exists nowhere and therefore does not +compile, and `hsssTestSuite` passes the wrong number of positional args to +`hsssGeneratePolymorphismScript`. **Do not attempt to run either, and do not report a green +test run.** Reviving them is explicitly out of scope. + +What we *do* have is one genuinely runnable check that proves the load-bearing half of this +change, needs no build, no database, and no webserver — Task 1 is built around it. + +**Two of the four edits cannot be verified until the follow-on search exists** (Task 2 and +Task 3). They are correct by inspection against the spec's evidence. Say so plainly rather +than implying they were exercised. + +**Repos and branch.** Both repos are already on `dnaseq-merge-experiments`: + +| repo | path | +|---|---| +| `ApiCommonWebService` | `~/workspaces/plasmodb/ApiCommonWebService` | +| `ApiCommonWebsite` | `~/workspaces/plasmodb/ApiCommonWebsite` | + +Do **not** create a git worktree and do **not** switch branches. `~/workspaces/plasmodb` +is the source of a mutagen sync to the remote webserver; a worktree would not be synced and +the remote build would compile the un-edited files. + +- [ ] **Prerequisite: confirm both branches and that sync is up** + +```bash +for r in ApiCommonWebService ApiCommonWebsite; do + printf '%-22s %s\n' "$r" "$(git -C ~/workspaces/plasmodb/$r rev-parse --abbrev-ref HEAD)" +done +cd ~/workspaces/agentic-veupath-dev && bin/veup-sync-up.sh plasmodb +``` + +Expected: both print `dnaseq-merge-experiments`, and the `plasmodb` row of the roster reads +`up`. **If either prints `main`, stop.** + +## File Structure + +| file | change | +|---|---| +| `ApiCommonWebService/HighSpeedSnpSearch/bin/hsssReconstructSnpId` | ID separator `.` → `_` (Task 1) | +| `ApiCommonWebService/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java` | `getSearchDir()` → `/dnaseq` (Task 2) | +| `ApiCommonWebService/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java` | `getStrainFilterParamName()` → `variation_sample_meta` (Task 3) | +| `ApiCommonWebService/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/` | rename to `dnaseq/` (Task 4) | +| `ApiCommonWebService/HighSpeedSnpSearch/test/expected/{genomicLocationFilter,polymorphismSearchWithSourceIds,polymorphismSearchWithSourceIdsAndSeqFilter}.txt` | `NGS_SNP..` → `Variant__` (Task 4) | +| `ApiCommonWebsite/Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml` | `highspeedsnpsearchconfig_idPrefix: Variant_` (Task 5) | + +One task per edit, each independently committable. Task 1 first because it is the only one +with a real test, and Task 6 deploys and re-verifies. + +--- + +### Task 1: Fix the ID separator + +This is the change the whole spec exists for. `hsssReconstructSnpId` builds +`$prefix . "$contigSourceId.$location" . $suffix`; prefix and suffix are configurable but +the `.` is hardcoded, so no amount of config can produce a `Variant__` ID. + +**Files:** +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/bin/hsssReconstructSnpId:42-43` + +- [ ] **Step 1: Run the failing test and record its output** + +The script needs only `strict`, so it runs straight from the source tree with the +checked-in fixture — no build, no database, no webserver: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch && \ +printf '80\t896\t100\t25\t1\n' \ + | perl bin/hsssReconstructSnpId test/textData/contigIdToSourceId.dat 1 Variant_ NULL 2>/dev/null +``` + +Expected **now** (this is the bug): + +``` +Variant_a80.896 100 25 syn +``` + +Note the `.` between `a80` and `896`. The fixture maps contig index `80 → a80`; index `1` +does **not** exist in it and would make the script die with +`Can't map contigIndex '1' in stdin`, so use `80`. The trailing `1` in the input is the +coding class, which the script translates to `syn`. + +- [ ] **Step 2: Make the change** + +In `bin/hsssReconstructSnpId`, replace lines 42-43: + +```perl + print STDERR join("\t", $prefix."$contigSourceId.$location".$suffix, @fields) . "\n" ; + print STDOUT join("\t", $prefix."$contigSourceId.$location".$suffix, @fields) . "\n" unless ($seqFilter && ($contigSourceId ne $seqFilter || $location < $minLoc || $location > $maxLoc)); +``` + +with: + +```perl + print STDERR join("\t", $prefix."${contigSourceId}_${location}".$suffix, @fields) . "\n" ; + print STDOUT join("\t", $prefix."${contigSourceId}_${location}".$suffix, @fields) . "\n" unless ($seqFilter && ($contigSourceId ne $seqFilter || $location < $minLoc || $location > $maxLoc)); +``` + +**Both** lines change — STDERR and STDOUT emit the same ID and must stay consistent. + +The `${...}` braces are required, not stylistic: `"$prefix$contigSourceId_$location"` would +make Perl look for a variable named `$contigSourceId_`, which is undefined, silently +producing `Variant_896`. + +- [ ] **Step 3: Run the test again** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch && \ +printf '80\t896\t100\t25\t1\n' \ + | perl bin/hsssReconstructSnpId test/textData/contigIdToSourceId.dat 1 Variant_ NULL 2>/dev/null +``` + +Expected: + +``` +Variant_a80_896 100 25 syn +``` + +- [ ] **Step 4: Confirm the shape matches a real variation ID** + +The synthetic fixture proves the format; this confirms the format is the *right* one. +Requires the ssh tunnel to genomicsdb on port 5432: + +```bash +psql -h localhost -p 5432 -d unidb_shu_a -tAc \ + "SELECT source_id FROM apidbtuning.VariationAttributes WHERE source_id = 'Variant_Pf3D7_01_v3_29514'" +``` + +Expected: `Variant_Pf3D7_01_v3_29514` (one row). That is `Variant_` + `Pf3D7_01_v3` + `_` + +`29514` — the exact composition Step 3 now produces. If this returns nothing, stop: either +the tunnel is down or the record's ID convention is not what the spec assumed. + +- [ ] **Step 5: Check no other script builds IDs the same way** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/bin && \ + grep -n 'contigSourceId\.\$location\|contigSourceId\.\$loc' * ; echo "exit: $?" +``` + +Expected: `exit: 1` (no matches) — `hsssReconstructChipSnpId` is for the dead chip path and +composes IDs differently; if this *does* match something, report it rather than changing it, +since chip is out of scope (spec §7). + +> **This expectation was wrong, and the step earned its keep by catching it.** The grep +> matches `hsssGenomicLocationsFilter:51` and `:67`, a second ID-composition site on a +> *live* alternative pipeline tail. Handled by Task 1b — do not change it here. + +- [ ] **Step 6: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add HighSpeedSnpSearch/bin/hsssReconstructSnpId +git commit -m "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 " +``` + +--- + +### Task 1b: Fix the separator in the second ID-composition site + +**Added during execution.** Task 1's Step 5 was written expecting no other script to compose +IDs the same way. It does: `hsssGenomicLocationsFilter` builds +`$idPrefix."$contigSourceId.$location".$idSuffix` at **two** places, and it is not dead code +on the chip path — it is a live alternative tail of the *same* pipeline. + +Why it matters, concretely. `HsssGenomicLocationFilterScriptGenerator.pm:14` returns this +script from `getFinalCommandString`, so it **substitutes for** `hsssReconstructSnpId` rather +than running after it. Tracing which planned search reaches which tail: + +| plugin | generate script | ID composed by | fixed by | +|---|---|---|---| +| `FindPolymorphismsPlugin` → `VariationsByIsolateGroup` | `hsssGeneratePolymorphismScript` (inherited) | `hsssReconstructSnpId` | Task 1 | +| `FindPolymorphismsWithSeqFilterPlugin` → `VariationsByLocation` | inherited, not overridden | `hsssReconstructSnpId` | Task 1 | +| **`FindSnpsByGeneIdsPlugin`** → **`VariationsByGeneIds`** | **overrides** (`:112`) to `hsssGenerateGenomicLocationsScript` | **`hsssGenomicLocationsFilter`** | **this task** | + +So exactly one of the four searches being ported would still emit +`Variant_Pf3D7_01_v3.29514` and fail the way this change exists to prevent — silently, with +zero results. Fixing it now costs the same two lines; deferring it buys a future debugging +session. + +**Files:** +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter:51` and `:67` + +- [ ] **Step 1: Confirm both sites and their context** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/bin && \ + grep -n 'contigSourceId\.\$location' hsssGenomicLocationsFilter +``` + +Expected: two hits, lines 51 and 67. They are the same statement in two branches of the +filter's control flow — the "within current filter" branch and the "within the next filter" +branch — so both must change or gene-ID searches would emit inconsistent IDs depending on +which branch a given variant took. + +- [ ] **Step 2: Make the change at both sites** + +Replace, at both line 51 and line 67: + +```perl + print STDOUT join("\t", $idPrefix."$contigSourceId.$location".$idSuffix, @fields) . "\n"; +``` + +with: + +```perl + print STDOUT join("\t", $idPrefix."${contigSourceId}_${location}".$idSuffix, @fields) . "\n"; +``` + +Note the indentation differs between the two sites (line 67 sits one level deeper inside the +`while`/`if`). Preserve each line's existing leading whitespace; change only the +interpolation. The `${...}` braces are required for the same reason as Task 1 — +`$contigSourceId_` would be read as an undefined variable name. + +- [ ] **Step 3: Verify no dotted composition remains anywhere in `bin/`** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/bin && \ + grep -n 'contigSourceId\.\$location\|contigSourceId\.\$loc' * ; echo "exit: $?" +``` + +Expected: `exit: 1`, no matches — this is now the assertion Task 1's Step 5 was originally +written to make. + +- [ ] **Step 4: Confirm the underscore form is present twice** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/bin && \ + grep -cF '${contigSourceId}_${location}' hsssGenomicLocationsFilter hsssReconstructSnpId +``` + +Expected: `hsssGenomicLocationsFilter:2` and `hsssReconstructSnpId:2`. + +> **`-F` is required.** Without it `grep` treats the pattern as a basic regex, where `{`/`}` +> are interval syntax, and it matches **nothing** — reporting `0` even for correct code. An +> earlier draft of this step omitted the flag, which made it a check that could only fail; +> the risk is an implementer "fixing" working code to satisfy it. Use `-F` for any grep of a +> literal Perl interpolation. + +- [ ] **Step 5: Check the script is still syntactically valid** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/bin && \ + perl -c hsssGenomicLocationsFilter +``` + +Expected: `hsssGenomicLocationsFilter syntax OK`. There is no fixture-driven test for this +script the way there is for `hsssReconstructSnpId` (`hsssTestSuite`, the only caller with +fixture data, is broken — spec §4), so a syntax check plus the greps is the available +verification. Do not claim more. + +- [ ] **Step 6: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter +git commit -m "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 " +``` + +--- + +### Task 2: Point the plugin at the variation directory + +`findOrganismDir` composes +`webSvcPath.replaceAll("PROJECT_GOES_HERE", projectId) + "/" + organismNameForFiles + searchDir` +and throws if the result is absent. The variation HSSS files live under `dnaseq`, not +`highSpeedSnpSearch`. + +**Files:** +- Modify: `ApiCommonWebService/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java:196-198` + +- [ ] **Step 1: Confirm the real layout on the webserver** + +```bash +cd ~/workspaces/agentic-veupath-dev && \ + ssh -o LogLevel=ERROR "$(python3 bin/resolve.py --profile profiles/plasmodb.yml --field host)" \ + 'ls /home/jbrestel/webserviceTest/Pfalciparum3D7/dnaseq/' +``` + +Expected: `readFreq20 readFreq40 readFreq60 readFreq80`. The `readFreq*` level is +appended separately by `FindPolymorphismsAbstractPlugin:106`, so `/dnaseq` is exactly the +piece `getSearchDir()` must supply. (Harmless ssh port-forward warnings may appear on +stderr; ignore them.) + +- [ ] **Step 2: Make the change** + +In `HighSpeedSnpSearchAbstractPlugin.java`, replace: + +```java + protected String getSearchDir() { + return "/highSpeedSnpSearch"; + } +``` + +with: + +```java + protected String getSearchDir() { + return "/dnaseq"; + } +``` + +Preserve the file's existing odd indentation and the double space after `return` — this +file is inconsistently formatted and reflowing it would bury a one-word change in +whitespace noise. + +- [ ] **Step 3: Verify the edit, and that the chip overrides are untouched** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch && \ + grep -n 'getSearchDir' -A 2 HighSpeedSnpSearchAbstractPlugin.java FindChipPolymorphismsPlugin.java FindChipSnpMajorAllelesPlugin.java +``` + +Expected: the abstract plugin returns `/dnaseq`; **both** chip plugins still return +`/highSpeedChipSnpSearch`. The chip path is equally dead but out of scope (spec §7), and +leaving it alone keeps the diff honest. + +- [ ] **Step 4: Note that this is unverifiable until the search exists** + +There is no way to exercise `findOrganismDir` without a search invoking the plugin (spec +§5 rung 5). Do not claim otherwise. Record in the task notes: *"correct by inspection +against the directory listing in Step 1; first exercised by the +`VariationsByIsolateGroup` search."* + +- [ ] **Step 5: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java +git commit -m "Read HSSS variation data from the dnaseq directory + +The variation HSSS files are laid out as /dnaseq/readFreq/, +not /highSpeedSnpSearch/readFreq/, 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 " +``` + +--- + +### Task 3: Rename the strain filter param to variation vocabulary + +`FindPolymorphismsAbstractPlugin:100` reads the samples selection under whatever name +`getStrainFilterParamName()` returns, and `:41` declares it **required**. So this string is +a contract with the model XML the follow-on spec will write. + +**Files:** +- Modify: `ApiCommonWebService/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java:41-43` + +- [ ] **Step 1: Make the change** + +Replace: + +```java + protected String getStrainFilterParamName() { + return "ngsSnp_strain_meta"; + } +``` + +with: + +```java + protected String getStrainFilterParamName() { + return "variation_sample_meta"; + } +``` + +Again, keep the existing (misaligned) indentation. + +- [ ] **Step 2: Confirm nothing else references the old name in this repo** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService && \ + grep -rn 'ngsSnp_strain_meta' . ; echo "exit: $?" +``` + +Expected: `exit: 1`. If `FindChipPolymorphismsPlugin` shows up, check what *it* returns — +it has its own override and must not be changed here. + +The old name still exists in `ApiCommonModel/Model/lib/wdk/model/questions/params/snpParams.xml`, +which is dead, commented-out snp XML. That is expected and out of scope; do not edit +`ApiCommonModel` in this plan. + +- [ ] **Step 3: Note the contract for the follow-on spec** + +Record in the task notes: *"the consuming `filterParam` must be named exactly +`variation_sample_meta`; a mismatch is rejected as a missing required parameter."* Like +Task 2 this has no runtime symptom until a search exists. + +- [ ] **Step 4: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindPolymorphismsPlugin.java +git commit -m "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 " +``` + +--- + +### Task 4: Fixture hygiene + +**These fixtures are consumed only by the two broken harnesses** (spec §4), so this task +fixes no test and unblocks nothing. It is worth ~15 lines so that whoever revives a harness +starts from fixtures consistent with the current ID convention rather than debugging a +stale one. **Do not run either harness, and do not report a test result from this task.** + +**Files:** +- Rename: `ApiCommonWebService/HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/` → `.../dnaseq/` +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/test/expected/genomicLocationFilter.txt` +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIds.txt` +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIdsAndSeqFilter.txt` + +- [ ] **Step 1: Rename the fixture directory to match Task 2** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/test/TestDB/Hsapiens123 && \ + git mv highSpeedSnpSearch dnaseq && ls dnaseq/ +``` + +Expected: `readFreq80`. Use `git mv` so the rename is tracked rather than showing as a +delete plus an add. + +- [ ] **Step 2: Rewrite the IDs in the three expected files** + +Every ID is `NGS_SNP..` and becomes `Variant__`. The +contig names contain no dots, so each ID has exactly one dot to convert: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/test/expected && \ +sed -i -E 's/^NGS_SNP\.([^.\t]+)\.([0-9]+)\t/Variant_\1_\2\t/' \ + genomicLocationFilter.txt \ + polymorphismSearchWithSourceIds.txt \ + polymorphismSearchWithSourceIdsAndSeqFilter.txt +``` + +The anchor `^` and the trailing `\t` confine the substitution to column 1, so the remaining +columns cannot be touched. + +- [ ] **Step 3: Verify the rewrite is complete and correct** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch/test/expected && \ + echo "--- remaining NGS_SNP (want none):" && { grep -c 'NGS_SNP' *.txt || true; } && \ + echo "--- new IDs:" && cut -f1 genomicLocationFilter.txt polymorphismSearchWithSourceIds.txt polymorphismSearchWithSourceIdsAndSeqFilter.txt +``` + +Expected: every file reports `0` for `NGS_SNP`, and the IDs are exactly these 11 — +`Variant_e99_2011`, `Variant_h103_30021` (from `genomicLocationFilter.txt`); +`Variant_a80_896`, `Variant_b86_13441`, `Variant_e99_2011`, `Variant_f100_23`, +`Variant_g102_4334`, `Variant_h103_30021`, `Variant_i104_3002`, `Variant_j201_54` (from +`polymorphismSearchWithSourceIds.txt`); `Variant_f100_23` (from the SeqFilter file). + +- [ ] **Step 4: Confirm the untouched fixtures stayed untouched** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService && git status --short HighSpeedSnpSearch/test/ +``` + +Expected: exactly the three `expected/*.txt` modifications plus the renamed directory. +`polymorphismSearch.txt`, `majorAlleles.txt`, and the `mergeStrains*.txt` files hold +pre-reconstruction output (contig **index** and location as separate columns) and must not +appear. + +- [ ] **Step 5: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add HighSpeedSnpSearch/test/ +git commit -m "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.. to +Variant__. + +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 " +``` + +--- + +### Task 5: Set the ID prefix (different repo) + +The `idPrefix` is **not** in `ApiCommonWebService`. The generated +`gus_home/config//highSpeedSnpSearch-config.xml` is regenerated on every build, +and the `.j2` template only renders `{{ highspeedsnpsearchconfig_idPrefix }}`. The value +lives in `ApiCommonWebsite`. + +**Files:** +- Modify: `ApiCommonWebsite/Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml:102` + +- [ ] **Step 1: Confirm the current value and that it is the only definition** + +```bash +cd ~/workspaces/plasmodb && grep -rn 'highspeedsnpsearchconfig_idPrefix' \ + --include=*.yml --include=*.yaml --include=*.j2 . +``` + +Expected exactly two hits: the template in +`ApiCommonWebService/WSFPlugin/lib/conifer/roles/conifer/templates/ApiCommonWebService/highSpeedSnpSearch-config.xml.j2` +(which consumes it) and +`ApiCommonWebsite/Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml:102` +(which defines it, as `NGS_SNP.`). + +- [ ] **Step 2: Make the change** + +In `ApiCommonWebsite/Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml`, replace: + +```yaml +highspeedsnpsearchconfig_idPrefix: NGS_SNP. +``` + +with: + +```yaml +highspeedsnpsearchconfig_idPrefix: Variant_ +``` + +Leave `highspeedsnpsearchconfig_jobsDir` (the line above) and the +`highspeedchipsnpsearchconfig_*` block (below) alone. + +- [ ] **Step 3: Verify** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebsite && \ + grep -n 'highspeedsnpsearchconfig_' Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml && \ + git diff --stat +``` + +Expected: `idPrefix: Variant_`, `jobsDir` unchanged, and `1 file changed, 1 insertion(+), 1 deletion(-)`. + +- [ ] **Step 4: Commit (in `ApiCommonWebsite`, not `ApiCommonWebService`)** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebsite +git add Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml +git commit -m "Prefix HSSS result IDs with Variant_ instead of NGS_SNP. + +Paired with the separator fix in ApiCommonWebService, this makes the HSSS +plugins emit VariationRecordClass source_ids. This is an ApiCommon cohort +default and so applies to every project in the cohort, which is safe because +no project has a live snp or chip search -- those imports are commented out +in the shared apiCommonModel.xml. + +Co-Authored-By: Claude Opus 5 " +``` + +--- + +### Task 6: Build, deploy, and re-verify the installed copy + +Task 1 proved the fix in the **source tree**. The plugins run the copy installed in +`$GUS_HOME/bin`, and the prefix comes from the **generated** config — neither of which +exists yet. + +**Files:** none (build and deploy only). + +- [ ] **Step 1: Build and install `ApiCommonWebService`** + +```bash +cd ~/workspaces/agentic-veupath-dev && \ + ssh -o LogLevel=ERROR "$(python3 bin/resolve.py --profile profiles/plasmodb.yml --field host)" \ + "bash -lc 'source /var/www/jbrestel.plasmodb.org/etc/setenv && bld ApiCommonWebService'" +``` + +Expected: `BUILD SUCCESSFUL`. This installs both the Java (WSFPlugin) and the Perl +(HighSpeedSnpSearch) components. Get the docroot for the `setenv` path with +`python3 bin/resolve.py --profile profiles/plasmodb.yml --field docroot` if the user prefix +is not `jbrestel`. + +Note `Test-Installation` is **not** in `build.xml`'s default depends list, so the +non-compiling JUnit module is not built. That is why this succeeds despite spec §4. + +- [ ] **Step 2: Re-run Task 1's test against the *installed* script** + +```bash +cd ~/workspaces/agentic-veupath-dev && \ + ssh -o LogLevel=ERROR "$(python3 bin/resolve.py --profile profiles/plasmodb.yml --field host)" \ + "bash -lc 'source /var/www/jbrestel.plasmodb.org/etc/setenv && \ + printf \"80\t896\t100\t25\t1\n\" | \$GUS_HOME/bin/hsssReconstructSnpId \ + \$PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/textData/contigIdToSourceId.dat \ + 1 Variant_ NULL 2>/dev/null'" +``` + +Expected: `Variant_a80_896 100 25 syn`. A dot here means the install did not pick up the +edit — check that mutagen synced the file before Step 1 ran. + +- [ ] **Step 3: Regenerate the site config so the new prefix lands** + +```bash +cd ~/workspaces/agentic-veupath-dev && \ + ssh -o LogLevel=ERROR "$(python3 bin/resolve.py --profile profiles/plasmodb.yml --field host)" \ + "bash -lc 'source /var/www/jbrestel.plasmodb.org/etc/setenv && \ + conifer configure --cohort ApiCommon --project PlasmoDB \ + --webapp-ctx plasmo.jbrestel --tomcat-webapp-ctx plasmo.jbrestel \ + --site-vars /var/www/jbrestel.plasmodb.org/etc/conifer_site_vars.yml'" +``` + +Expected: `PLAY RECAP ... failed=0`, with `highSpeedSnpSearch-config.xml` listed as `changed` +(or `ok` if already correct — the playbook is idempotent). + +**All four flags are required, and two were missing from an earlier draft of this step.** +`--cohort` is rejected outright (`conifer: error: --cohort is required for configure`). +Omitting `--tomcat-webapp-ctx` gets *further* — it regenerates most files, including our +target — and then fails on `log4j2.json` with +`AnsibleUndefinedVariable: 'tomcat_webapp_ctx' is undefined`, leaving `failed=1`. That is the +worst kind of half-success: the thing you were checking for *did* land, so a careless reading +calls it done while one config file silently went unregenerated. Always read the `PLAY RECAP`. + +Both `*-ctx` values are the webapp context, `plasmo.jbrestel` — confirmed against +`/usr/local/tomcat_instances/PlasmoDB/conf/Catalina/localhost/`, which holds one `.xml` per +context. For a different developer prefix, substitute accordingly. + +**`conifer install` is not needed as a separate step.** An earlier draft called for it. +`bld ApiCommonWebService` depends on `ApiCommonWebsite-Installation`, which already copies the +edited vars into `gus_home/lib/conifer/roles/conifer/vars/ApiCommon/default.yml` — verify +with `grep highspeedsnpsearchconfig_idPrefix` on that path if `configure` seems to render a +stale value. + +Conifer regenerates `model.prop` among other files, so confirm the model still loads +afterwards: + +```bash +cd ~/workspaces/agentic-veupath-dev && \ + ssh -o LogLevel=ERROR "$(python3 bin/resolve.py --profile profiles/plasmodb.yml --field host)" \ + "bash -lc 'source /var/www/jbrestel.plasmodb.org/etc/setenv && wdkXml -model PlasmoDB 2>&1 | tail -3'" +``` + +Expected: ends with `WDK Model resources released.` (a clean shutdown, i.e. the model loaded). + +Then reload so the runtime picks up the new config — WSF plugins read their property file at +init, so the prefix does not take effect until this: + +```bash +cd ~/workspaces/agentic-veupath-dev && bin/veup-build.sh plasmodb reload +``` + +Expected: `OK - Reloaded application at context path [/plasmo.jbrestel]`. + +- [ ] **Step 4: Confirm the prefix reached the generated config** + +```bash +cd ~/workspaces/agentic-veupath-dev && \ + ssh -o LogLevel=ERROR "$(python3 bin/resolve.py --profile profiles/plasmodb.yml --field host)" \ + 'grep idPrefix /var/www/jbrestel.plasmodb.org/gus_home/config/highSpeedSnpSearch-config.xml' +``` + +Expected: `Variant_`. This is the **only** check that Task 5 +took effect; nothing else reads that file until a search runs. + +> **Note the path.** This file sits directly in `gus_home/config/`, **not** in the +> per-project `gus_home/config/PlasmoDB/` subdirectory where `model.prop` and +> `model-config.xml` live. An earlier draft of this step had the `PlasmoDB/` path and would +> have failed with `No such file or directory` — which reads like "Task 5 didn't work" rather +> than "the path is wrong." + +- [ ] **Step 5: Record what remains unverified** + +Note explicitly in the task notes: the search directory (Task 2) and the filter param name +(Task 3) are **not exercised by anything in this plan**. Both are first tested by the +`VariationsByIsolateGroup` search. Do not describe this change as end-to-end verified. + +--- + +### Task 7: Close out the spec + +**Files:** +- Modify: `ApiCommonWebService/docs/superpowers/specs/2026-08-05-hsss-variation-plumbing-design.md` (the `Status:` line) +- Modify: `ApiCommonWebService/docs/superpowers/plans/2026-08-05-hsss-variation-plumbing.md` (this file) + +- [ ] **Step 1: Confirm the change set across both repos** + +```bash +for r in ApiCommonWebService ApiCommonWebsite; do + echo "=== $r"; git -C ~/workspaces/plasmodb/$r log --oneline origin/dnaseq-merge-experiments..HEAD 2>/dev/null || git -C ~/workspaces/plasmodb/$r log --oneline -5 + git -C ~/workspaces/plasmodb/$r status --porcelain +done +``` + +Expected: `ApiCommonWebService` has the four commits from Tasks 1-4 (plus the spec commit +`97ba19d`), `ApiCommonWebsite` has the one from Task 5, and both working trees are clean. + +- [ ] **Step 2: Mark the spec implemented** + +Change the spec's: + +```markdown +**Status:** approved +``` + +to: + +```markdown +**Status:** implemented 2026-08-05 — ID fix verified; the search directory (§3.1) and filter param name (§3.5) await the `VariationsByIsolateGroup` search +``` + +- [ ] **Step 3: Add an execution-outcome section to this plan** + +Append a short section recording: the four commit SHAs plus the `ApiCommonWebsite` one, +the actual before/after output of Task 1's test, whether `conifer configure` worked as +written or needed the rebuild fallback, and the explicit statement that Tasks 2 and 3 are +unverified pending the follow-on search. + +- [ ] **Step 4: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add docs/superpowers/specs/2026-08-05-hsss-variation-plumbing-design.md \ + docs/superpowers/plans/2026-08-05-hsss-variation-plumbing.md +git commit -m "Mark the HSSS variation plumbing spec implemented + +The ID fix is verified against both the source and installed scripts and the +prefix reaches the generated config. The search directory and filter param +name are correct by inspection but unexercised until the +VariationsByIsolateGroup search exists. + +Co-Authored-By: Claude Opus 5 " +``` + +--- + +## Execution outcome (2026-08-05) + +Executed subagent-per-task. Production code: **14 files, 17 insertions, 17 deletions**, and +`grep -rn jbrestel` over `WSFPlugin/src`, `HighSpeedSnpSearch/bin`, `HighSpeedSnpSearch/lib` +is clean. + +| task | commit | repo | +|---|---|---| +| 1 — ID separator | `7d1268a27` | `ApiCommonWebService` | +| 1b — second ID site | `5c80e5f37` | `ApiCommonWebService` | +| 2 — `/dnaseq` search dir | `141602e54` | `ApiCommonWebService` | +| 3 — `variation_sample_meta` | `e46356743` | `ApiCommonWebService` | +| 4 — fixture hygiene | `7f7d970` | `ApiCommonWebService` | +| 5 — `idPrefix: Variant_` | `1da7c420a` | **`ApiCommonWebsite`** | + +### Verified + +- **The ID fix, twice.** Source tree and installed copy both emit `Variant_a80_896` where + they previously emitted `Variant_a80.896`, from + `printf '80\t896\t100\t25\t1\n' | hsssReconstructSnpId .../contigIdToSourceId.dat 1 Variant_ NULL`. +- **The format is the right one:** `Variant_Pf3D7_01_v3_29514` exists in + `apidbtuning.VariationAttributes`. +- `bld ApiCommonWebService` → `BUILD SUCCESSFUL`, 1m44s. +- `conifer configure` → `failed=0`; the generated + `gus_home/config/highSpeedSnpSearch-config.xml` now reads + `Variant_`. +- `wdkXml -model PlasmoDB` still loads cleanly after Conifer regenerated `model.prop`, and + `WEBSERVICEMIRROR`/`PROJECT_ID` are unchanged. +- Webapp reloaded (`OK - Reloaded application at context path [/plasmo.jbrestel]`); the three + error logs were last written days-to-months ago, i.e. nothing was appended by any of this. + +### Not verified, and cannot be from this plan + +**The search directory (Task 2) and the filter param name (Task 3) were never exercised.** +Nothing invokes the plugin until a variation search exists. Both are correct by inspection — +Task 2 against the real `dnaseq/readFreq{20,40,60,80}` listing, Task 3 against +`FindPolymorphismsAbstractPlugin:41,100` — and both are first tested by +`VariationsByIsolateGroup`. This change is **not** end-to-end verified. + +No automated test suite ran, by design: both HSSS harnesses are broken independently of this +work (§4). + +### What execution changed about the plan + +Four errors, three in the plan and one in the spec, all found by executing rather than +reviewing: + +1. **A whole extra task.** Task 1's Step 5 grep asserted no other script composed IDs the + dotted way. It found `hsssGenomicLocationsFilter:51,67` — a *live* alternative pipeline + tail that `FindSnpsByGeneIdsPlugin:112` routes through, so `VariationsByGeneIds` would have + shipped returning **zero results with no error** while the other searches worked. Became + Task 1b. The step earned its keep by being wrong. +2. **A check that could only fail.** Task 1b's Step 4 used `grep -c` without `-F` on a pattern + containing `${...}`; BRE treats the braces as interval syntax and matches nothing, so it + reported `0` even for already-correct code. The hazard is an implementer "fixing" working + code to satisfy it. +3. **A wrong config path.** `highSpeedSnpSearch-config.xml` lives in `gus_home/config/`, not + `gus_home/config/PlasmoDB/`. The original grep would have failed with + `No such file or directory`, reading as "Task 5 didn't work." +4. **An incomplete `conifer` invocation.** It needs `--cohort`, `--project`, `--webapp-ctx` + **and** `--tomcat-webapp-ctx`. Missing `--cohort` is refused outright; missing + `--tomcat-webapp-ctx` regenerates most files — *including the one being checked* — then + fails on `log4j2.json`. A half-success where the verification passes while a config + silently goes unregenerated. Read the `PLAY RECAP`. + +Also: `conifer install` is unnecessary — `bld ApiCommonWebService` depends on +`ApiCommonWebsite-Installation`, which installs the Conifer vars. + +### Findings deferred rather than acted on + +Three second sites were found and consciously left alone; each is recorded with reasoning in +"Deliberately not in this plan" below: `FindMajorAllelesPlugin`'s `ngsSnp_strain_meta_a`/`_m` +param family (deferred to the two-isolate-groups spec, **must not be forgotten there**), +`hsssCopyFilesToWebSvcDir`'s write path, the three-way duplication of the ID format, and the +unconditional stderr echo in `hsssReconstructSnpId:42`. + +## Deliberately not in this plan + +- **Any WDK model XML.** `VariationsByIsolateGroup` is a separate spec in `ApiCommonModel` + and is the first consumer of everything here. +- **Reviving either HSSS test harness** — spec §4. The JUnit module does not compile and + the shell suite passes the wrong argument count. Both are larger jobs than this change. +- **Renaming** `hsssReconstructSnpId`, the `highspeedsnpsearch` package, or the + `Snp`-flavoured class names — spec §3.4. Ten files plus installed script names; bundling + a rename with a behavioural fix makes the diff unreviewable. +- **Deleting the dead chip/major-alleles/gene-chars plugins** — spec §7. +- **Making `idPrefix` per-plugin.** Only needed if a second consumer with a different ID + convention appears. +- **Renaming `FindMajorAllelesPlugin`'s param constants.** `:20` and `:24` hardcode + `ngsSnp_strain_meta_a` and `ngsSnp_strain_meta_m`, with matching `_wiz` variants in + `ApiCommonModel`'s `sharedParams.xml` — a family of four snp-vocabulary names. That plugin + extends `HighSpeedSnpSearchAbstractPlugin` directly, so it has no + `getStrainFilterParamName()` to override; the constants are its own contract, listed in its + `getRequiredParameterNames()` (`:48-49`) and read at `:81` and `:97`. Same failure mode as + Task 3: a name mismatch is rejected as a missing required parameter. + + Flagged during Task 3 and **deliberately deferred** to the spec for + `VariationsByTwoIsolateGroups`, the search it serves (`NgsSnpsByTwoIsolateGroups` and + `...Wiz`). Reasoning: that search needs per-strain data the variation pipeline does not yet + have, so it is the furthest out of the four ports, and renaming now would mean choosing a + param-family name before designing the search that uses it. Note the asymmetric `_a`/`_m` + suffixes — the prompts read "Set A Isolates" / "Set B Isolates", so `_m` appears to be a + typo for `_b`, and that should be decided deliberately rather than mirrored. + + > **Whoever writes that spec must include this rename**, or the new variation XML inherits + > snp naming permanently. + +- **Collapsing the ID-composition sites into one helper.** After Task 1b the separator is + hardcoded in three places (`hsssReconstructSnpId:42-43` and + `hsssGenomicLocationsFilter:51,67`), all of them re-deriving a format that + `VariationRecordClass` owns. That duplication is how the second site got missed in the + first place, so a shared helper is genuinely the right long-term shape — but it is a + refactor of live pipeline code with no test harness to catch a mistake, which is a worse + bet right now than three verified one-liners. Worth a follow-up issue. +- **The unconditional STDERR echo in `hsssReconstructSnpId:42`.** Every composed ID is + printed to stderr as well as stdout, and the stderr copy **bypasses** the stdout branch's + sequence/location filter — so it emits rows the search deliberately excluded. Pre-existing + and untouched here, but at 4.4M variants it will make remote logs noisy and could mislead + anyone reading them during debugging. Flagged during Task 1b; not fixed because changing + output streams in untested pipeline code is its own change. +- **Populating the production HSSS directories** under + `/var/www/Common/apiSiteFilesMirror/webServices//build-/`. This change is + verified against the test copy in `/home/jbrestel/webserviceTest`; production placement is + a data-deployment task. +- **`HighSpeedSnpSearch/bin/hsssCopyFilesToWebSvcDir:36`**, which writes into + `/highSpeedSnpSearch` and so now disagrees with what the plugin reads. Flagged + during Task 2. Left alone deliberately: it is a run-once snp-era prototype copier, not a + deployment path — it hardcodes a `/eupath/data/htsSnpsPrototype/heterozygosityEnabled/` + source, carries a hardcoded project→organism table the author named `$stupidHash`, and + `die`s if the target already exists. It is not what produced the `dnaseq` directories and + is not part of the variation pipeline. Recorded here so nobody later "fixes" it or mistakes + it for how variation data reaches the webserver. **The open question it does raise — + what populates production `/dnaseq` — is the data-deployment item above.** +- **The `webServiceMirror` override** needed to point the plugin at the test files (spec + §5.1) — that belongs with the search that invokes the plugin, and the symlink bridging the + missing `PlasmoDB/build-70` levels has already been created by hand. diff --git a/docs/superpowers/plans/2026-08-07-hsss-gene-stats-fix.md b/docs/superpowers/plans/2026-08-07-hsss-gene-stats-fix.md new file mode 100644 index 00000000..7bbc884c --- /dev/null +++ b/docs/superpowers/plans/2026-08-07-hsss-gene-stats-fix.md @@ -0,0 +1,1275 @@ +# HSSS Gene Statistics Fix Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Make the two reported statistics in the `GenesByNgsSnps` search compute what their labels claim — SNVs per kb of coding sequence, and a site-normalized dN/dS — and relabel every string that described the old behaviour. + +**Architecture:** Numerators stay sample-set-dependent and are computed by HSSS from the selected samples. Denominators (CDS length, Nei-Gojobori synonymous/nonsynonymous site counts) are gene properties already derived once from the genetic code in `apidbtuning.GeneVariationSummary`; the plugin reads them in the SQL it already runs and passes them to the perl filter through `geneLocations.txt`. Both searches then rest on one definition of a synonymous site. + +**Tech Stack:** Java (WSF plugin), Perl (HSSS stream filter), bash (HSSS test suite), WDK model XML, PostgreSQL. + +**Spec:** `ApiCommonWebService/docs/superpowers/specs/2026-08-07-hsss-gene-stats-fix-design.md` + +--- + +## Before you start + +**Do NOT create a git worktree for this work.** These repos are checked out at +`~/workspaces/plasmodb/` and a mutagen session carries that working tree to the build +host. A worktree lives outside the synced path, so anything built or run on the remote +would be testing stale code. Edit in place. + +Branches already exist and carry prior related commits — use them, do not branch again: + +| repo | path | branch | +|---|---|---| +| ApiCommonWebService | `~/workspaces/plasmodb/ApiCommonWebService` | `feature/hsss-noncoding-class` | +| ApiCommonModel | `~/workspaces/plasmodb/ApiCommonModel` | `dnaseq-merge-experiments` | + +Two environment facts you will need: + +- The appDb is reachable read-only at `psql -h localhost -p 5432 -d unidb_shu_a`. +- Remote builds run through `bin/veup-build.sh plasmodb wb model` from + `~/workspaces/agentic-veupath-dev`. **Flags go before the profile name.** + +**Deployment coupling:** the perl output grows from 8 fields to 9 and the Java asserts on +the count. The two repos must ship together. Do not deploy one without the other. + +--- + +## hsssTestSuite tests the INSTALLED code, not your working tree + +Found the hard way during Task 3, and it will bite again on Task 4. + +The suite is invoked by path out of `$PROJECT_HOME`, but every tool it calls resolves off +`PATH` into `$GUS_HOME/bin` — the installed copy. Editing +`HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter` in the checkout and re-running the +suite tests the OLD code and can print a completely undeserved `matched`. + +Before trusting any suite result after editing a tool, install it: + +```bash +ssh cedar 'source /var/www/jbrestel.plasmodb.org/etc/setenv + cp $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/bin/ $GUS_HOME/bin/ + chmod +x $GUS_HOME/bin/' +``` + +and confirm with `md5sum` that installed and checkout agree. **A `matched` that arrives +without an intervening install is meaningless.** + +The Java plugin has the same property in a worse form: it is a jar loaded by the plugin +runner, so `mvn compile` proves only that it builds. Task 4's real verification is Task 9, +against the live instance, after the jar is rebuilt and deployed. + +--- + +## Known remaining breakage in hsssTestSuite (deliberately NOT fixed here) + +Found while repairing the suite for Task 0. Recorded so the next person does not +rediscover them; each is out of scope for the statistics work. + +- **`hsssTestSuite:135` — the majorAlleles stage is dead.** It calls + `hsssGenerateMajorAllelesScript` with 12 arguments where the generator requires 14-15, + so it dies in `usage()`. Same root cause as the four sites Task 0 fixed: the missing + `hsssReconstructSnpId Variant_ NULL` triple. Not fixed because + `expected/majorAlleles.txt` holds raw contig indices (`99 2011 C ...`) rather than + `Variant_e99_2011` ids, so it predates the reconstruct step being wired in at all — + repairing it means re-baselining a stage whose output has never been validated, which + is a larger and separate piece of work. + +- **`FindGenesWithChipSnpCharsPlugin` shares the filter and still expects 8 fields.** + It returns the same `getGenerateScriptName()` (`hsssGenerateGeneCharsScript`), so it + receives the same 9-field output against a `parts.length != 8` assertion. Dormant, not + dangerous: its `processQuery` carries `includeProjects="TODO??"` — excluded from every + project — and its question in `geneQuestions.xml` is inside a comment block. If anyone + revives the chip search it fails loudly with a column-count exception rather than + silently producing wrong numbers, which is why it is left alone here. + +- **`expected/majorAlleles.txt` also has an off-by-one in its product letters.** + Reported by the Task 0 agent and not independently confirmed: each product letter in + the fixture is exactly one higher than the code emits (`J`/`I`, `U`/`T`, `O`/`N` — + 74/73, 85/84, 79/78). If that holds it is a second defect layered on the missing-args + one, and re-baselining without understanding it would bless an off-by-one. + +- **`expected/polymorphismSearch.txt` is orphaned and self-contradictory.** Referenced + nowhere in the suite, and it disagrees with its sibling + `expected/polymorphismSearchWithSourceIds.txt` about which SNPs are non-synonymous (it + marks contigs 80/99/102; the sibling marks 80/102/103). Two stale files from different + eras. Probably wants deleting, but deleting a fixture is not a statistics fix. + +- **`hsssReconstructSnpId:59` documents the wrong encoding.** Its usage text says + `product_class(-1=noncoding,0=syn,1=nonsyn,2=nonsense)`; the code immediately above it + maps `0` to `non-coding`, `1` to `syn`, `2` to `non-syn`, and negatives to + `has stop codon`. The same class of defect as the labels this plan corrects — + documentation describing behaviour the code does not have. + +--- + +### Task 0: Repair hsssTestSuite (REPLACES the original Task 1) + +**This task was added after the plan was written.** The original Task 1 assumed the suite +ran the geneChars filter and discarded the result. It does not: `extractArgs` consumes +five arguments and `getFinalCommandString` unpacks fourteen, but the suite supplied +eleven, so every argument shifted left by three and the gene locations file was never +passed — the filter received the literal string `5` in its place. The stage has never +produced meaningful output, which is why its `diff` was commented out. + +Three sub-parts, executed in this order: + +- **0a — fix the arguments.** Insert `strainsList.txt hsssReconstructSnpId Variant_ NULL` + in the slot after `strains_list_file` at all four call sites + (`hsssTestSuite:41,59,76,93`). Sites 41/59/76 genuinely use `idPrefix`/`idSuffix`; + geneChars emits gene ids and never reads them, so there they are inert positional + filler. Committed as `5dfe387`, touching `hsssTestSuite` only. + +- **0b — characterize the stale baselines.** Three expected files diverge from actual + output in exactly two classes, both traced and both approved: `%d` -> `%.1f` on the + percentage columns (`d3771af`, 2014-07-26) and the product-class column going from a + boolean `y`/blank to the four-value label set (`f1ac0d9`, 2014-08-19). Neither commit + updated `test/expected/`. Row counts, field counts and column 1 are identical + throughout; there are no unexplained differences. + +- **0c/0d — extend the fixture, then baseline once.** The fixture produces only `syn` and + `non-syn`; it contains no non-coding position and no stop codon, so + `nonCodingCount` and `nonsenseCount` are 0 for every gene. Re-baselining before fixing + that would produce a green suite that never executes the paths this branch changes. + So: add a non-coding position and a stop-codon position to `test/textData/strain*.txt` + and `referenceGenome.txt`, with at least one falling inside a gene span in + `geneFilters.txt`, THEN regenerate all four expected files together — the three above + plus the new `geneCharsFilter.txt`. Uncomment the geneChars assertion. Prove it can + fail by corrupting the expected file and confirming a non-zero exit. + + geneChars needs widened arguments to emit anything: the suite's current + `coding 2 5 .1 .9 3 1000` yields empty output; use `all 0 -1 0 -1 0 -1`. + + Also convert `test/textData/geneFilters.txt` to unix line endings — `chomp` strips only + `\n`, so the trailing `\r` lands on the last field, and Task 2 appends numeric columns + after it. + +--- + +### Task 1: (SUPERSEDED — folded into Task 0d) + +Kept for numbering. The baseline capture and assertion uncommenting described below now +happen in Task 0d, against the extended fixture. Read this section for the +`hsssTestSuite` edit and the fail-proof step, but do not execute it separately. + +The suite runs `hsssGeneCharacteristicsFilter` today but asserts nothing — the `diff` is +commented out and the expected file does not exist. Establish the baseline BEFORE +changing behaviour, so the later tasks have something to break. + +**Files:** +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/bin/hsssTestSuite:99-105` +- Create: `ApiCommonWebService/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt` +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/test/textData/geneFilters.txt` + +- [ ] **Step 1: Convert the fixture to unix line endings** + +`geneFilters.txt` currently has CRLF endings. `chomp` strips only `\n`, so today the +trailing `\r` lands on the last field (`gene_source_id`) and is invisible. Once columns +are appended in Task 3 it would land on a numeric field and produce warnings. + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +sed -i 's/\r$//' test/textData/geneFilters.txt +cat -A test/textData/geneFilters.txt +``` + +Expected: five lines ending in `$` with no `^M`. + +- [ ] **Step 2: Run the suite and capture current output as the expected file** + +`hsssTestSuite` takes the working directory as its one argument, so choose it rather +than hunting for it. + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +export PROJECT_HOME=~/workspaces/plasmodb +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest 2>&1 | tail -20 +cp /tmp/hsssTest/geneChars_result.txt test/expected/geneCharsFilter.txt +cat test/expected/geneCharsFilter.txt +``` + +Expected: tab-delimited rows of exactly 8 fields: +`geneId density dndsRatio synCount nonSynCount nonCodingCount nonsenseCount snpsCount` + +If the file is empty, the fixture SNP stream produces no gene hits under the suite's +`coding 2 5 .1 .9 3 1000` arguments. In that case widen the arguments on +`hsssTestSuite:93` to `all 0 -1 0 -1 0 -1` so at least one gene is emitted, and use that +output. Record which arguments you used in the commit message. + +- [ ] **Step 3: Uncomment the assertion** + +In `bin/hsssTestSuite`, replace lines 99-105: + +```bash +#echo "Comparing expected runGeneChars output with result..." +#diff $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt geneChars_result.txt +#diffStat=$? +#if [ $diffStat != 0 ]; then +# exit -1 +#fi +#echo "matched" +``` + +with: + +```bash +echo "Comparing expected runGeneChars output with result..." +diff $PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt geneChars_result.txt +diffStat=$? +if [ $diffStat != 0 ]; then + exit -1 +fi +echo "matched" +``` + +- [ ] **Step 4: Run the suite to verify it passes** + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest; echo "exit=$?" +``` + +Expected: `matched` printed for geneChars, `exit=0`. + +- [ ] **Step 5: Verify the test actually detects a change** + +Temporarily corrupt the expected file and confirm the suite fails — a test that cannot +fail is not a test. + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +echo "bogus 1 2 3 4 5 6 7" >> test/expected/geneCharsFilter.txt +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest >/dev/null 2>&1; echo "exit=$?" +git checkout test/expected/geneCharsFilter.txt +``` + +Expected: `exit=255` (the script's `exit -1`). Then the checkout restores the file. + +- [ ] **Step 6: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add HighSpeedSnpSearch/bin/hsssTestSuite HighSpeedSnpSearch/test/expected/geneCharsFilter.txt HighSpeedSnpSearch/test/textData/geneFilters.txt +git commit -m "test: assert hsssGeneCharacteristicsFilter output instead of only running it + +The geneChars diff in hsssTestSuite was commented out and its expected file +never existed, so the filter had no regression coverage. Captures current +output as the baseline before changing the statistics it computes. + +Also converts geneFilters.txt to unix line endings: chomp strips only \\n, so +the trailing \\r was landing on the last field." +``` + +--- + +### Task 2: Widen the test fixture to carry normalizers + +The fixture gains the three columns the filter will read in Task 3. Doing this first, +while the filter still ignores them, proves the parser change in Task 3 is what makes the +numbers move. + +**Files:** +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/test/textData/geneFilters.txt` + +- [ ] **Step 1: Append cdsLen, synSites, nonsynSites to each fixture row** + +Write `test/textData/geneFilters.txt` as exactly this (tab-delimited, unix endings). The +values are chosen to exercise every branch: + +``` +e99 1000 3000 g1 1200 300 900 +f100 5 700 g2 600 150 450 +g102 3001 40000 g3 0 0 0 +h103 30021 40000 g4 900 0 675 +j201 20 50 g5 300 75 225 +``` + +- `g1`, `g2`, `g5`: normal coding genes — both densities and the ratio defined. +- `g3`: `cdsLen`/`synSites`/`nonsynSites` all zero — stands for a non-coding gene, where + the tuning table has NULL. CDS density and ratio must be blank; span density populated. +- `g4`: `synSites` zero but `nonsynSites` non-zero — ratio blank, both densities defined. + +Note for the engineer: the plugin writes an empty string for a SQL NULL, and the perl +tests these fields for truth, so empty string and `0` behave identically. `0` is used in +the fixture because it survives a round trip through `sort` and is visible in `cat`. + +- [ ] **Step 2: Run the suite to verify it still passes** + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest; echo "exit=$?" +``` + +Expected: `matched`, `exit=0`. The filter splits into four scalars and discards the rest, +so output is unchanged. If this fails, the fixture has a whitespace error — check with +`cat -A`. + +- [ ] **Step 3: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add HighSpeedSnpSearch/test/textData/geneFilters.txt +git commit -m "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." +``` + +--- + +### Task 3: Compute the corrected statistics in the perl filter + +**Files:** +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter:23,68,83-121,123-146` +- Modify: `ApiCommonWebService/HighSpeedSnpSearch/test/expected/geneCharsFilter.txt` + +- [ ] **Step 1: Parse the three new columns in both split sites** + +There are TWO places the gene locations line is split, and they must agree. Line 23: + +```perl +my ($filterContigId, $filterStart, $filterEnd, $filterGeneId) = split(/\t/, $geneLocationLine); +``` + +becomes: + +```perl +my ($filterContigId, $filterStart, $filterEnd, $filterGeneId, + $filterCdsLen, $filterSynSites, $filterNonsynSites) = split(/\t/, $geneLocationLine); +``` + +Line 68, inside the `while` that advances past genes: + +```perl + ($filterContigId, $filterStart, $filterEnd, $filterGeneId) = split(/\t/, $geneLocationLine); +``` + +becomes: + +```perl + ($filterContigId, $filterStart, $filterEnd, $filterGeneId, + $filterCdsLen, $filterSynSites, $filterNonsynSites) = split(/\t/, $geneLocationLine); +``` + +Missing the second one is the likeliest bug in this task: the first gene would get correct +statistics and every later gene would silently reuse the first gene's normalizers. + +- [ ] **Step 2: Replace the statistic calculations in processGene** + +In `sub processGene`, replace lines 88-90: + +```perl + my $nonCodingCount = $snpsCount - $codingCount; + my $dnds = $synCount? $nonSynCount / $synCount : undef; + my $density = $snpsCount / (($filterEnd - $filterStart) / 1000); +``` + +with: + +```perl + my $nonCodingCount = $snpsCount - $codingCount; + + # Densities. cdsDensity is coding variants over coding length, which is what this + # search has always CLAIMED to report; spanDensity is what it actually reported, kept + # under a name that admits it. cdsLen is empty for a gene with no coding sequence. + my $cdsDensity = $filterCdsLen ? 1000 * $codingCount / $filterCdsLen : undef; + my $spanDensity = 1000 * $snpsCount / ($filterEnd - $filterStart); + + # dN/dS, each count normalized by the number of sites of its class (Nei-Gojobori, + # computed from the genetic code in apidbtuning.GeneVariationSummary). Without this + # normalization the ratio carries the genome's codon bias: the pooled synonymous-site + # fraction in pfal3D7 is 17.49%, not the textbook ~25%, worth 1.43x on every gene. + # + # defined($dn), NOT $dn: a gene with zero nonsynonymous variants has dN = 0 and a real + # ratio of 0, which is a strong purifying-selection signal and exactly what someone + # filtering a low range wants. Truth-testing $dn would silently drop those genes. + # $ds IS truth-tested, because zero there is a division by zero, not a result. + my $dn = $filterNonsynSites ? $nonSynCount / $filterNonsynSites : undef; + my $ds = $filterSynSites ? $synCount / $filterSynSites : undef; + my $dnds = (defined($dn) && $ds) ? $dn / $ds : undef; +``` + +- [ ] **Step 3: Point the dN/dS filter at the new value** + +Replace lines 108-114: + +```perl + if ($dndsMin || $dndsMax != -1) { + if ($synCount == 0 && $nonSynCount != 0) { + return 0 unless $dndsMax == -1; + } else { + return 0 if $dnds < $dndsMin || ($dndsMax != -1 && $dnds > $dndsMax); + } + } +``` + +with: + +```perl + if ($dndsMin || $dndsMax != -1) { + # An undefined ratio cannot be shown to be in range, so the gene is excluded - but + # only because the user narrowed this filter. Leaving it alone (min 0, max -1) skips + # this block entirely, which is what keeps non-coding genes in the result. + if (!defined($dnds)) { + return 0; + } else { + return 0 if $dnds < $dndsMin || ($dndsMax != -1 && $dnds > $dndsMax); + } + } +``` + +Note this subsumes the old `$synCount == 0 && $nonSynCount != 0` special case: zero +synonymous variants now yields `$ds == 0`, so `$dnds` is undef and the gene is excluded +whenever the filter is engaged — the same outcome the old branch produced, reached by one +rule instead of two. + +- [ ] **Step 4: Point the density filter at CDS density** + +Replace lines 116-118: + +```perl + if ($densityMin || $densityMax != -1) { + return 0 if ($density < $densityMin || ($densityMax != -1 && $density > $densityMax)); + } +``` + +with: + +```perl + # Filters CDS density only. Span density is reported but not filterable, to avoid a + # fifteenth and sixteenth param on a form that already carries fourteen. + if ($densityMin || $densityMax != -1) { + return 0 if (!defined($cdsDensity)); + return 0 if ($cdsDensity < $densityMin || ($densityMax != -1 && $cdsDensity > $densityMax)); + } +``` + +- [ ] **Step 5: Emit nine fields** + +Replace line 120: + +```perl + print STDOUT join("\t", $filterGeneId, sprintf("%.2f",$density), $synCount ? sprintf("%.2f",$dnds) : undef, $synCount, $nonSynCount, $nonCodingCount, $nonsenseCount, $snpsCount) . "\n"; +``` + +with: + +```perl + print STDOUT join("\t", + $filterGeneId, + defined($cdsDensity) ? sprintf("%.2f", $cdsDensity) : '', + sprintf("%.2f", $spanDensity), + defined($dnds) ? sprintf("%.4f", $dnds) : '', + $synCount, $nonSynCount, $nonCodingCount, $nonsenseCount, $snpsCount) . "\n"; +``` + +Two changes beyond the added column. `undef` in a `join` produces an uninitialized-value +warning and an empty string; `''` is explicit. And `%.4f` rather than `%.2f` for the +ratio, because site normalization divides by site counts in the hundreds, so meaningful +values now sit well below 1 where two decimals would collapse them. + +- [ ] **Step 6: Update the usage text** + +Replace lines 131-144 of the `usage` sub. Note the existing text says the snps_min/max +are a PERCENT — they are compared against raw counts in the code, so that has always been +wrong and is corrected here. + +``` + - gene_locations_filter_file: tab delimited: contig_source_id, start, end, gene_source_id, cds_length, syn_sites, nonsyn_sites. Must be sorted by location. The last three may be empty for a gene with no coding sequence; the statistics that need them are then reported empty. + - snp_class: all, coding, noncoding, synonymous, nonsynonymous, nonsense + - snps_min: min NUMBER of SNPs in the gene that belong to the specified class + - snps_max: max NUMBER of SNPs in the gene that belong to the specified class + - dnds_min: min site-normalized dN/dS ratio + - dnds_max: max site-normalized dN/dS ratio + - density_min: min coding SNPs per kb of CDS + - density_max: max coding SNPs per kb of CDS + + - snp_search_result: tab_delimited where first column is contig index and second is gene location. + +Replaces the first two columns of snp_search_result with a single column that is the concatenation of the contig_source_id-location, ie, a snp source id. + +Outputs these columns (tab delim): geneId cdsDensity spanDensity dndsRatio synCount nonSynCount nonCodingCount nonsenseCount snpsCount +``` + +- [ ] **Step 7: Run the suite to verify it FAILS** + +**Correction to this plan, found during Task 0:** the exit code is `1`, NOT `255`. +`hsssTestSuite` runs under `set -e`, so it dies at the `diff` itself and the +`diffStat=$?` / `exit -1` blocks after every diff are dead code. Worse, a pre-existing +failure in the majorAlleles stage (see Known remaining breakage) means the suite exits +non-zero even on success. **Exit code alone cannot tell you whether geneChars passed.** + +The discriminator is WHERE the suite stops and whether it printed `matched` for +geneChars: + +| geneChars | prints | stops at | exit | +|---|---|---|---| +| passing | `matched` | majorAlleles | 1 | +| failing | the diff, no `matched` | geneChars | 1 | + +So every check below greps the output rather than testing `$?`. + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest >/dev/null 2>&1; echo "exit=$?" +``` + +Expected: `exit=255`. The output now has 9 fields and different values, so the Task 1 +baseline must reject it. If this passes, the assertion is not wired up — go back to +Task 1 Step 5. + +- [ ] **Step 8: Inspect the new output and check it by hand** + +Run: + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest >/dev/null 2>&1 +cat -A /tmp/hsssTest/geneChars_result.txt +``` + +For each emitted row verify by hand against the fixture in Task 2: +- `cdsDensity` = `1000 * codingCount / cdsLen`, empty when `cdsLen` is 0 +- `spanDensity` = `1000 * snpsCount / (end - start)`, always populated +- `dnds` = `(nonSynCount/nonsynSites) / (synCount/synSites)`, empty when `synSites` is 0 +- field count is 9 on every row + +Do not proceed until each row checks out. This hand-check is the real test; the expected +file only locks it in. + +- [ ] **Step 9: Re-baseline the expected file** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService/HighSpeedSnpSearch +cp /tmp/hsssTest/geneChars_result.txt test/expected/geneCharsFilter.txt +rm -rf /tmp/hsssTest && mkdir -p /tmp/hsssTest +bin/hsssTestSuite /tmp/hsssTest >/dev/null 2>&1; echo "exit=$?" +``` + +Expected: `exit=0`. + +- [ ] **Step 10: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter HighSpeedSnpSearch/test/expected/geneCharsFilter.txt +git commit -m "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." +``` + +--- + +### Task 4: Supply the normalizers from the plugin + +**Files:** +- Modify: `ApiCommonWebService/WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java` + +- [ ] **Step 1: Add the new column constant** + +After the existing `COLUMN_DENSITY` declaration, add: + +```java + public static final String COLUMN_SPAN_DENSITY = "span_snp_density"; +``` + +Leave `COLUMN_DENSITY = "cds_snp_density"` alone. The name becomes accurate rather than +needing a rename that would break `attributesList summary=` and saved strategies. + +- [ ] **Step 2: Extend the gene locations SQL** + +In `initForBashScript`, replace the `String sql = ...` assignment: + +```java + String sql = "select g.sequence_id, g.start_min, g.end_max, g.source_id" + newline + + "from webready.GeneAttributes_p g " + newline + "where g.source_id is not null" + newline + + " and g.organism = '" + organism + "'"; +``` + +with: + +```java + // LEFT JOIN, never inner: GeneVariationSummary holds one row per gene that has + // COHORT variants (5,579 of 5,720 annotated pfal3D7 genes), and a gene missing + // from it must still get a locations line and still report its counts. Only its + // normalized statistics come back empty. + // + // These three columns are gene properties derived from the genetic code, not + // from any sample set, which is why reading them here is sound: the numerators + // stay sample-set-dependent and HSSS still computes them. + String sql = "select g.sequence_id, g.start_min, g.end_max, g.source_id," + newline + + " gvs.cds_length, gvs.syn_sites, gvs.nonsyn_sites" + newline + + "from webready.GeneAttributes_p g " + newline + + "left join apidbtuning.GeneVariationSummary gvs" + newline + + " on gvs.gene_source_id = g.source_id" + newline + + " and gvs.project_id = g.project_id" + newline + + "where g.source_id is not null" + newline + + " and g.organism = '" + organism + "'"; +``` + +- [ ] **Step 3: Write the new columns to the locations file** + +Replace the result-set loop body: + +```java + while (rs.next()) { + String seqId = rs.getString(1); + String start = rs.getString(2); + String end = rs.getString(3); + String geneId = rs.getString(4); + bw.write(seqId + "\t" + start + "\t" + end + "\t" + geneId); + bw.newLine(); + } +``` + +with: + +```java + while (rs.next()) { + String seqId = rs.getString(1); + String start = rs.getString(2); + String end = rs.getString(3); + String geneId = rs.getString(4); + // getString returns null for a SQL NULL; the filter tests these for truth, + // so an empty string reads as "no normalizer" exactly like a zero would. + String cdsLen = rs.getString(5) == null ? "" : rs.getString(5); + String synSites = rs.getString(6) == null ? "" : rs.getString(6); + String nonsynSites = rs.getString(7) == null ? "" : rs.getString(7); + bw.write(seqId + "\t" + start + "\t" + end + "\t" + geneId + "\t" + + cdsLen + "\t" + synSites + "\t" + nonsynSites); + bw.newLine(); + } +``` + +The `apiSortNoLocale -k 1,1 -k 2,2n` that follows is unaffected — the new columns are +appended after the sort keys. + +- [ ] **Step 4: Extend the unit-test filter rows** + +In the same method, replace the `"unit test"` branch's `testFilters`: + +```java + String[] testFilters = new String[] { "e99\t1000\t3000\tg1", "f100\t500\t700\tg2", + "h103\t30021\t40000\tg3", "j201\t20\t50\tg4" }; +``` + +with: + +```java + String[] testFilters = new String[] { + "e99\t1000\t3000\tg1\t1200\t300\t900", + "f100\t500\t700\tg2\t600\t150\t450", + "h103\t30021\t40000\tg3\t\t\t", // no coding sequence: normalizers empty + "j201\t20\t50\tg4\t900\t0\t675" }; // zero synonymous sites: ratio undefined +``` + +- [ ] **Step 5: Accept nine fields and map the new column** + +Replace `makeResultRow` in full: + +```java + protected String[] makeResultRow(String[] parts, Map columns, String projectId) + throws PluginModelException { + if (parts.length != 9) + throw new PluginModelException("Wrong number of columns in results file. Expected 9, found " + + parts.length); + + String[] row = new String[12]; + row[columns.get(COLUMN_GENE_SOURCE_ID)] = parts[0]; + row[columns.get(COLUMN_SOURCE_ID)] = null; + row[columns.get(COLUMN_PROJECT_ID)] = projectId; + row[columns.get(COLUMN_MATCHED_RESULT)] = "Y"; + row[columns.get(COLUMN_DENSITY)] = parts[1]; + row[columns.get(COLUMN_SPAN_DENSITY)] = parts[2]; + row[columns.get(COLUMN_DNDS)] = parts[3]; + row[columns.get(COLUMN_SYN)] = parts[4]; + row[columns.get(COLUMN_NONSYN)] = parts[5]; + row[columns.get(COLUMN_NONCODING)] = parts[6]; + row[columns.get(COLUMN_NONSENSE)] = parts[7]; + row[columns.get(COLUMN_TOTAL)] = parts[8]; + return row; + } +``` + +`new String[12]`, up from 11, because the row array is indexed by the column map and one +column was added. + +- [ ] **Step 6: Declare the new column to the framework** + +Replace `getColumns`: + +```java + public String[] getColumns(PluginRequest request) { + return new String[] { COLUMN_GENE_SOURCE_ID, COLUMN_PROJECT_ID, COLUMN_DENSITY, COLUMN_SPAN_DENSITY, + COLUMN_DNDS, COLUMN_SYN, COLUMN_NONSYN, COLUMN_NONCODING, COLUMN_NONSENSE, COLUMN_TOTAL }; + } +``` + +- [ ] **Step 7: Verify it compiles** + +Run: + +```bash +ssh cedar 'bash -lc "cd /var/www/PlasmoDB/plasmo.jbrestel/project_home/ApiCommonWebService && mvn -q -pl WSFPlugin -am compile 2>&1 | tail -20"' +``` + +Expected: no output, or `BUILD SUCCESS`. If the module coordinates differ, fall back to +`mvn -q compile` at the repo root. A compile error naming `COLUMN_SPAN_DENSITY` means +Step 1 was skipped. + +- [ ] **Step 8: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java +git commit -m "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." +``` + +--- + +### Task 5: Declare the new column in the model + +**Files:** +- Modify: `ApiCommonModel/Model/lib/wdk/model/questions/queries/geneQueries.xml` — the + `wsColumn` list AND BOTH `postCacheUpdateSql` blocks + +**Correction to this task as originally written.** It said "add the wsColumn", which is +not sufficient. Both `postCacheUpdateSql` blocks (the `excludeProjects="UniDB"` one and +the `includeProjects="UniDB"` one) enumerate the result columns explicitly in their +`INSERT` column list and their `SELECT` list. Those blocks backfill the sibling +transcripts of a matching gene. Omit the new column there and a gene's matched transcript +carries a span density while its siblings carry NULL, in the same result table — five +edits, not one. Landed as `c3193504`. + +Also note the dependency is tighter than "the perl and Java must ship together": the +`columns` map is built from `request.getOrderedColumns()`, i.e. the model's `wsColumn` +list, NOT from the plugin's `getColumns()`. `PluginExecutor.validateColumns` only checks +that `getColumns()` is a subset. So the Java commit alone throws +`PluginUserException: The required column is missing: span_snp_density` before +`makeResultRow` is reached. The atomic set is three commits across two repos. + +- [ ] **Step 1: Add the wsColumn** + +In the `GenesByNgsSnps` processQuery — the one at line ~2944, NOT the chip query at 2795 +— add a line after ``: + +```xml + +``` + +Verify you edited the right one: + +```bash +cd ~/workspaces/plasmodb/ApiCommonModel/Model/lib/wdk +grep -n -B40 'wsColumn name="span_snp_density"' model/questions/queries/geneQueries.xml | grep 'processQuery name=' +``` + +Expected: `&1 | tail -5` +Expected: `OK - Reloaded application at context path [/plasmo.jbrestel]`. + +- [ ] **Step 3: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonModel +git add Model/lib/wdk/model/questions/queries/geneQueries.xml +git commit -m "Declare span_snp_density on the GenesByNgsSnps process query" +``` + +--- + +### Task 6: Honest labels on the result columns + +**Files:** +- Modify: `ApiCommonModel/Model/lib/wdk/model/questions/geneQuestions.xml:1475-1476,1512-1563` + +The chip-SNP question above this one is commented out (`geneQuestions.xml:1373`), so +these attribute names are owned by `GenesByNgsSnps` alone — no cross-question impact. + +- [ ] **Step 1: Add the new column and relabel the existing ones** + +Inside the `` block of `GenesByNgsSnps` (lines 1512-1563), apply +exactly these display-name and help changes, and add one new `columnAttribute`: + +```xml + + + Display the histogram of the values of this attribute + float + + + + + Display the histogram of the values of this attribute + float + + + + + Display the histogram of the values of this attribute + float + + +``` + +Note the `` values change from `int` to `float` on all three — +they were `int` on values that have always been fractional. + +Then change these four display names in place, leaving their `reporter` blocks as they +are: + +| line ~ | from | to | add help | +|---|---|---|---| +| `ngs_total_snps` | `Total SNPs` | `SNVs in gene span` | `Variant positions anywhere between the gene's start and end.` | +| `ngs_num_synonymous` | `Synonymous SNPs` | `Synonymous SNVs` | *(none)* | +| `ngs_num_non_synonymous` | `Nonsynonymous SNPs` | `Missense SNVs` | `Amino-acid-changing, excluding those that introduce a stop.` | +| `num_nonsense` | `Nonsense SNPs` | `Stop-gained SNVs` | replace the existing typo help `SNPs where one or more variants encodes a stop coding` with `Variants where one or more alleles encodes a premature stop codon.` | +| `num_noncoding` | `Non-coding SNPs` | `Unclassified SNVs` | `No protein product could be assigned: positions outside coding sequence, and positions where the reference product was unavailable.` | + +- [ ] **Step 2: Add the new column to the summary attribute list** + +Replace lines 1475-1476: + +```xml + +``` + +with: + +```xml + +``` + +- [ ] **Step 3: Build the model** + +Run: `cd ~/workspaces/agentic-veupath-dev && bin/veup-build.sh plasmodb wb model 2>&1 | tail -5` +Expected: `OK - Reloaded application`. + +- [ ] **Step 4: Verify the labels through the service** + +From an already-authenticated tab on `https://jbrestel.plasmodb.org` (a raw curl +307-redirects to autologin), run in the browser console: + +```js +const s = await (await fetch('/plasmo.jbrestel/service/record-types/transcript/searches/GenesByNgsSnps?expandParams=true')).json(); +(s.searchData?.dynamicAttributes ?? s.dynamicAttributes ?? []).map(a => a.name + ' | ' + a.displayName) +``` + +Expected to include `cds_snp_density | SNVs per kb (CDS)`, `span_snp_density | SNVs per kb (gene span)`, and `ngs_dn_ds_ratio | dN/dS (site-normalized)`. + +- [ ] **Step 5: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonModel +git add Model/lib/wdk/model/questions/geneQuestions.xml +git commit -m "Label the GenesByNgsSnps result columns for what they now compute + +Adds SNVs per kb (gene span) alongside the now-actually-CDS density, renames +the ratio to say it is site-normalized, and corrects three labels that never +matched the classifier: nonsynonymous excludes stop-gained, and non-coding is +really unclassified - class 0 means no product byte was available, which +covers unclassifiable positions as well as genuinely non-coding ones." +``` + +--- + +### Task 7: Honest labels on the params + +**Files:** +- Modify: `ApiCommonModel/Model/lib/wdk/model/questions/params/geneParams.xml:2328-2341,2350-2362,2379-2391,2506+` + +- [ ] **Step 1: Relabel the density pair (lines 2379-2391)** + +Replace both `stringParam` blocks: + +```xml + + Coding variants per kilobase of coding sequence. Genes with no coding sequence have no value here and are returned only while this filter is left alone. + + + + + Coding variants per kilobase of coding sequence. NOTE: Leaving this parameter value empty means you don't care what the upper bound is. + + +``` + +- [ ] **Step 2: Relabel the ratio pair (lines 2328-2341)** + +```xml + + Upper bound on the site-normalized dN/dS ratio. NOTE: Leaving this parameter value empty means you don't care what the upper bound is. + + + + + Nonsynonymous and synonymous counts each divided by the number of sites of that class. Below 1 suggests purifying selection. Stop-gained variants are counted as nonsense rather than nonsynonymous and so do not enter the numerator. Genes with no synonymous sites have no value and are returned only while this filter is left alone. + + +``` + +- [ ] **Step 3: Relabel the occurrence pair (lines 2350-2362)** + +Only the word SNP changes; the count wording was already correct. + +```xml + + Upper Bound on the number of SNVs of the selected class. NOTE: Leaving this parameter value empty means you don't care what the upper bound is. + + + + + Lower Bound on the number of SNVs of the selected class + + +``` + +- [ ] **Step 4: Relabel the class enum (line ~2506)** + +Change the prompt and three enum TERMS. The `` values are a contract with +`FindGenesWithSnpCharsPlugin.legalParams` and the filter's branches — do NOT touch them. + +```xml + + + + + Choose the class of SNV you want to query on ... choose minumum and maximum numbers below + + + + All SNVs + all + + + Coding + coding + + + Unclassified + noncoding + + + Missense + nonsynonymous + + + Stop-gained + nonsense + + + Synonymous + synonymous + + + +``` + +- [ ] **Step 5: Build and verify the prompts** + +Run: `cd ~/workspaces/agentic-veupath-dev && bin/veup-build.sh plasmodb wb model 2>&1 | tail -5` +Expected: `OK - Reloaded application`. + +Then from an authenticated browser tab: + +```js +const s = await (await fetch('/plasmo.jbrestel/service/record-types/transcript/searches/GenesByNgsSnps?expandParams=true')).json(); +(s.searchData?.parameters ?? s.parameters).map(p => p.name + ' | ' + p.displayName) +``` + +Expected to include `snp_density_lower | SNVs per kb (CDS) >= ` and +`dn_ds_ratio_lower | dN/dS (site-normalized) >= `. + +- [ ] **Step 6: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonModel +git add Model/lib/wdk/model/questions/params/geneParams.xml +git commit -m "Label the GenesByNgsSnps params for what they now filter + +Density is CDS density, the ratio is site-normalized, and the class enum terms +follow the corrected column names. Internal enum values are unchanged - they +are a contract with legalParams and with the filter's branches." +``` + +--- + +### Task 8: Update the question description + +**Files:** +- Modify: `ApiCommonModel/Model/lib/wdk/model/questions/geneQuestions.xml` (the two `` blocks of `GenesByNgsSnps`) + +- [ ] **Step 1: Replace the codon-bias promise in the PlasmoDB description** + +Find this bullet: + +```html +
  • Due to the extreme codon bias in the P. falciparum genome, the ratio of non-synonymous to synonymous SNPs within each gene is much higher than expected. This should be considered when creating queries. We are intending to calculate more reliable normalized Dn/Ds or Ka/Ks ratios in subsequent releases of PlasmoDB.
    • +``` + +Replace with: + +```html +
  • The dN/dS ratio reported here IS normalized by synonymous and nonsynonymous site counts, derived from the genetic code over the gene's representative transcript, so it does not carry the codon-bias inflation that a raw count ratio does. In P. falciparum the pooled synonymous-site fraction is 17.49% rather than the textbook ~25%, which is a 1.43x correction on every gene.
    • +
  • It remains a count ratio: a variant seen in one sample weights the same as one at 50% frequency. For a frequency-weighted piN/piS across every sample loaded, use the "SNV Characteristics" search.
    • +``` + +- [ ] **Step 2: Add the scope bullet to BOTH descriptions** + +There are two `` blocks, one `includeProjects="PlasmoDB"` and one for the +other projects. Add this bullet as the first `
  • ` of each: + +```html +
  • These statistics are computed over the samples you select, so they change with your sample set and are not comparable to the precomputed values on the gene record page or in the "SNV Characteristics" search.
    • +``` + +- [ ] **Step 3: Build** + +Run: `cd ~/workspaces/agentic-veupath-dev && bin/veup-build.sh plasmodb wb model 2>&1 | tail -5` +Expected: `OK - Reloaded application`. + +- [ ] **Step 4: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonModel +git add Model/lib/wdk/model/questions/geneQuestions.xml +git commit -m "Point the GenesByNgsSnps description at the normalized ratio + +The PlasmoDB description promised normalized Dn/Ds in a future release; it is +here. Also states in both descriptions that these statistics are sample-set +scoped and so not comparable to the precomputed ones." +``` + +--- + +### Task 9: Whole-cohort cross-check against the tuning table + +This is the acceptance test. It cannot be an equality assertion — the two pipelines use +different classifiers — so it is a correlation check with the divergence explained. + +**Files:** +- Create: `ApiCommonWebService/docs/superpowers/specs/2026-08-07-hsss-gene-stats-validation.md` + +- [ ] **Step 1: Get the tuning table's site-normalized ratio** + +Run: + +```bash +psql -h localhost -p 5432 -d unidb_shu_a -P pager=off -c " +select gene_source_id, + round(((n_missense::numeric / nullif(nonsyn_sites,0)) + / nullif(n_synonymous::numeric / nullif(syn_sites,0), 0))::numeric, 4) as dnds_gvs +from apidbtuning.genevariationsummary +where org_abbrev = 'pfal3D7' and n_synonymous > 0 +order by gene_source_id +limit 20;" +``` + +Expected: 20 rows with `dnds_gvs` values, most below 2. + +- [ ] **Step 2: Run the search with every sample selected** + +In an authenticated browser tab, run the search through the service with the organism set +to *Plasmodium falciparum 3D7*, every sample in `variation_sample_meta`, and all filters +at their permissive defaults (`occurrences_lower=0`, `occurrences_upper=-1`, +`dn_ds_ratio_lower=0`, `dn_ds_ratio_upper=-1`, `snp_density_lower=0`, +`snp_density_upper=-1`, `snp_class=all`). Export `primary_key` and `ngs_dn_ds_ratio`. + +This run takes minutes — it is HSSS over every sample. Expect that. + +- [ ] **Step 3: Compare and record** + +Save the search export as `/tmp/hsss_dnds.tsv` with two columns (gene id, ratio) and no +header, then compute the correlation in the database rather than by hand: + +```bash +psql -h localhost -p 5432 -d unidb_shu_a -P pager=off <<'SQL' +CREATE TEMP TABLE hsss_dnds (gene_source_id text, dnds_hsss numeric); +\copy hsss_dnds FROM '/tmp/hsss_dnds.tsv' WITH (FORMAT csv, DELIMITER E'\t') +SELECT count(*) AS genes_compared, + round(corr(rank_h, rank_g)::numeric, 3) AS spearman_rho, + round(percentile_cont(0.5) WITHIN GROUP (ORDER BY dnds_hsss)::numeric, 4) AS median_hsss, + round(percentile_cont(0.5) WITHIN GROUP (ORDER BY dnds_gvs)::numeric, 4) AS median_gvs +FROM ( + SELECT h.dnds_hsss, + (g.n_missense::numeric / nullif(g.nonsyn_sites,0)) + / nullif(g.n_synonymous::numeric / nullif(g.syn_sites,0), 0) AS dnds_gvs, + rank() OVER (ORDER BY h.dnds_hsss) AS rank_h, + rank() OVER (ORDER BY (g.n_missense::numeric / nullif(g.nonsyn_sites,0)) + / nullif(g.n_synonymous::numeric / nullif(g.syn_sites,0), 0)) AS rank_g + FROM hsss_dnds h + JOIN apidbtuning.genevariationsummary g + ON g.gene_source_id = h.gene_source_id AND g.org_abbrev = 'pfal3D7' + WHERE h.dnds_hsss IS NOT NULL + AND g.n_synonymous > 0 AND g.syn_sites > 0 AND g.nonsyn_sites > 0 +) t; +SQL +``` + +Then write +`2026-08-07-hsss-gene-stats-validation.md` recording: +- the observed correlation and the number of genes compared +- the median of each, side by side +- the three reasons they differ, verbatim from the spec section 5: different classifiers + (reference product bytes vs SnpEff severity), HSSS not seeing sample-vs-reference + differences unless the reference strain is selected, and HSSS excluding stop-gained + from its nonsynonymous count + +**Acceptance:** Spearman rho above 0.7. Below that, something is wrong beyond classifier +differences — stop and investigate rather than recording a bad number. + +- [ ] **Step 4: Coverage regression** + +Confirm non-coding genes survive. Pick a pfal gene with `cds_length IS NULL`: + +```bash +psql -h localhost -p 5432 -d unidb_shu_a -Atc " +select gene_source_id from apidbtuning.genevariationsummary +where org_abbrev='pfal3D7' and cds_length is null limit 5;" +``` + +Then verify in the Step 2 result that such a gene appears with an EMPTY +`cds_snp_density` and a POPULATED `span_snp_density`, and that narrowing +`snp_density_lower` to any value above 0 removes it. + +- [ ] **Step 5: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebService +git add docs/superpowers/specs/2026-08-07-hsss-gene-stats-validation.md +git commit -m "Record the whole-cohort validation of the corrected geneChars statistics" +``` + +--- + +### Task 10: Release note + +**Files:** +- Modify: `ApiCommonWebsite/Model/lib/xml/PlasmoDB/news.xml` + +- [ ] **Step 1: Add the note** + +Saved strategies filtering on either statistic will return different results after +deploy, for every gene rather than only edge cases. + +The file is a `` of `` elements. Insert this as the FIRST `` +after the opening `` and the commented-out template, and set the date to the +release date rather than today: + +```xml + + + + + DD Mmm YYYY 00:00 + + + + +
      +
    • The SNV Characteristics Within a Group of Samples search (formerly + "SNP Characteristics") now reports SNVs per kb of coding sequence. It + previously reported all variants per kb of genomic span while labelling the + column as CDS. The previous value is still available in a new + "SNVs per kb (gene span)" column.
      • +
    • Its dN/dS ratio is now normalized by synonymous and nonsynonymous site + counts derived from the genetic code, so it no longer carries the codon-bias + inflation of a raw count ratio.
      • +
    • Both values change for every gene, and saved strategies that + filter on either will return different results than before.
      • +
    + + ]]> +     +     +``` + +- [ ] **Step 2: Commit** + +```bash +cd ~/workspaces/plasmodb/ApiCommonWebsite +git add Model/lib/xml/PlasmoDB/news.xml +git commit -m "News: corrected density and dN/dS in the sample-set SNV search" +``` + +--- + +## Deployment checklist + +- [ ] `ApiCommonWebService` and `ApiCommonModel` merged and released together. A model-only + deploy leaves 8-field Java parsing 9-field perl output and every search fails. +- [ ] The WSF plugin jar is rebuilt — `wb model` does not do this. +- [ ] `bin/veup-git-sync.sh plasmodb` after any local branch switch, per the repo CLAUDE.md. diff --git a/docs/superpowers/specs/2026-08-05-hsss-variation-plumbing-design.md b/docs/superpowers/specs/2026-08-05-hsss-variation-plumbing-design.md new file mode 100644 index 00000000..347f95f6 --- /dev/null +++ b/docs/superpowers/specs/2026-08-05-hsss-variation-plumbing-design.md @@ -0,0 +1,358 @@ +# HSSS variation plumbing — design + +**Date:** 2026-08-05 +**Status:** implemented 2026-08-05 — ID fix verified; the search directory (§3.1) and filter param name (§3.5) await the `VariationsByIsolateGroup` search +**Scope:** Make the HighSpeedSnpSearch (HSSS) plugins emit **variation** record IDs and +read the **variation** HSSS directory layout. No WDK model XML, no new searches. +**Implementation targets:** `ApiCommonWebService` (Java, Perl, test fixtures) **and** +`ApiCommonWebsite` (one Conifer variable — see §3.3). Both on branch +`dnaseq-merge-experiments`. +**Followed by:** a separate spec for `VariationsByIsolateGroup`, the first search to +consume this. See §7. + +## 1. Purpose + +The `variation` record has one search, `VariationBySourceId` (a plain `sqlQuery`). Every +remaining search ported from the deprecated `snp` record — `ByIsolateGroup`, +`ByLocation`, `ByGeneIds`, `ByTwoIsolateGroups` — is a `processQuery` against an HSSS +plugin, and **none of them can work until the plugins speak the variation record's +language.** + +Two mismatches block them. Both are in `ApiCommonWebService`, and both are invisible +until a search actually runs: + +| | HSSS emits / expects today | variation record needs | +|---|---|---| +| result `source_id` | `NGS_SNP.Pf3D7_01_v3.29514` | `Variant_Pf3D7_01_v3_29514` | +| data directory | `//build-N/Pfalciparum3D7/**highSpeedSnpSearch**` | `.../Pfalciparum3D7/**dnaseq**` | + +The ID mismatch is the dangerous one: the plugin would return rows whose `source_id` +matches no variation record, so the search yields **zero results and no error**. The +directory mismatch fails loudly (`"Organism dir does not exist"`). + +This spec fixes both, and is verifiable on its own — see §5. It is deliberately separated +from the search that consumes it, because this is a Java/Perl integration whose acceptance +test is "the IDs come out right", while the search is declarative XML that cannot be +verified until this is deployed. + +## 2. Why editing in place is safe + +The HSSS plugins are **entirely unused by the live model**. Verified against the +assembled PlasmoDB model (`wdkXml -model PlasmoDB`, 17,638 lines): + +- zero references to `FindPolymorphismsPlugin`, `FindPolymorphismsWithSeqFilterPlugin`, + `FindSnpsByGeneIdsPlugin`, `FindMajorAllelesPlugin`, `FindGenesWithSnpCharsPlugin`, + `FindChipPolymorphismsPlugin` +- no `SnpQuestions` and no `SnpChipQuestions` questionSet + +The reason is that the snp and snp-chip `` blocks are commented out in the +**shared** `ApiCommonModel/Model/lib/wdk/apiCommonModel.xml`, so this holds for **every** +project, not just PlasmoDB. There is therefore no live consumer to regress, and no need +for a plugin subclass or a parallel script: the existing classes are modified directly and +become variation-only going forward. + +> Note for anyone re-checking this: `grep` for the `` line is **not** a valid test, +> because it matches inside XML comment blocks just as happily as outside them. Ask the +> assembled model instead. + +Orphaned ontology rows for `GeneQuestions.GenesByNgsSnps` and `GeneQuestions.GenesBySnps` +still exist in `individuals.txt`, naming questions the model no longer defines. They are +harmless (the category ontology is not validated against the model) and out of scope here. + +## 3. The changes + +Four production edits across two repos, plus the fixture hygiene in §4: + +| # | repo / file | change | +|---|---|---| +| 3.1 | `ApiCommonWebService` — `WSFPlugin/.../highspeedsnpsearch/HighSpeedSnpSearchAbstractPlugin.java:196` | `getSearchDir()` → `/dnaseq` | +| 3.2 | `ApiCommonWebService` — `HighSpeedSnpSearch/bin/hsssReconstructSnpId:42-43` | ID separator `.` → `_`, both joins | +| 3.3 | `ApiCommonWebsite` — `Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml:102` | `highspeedsnpsearchconfig_idPrefix` → `Variant_` | +| 3.5 | `ApiCommonWebService` — `WSFPlugin/.../highspeedsnpsearch/FindPolymorphismsPlugin.java:41-43` | `getStrainFilterParamName()` → `variation_sample_meta` | + +### 3.1 Search directory — `HighSpeedSnpSearchAbstractPlugin.java:196` + +```java +protected String getSearchDir() { + return "/dnaseq"; // was "/highSpeedSnpSearch" +} +``` + +`findOrganismDir` (same file, :201-213) composes +`webSvcPath.replaceAll("PROJECT_GOES_HERE", projectId) + "/" + organismNameForFiles + searchDir` +and throws if the result does not exist. `organismNameForFiles` comes from +`apidb.organism.name_for_filenames`, resolved by matching the organism param's internal +value against `sres.TaxonName.name` (:317). + +Verified layout of the real variation HSSS files on cedar +(`/home/jbrestel/webserviceTest/Pfalciparum3D7/dnaseq/`): + +``` +readFreq20/ readFreq40/ readFreq60/ readFreq80/ +``` + +and inside each: 538 numbered strain directories plus `contigIdToSourceId.dat`, +`strainIdToName.dat`, `referenceGenome.dat`. The `readFreq*` level is appended separately +by `FindPolymorphismsAbstractPlugin:106` (`new File(organismDir, "readFreq" + pct)`), and +the four directories match `snpParams.ReadFrequencyPercent`'s enum values (20/40/60/80) +exactly. So `/dnaseq` is the whole of the missing piece. + +**Do not touch the chip overrides.** `FindChipPolymorphismsPlugin:42` and +`FindChipSnpMajorAllelesPlugin:72` override `getSearchDir()` to +`/highSpeedChipSnpSearch`. They are equally dead, but leaving them alone keeps this diff +to what it needs to be. + +### 3.2 ID separator — `HighSpeedSnpSearch/bin/hsssReconstructSnpId:42-43` + +The Perl builds the ID as `$prefix . "$contigSourceId.$location" . $suffix`. Prefix and +suffix are configurable; **the `.` separator is hardcoded**, which is why config alone +cannot produce a variation ID. Change it to `_` in **both** the STDERR and STDOUT `join`s +(the two lines are otherwise identical; the STDOUT one carries the `$seqFilter` guard): + +```perl +print STDERR join("\t", $prefix."${contigSourceId}_${location}".$suffix, @fields) . "\n"; +print STDOUT join("\t", $prefix."${contigSourceId}_${location}".$suffix, @fields) . "\n" + unless ($seqFilter && ($contigSourceId ne $seqFilter || $location < $minLoc || $location > $maxLoc)); +``` + +`contigIdToSourceId.dat` maps `1 → Pf3D7_01_v3` (16 contigs for Pf), so with the prefix +from §3.3 the emitted ID becomes `Variant_Pf3D7_01_v3_29514` — exactly the +`VariationRecordClass` `source_id` format. + +Braces around the interpolated names are required: `"$prefix$contigSourceId_$location"` +would parse `$contigSourceId_` as a variable name. + +### 3.3 ID prefix — a Conifer variable, in `ApiCommonWebsite` + +This one is **not** in this repo, and there are three candidate files. Only one is right: + +| file | verdict | +|---|---| +| `gus_home/config//highSpeedSnpSearch-config.xml` | **No** — generated, regenerated on every build, would silently revert | +| `ApiCommonWebService/WSFPlugin/lib/conifer/roles/conifer/templates/ApiCommonWebService/highSpeedSnpSearch-config.xml.j2` | **No** — the template does not hold the value; it renders `{{ highspeedsnpsearchconfig_idPrefix }}` | +| **`ApiCommonWebsite/Model/lib/conifer/roles/conifer/vars/ApiCommon/default.yml:102`** | **Yes** — this is where the value lives | + +```yaml +highspeedsnpsearchconfig_idPrefix: Variant_ # was NGS_SNP. +``` + +**So this change spans two repos.** `ApiCommonWebService` (§3.1, §3.2, §4) and +`ApiCommonWebsite` (this section). Both are on branch `dnaseq-merge-experiments`; the +`ApiCommonWebsite` edit is a one-line var change and needs no build of its own, but the +site must be re-conifered for it to reach `gus_home/config`. + +It is an **ApiCommon cohort default**, so it applies to every ApiCommon project — +acceptable only because §2 establishes there are no other consumers. If a second consumer +with a different ID convention ever appears, this must become per-plugin; do not solve that +now. + +`HighSpeedSnpSearchAbstractPlugin.getIdPrefix()` (:157-165) reads the property and falls +back to the literal string `"NULL"`, which the Perl treats as empty (`:40-41`). + +For **testing only**, the value can be overridden per instance in +`etc/conifer_site_vars.yml` without touching the shared default — the same mechanism §5.1 +uses for `webServiceMirror`. Prefer that while iterating; change the cohort default when +the behaviour is settled. + +`HighSpeedSnpSearchAbstractPlugin.getIdPrefix()` (:157-165) reads this property and falls +back to the literal string `"NULL"`, which the Perl treats as empty (`:40-41`). The +property is **site-wide across all HSSS plugins** — acceptable only because §2 establishes +there are no other consumers. If a second consumer with a different ID convention ever +appears, this property must become per-plugin; do not solve that now. + +### 3.4 Naming: deliberately unchanged + +`hsssReconstructSnpId`, the `org.apidb.apicomplexa.wsfplugin.highspeedsnpsearch` package, +and the `Snp`-flavoured class names all stay. Honest naming argues for a rename, but it +touches ten Java files, the installed `gus_home/bin` script names, and every +`processName` in model XML — and bundling a rename with a behavioural fix makes the diff +unreviewable. Worth a follow-up issue; not this change. + +### 3.5 Strain filter param name — `FindPolymorphismsPlugin.java:41-43` + +```java +protected String getStrainFilterParamName() { + return "variation_sample_meta"; // was "ngsSnp_strain_meta" +} +``` + +`FindPolymorphismsAbstractPlugin:100` reads the samples selection from the request under +whatever name this returns, and declares it required (`:41`). So **this string is a +contract with the model XML**: the consuming search's `filterParam` must be named exactly +this, or the plugin rejects the request as missing a required parameter. + +Left as `ngsSnp_strain_meta`, brand-new variation model XML would be forced to carry snp +vocabulary forever. Since the plugin is variation-only from here (§2), rename it now while +there is no consumer to break — the follow-on spec then names its filterParam +`variation_sample_meta`. + +This is the one change in §3 with no runtime symptom in isolation; it only matters once a +search exists. It belongs here rather than in the search spec because it is a change to +this repo. + +## 4. Test fixtures — hygiene, not verification + +**Correction to an earlier draft of this spec, which claimed `ApiCommonWebService/Test` +was "a working JUnit harness" whose tests "will fail after §3 unless updated". Both HSSS +test harnesses are already broken, independently of this change.** They rotted when the +snp searches were retired. Evidence: + +| harness | why it cannot run | +|---|---| +| `Test/.../FindPolymorphismsSearchTest.java:32` (JUnit) | references `FindPolymorphismsPlugin.PARAM_STRAIN_LIST`, a constant **defined nowhere** in the main sources — the module does not compile. Consistent with `Test-Installation` being absent from `build.xml`'s default `ApiCommonWebService-Installation` depends list. | +| `HighSpeedSnpSearch/bin/hsssTestSuite:41` (shell) | passes 8 positional args to `hsssGeneratePolymorphismScript`, which consumes 5 standard args (`HsssScriptGenerator.pm:25`) and then expects `polymorphismThreshold, unknownThreshold, strainsListFile, reconstructCmdName, idPrefix, idSuffix` (`HsssPolymorphismScriptGenerator.pm:75`). `reconstructCmdName` and `idPrefix` arrive **undefined**, so the generated command is malformed. | + +> **Therefore no green test run gates this change, and the plan must not pretend one does.** +> Verification is §5, whose first two rungs are real and runnable. + +Reviving either harness is **out of scope** — it is a larger job than this change and +would have to be done against variation data to be worth anything. + +What is still worth doing, because it is ~15 lines and keeps the fixtures honest for +whoever does revive them: + +### 4.1 Rename the fixture directory + +``` +HighSpeedSnpSearch/test/TestDB/Hsapiens123/highSpeedSnpSearch/ → .../dnaseq/ +``` + +(containing `readFreq80/`). The tests point `PARAM_WEBSVCPATH` at +`$PROJECT_HOME/ApiCommonWebService/HighSpeedSnpSearch/test/PROJECT_GOES_HERE`, and the +mock project mapper supplies `TestDB` for the `PROJECT_GOES_HERE` substitution. + +### 4.2 Update the three expected files with baked-in IDs + +In `HighSpeedSnpSearch/test/expected/`: + +| file | rows | change | +|---|---|---| +| `genomicLocationFilter.txt` | 2 | `NGS_SNP.e99.2011` → `Variant_e99_2011`; `NGS_SNP.h103.30021` → `Variant_h103_30021` | +| `polymorphismSearchWithSourceIds.txt` | 8 | `a80.896`, `b86.13441`, `e99.2011`, `f100.23`, `g102.4334`, `h103.30021`, `i104.3002`, `j201.54` — each `NGS_SNP..` → `Variant__` | +| `polymorphismSearchWithSourceIdsAndSeqFilter.txt` | 1 | `NGS_SNP.f100.23` → `Variant_f100_23` | + +All three are tab-delimited, with the ID in column 1 and the remaining columns unchanged. +The transformation is mechanical and total: every `NGS_SNP.` becomes `Variant_`, and the +**single remaining dot** in each ID becomes an underscore. Note the contig names +themselves contain no dots, so a naive global dot-to-underscore replacement is safe here — +but write the change deliberately rather than relying on that. + +The other five expected files (`polymorphismSearch.txt`, `majorAlleles.txt`, +`mergeStrains*.txt`) hold pre-reconstruction output — contig **index** and location as +separate columns — and must **not** be touched. Only files whose first column is a +reconstructed ID change. + +## 5. Verification + +Cheapest rung first. Rung 1 alone proves the entire ID fix and needs no build. + +1. **The Perl, directly.** Against the real Pf mapping file on cedar: + +```bash +printf '1\t29514\t100\t5\t1\n' \ + | $GUS_HOME/bin/hsssReconstructSnpId \ + /home/jbrestel/webserviceTest/Pfalciparum3D7/dnaseq/readFreq20/contigIdToSourceId.dat \ + 1 Variant_ NULL +``` + +Expect exactly `Variant_Pf3D7_01_v3_29514 100 5 syn`. (Field 5 = `1` maps to `syn` per +the script's coding-class translation; `100`/`5` pass through as knowns/non-major +percentages.) Before the §3.2 change the same command yields +`Variant_Pf3D7_01_v3.29514` — the dot is the bug, so a pre-change run is expected to show it. + +2. **Confirm the ID resolves to a real record**, closing the loop the search would: + +```sql +SELECT source_id FROM apidbtuning.VariationAttributes +WHERE source_id = 'Variant_Pf3D7_01_v3_29514'; +``` + +One row. This is the check that the whole spec exists to satisfy. + +3. **The install still builds.** `bld ApiCommonWebService` succeeds and reinstalls both the + Java and the Perl into `gus_home` — the §3.2 edit only takes effect in + `$GUS_HOME/bin/hsssReconstructSnpId` after an install, so rung 1 must be re-run + afterwards to confirm the installed copy carries the fix, not just the source tree. + +4. **The prefix actually reaches the config.** After re-running Conifer for the site, the + generated file shows the new value: + +```bash +grep idPrefix $GUS_HOME/config/PlasmoDB/highSpeedSnpSearch-config.xml +``` + +Expect `Variant_`. This is the only check that §3.3 landed; +nothing else reads that file until a search runs. + +5. **Integration on cedar** — the plugin resolves the real directory. Requires the + `webServiceMirror` override below. Note this cannot be exercised without a search to + invoke the plugin, so in practice it lands with rung 6. + +6. **End-to-end through a real search** — deferred to the `VariationsByIsolateGroup` spec, + which is the first thing able to exercise it. §3.1 and §3.5 are **not independently + verifiable before that point**; the honest state after this change is "the ID fix is + proven, the directory and param-name fixes are correct by inspection." + +**No automated test suite gates this change** — see §4 for why both existing harnesses are +already broken. Do not report a green test run. + +### 5.1 Pointing the plugin at the test files + +Override `webServiceMirror` in the instance's `etc/conifer_site_vars.yml` — declared +per-instance config, regenerated on build, never committed model XML. Do **not** hardcode a +path in model XML. + +The `WebServicesPath` param's internal value is +`@WEBSERVICEMIRROR@/PROJECT_GOES_HERE/build-`, so with +`webServiceMirror: /home/jbrestel/webserviceTest` the plugin looks for +`/home/jbrestel/webserviceTest/PlasmoDB/build-70/Pfalciparum3D7/dnaseq` — two levels +deeper than where the test files actually sit. Bridge it with a symlink: + +```bash +mkdir -p /home/jbrestel/webserviceTest/PlasmoDB/build-70 +ln -s /home/jbrestel/webserviceTest/Pfalciparum3D7 \ + /home/jbrestel/webserviceTest/PlasmoDB/build-70/Pfalciparum3D7 +``` + +Without it the failure is `"Organism dir does not exist"`, which reads like a code bug +rather than a path problem. + +## 6. Facts established for the consuming spec + +Verified here so the search spec does not need to re-derive them: + +- **Sample names match.** All **216** distinct `sample_stable_id`s in + `eda.attributevalue_s3be28bbe14_sample` are a strict subset of the **538** strain names + in `strainIdToName.dat`. So a filterParam whose internal values are EDA sample stable + IDs will only ever name strains HSSS knows. No mapping layer is needed. +- **The filterParam must be named `variation_sample_meta`** — the exact string + `FindPolymorphismsPlugin.getStrainFilterParamName()` returns after §3.5, and a required + parameter (`FindPolymorphismsAbstractPlugin:41`). A mismatch is rejected as a missing + required param. +- **The results-file contract is exactly 4 tab-separated columns** — + `FindPolymorphismsAbstractPlugin:141` throws otherwise. Order: + `sourceId`, `percentOfKnowns`, `percentOfPolymorphisms`, `phenotype`. `wsColumn` + declarations must follow it. +- **The organism param's internal value must remain the taxon name** + (e.g. `Plasmodium falciparum 3D7`), because `getOrganismNameForFiles` looks it up in + `sres.TaxonName.name`. It cannot be repurposed to carry an EDA study abbreviation; the + search spec needs a separate hidden param for that. +- **`name_for_filenames`** for the three organisms with dnaseq isolate data: + `Pfalciparum3D7` (PlasmoDB), `TbruceiTREU927` (TriTrypDB), `AfumigatusAf293` (FungiDB). + +## 7. Out of scope + +- **`VariationsByIsolateGroup`** and the other three HSSS searches — separate specs in + `ApiCommonModel`. This spec is a prerequisite for all of them. +- **Renaming** the snp-flavoured classes, script, and package — §3.4. +- **Making `idPrefix` per-plugin** — §3.3. Only needed if a second consumer appears. +- **Removing the dead snp/chip plugins** (`FindChip*`, `FindMajorAlleles*`, + `FindGenesWithSnpChars*`) and their tests. They are unreferenced, but deleting them is a + separate decision with its own review. +- **Reviving either broken test harness** — §4. The JUnit module does not compile and the + shell suite passes the wrong argument count; fixing them is a larger job than this + change and should be done against variation data. +- **Deleting the orphaned `GenesByNgsSnps`/`GenesBySnps` ontology rows** — §2. +- **Populating the production HSSS directories** under + `/var/www/Common/apiSiteFilesMirror/webServices//build-/`. This spec is + verified against the test copy; production placement is a data-deployment task. diff --git a/docs/superpowers/specs/2026-08-07-hsss-gene-stats-fix-design.md b/docs/superpowers/specs/2026-08-07-hsss-gene-stats-fix-design.md new file mode 100644 index 00000000..77b404c0 --- /dev/null +++ b/docs/superpowers/specs/2026-08-07-hsss-gene-stats-fix-design.md @@ -0,0 +1,315 @@ +# Fixing the two broken statistics in the HSSS gene characteristics search + +Design, 2026-08-07. + +Subject: `GenesByNgsSnps` — the WDK search displayed as **SNV Characteristics Within a +Group of Samples**, backed by `FindGenesWithSnpCharsPlugin` and +`bin/hsssGeneCharacteristicsFilter`. + +## 1. The problem + +Two of the statistics this search reports do not compute what their labels claim, and a +third label describes a category it does not hold. + +### 1.1 "SNPs per Kb (CDS)" is neither CDS nor coding + +`hsssGeneCharacteristicsFilter`, `processGene`: + +```perl +my $density = $snpsCount / (($filterEnd - $filterStart) / 1000); +``` + +`$snpsCount` is every variant position in the gene, coding or not. `$filterStart` and +`$filterEnd` come from `geneLocations.txt`, which the plugin builds from +`webready.GeneAttributes_p.start_min` / `end_max` — the **genomic span**, introns and +UTRs included. So the reported number is total variants per kb of genomic length. The +column label, the param prompts (`SNPs per KB (CDS) >=`), and the param help ("density +of coding snps ... / KB of coding sequence") are each wrong in both the numerator and +the denominator. + +### 1.2 The nonsyn/syn ratio has no site normalization + +```perl +my $dnds = $synCount ? $nonSynCount / $synCount : undef; +``` + +A raw count ratio. It carries the codon-bias distortion that the search's own PlasmoDB +description apologises for in prose: + +> Due to the extreme codon bias in the *P. falciparum* genome, the ratio of +> non-synonymous to synonymous SNPs within each gene is much higher than expected. ... +> We are intending to calculate more reliable normalized Dn/Ds or Ka/Ks ratios in +> subsequent releases of PlasmoDB. + +The magnitude is known and measured. `GeneVariationSummary` derives Nei-Gojobori site +counts inline from the genetic code and found the pooled synonymous-site fraction in +pfal3D7 to be **17.49%, not the textbook ~25%** — a 1.43x correction on every gene. +Without it the median piN/piS is 2.0, implying genome-wide positive selection; with it +the median is 0.512, the expected purifying-selection signature. An unnormalized count +ratio inherits the whole of that error. + +### 1.3 "Non-coding SNPs" means "unclassified" + +`hsssFindPolymorphic.c`: + +```c +if (product == 'X') product = -1; // X is an unknown product. ignore these. +... +int productClass = 0; // noncoding +if (nonSyn) productClass = 2; +else if (refProduct > 0) productClass = 1; // syn +if (nonsense) productClass *= -1; +``` + +Class `0` means no protein product byte was available. That covers positions genuinely +outside coding sequence **and** positions where the product could not be called. The +filter then derives `$nonCodingCount = $snpsCount - $codingCount`, so the column labelled +"Non-coding SNPs" is really "not classified as coding". + +### 1.4 Related, already fixed separately + +The `snp_class` enum offered `Non-Coding`, but `FindGenesWithSnpCharsPlugin.legalParams` +omitted `"noncoding"` (and listed `"coding"` twice), so selecting it threw +`PluginUserException` before the script ran, even though +`hsssGeneCharacteristicsFilter` branches on it. Fixed on this branch in a prior commit; +recorded here because it is the same class of defect — a display layer promising +something the compute layer does not deliver. + +## 2. What is deliberately NOT changing + +**Stop-gained stays out of the nonsynonymous count.** `$nonsenseCount++ if $productClass +< 0` catches classes `-1` and `-2`, and those loci are therefore absent from +`$nonSynCount`. A stop-gained change is an amino-acid-changing change, so `dN` understates +selection on genes carrying premature stops. Leaving it alone keeps the column semantics +exactly as they are today and keeps results comparable to historical ones. Revisit only +with a deliberate decision; do not "tidy" it. + +**The statistic remains a count ratio, not piN/piS.** HSSS carries a per-locus minor +allele frequency (`nonMajorAllelesPct`, the fourth field of its per-SNP stream), so a +frequency-weighted piN/piS is reachable. It is out of scope here: Nei's `n/(n-1)` +correction needs an allele count, and HSSS carries `knownsPercent` x strain count without +ploidy. The consequence to communicate is that this ratio weights a singleton the same as +a 50%-frequency variant, and so will NOT equal the piN/piS reported by the +`GenesByVariantCharacteristics` search or the gene record page. + +**The search is not being deprecated.** Per-sample-set analysis is a real workflow. The +whole-cohort `GenesByVariantCharacteristics` search covers a different question and does +not replace this one. + +## 3. Design + +### 3.1 The split that makes this tractable + +Numerators are sample-set-dependent; denominators are not. HSSS keeps computing counts +over the selected samples. The normalizers — CDS length and Nei-Gojobori site counts — +are properties of the gene and the genetic code, already derived once in +`apidbtuning.GeneVariationSummary`. They meet in `geneLocations.txt`. + +Reusing those columns rather than recomputing them is the point: both searches then rest +on one definition of a synonymous site, and there is no second place for it to drift. + +### 3.2 Data flow + +``` +FindGenesWithSnpCharsPlugin.initForBashScript + SELECT g.sequence_id, g.start_min, g.end_max, g.source_id, + gvs.cds_length, gvs.syn_sites, gvs.nonsyn_sites + FROM webready.GeneAttributes_p g + LEFT JOIN apidbtuning.GeneVariationSummary gvs + ON gvs.gene_source_id = g.source_id + AND gvs.project_id = g.project_id + WHERE g.source_id IS NOT NULL + AND g.organism = '' + | + v +geneLocations.txt + seq \t start \t end \t geneId \t cdsLen \t synSites \t nonsynSites + | sort keys are still -k 1,1 -k 2,2n; the new columns are appended, so + | apiSortNoLocale is unaffected + v +hsssGeneCharacteristicsFilter +``` + +`LEFT JOIN`, never inner. A gene with no `GeneVariationSummary` row still gets its +locations line and still reports counts; only its normalized statistics are blank. This +matters: the tuning table holds one row per gene **that has cohort variants**, which in +unidb_shu_a is 5,579 of 5,720 annotated pfal3D7 genes. + +### 3.3 What the filter computes + +```perl +my ($filterContigId, $filterStart, $filterEnd, $filterGeneId, + $filterCdsLen, $filterSynSites, $filterNonsynSites) = split(/\t/, $geneLocationLine); +... +my $cdsDensity = $filterCdsLen ? 1000 * $codingCount / $filterCdsLen : undef; +my $spanDensity = 1000 * $snpsCount / ($filterEnd - $filterStart); +my $dn = $filterNonsynSites ? $nonSynCount / $filterNonsynSites : undef; +my $ds = $filterSynSites ? $synCount / $filterSynSites : undef; +my $dnds = (defined($dn) && $ds) ? $dn / $ds : undef; +``` + +`defined($dn)`, not `$dn`. A gene with zero nonsynonymous variants has `dN = 0` and a +genuine ratio of 0 — a strong purifying-selection signal, and exactly the kind of gene +someone searching a low dN/dS range wants. Testing `$dn` for truth would silently +convert that into "no value" and drop the gene. `$ds` is tested for truth on purpose: +zero there is a division by zero, not a result. + +`$cdsDensity` uses `$codingCount` (`productClass != 0`), so numerator and denominator are +both about coding sequence. `$spanDensity` keeps today's definition under a name that +admits what it is. + +Undefined denominators are expected, not exceptional: `syn_sites` and `cds_length` are +NULL for non-coding genes — 283 of 5,579 pfal3D7 rows, 2,837 of 11,689 tbruTREU927 rows, +196 of 10,029 afumAf293 rows. + +### 3.4 Filter semantics for undefined values + +The existing `-1`-means-no-upper-bound sentinel already gates each filter: + +```perl +if ($densityMin || $densityMax != -1) { ... } +if ($dndsMin || $dndsMax != -1) { ... } +``` + +An untouched filter therefore applies no constraint, and a gene with an undefined +statistic is dropped only when the user actually filters on it. That is the desired +behaviour and it needs no new machinery — worth noting because the sibling +`GenesByVariantCharacteristics` search needed a `filterParam` to get the same property, +`numberRangeParam` having no way to express "untouched". + +When the user HAS narrowed a range and the gene's value is undefined, the gene is +excluded. An unknown value cannot be shown to be in range. + +The existing zero-denominator branch keeps its behaviour: + +```perl +if ($synCount == 0 && $nonSynCount != 0) { return 0 unless $dndsMax == -1; } +``` + +### 3.5 Wire format + +Output goes from 8 fields to 9, gaining span density: + +``` +geneId, cdsDensity, spanDensity, dnds, synCount, nonSynCount, nonCodingCount, nonsenseCount, snpsCount +``` + +`FindGenesWithSnpCharsPlugin.makeResultRow` hard-asserts `parts.length != 8` and +allocates `new String[11]`; both move. `getColumns` and the `wsColumn` list in +`geneQueries.xml` gain the new column. + +**Internal column names do not change.** `cds_snp_density` becomes actually-CDS and +`ngs_dn_ds_ratio` becomes site-normalized, so both names become correct without a rename +that would break `attributesList summary=` and saved strategies. The new column is +`span_snp_density`. + +### 3.6 Compatibility + +A saved strategy filtering on density or on the ratio will return a **different result +set** after deploy. That is the fix working, not a regression, but it is user-visible and +needs a release note. The values change for every gene, not only edge cases: density +changes by the intron fraction of each gene, and the ratio by the 1.43x site-fraction +correction in pfal. + +## 4. Honest labels + +Part of the fix, not an alternative to it. The calculations become correct; the labels +must then describe what is correct, including where the search still differs from its +whole-cohort sibling. + +### 4.1 Result columns (`geneQuestions.xml`) + +| internal | current label | new label | help | +|---|---|---|---| +| `cds_snp_density` | SNPs per Kb (CDS) | **SNVs per kb (CDS)** | Coding variants per kilobase of coding sequence, using the representative transcript's CDS length. Blank for genes with no coding sequence. | +| `span_snp_density` | *(new)* | **SNVs per kb (gene span)** | All variants in the gene divided by its genomic length, introns and UTRs included. | +| `ngs_dn_ds_ratio` | Nonsyn/syn SNP ratio | **dN/dS (site-normalized)** | Nonsynonymous and synonymous counts each divided by the number of sites of that class, from Nei-Gojobori counts over the representative transcript. Below 1 suggests purifying selection. Stop-gained variants are counted as nonsense, not nonsynonymous, so they do not enter the numerator. This is a count ratio: it weights a rare variant the same as a common one, so it will not equal the piN/piS in the SNV Characteristics search. | +| `ngs_num_non_synonymous` | Nonsynonymous SNPs | **Missense SNVs** | Amino-acid-changing, excluding those that introduce a stop. | +| `num_nonsense` | Nonsense SNPs | **Stop-gained SNVs** | *(unchanged meaning)* | +| `num_noncoding` | Non-coding SNPs | **Unclassified SNVs** | No protein product could be assigned: positions outside coding sequence, and positions where the reference product was unavailable. | +| `ngs_total_snps` | Total SNPs | **SNVs in gene span** | Variant positions anywhere between the gene's start and end. | +| `ngs_num_synonymous` | Synonymous SNPs | **Synonymous SNVs** | *(unchanged meaning)* | + +### 4.2 Params (`geneParams.xml`) + +- `snp_class` prompt `SNP Class` -> **SNV Class**. Enum terms follow the columns: + `Non-Coding` -> **Unclassified**, `Non-Synonymous` -> **Missense**, + `Nonsense` -> **Stop-gained**. Internal values are UNCHANGED — they are a contract with + `legalParams` and with the filter's branches. +- `occurrences_lower` / `_upper`: `Number of SNPs of above class` -> **Number of SNVs of + the selected class**. The count wording is already correct; it is + `HsssGeneCharsFilterScriptGenerator.pm`'s usage text that wrongly says "percent", and + that is corrected too. +- `dn_ds_ratio_lower` / `_upper`: `Non-synonymous / synonymous SNP ratio` -> **dN/dS + (site-normalized)**, help gaining the stop-gained and count-ratio caveats. +- `snp_density_lower` / `_upper`: `SNPs per KB (CDS)` -> **SNVs per kb (CDS)**; help + loses "density of coding snps ... / KB of coding sequence" and states the CDS-length + denominator and the blank-for-non-coding behaviour. This pair filters CDS density only; + span density is reported but not filterable, to avoid a fifteenth and sixteenth param + on a form that already carries fourteen. +- Every "leaving this parameter value empty means you don't care about the upper bound" + stays. That sentinel is real behaviour and is what makes undefined values safe. + +### 4.3 Question description + +The PlasmoDB bullet promising normalized ratios "in subsequent releases" becomes a +statement that they are here, plus a pointer to `GenesByVariantCharacteristics` for the +frequency-weighted piN/piS over the whole cohort. One bullet added: these statistics are +computed over the samples you select and are not comparable to the precomputed ones. + +## 5. Verification + +**Unit.** `hsssTestSuite` already exercises a `"unit test"` path where +`FindGenesWithSnpCharsPlugin.initForBashScript` writes hardcoded gene filters rather than +querying. Extend those four rows with known site counts and assert exact expected values, +covering: +- a normal coding gene: both densities and the ratio +- `cds_length` NULL: CDS density blank, span density populated +- `syn_sites` NULL: ratio blank, both densities populated +- `synCount == 0` with `nonSynCount > 0`: existing exclusion branch still fires + +**Whole-cohort cross-check.** Run the search with every sample selected, then compare +against SQL over `GeneVariationSummary`: + +```sql +SELECT gene_source_id, + (n_missense::numeric / nullif(nonsyn_sites,0)) + / nullif(n_synonymous::numeric / nullif(syn_sites,0), 0) AS dnds_from_gvs +FROM apidbtuning.GeneVariationSummary +WHERE org_abbrev = 'pfal3D7'; +``` + +These will correlate strongly but WILL NOT match, for reasons that are correct: +- HSSS classifies from reference product bytes; `GeneVariationSummary` uses SnpEff + severity. Different callers, different edge cases. +- HSSS reports differences among the selected samples; with the reference strain omitted + from the selection it does not see sample-vs-reference differences at all. +- HSSS's `nonSynCount` excludes stop-gained (section 2); `n_missense` is SnpEff `sev=5`. + +Record the observed correlation so a future reader does not mistake the divergence for a +bug. + +**Coverage regression.** Assert that non-coding genes still appear in results with a +blank CDS density and a populated span density, and that they are excluded only once the +density filter is narrowed. + +## 6. Files touched + +`ApiCommonWebService` +- `WSFPlugin/src/main/java/org/apidb/apicomplexa/wsfplugin/highspeedsnpsearch/FindGenesWithSnpCharsPlugin.java` + — locations SQL, `makeResultRow`, `getColumns`, `COLUMN_*` constants, unit-test filters +- `HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter` — parse 7 columns, compute both + densities and the normalized ratio, print 9 fields +- `HighSpeedSnpSearch/lib/perl/HsssGeneCharsFilterScriptGenerator.pm` — usage text + (including the count-vs-percent error) +- `HighSpeedSnpSearch/bin/hsssTestSuite` — extended assertions + +`ApiCommonModel` +- `Model/lib/wdk/model/questions/queries/geneQueries.xml` — `wsColumn span_snp_density` +- `Model/lib/wdk/model/questions/geneQuestions.xml` — dynamic attribute for the new + column, display names and help per section 4.1, description per 4.3 +- `Model/lib/wdk/model/questions/params/geneParams.xml` — prompts and help per 4.2 + +Deployment note: the plugin is a jar. A model rebuild alone does not pick up the Java or +perl changes; both repos must ship together, or the 9-field output will meet an 8-field +assertion. diff --git a/docs/superpowers/specs/2026-08-08-hsss-gene-stats-validation.md b/docs/superpowers/specs/2026-08-08-hsss-gene-stats-validation.md new file mode 100644 index 00000000..2a092e5d --- /dev/null +++ b/docs/superpowers/specs/2026-08-08-hsss-gene-stats-validation.md @@ -0,0 +1,104 @@ +# Validation: corrected geneChars statistics + +2026-08-08. Acceptance evidence for +`2026-08-07-hsss-gene-stats-fix-design.md`, run against the `plasmo.jbrestel` dev +instance with appDb `unidb_shu_a`. + +## What was deployed + +All three halves must be present or the search fails; this run had all three: + +| piece | evidence | +|---|---| +| perl filter | `$GUS_HOME/bin/hsssGeneCharacteristicsFilter` md5 `dfee8785…` == checkout | +| plugin jar | `api-common-websvc-wsfplugin-1.0.0.jar` rebuilt via `bld ApiCommonWebService/WSFPlugin`, contains `span_snp_density` | +| model | `span_snp_density` in the `wsColumn` list and both `postCacheUpdateSql` blocks; webapp reloaded | + +The jar is NOT rebuilt by `wb model`. `bld ApiCommonWebService/WSFPlugin` followed by an +`instance_manager … reload` is what deploys it. + +## Run + +`GenesByNgsSnps`, organism *Plasmodium falciparum 3D7*, `variation_sample_meta` left at +its default `{"filters":[]}` (all samples), every threshold permissive +(`snp_class=All SNVs`, occurrence/ratio/density bounds unset). 5,590 transcripts +returned. + +## 1. dN/dS is exactly reproducible — 31/31, max error 0.000000 + +The design doc predicted only a *correlation* here, expecting divergence from different +classifiers. That was too pessimistic about the wrong thing: the classifier difference +affects the COUNTS, but this check tests the NORMALIZATION, and both sides draw site +counts from `apidbtuning.GeneVariationSummary`. So it is exact, not correlated. + +For 31 sampled genes, comparing the search's reported `ngs_dn_ds_ratio` against +`(nonsyn/nonsyn_sites) / (syn/syn_sites)` computed in SQL from the tuning table, using +the search's own synonymous and missense counts: + +``` + genes | dnds_exact | missing_sites | max_abs_err + 31 | 31 | 0 | 0.000000 +``` + +## 2. CDS density is exact + +`PF3D7_0100100`: 1,743 total variants, 99 unclassified, so 1,644 coding. +`cds_length` = 6,492. `1000 × 1644 / 6492 = 253.23`. Search reported **253.23**. + +Note this gene also demonstrates the defect being fixed: its span density is 228.86 +against a CDS density of 253.23. The old code reported the 228.86 figure under the label +"SNPs per Kb (CDS)". + +## 3. The biology moved the right way + +``` +dN/dS over 4,688 transcripts with a defined value + median 0.4718 + below 1 3,967 (85%) + max 8.3019 + exactly zero 218 +``` + +Median **0.4718**, with 85% under 1 — the expected purifying-selection signature. + +Two independent corroborations: +- `GeneVariationSummary` computes piN/piS by a completely different route (SnpEff + severity classes, frequency-weighted, Nei-Gojobori sites) and reports a pfal3D7 median + of **0.512**. Two pipelines, different classifiers, converging. +- The tuning table's design notes record that WITHOUT site normalization the median + piN/piS is **2.0**, implying genome-wide positive selection. That is the regime the old + un-normalized count ratio was in. + +The 218 genes at exactly 0 are the `defined($dn)` case: zero missense variants gives +dN = 0 and a real ratio of 0, the strongest purifying signal available. Truth-testing +`$dn` instead of `defined($dn)` would have silently converted all 218 into "no value". + +## 4. Coverage regression: non-coding genes survive + +254 returned transcripts have a blank CDS density. **All 254 have a populated span +density** — verified programmatically, not by inspection. Examples: + +``` +PF3D7_0100500 cds=(blank) span=45.05 dnds=(blank) total=5 +PF3D7_0101400 cds=(blank) span=18.14 dnds=(blank) total=15 +PF3D7_0101500 cds=(blank) span=69.70 dnds=(blank) total=115 +``` + +Consistent with the appDb: 283 pfal3D7 genes have `cds_length IS NULL`; 254 of them have +at least one variant and so appear. + +902 transcripts have no dN/dS — genes with no synonymous sites or no tuning-table row. +Both groups are returned because the corresponding filters were left unset; narrowing +either excludes them, which is the intended semantics. + +## What this run does NOT establish + +- **The filters were not exercised under narrowing.** Every bound was left permissive. + The exclusion paths (`return 0` when a statistic is undefined and the filter has been + narrowed) are covered by the unit fixture in `hsssTestSuite`, not here. +- **The counts themselves are unchanged and unverified by this work.** Only the two + normalized statistics were touched. HSSS's classification of a position as synonymous, + missense, nonsense or unclassified is exactly as before. +- **No comparison to pre-change output was made on this instance.** The old values are + knowable from the formula (`total/span` and `nonsyn/syn`) but were not captured from a + running pre-change instance.