From ce2067534fc6049028f0441878e624eff20ab69b Mon Sep 17 00:00:00 2001
From: marco cominelli <cominelli@ingm.org>
Date: Wed, 21 Jan 2026 16:26:23 +0000
Subject: [PATCH 01/11] Implement fragments aggregation into instances:
 modified main.py and map.py

---
 irescue/main.py | 17 ++++++-----
 irescue/map.py  | 76 ++++++++++++++++++++++++++++++++++++++++++-------
 2 files changed, 76 insertions(+), 17 deletions(-)

diff --git a/irescue/main.py b/irescue/main.py
index dc7d5e5..8fff610 100644
--- a/irescue/main.py
+++ b/irescue/main.py
@@ -100,12 +100,15 @@ def parseArguments():
         help="BAM tag containing the UMI sequence (default: %(default)s).",
     )
     parser.add_argument(
-        "-l",
-        "--locus",
-        action="store_true",
+        "--locus-level",
+        type=str,
+        metavar="STR",
+        choices=["disabled", "fragment", "instance"],
+        default="disabled", 
         help=(
-            "Perform locus-level quantification, instead of subfamily-level"
-            " (default: %(default)s)."
+            "Perform locus-level quantification. "
+            "One of: disabled, fragment, instance (default: %(default)s). "
+            "disabled means subfamily-level quantification."
         ),
     )
     parser.add_argument(
@@ -306,7 +309,7 @@ def main():
         genome=args.genome,
         genomes=__genomes__,
         outdir=dirs["out"],
-        locus=args.locus,
+        locus=args.locus_level,
         outname="rmsk.bed.gz",
     )
 
@@ -353,7 +356,7 @@ def main():
         threads=args.threads,
         outdir=dirs["mex"],
         tmpdir=dirs["tmp"],
-        locus=args.locus,
+        locus=args.locus_level,
         bedtools=args.bedtools,
         verbose=args.verbose,
     )
diff --git a/irescue/map.py b/irescue/map.py
index 06ddecf..5d97eb7 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -34,9 +34,8 @@ def checkIndex(bamFile, verbose):
                     send=verbose,
                 )
 
