feat(pipeline): organize step outputs and emit protospacer BED spans

Author: emattei

docs/python/tutorial-guide-alignment-qc.md

@@ -47,17 +47,29 @@ python scripts/run_guide_alignment_qc.py \
 
 In `results/guide_alignment_qc/`:
 
-- `guides_input.fastq`
+- `logs/pipeline.log`
+
+In `results/guide_alignment_qc/gem_index/`:
+
 - `genome_index.gem` (+ index sidecar files)
+- `gem_index.log`
+- `index_command.sh`
+- `index_inputs.txt`
+
+In `results/guide_alignment_qc/alignment_outputs/`:
+
+- `guides_input.fastq`
 - `guides_mapped.sam`
 - `guides_mapped.log`
-- `guides_valid_unique.sam`
-- `guides_valid_multi.sam`
-- `valid_alignments.bed`
+
+In `results/guide_alignment_qc/guide_alignments_outputs/`:
+
+- `valid_alignments.bed` (protospacer coordinates; excludes PAM)
 - `discarded_alignments.tsv`
 - `unmapped.tsv`
 - `guide_alignment_log.tsv`
 - `invalid_alignments.tsv`
+- `alignment_summary.tsv`
 
 ## 6. Run on SLURM (one Python script)
 

scripts/run_guide_alignment_qc.py

@@ -14,6 +14,7 @@
 from datetime import datetime
 import logging
 from pathlib import Path
+import shlex
 import time
 
 import crispr_analysis_utils as cau
@@ -46,6 +47,13 @@ def main() -> None:
 
     outdir = Path(args.outdir)
     outdir.mkdir(parents=True, exist_ok=True)
+
+    gem_index_dir = outdir / "gem_index"
+    alignment_dir = outdir / "alignment_outputs"
+    filtering_dir = outdir / "guide_alignments_outputs"
+    for p in (gem_index_dir, alignment_dir, filtering_dir):
+        p.mkdir(parents=True, exist_ok=True)
+
     log_dir = outdir / "logs"
     log_dir.mkdir(parents=True, exist_ok=True)
     log_path = log_dir / "pipeline.log"
@@ -65,12 +73,10 @@ def main() -> None:
     logging.info("Output directory: %s", outdir)
     logging.info("Log file: %s", log_path)
 
-    guides_fastq = outdir / "guides_input.fastq"
-    gem_index_prefix = outdir / "genome_index"
-    gem_index = outdir / "genome_index.gem"
-    mapped_sam = outdir / "guides_mapped.sam"
-    valid_unique_sam = outdir / "guides_valid_unique.sam"
-    valid_multi_sam = outdir / "guides_valid_multi.sam"
+    guides_fastq = alignment_dir / "guides_input.fastq"
+    gem_index_prefix = gem_index_dir / "genome_index"
+    gem_index = gem_index_dir / "genome_index.gem"
+    mapped_sam = alignment_dir / "guides_mapped.sam"
 
     logging.info("Step 1/4: Build guide FASTQ")
     cau.guide_qc.guides_to_fastq(
@@ -81,10 +87,26 @@ def main() -> None:
     )
 
     logging.info("Step 2/4: Build GEM index")
+    index_cmd_parts = [
+        "gem-indexer",
+        "-i",
+        str(args.reference_fasta),
+        "-o",
+        str(gem_index_prefix),
+        "-t",
+        str(args.threads),
+    ]
+    index_cmd = " ".join(shlex.quote(x) for x in index_cmd_parts)
+    (gem_index_dir / "index_command.sh").write_text(index_cmd + "\n", encoding="utf-8")
+    (gem_index_dir / "index_inputs.txt").write_text(
+        f"reference_fasta={Path(args.reference_fasta).resolve()}\n",
+        encoding="utf-8",
+    )
     cau.gem_mapper.build_gem_index(
         args.reference_fasta,
         gem_index_prefix,
         threads=args.threads,
+        log_path=gem_index_dir / "gem_index.log",
     )
 
     logging.info("Step 3/4: Map guides with GEM")
@@ -95,16 +117,23 @@ def main() -> None:
         threads=args.threads,
         mapping_mode=args.mapping_mode,
         sam_compact=False,
