Fix SNV density and dN/dS in the sample-set SNV search, and revive its test suite

HighSpeedSnpSearch/bin/hsssGeneCharacteristicsFilter

@@ -20,7 +20,8 @@ close(C);
 open(F, $geneLocationsFile) || die "Can't open gene locations file '$geneLocationsFile'\n";
 my $geneLocationLine = <F>;
 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(<STDIN>) {
       $geneLocationLine = <F>;
       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
 ";
 }

HighSpeedSnpSearch/bin/hsssGenomicLocationsFilter

@@ -48,7 +48,7 @@ while(<STDIN>) {
   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(<STDIN>) {
     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";
     }
   }
 }

HighSpeedSnpSearch/bin/hsssReconstructSnpId

@@ -39,8 +39,8 @@ while(<STDIN>) {
   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 {

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

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

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

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

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

HighSpeedSnpSearch/test/expected/polymorphismSearchWithSourceIdsAndSeqFilter.txt

@@ -1 +1 @@
-NGS_SNP.f100.23	100	20	
+Variant_f100_23	100.0	20.0	syn

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

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

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

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