-
 def makeRmsk(
-    regions, genome, genomes, outdir, locus=False, outname="rmsk.bed.gz"
+    regions, genome, genomes, outdir, locus="disabled", outname="rmsk.bed.gz"
 ):
     """Format and/or download RepeatMasker annotation.
 
@@ -49,7 +48,7 @@ def makeRmsk(
         genome (str): Genome assembly name.
         genomes (dict): Dictionary of genome assembly names and URLs.
         outdir (str): Path to output directory.
-        locus (bool): If True, prepare for locus-level quantification.
+        locus (str): If not disabled, prepare for locus-level quantification.
         outname (str): Name of the output repeatmasker bed file.
 
     Returns:
@@ -101,7 +100,7 @@ def rl(x, decode=False):
         out = os.path.join(outdir, outname)
         with gzip.GzipFile(out, "wb", mtime=0) as f:
             # print header
-            h = ["#chr", "start", "end", "name", "locus_index", "strand"]
+            h = ["#chr", "start", "end", "name", "repInstance", "strand"]
             h = "\t".join(h) + "\n"
             f.write(h.encode())
             # skip rmsk header
@@ -120,13 +119,14 @@ def rl(x, decode=False):
             for line in rmsk:
                 lst = line.decode("utf-8").strip().split()
                 strand, repname, famclass = lst[8:11]
+                repInstance = lst[-1]
                 if famclass.split("/")[0] in fams_to_skip:
                     continue
                 # concatenate family and class with subfamily
                 repname += "#" + famclass
                 subfamilies[repname] += 1
                 locus_index = subfamilies[repname]
-                if locus:
+                if locus != "disabled":
                     # make unique locus names
                     repname += f"~{locus_index}"
                 chr, start, end = lst[4:7]
@@ -135,20 +135,76 @@ def rl(x, decode=False):
                 if strand != "+":
                     strand = "-"
                 outl = "\t".join(
-                    [chr, start, end, repname, str(locus_index), strand]
+                    [chr, start, end, repname, str(repInstance), strand]
                 )
                 outl += "\n"
                 f.write(outl.encode())
-        writerr(f"Wrote RepeatMasker annotation to {out}.")
+        
+        if locus == "instance":
+            outname = "rmsk_instance.bed.gz"
+            out_instance = os.path.join(outdir, outname)
+            rmsk_dict = dict()
+            with gzip.open(out, "rb") as rmsk:
+                # skip header row
+                next(rmsk)
+                for line in rmsk:
+                    chr,start,end,name,repInstance,strand = line.decode("utf-8").strip().split("\t")
+                    repInstance = str(repInstance)
+                    repSubfam, class_fam_locus = name.split("#")
+                    class_fam = class_fam_locus.split("~")[0]
+                    class_fam = class_fam.split("/")
+                    # loci without repFamily annotated
+                    if len(class_fam) == 1:
+                        repClass = class_fam[0]
+                        repFam = ""
+                    else:
+                        repClass, repFam = class_fam
+                    
+                    if repInstance not in rmsk_dict.keys():
+                        rmsk_dict[repInstance] = dict({
+                            "repSubfam":set([repSubfam]),
+                            "repFam":set([repFam]),
+                            "repClass":set([repClass])
+                        })
+                    else:
+                        rmsk_dict[repInstance]["repSubfam"].add(repSubfam)
+                        rmsk_dict[repInstance]["repFam"].add(repFam)
+                        rmsk_dict[repInstance]["repClass"].add(repClass)
+        
+            rmsk = gzip.open(out, "rb")
+            with gzip.GzipFile(out_instance, "wb", mtime=0) as f:
+                next(rmsk)
+                # print header
+                h = ["#chr", "start", "end", "name", "repInstance", "strand"]
+                h = "\t".join(h) + "\n"
+                f.write(h.encode())
+                for line in rmsk:
+                    chr,start,end,name,repInstance,strand = line.decode("utf-8").strip().split("\t")
+                    repInstance = str(repInstance)
+                    repSubfam = "|".join(rmsk_dict[repInstance]["repSubfam"])
+                    repFam = "|".join(rmsk_dict[repInstance]["repFam"])
+                    repFam = "" if not repFam else "/"+repFam # to take into account cases in which repFamily is not annotated
+                    repClass = "|".join(rmsk_dict[repInstance]["repClass"])
+                    
+                    repname = repSubfam + "#" + repClass + repFam + "~" + str(repInstance)
+                    outl = "\t".join(
+                        [chr, start, end, repname, str(repInstance), strand]
+                    )
+                    outl += "\n"
+                    f.write(outl.encode())
+
+        writerr(f"Wrote RepeatMasker {"subfamily" if locus == "disabled" else locus}-level annotation to {out if locus != "instance" else out_instance}.")
     else:
         writerr(
             "Error: it is mandatory to define either --regions OR "
             "--genome parameter.",
             error=True,
         )
+    
+    if locus == "instance":
+        return out_instance
     return out
 
-
 def prepare_whitelist(whitelist, tmpdir):
     """Uncompress the whitelist file if compressed.
     Return the whitelist path, or False if not using a whitelist.
@@ -300,7 +356,7 @@ def chrcat(
     threads,
     outdir,
     tmpdir,
-    locus=False,
+    locus="disabled",
     bedtools="bedtools",
     verbose=0,
 ):
