From 55ac52e8942116235a4311ddea48249dbf31b09f Mon Sep 17 00:00:00 2001
From: "Mitchell R. Vollger" <mvollger@gmail.com>
Date: Fri, 14 Aug 2026 14:08:32 -0700
Subject: [PATCH] refactor: remove the unused leviosam2 chain-lift rules

Nothing consumes the leviosam2 outputs, chain mode has been broken since
the fire rule stopped producing a bam (its input referenced a bam output
that no longer exists), and the rules carried further latent bugs (index
written to a hardcoded shared path instead of the declared output).
Delete the rules and the chain/levio_exe config options.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 workflow/Snakefile       |  10 ----
 workflow/rules/levio.smk | 102 ---------------------------------------
 2 files changed, 112 deletions(-)
 delete mode 100644 workflow/rules/levio.smk

diff --git a/workflow/Snakefile b/workflow/Snakefile
index 967e69371..4cb4ca372 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -74,11 +74,6 @@ bw_types = ["log_FDR"]  # "score", "FDR",
 el_types = ["fire", "linker", "nucleosome"]
 
 
-# DSA options
-DSA_CHAIN = config.get("chain", None)
-DSA = DSA_CHAIN is not None
-LEVIO_EXE = config.get("levio_exe", "leviosam2")
-
 # developer options
 FT_EXE = config.get("ft_exe", "ft")
 if FT_EXE != "ft":
@@ -104,11 +99,6 @@ include: "rules/decorated-reads.smk"
 include: "rules/track-hub.smk"
 
 
-if DSA:
-
-    include: "rules/levio.smk"
-
-
 wildcard_constraints:
     chrom="|".join(get_chroms()),
     call="|".join(["msp", "m6a"]),
diff --git a/workflow/rules/levio.smk b/workflow/rules/levio.smk
deleted file mode 100644
index 409303abc..000000000
--- a/workflow/rules/levio.smk
+++ /dev/null
@@ -1,102 +0,0 @@
-#
-# Index the chain file for leviosam2.
-#
-# Input here is a chain file that defines the alignment between the DSA and the reference genome
-# at a contig level (>100 kbp of alignment).
-#
-# The output is a special leviosam2 index file that is used to lift over the alignments from the DSA to the reference genome.
-#
-rule leviosam2_index:
-    input:
-        chain=DSA_CHAIN,
-        fai=FAI,
-    output:
-        index=temp("temp/{sm}/leviosam2-index/index.clft"),
-    conda:
-        DEFAULT_ENV
-    threads: 1
-    resources:
-        mem_mb=64 * 1024,
-        runtime=16 * 60,
-    shell:
-        """
-        {LEVIO_EXE} index \
-            -p results/leviosam2-index/index \
-            -c {input.chain} \
-            -F {input.fai}
-        """
-
-
-#
-# Lift over the alignments from the DSA to the reference genome using the chain file / leviosam2 index.
-#
-# This is not a realignment, but a lift over of the reads from the DSA to the reference genome.
-#
-rule leviosam2:
-    input:
-        bam=rules.fire.output.bam,
-        levio_index=rules.leviosam2_index.output.index,
-        ref=REF,
-    output:
-        lifted=temp("temp/{sm}/leviosam2/{sm}-{chrom}-committed.bam"),
-        deferred=temp("temp/{sm}/leviosam2/{sm}-{chrom}-deferred.bam"),
-        unliftable=temp("temp/{sm}/leviosam2/{sm}-{chrom}-unliftable.bam"),
-    conda:
-        DEFAULT_ENV
-    threads: MAX_THREADS
-    resources:
-        mem_mb=MAX_THREADS * 4 * 1024,
-        runtime=16 * 60,
-    params:
-        # maximum number of CIGAR opts to change, also the max gap size that can be spanned
-        G=config.get("levio_G", 100_000),
-        # Using -S clipped_frac 0.05 means when a read has >5% clipped bases, it is deferred. A lower value is more stringent (by deferring more reads).
-        # aln_score is the minumum score before the alignment is lifted over
-        S=config.get(
-            "levio_S",
-            f"-S mapq:0 -S hdist:{100_000} -S isize:{100_000} -S clipped_frac:0.95 -S aln_score:100",
-        ),
-        # number of reads per thread
-        T=config.get("levio_T", 4 * 256),
-    shell:
-        """
-        PRE="temp/{wildcards.sm}/leviosam2/{wildcards.sm}-{wildcards.chrom}"
-        {LEVIO_EXE} lift -t {threads} -a {input.cram} \
-            -T {params.T} -G {params.G} {params.S} \
-            -C {input.levio_index} -p $PRE -f {input.ref} -m -O bam
-        """
-
-
-#
-# This step sorted the leviosam2 output and fixes some tags in the CRAM file.
-#
-# Specifically, the MAPQ is reset to 60 for all reads that were previously aligned to the DSA.
-# And the XS tag is set to zero for all reads that were aligned to the DSA.
-# This is a hueristic that we may need to return to in the future.
-#
-# Other tags and fields like CIGAR, bitflags, and MD are correctly updated by
-# leviosam2 during liftover.
-#
-rule leviosam2_sorted:
-    input:
-        lifted=rules.leviosam2.output.lifted,
-        ref=REF,
-    output:
-        bam=temp("temp/{sm}/leviosam2/{sm}-{chrom}-sorted.bam"),
-    conda:
-        DEFAULT_ENV
-    threads: SORT_THREADS
-    resources:
-        mem_mb=SORT_THREADS * 4 * 1024,
-        runtime=16 * 60,
-    shell:
-        """
-        samtools sort {input.lifted} \
-            -@ {threads} -m 3G \
-            -o {output.bam}
-        """
-
-
-# params:
-# reset_mapq=workflow.source_path("../scripts/reset-mapq.py"),
-# python {params.reset_mapq} -t {threads} {input.lifted} \