Add aggressive odds ratio pre-filtering by BH-adjusted p-value

Add filter_odds_ratios.py script and Snakemake rule to pre-filter
odds ratio files using total_obs >= 100 AND p_adjusted < 0.01 (BH FDR).
This reduces filtered file sizes ~20x (1.3GB -> 63MB across 71 samples)
and speeds up downstream R analysis. Also bumps modkit memory to 192GB
and reduces modkit_extract_full threads to 4.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: jayhesselberth

cluster/slurm/config.yaml

@@ -62,14 +62,14 @@ set-resources:
   # Modkit modification extraction - high memory
   modkit_extract_calls:
     runtime: 360               # 6 hours
-    mem_mb: 96000              # 96GB
+    mem_mb: 192000             # 192GB
     cpus_per_task: 4
 
   # Modkit full modification extraction
   modkit_extract_full:
     runtime: 360               # 6 hours
-    mem_mb: 96000              # 96GB
-    cpus_per_task: 12
+    mem_mb: 192000             # 192GB
+    cpus_per_task: 4
 
   # WarpDemuX demultiplexing rules
   warpdemux:

workflow/rules/aatrnaseq-modifications.smk

@@ -150,7 +150,7 @@ rule modkit_extract_full:
         tsv=os.path.join(
             outdir, "summary", "modkit", "{sample}", "{sample}.mod_full.tsv.gz"
         ),
-    threads: 12
+    threads: 4
     log:
         os.path.join(outdir, "logs", "modkit", "extract_full", "{sample}"),
     params:

workflow/rules/aatrnaseq-odds-ratios.smk

@@ -34,3 +34,38 @@ rule compute_odds_ratios:
             --min-coverage {params.min_cov} \
             2>&1 | tee {log}
         """
+
+
+rule filter_odds_ratios:
+    """
+    Pre-filter odds ratios to keep only well-observed position pairs.
+
+    Keeps rows where total_obs >= min_obs (default 100), matching the
+    downstream clover::filter_linkages() threshold. This reduces file
+    sizes ~10x and speeds up R data loading.
+    """
+    input:
+        tsv=rules.compute_odds_ratios.output.tsv,
+    output:
+        tsv=os.path.join(
+            outdir,
+            "summary",
+            "tables",
+            "{sample}",
+            "{sample}.odds_ratios_filtered.tsv.gz",
+        ),
+    log:
+        os.path.join(outdir, "logs", "odds_ratios", "{sample}.filter"),
+    params:
+        src=SCRIPT_DIR,
+        min_obs=config.get("odds_ratios", {}).get("min_obs_filter", 100),
+        max_p=config.get("odds_ratios", {}).get("max_p_filter", 0.01),
+    shell:
+        """
+        python {params.src}/filter_odds_ratios.py \
+            --input {input.tsv} \
+            --output {output.tsv} \
+            --min-obs {params.min_obs} \
+            --max-p {params.max_p} \
+            2>&1 | tee {log}
+        """

workflow/rules/common.smk

@@ -268,12 +268,12 @@ def pipeline_outputs():
         sample=samples.keys(),
     )
 
-    outs += expand(
-        os.path.join(
-            outdir, "summary", "modkit", "{sample}", "{sample}.mod_full.tsv.gz"
-        ),
-        sample=samples.keys(),
-    )
+    # outs += expand(
+    #     os.path.join(
+    #         outdir, "summary", "modkit", "{sample}", "{sample}.mod_full.tsv.gz"
+    #     ),
+    #     sample=samples.keys(),
+    # )
 
     if (
         "remora_kmer_table" in config
@@ -294,6 +294,17 @@ def pipeline_outputs():
         sample=samples.keys(),
     )
 
+    outs += expand(
+        os.path.join(
+            outdir,
+            "summary",
+            "tables",
+            "{sample}",
+            "{sample}.odds_ratios_filtered.tsv.gz",
+        ),
+        sample=samples.keys(),
+    )
+
     # tRNA-only reference FASTA (adapters stripped)
     outs.append(get_trna_fasta())
 

workflow/scripts/filter_odds_ratios.py

@@ -0,0 +1,67 @@
+"""
+Filter per-sample odds ratio files to keep only well-observed position pairs.
+
+Reads a (potentially large) odds_ratios.tsv.gz in chunks and writes only rows
+where total_obs >= min_obs AND p_adjusted < max_p.  This pre-filter matches the
+downstream threshold used by clover::filter_linkages() and drastically reduces
+file size / R load time.
+"""
+
+import argparse
+import sys
+
+import pandas as pd
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Filter odds ratios by minimum observation count."
+    )
+    parser.add_argument("--input", required=True, help="Input .odds_ratios.tsv.gz")
+    parser.add_argument("--output", required=True, help="Output filtered .tsv.gz")
+    parser.add_argument(
+        "--min-obs",
+        type=int,
+        default=100,
+        help="Minimum total_obs to keep a row (default: 100)",
+    )
+    parser.add_argument(
+        "--max-p",
+        type=float,
+        default=0.01,
+        help="Maximum p_adjusted (BH FDR) to keep a row (default: 0.01)",
+    )
+    args = parser.parse_args()
+
+    kept = 0
+    total = 0
+    first_chunk = True
+
+    for chunk in pd.read_csv(
+        args.input, sep="\t", compression="gzip", chunksize=100_000
+    ):
+        total += len(chunk)
+        filtered = chunk[
+            (chunk["total_obs"] >= args.min_obs)
+            & (chunk["p_adjusted"] < args.max_p)
+        ]
+        kept += len(filtered)
+        filtered.to_csv(
+            args.output,
+            sep="\t",
+            index=False,
+            compression="gzip",
+            mode="w" if first_chunk else "a",
+            header=first_chunk,
+        )
+        first_chunk = False
+
+    print(
+        f"Filtered odds ratios: kept {kept:,} / {total:,} rows "
+        f"(total_obs >= {args.min_obs}, p_adjusted < {args.max_p})",
+        file=sys.stderr,
+    )
+
+
+if __name__ == "__main__":
+    main()