Add long-read barcode assignment via pbmm2

Author: BennyKrup

config/example_assignment_pbmm2.yaml

@@ -0,0 +1,23 @@
+version: "0.6"
+
+assignments:
+  exampleLongRead:
+    bc_length: 15
+    long_read_input: "path/to/reads.bam"
+    design_file: "path/to/design.fa"
+    alignment_tool:
+      split_number: 1
+      tool: "pbmm2"
+      configs:
+        preset: "SUBREAD"
+        min_concordance: 100.0
+        pattern: "GCAAAGTGAACACATCGCTAAGCGAAAGCTAAG"
+    design_check:
+      fast: true
+      sequence_collisions: false
+    strand_sensitive:
+      enable: false
+    configs:
+      default:
+        min_support: 1
+        fraction: 0.0

workflow/envs/pbmm2_pysam.yaml

@@ -0,0 +1,8 @@
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - pbmm2
+  - pysam
+  - biopython
+  - python>=3.10

workflow/rules/assignment/mapping_pbmm2.smk

@@ -0,0 +1,76 @@
+rule assignment_mapping_pbmm2_index:
+    """
+    Create pbmm2 index from design reference.
+    """
+    conda:
+        getCondaEnv("pbmm2_pysam.yaml")
+    input:
+        ref="results/assignment/{assignment}/reference/reference.fa",
+        check="results/assignment/{assignment}/design_check.done",
+    output:
+        "results/assignment/{assignment}/reference/reference.fa.mmi",
+    log:
+        temp("results/logs/assignment/mapping_pbmm2_index.{assignment}.log"),
+    shell:
+        """
+        pbmm2 index {input.ref} {output} &> {log}
+        """
+
+
+rule assignment_mapping_pbmm2_align:
+    """
+    Align long reads (BAM or FASTA) to reference using pbmm2.
+    """
+    conda:
+        getCondaEnv("pbmm2_pysam.yaml")
+    threads: 8
+    input:
+        reads=lambda wc: config["assignments"][wc.assignment]["long_read_input"],
+        index="results/assignment/{assignment}/reference/reference.fa.mmi",
+    output:
+        temp("results/assignment/{assignment}/pbmm2/aligned.bam"),
+    params:
+        preset=lambda wc: config["assignments"][wc.assignment]["alignment_tool"]["configs"].get("preset", "SUBREAD"),
+        min_concordance=lambda wc: config["assignments"][wc.assignment]["alignment_tool"]["configs"].get("min_concordance", 100.0),
+    log:
+        temp("results/logs/assignment/mapping_pbmm2_align.{assignment}.log"),
+    shell:
+        """
+        pbmm2 align {input.index} {input.reads} {output} \
+            --preset {params.preset} \
+            --sort \
+            --best-n 1 \
+            --min-concordance-perc {params.min_concordance} \
+            -j {threads} &> {log}
+        """
+
+
+rule assignment_mapping_pbmm2_getBCs:
+    """
+    Extract barcodes from aligned long reads. Produces the standard
+    barcode TSV for downstream collection and filtering.
+    """
+    conda:
+        getCondaEnv("pbmm2_pysam.yaml")
+    input:
+        bam="results/assignment/{assignment}/pbmm2/aligned.bam",
+        script=getScript("assignment/longread_extract.py"),
+    output:
+        temp("results/assignment/{assignment}/BCs/barcodes.pbmm2.0.tsv"),
+    params:
+        pattern=lambda wc: config["assignments"][wc.assignment]["alignment_tool"]["configs"].get(
+            "pattern", "GCAAAGTGAACACATCGCTAAGCGAAAGCTAAG"
+        ),
+        bc_length=lambda wc: config["assignments"][wc.assignment]["bc_length"],
+    log:
+        temp("results/logs/assignment/mapping_pbmm2_getBCs.{assignment}.log"),
+    shell:
+        """
+        export LC_ALL=C
+        python {input.script} \
+            --bam {input.bam} \
+            --pattern {params.pattern} \
+            --bc-length {params.bc_length} \
+            --output /dev/stdout 2> {log} | \
+        sort -k1,1 -k2,2 -k3,3 -S 7G > {output}
+        """

workflow/scripts/assignment/longread_extract.py

@@ -0,0 +1,74 @@
+"""Extract barcodes from aligned long-read BAM.
+
+Outputs a 3-column TSV (barcode, reference, metadata) matching the
+MPRAsnakeflow assignment format used by assignment_collectBCs and
+assignment_filter.
+"""
+
+import argparse
+import logging
+import sys
+
+import pysam
+
+logging.basicConfig(
+    level=logging.INFO,
+    format="%(asctime)s [%(levelname)s] %(message)s",
+    handlers=[logging.StreamHandler(sys.stderr)],
+)
+logger = logging.getLogger(__name__)
+
+
+def extract_barcode(seq, pattern, bc_length):
+    """Return barcode downstream of pattern, or None."""
+    idx = seq.find(pattern)
+    if idx == -1:
+        return None
+    start = idx + len(pattern)
+    if start + bc_length > len(seq):
+        return None
+    return seq[start : start + bc_length]
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Extract barcodes from aligned long-read BAM")
+    parser.add_argument("--bam", required=True, help="Aligned BAM from pbmm2")
+    parser.add_argument("--pattern", required=True, help="Fixed sequence upstream of barcode")
+    parser.add_argument("--bc-length", type=int, required=True, help="Barcode length")
+    parser.add_argument("--output", required=True, help="Output TSV")
+    args = parser.parse_args()
+
+    stats = {"total": 0, "mapped": 0, "unmapped": 0, "found": 0, "missing": 0}
+
+    out_fh = sys.stdout if args.output == "/dev/stdout" else open(args.output, "w")
+    with pysam.AlignmentFile(args.bam, "rb") as bam, out_fh as out:
+        for read in bam.fetch():
+            stats["total"] += 1
+
+            if read.is_unmapped:
+                stats["unmapped"] += 1
+                continue
+
+            stats["mapped"] += 1
+            seq = read.query_sequence
+            if seq is None:
+                continue
+
+            bc = extract_barcode(seq, args.pattern, args.bc_length)
+            if bc:
+                stats["found"] += 1
+                cigar = read.cigarstring or "NA"
+                mapq = read.mapping_quality
+                out.write(f"{bc}\t{read.reference_name}\t{cigar};{mapq}\n")
+            else:
+                stats["missing"] += 1
+
+    t, m = stats["total"], stats["mapped"]
+    logger.info("Total reads: %s", f"{t:,}")
+    logger.info("  Mapped: %s (%.1f%%)", f"{m:,}", m / t * 100 if t else 0)
+    logger.info("  Pattern found: %s (%.1f%%)", f"{stats['found']:,}", stats["found"] / m * 100 if m else 0)
+    logger.info("  Pattern missing: %s", f"{stats['missing']:,}")
+
+
+if __name__ == "__main__":
+    main()