@@ -337,7 +393,7 @@ def chrcat(
     cmd2 = f"zcat {mappings_file} "
     cmd2 += " | cut -f3 | sed 's/,/\\n/g' | gawk '!x[$1]++ { "
     cmd2 += ' print $1"\\t"gensub(/#'
-    cmd2 += "[^~]" if locus else "."
+    cmd2 += "[^~]" if locus != "disabled" else "."
     cmd2 += '+/,"",1,$1)"\\tGene Expression" }\' '
     cmd2 += f" | LC_ALL=C sort -u | gzip > {features_file} "
 

From b8b3c3c41da6816b3aef555207bd9c581ee1ecfd Mon Sep 17 00:00:00 2001
From: Marco Cominelli <161143279+marco-cominelli01@users.noreply.github.com>
Date: Wed, 21 Jan 2026 18:29:15 +0100
Subject: [PATCH 02/11] Refactor logging for RepeatMasker annotation

---
 irescue/map.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/irescue/map.py b/irescue/map.py
index 5d97eb7..3ff6b70 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -193,7 +193,10 @@ def rl(x, decode=False):
                     outl += "\n"
                     f.write(outl.encode())
 
-        writerr(f"Wrote RepeatMasker {"subfamily" if locus == "disabled" else locus}-level annotation to {out if locus != "instance" else out_instance}.")
+        level = "subfamily" if locus == "disabled" else locus
+        outfile = out if locus != "instance" else out_instance
+        writerr(f"Wrote RepeatMasker {level}-level annotation to {outfile}.")
+
     else:
         writerr(
             "Error: it is mandatory to define either --regions OR "

From a02ca099fda3416fdd2f038d87d808f042e737d2 Mon Sep 17 00:00:00 2001
From: marco cominelli <cominelli@ingm.org>
Date: Thu, 22 Jan 2026 16:01:49 +0000
Subject: [PATCH 03/11] Modified test.yml for compatibility with new
 locus-level parameter

---
 tests/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test.yml b/tests/test.yml
index a743008..18c4975 100644
--- a/tests/test.yml
+++ b/tests/test.yml
@@ -129,7 +129,7 @@
 - name: locus
   tags:
     - locus
-  command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus
+  command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus-level fragment
   files:
     - path: "irescue_out/counts/barcodes.tsv.gz"
       md5sum: 1a74fa12e65ac1703bbe61282854f151

From 1dbbc4848352972e13def7b7f9acc4c81175bea2 Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Fri, 23 Jan 2026 16:42:37 +0100
Subject: [PATCH 04/11] tiny typo

---
 irescue/main.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/irescue/main.py b/irescue/main.py
index 8fff610..184a1c6 100644
--- a/irescue/main.py
+++ b/irescue/main.py
@@ -108,7 +108,7 @@ def parseArguments():
         help=(
             "Perform locus-level quantification. "
             "One of: disabled, fragment, instance (default: %(default)s). "
-            "disabled means subfamily-level quantification."
+            "Disabled means subfamily-level quantification."
         ),
     )
     parser.add_argument(

From 696ead6873113187a114b21644962fbd5429e9f7 Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Fri, 23 Jan 2026 16:45:07 +0100
Subject: [PATCH 05/11] Refactor makeRmsk function for improved readability and
 performance

---
 irescue/map.py | 148 +++++++++++++++++++++++--------------------------
 1 file changed, 68 insertions(+), 80 deletions(-)

diff --git a/irescue/map.py b/irescue/map.py
index 3ff6b70..0ae273c 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -34,9 +34,8 @@ def checkIndex(bamFile, verbose):
                     send=verbose,
                 )
 
-def makeRmsk(
-    regions, genome, genomes, outdir, locus="disabled", outname="rmsk.bed.gz"
-):
+
+def makeRmsk(regions, genome, genomes, outdir, locus="disabled", outname="rmsk.bed.gz"):
     """Format and/or download RepeatMasker annotation.
 
     Check repeatmasker regions bed file format. Download if not provided.
@@ -96,18 +95,20 @@ def rl(x, decode=False):
                 f"Couldn't connect to host.\n\n{e}",
                 error=True,
             )
-        rmsk = gzip.open(io.BytesIO(response.content), "rb")
         out = os.path.join(outdir, outname)
-        with gzip.GzipFile(out, "wb", mtime=0) as f:
+        with (
+            gzip.open(io.BytesIO(response.content), "rb") as rmsk_in,
+            gzip.GzipFile(out, "wb", mtime=0) as rmsk_out,
+        ):
             # print header
             h = ["#chr", "start", "end", "name", "repInstance", "strand"]
             h = "\t".join(h) + "\n"
-            f.write(h.encode())
+            rmsk_out.write(h.encode())
             # skip rmsk header
             for _ in range(header_lines):
-                next(rmsk)
+                next(rmsk_in)
             # parse rmsk
-            fams_to_skip = [
+            repeats_to_skip = [
                 "Low_complexity",
                 "Simple_repeat",
                 "rRNA",
@@ -116,82 +117,72 @@ def rl(x, decode=False):
                 "tRNA",
             ]
             subfamilies = defaultdict(int)
-            for line in rmsk:
-                lst = line.decode("utf-8").strip().split()
-                strand, repname, famclass = lst[8:11]
-                repInstance = lst[-1]
-                if famclass.split("/")[0] in fams_to_skip:
+            rmsk_dict = dict()
+            for line in rmsk_in:
+                fields = line.decode("utf-8").strip().split()
+                chr, start, end = fields[4:7]
+                strand, subfamily, clasfam = fields[8:11]
+                instance = fields[14]
+                clas, family = clasfam.split("/") if "/" in clasfam else (clasfam, "")
+                # skip problematic repeats
+                if clas in repeats_to_skip:
                     continue
                 # concatenate family and class with subfamily
-                repname += "#" + famclass
-                subfamilies[repname] += 1
-                locus_index = subfamilies[repname]
+                name = f"{subfamily}#{clasfam}"
+
+                # repeat names for locus-level quantifications
                 if locus != "disabled":
-                    # make unique locus names
-                    repname += f"~{locus_index}"
-                chr, start, end = lst[4:7]
-                # make coordinates 0-based
-                start = str(int(start) - 1)
+                    subfamilies[name] += 1
+                    locus_index = subfamilies[name]
+                    name += f"~{locus_index}"
+                    if locus == "instance":
+                        if instance in rmsk_dict.keys():
+                            rmsk_dict[instance]["repSubfam"].add(subfamily)
+                            rmsk_dict[instance]["repFam"].add(family)
+                            rmsk_dict[instance]["repClass"].add(clas)
+                        else:
+                            rmsk_dict[instance] = {
+                                "repSubfam": {subfamily},
+                                "repFam": {family},
+                                "repClass": {clas},
+                            }
+
+                start = str(int(start) - 1)  # make coordinates 0-based
                 if strand != "+":
                     strand = "-"
-                outl = "\t".join(
-                    [chr, start, end, repname, str(repInstance), strand]
-                )
+                outl = "\t".join([chr, start, end, name, instance, strand])
                 outl += "\n"
-                f.write(outl.encode())
-        
+                rmsk_out.write(outl.encode())
+
+        # if locus-level quantification by instance is requested,
+        # prepare a separate file with instance-level annotation
         if locus == "instance":
             outname = "rmsk_instance.bed.gz"
             out_instance = os.path.join(outdir, outname)
-            rmsk_dict = dict()
-            with gzip.open(out, "rb") as rmsk:
-                # skip header row
-                next(rmsk)
-                for line in rmsk:
-                    chr,start,end,name,repInstance,strand = line.decode("utf-8").strip().split("\t")
-                    repInstance = str(repInstance)
-                    repSubfam, class_fam_locus = name.split("#")
-                    class_fam = class_fam_locus.split("~")[0]
-                    class_fam = class_fam.split("/")
-                    # loci without repFamily annotated
-                    if len(class_fam) == 1:
-                        repClass = class_fam[0]
-                        repFam = ""
-                    else:
-                        repClass, repFam = class_fam
-                    
-                    if repInstance not in rmsk_dict.keys():
-                        rmsk_dict[repInstance] = dict({
-                            "repSubfam":set([repSubfam]),
-                            "repFam":set([repFam]),
-                            "repClass":set([repClass])
-                        })
-                    else:
-                        rmsk_dict[repInstance]["repSubfam"].add(repSubfam)
-                        rmsk_dict[repInstance]["repFam"].add(repFam)
-                        rmsk_dict[repInstance]["repClass"].add(repClass)
-        
-            rmsk = gzip.open(out, "rb")
-            with gzip.GzipFile(out_instance, "wb", mtime=0) as f:
-                next(rmsk)
+            with (
+                gzip.open(out, "rb") as rmsk_in,
+                gzip.GzipFile(out_instance, "wb", mtime=0) as rmsk_out,
+            ):
+                next(rmsk_in)  # skip header row
                 # print header
                 h = ["#chr", "start", "end", "name", "repInstance", "strand"]
                 h = "\t".join(h) + "\n"
-                f.write(h.encode())
-                for line in rmsk:
-                    chr,start,end,name,repInstance,strand = line.decode("utf-8").strip().split("\t")
-                    repInstance = str(repInstance)
-                    repSubfam = "|".join(rmsk_dict[repInstance]["repSubfam"])
-                    repFam = "|".join(rmsk_dict[repInstance]["repFam"])
-                    repFam = "" if not repFam else "/"+repFam # to take into account cases in which repFamily is not annotated
-                    repClass = "|".join(rmsk_dict[repInstance]["repClass"])
-                    
-                    repname = repSubfam + "#" + repClass + repFam + "~" + str(repInstance)
-                    outl = "\t".join(
-                        [chr, start, end, repname, str(repInstance), strand]
+                rmsk_out.write(h.encode())
+                for line in rmsk_in:
+                    chr, start, end, name, instance, strand = (
+                        line.decode("utf-8").strip().split("\t")
                     )
+                    repSubfam = "|".join(sorted(rmsk_dict[instance]["repSubfam"]))
+                    repFam = "|".join(sorted(rmsk_dict[instance]["repFam"]))
+                    repFam = (
+                        "" if not repFam else "/" + repFam
+                    )  # to take into account cases in which repFamily is not annotated
+                    repClass = "|".join(sorted(rmsk_dict[instance]["repClass"]))
+
+                    repname = f"{repSubfam}#{repClass}{repFam}~{instance}"
+                    outl = "\t".join([chr, start, end, repname, instance, strand])
                     outl += "\n"
-                    f.write(outl.encode())
+                    rmsk_out.write(outl.encode())
 
         level = "subfamily" if locus == "disabled" else locus
         outfile = out if locus != "instance" else out_instance
@@ -199,14 +190,15 @@ def rl(x, decode=False):
 
     else:
         writerr(
-            "Error: it is mandatory to define either --regions OR "
-            "--genome parameter.",
+            "Error: it is mandatory to define either --regions OR --genome parameter.",
             error=True,
         )
-    
+
     if locus == "instance":
         return out_instance
-    return out
+    else:
+        return out
+
 
 def prepare_whitelist(whitelist, tmpdir):
     """Uncompress the whitelist file if compressed.
@@ -301,15 +293,11 @@ def isec(
     else:
         stream = f" <({samtools} view -h {bamFile} {chrom} | "
     stream += ' gawk \'!($1~/^@/) { split("", tags); '
-    stream += (
-        ' for (i=12;i<=NF;i++) {split($i,tag,":"); tags[tag[1]]=tag[3]}; '
-    )
+    stream += ' for (i=12;i<=NF;i++) {split($i,tag,":"); tags[tag[1]]=tag[3]}; '
     # Discard records without CB tag, unvalid STARSolo CBs, missing UMI tag,
     # UMIs with Ns and homopolymer UMIs
     if UMItag:
-        stream += (
-            f' if(tags["{CBtag}"]~/^(|-)$/ || tags["{UMItag}"]~/.*N.*/ || '
-        )
+        stream += f' if(tags["{CBtag}"]~/^(|-)$/ || tags["{UMItag}"]~/.*N.*/ || '
         stream += f' tags["{UMItag}"]~/^$|^(A+|G+|T+|C+)$/) {{next}}; '
         # Append CB and UMI to read name
         stream += f' $1=$1"/"tags["{CBtag}"]"/"tags["{UMItag}"]; '

From 77201199537eb39a2f383cf50c8694db0f200faf Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Fri, 23 Jan 2026 16:45:36 +0100
Subject: [PATCH 06/11] update md5sums

---
 tests/test.yml | 26 +++++++++++++++++++++++---
 1 file changed, 23 insertions(+), 3 deletions(-)

diff --git a/tests/test.yml b/tests/test.yml
index 18c4975..88cb7d3 100644
--- a/tests/test.yml
+++ b/tests/test.yml
@@ -28,7 +28,7 @@
     - path: "irescue_out/tmp/mappings.tsv.gz"
       md5sum: d404e6c3123f8cfde7689b5fb7763a89
     - path: "irescue_out/rmsk.bed.gz"
-      md5sum: 169571f538624496a00f189771be2f5e
+      md5sum: 663134e9bb5519d136201d2bd5d35d63
 
 - name: multi
   tags:
@@ -126,7 +126,7 @@
     - path: "irescue_out/tmp/mappings.tsv.gz"
       md5sum: f5aae354f59bef7a4bdb9cf3c5c8dafe
 
-- name: locus
+- name: fragment
   tags:
     - locus
   command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus-level fragment
@@ -142,4 +142,24 @@
     - path: "irescue_out/tmp/mappings.tsv.gz"
       md5sum: d7db121c92c04b36e3cc26ae0223426d
     - path: "irescue_out/rmsk.bed.gz"
-      md5sum: 75ed0333b049a02672da8b1379cfa6bc
+      md5sum: 050b0a484a04e5d466afd54c794e35fd
+
+- name: instance
+  tags:
+    - locus
+  command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus-level instance
+  files:
+    - path: "irescue_out/counts/barcodes.tsv.gz"
+      md5sum: 1a74fa12e65ac1703bbe61282854f151
+    - path: "irescue_out/counts/features.tsv.gz"
+      md5sum: 045273e35f85522ec70a7be127fd5d7f
+    - path: "irescue_out/counts/matrix.mtx.gz"
+      md5sum: a63c59d45b79bd4bfae1d7ffb26b6c6e
+    - path: "irescue_out/ec_dump.tsv.gz"
+      md5sum: dc7a39ce64959f0b83dd53c158041746
+    - path: "irescue_out/tmp/mappings.tsv.gz"
+      md5sum: e0986db4bea3929a8f14eb61f2596a41
+    - path: "irescue_out/rmsk.bed.gz"
+      md5sum: 050b0a484a04e5d466afd54c794e35fd
+    - path: "irescue_out/rmsk_instance.bed.gz"
+      md5sum: 7ba1ba91b19292a6d49a34e818470630

From 0f621af55ed03fe54d5b2f4aea4f0f0906627623 Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Mon, 26 Jan 2026 10:35:03 +0100
Subject: [PATCH 07/11] [fix

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 irescue/map.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/irescue/map.py b/irescue/map.py
index 0ae273c..2e586fc 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -173,7 +173,9 @@ def rl(x, decode=False):
                         line.decode("utf-8").strip().split("\t")
                     )
                     repSubfam = "|".join(sorted(rmsk_dict[instance]["repSubfam"]))
-                    repFam = "|".join(sorted(rmsk_dict[instance]["repFam"]))
+                    repFam = "|".join(
+                        sorted(f for f in rmsk_dict[instance]["repFam"] if f)
+                    )
                     repFam = (
                         "" if not repFam else "/" + repFam
                     )  # to take into account cases in which repFamily is not annotated

From d927b1bd4ab06a25de9019352c215d714f9c63fe Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Mon, 26 Jan 2026 10:49:53 +0100
Subject: [PATCH 08/11] Remove trailing whitespace

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 irescue/main.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/irescue/main.py b/irescue/main.py
index 184a1c6..226197e 100644
--- a/irescue/main.py
+++ b/irescue/main.py
@@ -104,7 +104,7 @@ def parseArguments():
         type=str,
         metavar="STR",
         choices=["disabled", "fragment", "instance"],
-        default="disabled", 
+        default="disabled",
         help=(
             "Perform locus-level quantification. "
             "One of: disabled, fragment, instance (default: %(default)s). "

From 33cacd63fde7586ffaa260ce9ad4897cb4b5c5f7 Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Mon, 26 Jan 2026 10:50:25 +0100
Subject: [PATCH 09/11] Remove redundant call

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 irescue/map.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/irescue/map.py b/irescue/map.py
index 2e586fc..79c7439 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -136,7 +136,7 @@ def rl(x, decode=False):
                     locus_index = subfamilies[name]
                     name += f"~{locus_index}"
                     if locus == "instance":
-                        if instance in rmsk_dict.keys():
+                        if instance in rmsk_dict:
                             rmsk_dict[instance]["repSubfam"].add(subfamily)
                             rmsk_dict[instance]["repFam"].add(family)
                             rmsk_dict[instance]["repClass"].add(clas)

From b1b539c435259d1648b8cd479231a34ab2a5921d Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Mon, 26 Jan 2026 11:13:03 +0100
Subject: [PATCH 10/11] Fix critical bug caused by --regions and --locus-level
 instance parameters in makeRmsk

---
 irescue/map.py | 31 ++++++++++++++-----------------
 1 file changed, 14 insertions(+), 17 deletions(-)

diff --git a/irescue/map.py b/irescue/map.py
index 79c7439..bd36e09 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -43,7 +43,7 @@ def makeRmsk(regions, genome, genomes, outdir, locus="disabled", outname="rmsk.b
 
     Args:
         regions (str): Path to repeatmasker bed file.
-            Takes priority over genome.
+            Ignores "genome", "genomes" and "locus" parameters.
         genome (str): Genome assembly name.
         genomes (dict): Dictionary of genome assembly names and URLs.
         outdir (str): Path to output directory.
@@ -57,8 +57,9 @@ def makeRmsk(regions, genome, genomes, outdir, locus="disabled", outname="rmsk.b
         SystemExit: If neither regions nor genome is provided, or if the
                     regions file is not properly formatted.
     """
-    # if a repeatmasker bed file is provided, use that
+
     if regions:
+    # if a repeatmasker bed file is provided, use that
         is_gz = testGz(regions)
         f = gzip.open(regions, "rb") if is_gz else open(regions, "r")
 
@@ -73,15 +74,18 @@ def rl(x, decode=False):
         if len(line.strip().split("\t")) < 4:
             writerr(
                 "Error: please provide a tab-separated BED file with at least"
-                " 4 columns and TE feature name (e.g. locus or subfamily)"
-                " in 4th column.",
+                " 4 columns with the 4th column being the TE feature name, "
+                "such as the TE subfamily or locus "
+                "(as in fragment or instance ID).",
                 error=True,
             )
         f.close()
-        out = regions
+        writerr(f"Using provided RepeatMasker annotation from {regions}.")
+        return regions
+
+    elif genome:
     # if no repeatmasker file is provided, and a genome assembly name is
     # provided, download and prepare a rmsk.bed file
-    elif genome:
         url, header_lines = genomes[genome]
         writerr(
             "Downloading and parsing RepeatMasker annotation for "
@@ -153,12 +157,11 @@ def rl(x, decode=False):
                 outl = "\t".join([chr, start, end, name, instance, strand])
                 outl += "\n"
                 rmsk_out.write(outl.encode())
-
+        if locus == "instance":
         # if locus-level quantification by instance is requested,
         # prepare a separate file with instance-level annotation
-        if locus == "instance":
-            outname = "rmsk_instance.bed.gz"
-            out_instance = os.path.join(outdir, outname)
+            outname_instance = "rmsk_instance.bed.gz"
+            out_instance = os.path.join(outdir, outname_instance)
             with (
                 gzip.open(out, "rb") as rmsk_in,
                 gzip.GzipFile(out_instance, "wb", mtime=0) as rmsk_out,
@@ -180,15 +183,14 @@ def rl(x, decode=False):
                         "" if not repFam else "/" + repFam
                     )  # to take into account cases in which repFamily is not annotated
                     repClass = "|".join(sorted(rmsk_dict[instance]["repClass"]))
-
                     repname = f"{repSubfam}#{repClass}{repFam}~{instance}"
                     outl = "\t".join([chr, start, end, repname, instance, strand])
                     outl += "\n"
                     rmsk_out.write(outl.encode())
-
         level = "subfamily" if locus == "disabled" else locus
         outfile = out if locus != "instance" else out_instance
         writerr(f"Wrote RepeatMasker {level}-level annotation to {outfile}.")
+        return outfile
 
     else:
         writerr(
@@ -196,11 +198,6 @@ def rl(x, decode=False):
             error=True,
         )
 
-    if locus == "instance":
-        return out_instance
-    else:
-        return out
-
 
 def prepare_whitelist(whitelist, tmpdir):
     """Uncompress the whitelist file if compressed.

From 19862a770b0ea771b9a6f07d1f997ca9582830f5 Mon Sep 17 00:00:00 2001
From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com>
Date: Mon, 26 Jan 2026 11:17:10 +0100
Subject: [PATCH 11/11] syntax

---
 irescue/map.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/irescue/map.py b/irescue/map.py
index bd36e09..bda5a93 100644
--- a/irescue/map.py
+++ b/irescue/map.py
@@ -59,7 +59,7 @@ def makeRmsk(regions, genome, genomes, outdir, locus="disabled", outname="rmsk.b
     """
 
     if regions:
-    # if a repeatmasker bed file is provided, use that
+        # if a repeatmasker bed file is provided, use that
         is_gz = testGz(regions)
         f = gzip.open(regions, "rb") if is_gz else open(regions, "r")
 
@@ -84,8 +84,8 @@ def rl(x, decode=False):
         return regions
 
     elif genome:
-    # if no repeatmasker file is provided, and a genome assembly name is
-    # provided, download and prepare a rmsk.bed file
+        # if no repeatmasker file is provided, and a genome assembly name is
+        # provided, download and prepare a rmsk.bed file
         url, header_lines = genomes[genome]
         writerr(
             "Downloading and parsing RepeatMasker annotation for "
@@ -158,8 +158,8 @@ def rl(x, decode=False):
                 outl += "\n"
                 rmsk_out.write(outl.encode())
         if locus == "instance":
-        # if locus-level quantification by instance is requested,
-        # prepare a separate file with instance-level annotation
+            # if locus-level quantification by instance is requested,
+            # prepare a separate file with instance-level annotation
             outname_instance = "rmsk_instance.bed.gz"
             out_instance = os.path.join(outdir, outname_instance)
             with (