add logSED

src/baskerville/dataset.py

@@ -325,16 +325,25 @@ def make_strand_transform(targets_df, targets_strand_df):
     # initialize sparse matrix
     strand_transform = dok_matrix((targets_df.shape[0], targets_strand_df.shape[0]))
 
+    # track which strand pairs we've seen
+    seen_pairs = set()
+
     # fill in matrix
     ti = 0
     sti = 0
     for _, target in targets_df.iterrows():
         strand_transform[ti, sti] = True
         if target.strand_pair == target.name:
+            # Unstranded target
             sti += 1
         else:
-            if target.identifier[-1] == "-":
+            # Stranded target - check if we've seen its pair
+            if target.strand_pair in seen_pairs:
+                # This is the second member of the pair, increment sti
                 sti += 1
+            else:
+                # This is the first member of the pair, mark it as seen
+                seen_pairs.add(target.name)
         ti += 1
 
     return strand_transform

src/baskerville/scripts/hound_snp.py

@@ -32,6 +32,36 @@
 Compute variant effect predictions for SNPs in a VCF file.
 """
 
+"""
+Available statistics for variant effect prediction (--stats parameter):
+
+Basic Sum-based Scores:
+- SUM: Direct subtraction of alternative from reference predictions, summed across sequence length.
+- logSUM: First applies log2(x+1) to bin, then sums across sequence length and averages across shifts.
+- logSED: First sums across sequence length, then applies log2(x+1), and finally averages across shifts.
+- sqrtSUM: First square-root each bin, then sum across sequence length and averages across shifts.
+
+Maximum Effect Scores:
+- SAX: Maximum Absolute Effect - finds position of maximum difference and reports that value
+
+Distance/Norm-based Scores:
+- D1: L1 (Manhattan) distance between alt and ref predictions
+- logD1: L1 distance between log-transformed predictions
+- sqrtD1: L1 distance between sqrt-transformed predictions
+- D2: L2 (Euclidean) distance between alt and ref predictions
+- logD2: L2 distance between log-transformed predictions  
+- sqrtD2: L2 distance between sqrt-transformed predictions
+
+Information Theory Scores:
+- JS: Jensen-Shannon divergence between normalized predictions
+- logJS: Jensen-Shannon divergence between log-transformed normalized predictions
+
+Notes:
+- All scores are averaged across prediction shifts if multiple shifts are used
+- Log transformations use log2(x+1) to handle zeros
+- Scores can be specified as comma-separated list via --stats parameter
+- Default statistic is "logSUM"
+"""
 
 def main():
     usage = "usage: %prog [options] <params_file> <model_file> <vcf_file>"

src/baskerville/scripts/hound_snpgene.py

@@ -68,7 +68,7 @@ def main():
     parser.add_option(
         "-g",
         dest="genes_gtf",
-        default="%s/genes/gencode41/gencode41_basic_nort.gtf" % os.environ["HG38"],
+        default="%s/genes/gencode41/gencode41_basic_nort.gtf" % os.environ.get("HG38", "/path/to/hg38"),
         help="GTF for gene definition [Default %default]",
     )
     parser.add_option(

src/baskerville/snps.py

@@ -624,6 +624,10 @@ def compute_scores(ref_preds, alt_preds, snp_stats, strand_transform=None):
     ref_preds_sqrt_sum = ref_preds_sqrt.sum(axis=(0, 1)) / num_shifts
     alt_preds_sqrt_sum = alt_preds_sqrt.sum(axis=(0, 1)) / num_shifts
 
+    # SED computation (sum first, then log)
+    ref_preds_log_sed = np.log2(ref_preds.sum(axis=1) + 1).mean(axis=0)
+    alt_preds_log_sed = np.log2(alt_preds.sum(axis=1) + 1).mean(axis=0)
+
     # difference
     altref_diff = alt_preds - ref_preds
     altref_log_diff = alt_preds_log - ref_preds_log
@@ -650,6 +654,9 @@ def strand_clip_save(key, score, d2=False):
     if "logSUM" in snp_stats:
         log_sad = alt_preds_log_sum - ref_preds_log_sum
         strand_clip_save("logSUM", log_sad)
+    if "logSED" in snp_stats:
+        log_sed = alt_preds_log_sed - ref_preds_log_sed
+        strand_clip_save("logSED", log_sed)
     if "sqrtSUM" in snp_stats:
         sqrt_sad = alt_preds_sqrt_sum - ref_preds_sqrt_sum
         strand_clip_save("sqrtSUM", sqrt_sad)