+        log_path=alignment_dir / "guides_mapped.log",
     )
 
     logging.info("Step 4/4: Filter alignments and compute QC outputs")
     summary = cau.guide_qc.filter_guide_alignments(
         mapped_sam,
-        valid_unique_sam,
-        valid_multi_sam,
+        None,
+        None,
         pam=args.pam,
         chromsizes=args.chromsizes,
         allow_leading_g_softclip=args.allow_leading_g_softclip,
+        output_valid_bed=filtering_dir / "valid_alignments.bed",
+        output_discarded_tsv=filtering_dir / "discarded_alignments.tsv",
+        output_unmapped_tsv=filtering_dir / "unmapped.tsv",
+        output_guide_log_tsv=filtering_dir / "guide_alignment_log.tsv",
+        output_invalid_tsv=filtering_dir / "invalid_alignments.tsv",
+        output_summary_tsv=filtering_dir / "alignment_summary.tsv",
     )
     finished = time.time()
     finish_iso = datetime.now().isoformat(timespec="seconds")

src/crispr_analysis_utils/gem_mapper.py

@@ -14,6 +14,7 @@ def build_gem_index(
     *,
     threads: int | None = None,
     gem_indexer_bin: str = "gem-indexer",
+    log_path: str | Path | None = None,
 ) -> str:
     """Build a GEM index from a reference FASTA."""
     cmd = [
@@ -25,7 +26,10 @@ def build_gem_index(
     ]
     if threads is not None:
         cmd.extend(["-t", str(threads)])
-    return run_shell_cmd(" ".join(shlex.quote(token) for token in cmd))
+    cmd_str = " ".join(shlex.quote(token) for token in cmd)
+    if log_path is not None:
+        cmd_str = f"{cmd_str} > {shlex.quote(str(log_path))} 2>&1"
+    return run_shell_cmd(cmd_str)
 
 
 def map_guides_with_gem(

src/crispr_analysis_utils/guide_qc.py

@@ -183,8 +183,8 @@ def _evaluate_alignment_layout(
 
 def filter_guide_alignments(
     input_sam: str | Path,
-    output_unique_sam: str | Path,
-    output_multi_sam: str | Path,
+    output_unique_sam: str | Path | None = None,
+    output_multi_sam: str | Path | None = None,
     *,
     output_invalid_tsv: str | Path = "auto",
     output_valid_bed: str | Path = "auto",
@@ -270,14 +270,20 @@ def filter_guide_alignments(
             if is_valid:
                 valid_by_guide.setdefault(guide_id, []).append(aln)
                 guide_stats[guide_id]["n_valid"] += 1
-                if aln.reference_end is not None:
-                    nm_tag = int(aln.get_tag("NM")) if aln.has_tag("NM") else -1
-                    as_tag = int(aln.get_tag("AS")) if aln.has_tag("AS") else -1
+                nm_tag = int(aln.get_tag("NM")) if aln.has_tag("NM") else -1
+                as_tag = int(aln.get_tag("AS")) if aln.has_tag("AS") else -1
+                bed_span = _protospacer_bed_span(
+                    aln,
+                    query_len=query_len,
+                    pam_len=len(pam),
+                    allow_leading_g_softclip=allow_leading_g_softclip,
+                )
+                if bed_span is not None:
                     valid_bed_rows.append(
                         (
                             contig,
-                            int(aln.reference_start),
-                            int(aln.reference_end),
+                            bed_span[0],
+                            bed_span[1],
                             guide_id,
                             int(aln.mapping_quality),
                             "-" if aln.is_reverse else "+",
@@ -312,22 +318,31 @@ def filter_guide_alignments(
                     )
                 )
 
-    unique_path = Path(output_unique_sam)
-    multi_path = Path(output_multi_sam)
-    unique_path.parent.mkdir(parents=True, exist_ok=True)
-    multi_path.parent.mkdir(parents=True, exist_ok=True)
-
-    with (
-        pysam.AlignmentFile(str(unique_path), "w", header=header) as unique_sam,
-        pysam.AlignmentFile(str(multi_path), "w", header=header) as multi_sam,
-    ):
-        for alignments in valid_by_guide.values():
-            target = unique_sam if len(alignments) == 1 else multi_sam
-            for aln in alignments:
-                target.write(aln)
+    base_output_dir = Path(".")
+    if output_unique_sam is not None:
+        unique_path = Path(output_unique_sam)
+        base_output_dir = unique_path.parent
+    elif output_multi_sam is not None:
+        multi_path = Path(output_multi_sam)
+        base_output_dir = multi_path.parent
+
+    if output_unique_sam is not None and output_multi_sam is not None:
+        unique_path = Path(output_unique_sam)
+        multi_path = Path(output_multi_sam)
+        unique_path.parent.mkdir(parents=True, exist_ok=True)
+        multi_path.parent.mkdir(parents=True, exist_ok=True)
+
+        with (
+            pysam.AlignmentFile(str(unique_path), "w", header=header) as unique_sam,
+            pysam.AlignmentFile(str(multi_path), "w", header=header) as multi_sam,
+        ):
+            for alignments in valid_by_guide.values():
+                target = unique_sam if len(alignments) == 1 else multi_sam
+                for aln in alignments:
+                    target.write(aln)
 
     if output_invalid_tsv == "auto":
-        invalid_path = unique_path.with_name("invalid_alignments.tsv")
+        invalid_path = base_output_dir / "invalid_alignments.tsv"
     else:
         invalid_path = Path(output_invalid_tsv)
     invalid_path.parent.mkdir(parents=True, exist_ok=True)
@@ -337,7 +352,7 @@ def filter_guide_alignments(
             handle.write("\t".join(row) + "\n")
 
     if output_valid_bed == "auto":
-        valid_bed_path = unique_path.with_name("valid_alignments.bed")
+        valid_bed_path = base_output_dir / "valid_alignments.bed"
     else:
         valid_bed_path = Path(output_valid_bed)
     valid_bed_path.parent.mkdir(parents=True, exist_ok=True)
@@ -346,7 +361,7 @@ def filter_guide_alignments(
             handle.write("\t".join(map(str, row)) + "\n")
 
     if output_discarded_tsv == "auto":
-        discarded_path = unique_path.with_name("discarded_alignments.tsv")
+        discarded_path = base_output_dir / "discarded_alignments.tsv"
     else:
         discarded_path = Path(output_discarded_tsv)
     discarded_path.parent.mkdir(parents=True, exist_ok=True)
@@ -358,7 +373,7 @@ def filter_guide_alignments(
             handle.write("\t".join(map(str, row)) + "\n")
 
     if output_unmapped_tsv == "auto":
-        unmapped_path = unique_path.with_name("unmapped.tsv")
+        unmapped_path = base_output_dir / "unmapped.tsv"
     else:
         unmapped_path = Path(output_unmapped_tsv)
     unmapped_path.parent.mkdir(parents=True, exist_ok=True)
@@ -368,7 +383,7 @@ def filter_guide_alignments(
             handle.write("\t".join(map(str, row)) + "\n")
 
     if output_guide_log_tsv == "auto":
-        guide_log_path = unique_path.with_name("guide_alignment_log.tsv")
+        guide_log_path = base_output_dir / "guide_alignment_log.tsv"
     else:
         guide_log_path = Path(output_guide_log_tsv)
     guide_log_path.parent.mkdir(parents=True, exist_ok=True)
@@ -402,7 +417,7 @@ def filter_guide_alignments(
             n_guides_multi_valid += 1
 
     if output_summary_tsv == "auto":
-        summary_path = unique_path.with_name("alignment_summary.tsv")
+        summary_path = base_output_dir / "alignment_summary.tsv"
     else:
         summary_path = Path(output_summary_tsv)
     summary_path.parent.mkdir(parents=True, exist_ok=True)
@@ -433,6 +448,40 @@ def filter_guide_alignments(
     }
 
 
+def _protospacer_bed_span(
+    aln,
+    *,
+    query_len: int,
+    pam_len: int,
+    allow_leading_g_softclip: bool,
+) -> tuple[int, int] | None:
+    """Map protospacer query region (no PAM) to genomic BED span."""
+    if query_len <= pam_len:
+        return None
+
+    query_ops = [None] * query_len
+    qpos = 0
+    for op, length in (aln.cigartuples or []):
+        if op in _CIGAR_MATCH | {_CIGAR_INS, _CIGAR_SOFT}:
+            for _ in range(length):
+                if qpos < query_len:
+                    query_ops[qpos] = op
+                qpos += 1
+
+    protospacer_start = 0
+    if allow_leading_g_softclip and query_ops and query_ops[0] == _CIGAR_SOFT:
+        protospacer_start = 1
+    protospacer_end = query_len - pam_len
+    if protospacer_start >= protospacer_end:
+        return None
+
+    q2r = {q: r for q, r in aln.get_aligned_pairs(matches_only=True)}
+    ref_positions = [q2r[i] for i in range(protospacer_start, protospacer_end) if i in q2r]
+    if not ref_positions:
+        return None
+    return min(ref_positions), max(ref_positions) + 1
+
+
 def _resolve_allowed_contigs(
     *,
     primary_contigs: set[str] | list[str] | tuple[str, ...] | None,