crispr_analysis_utils.guide_qc.filter_guide_alignments

crispr_analysis_utils.guide_qc.filter_guide_alignments

filter_guide_alignments(input_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', output_discarded_tsv: str | Path = 'auto', output_unmapped_tsv: str | Path = 'auto', output_guide_log_tsv: str | Path = 'auto', output_summary_tsv: str | Path = 'auto', alias_by_guide_id: dict[str, str] | None = None, pam: str = 'NGG', allow_leading_g_softclip: bool = True, primary_contigs: set[str] | list[str] | tuple[str, ...] | None = None, chromsizes: str | Path | DataFrame | None = None) -> dict[str, int]

Filter SAM/BAM guide alignments and split valid unique vs multi-mapping.

Rules: - must map to a primary contig - PAM must be aligned; mismatches/indels/soft clips not allowed in PAM except PAM's first base (N) - spacer cannot contain insertions or deletions - leading G soft-clipping is allowed by default

Source code in src/crispr_analysis_utils/guide_qc.py

def filter_guide_alignments(
    input_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",
    output_discarded_tsv: str | Path = "auto",
    output_unmapped_tsv: str | Path = "auto",
    output_guide_log_tsv: str | Path = "auto",
    output_summary_tsv: str | Path = "auto",
    alias_by_guide_id: dict[str, str] | None = None,
    pam: str = "NGG",
    allow_leading_g_softclip: bool = True,
    primary_contigs: set[str] | list[str] | tuple[str, ...] | None = None,
    chromsizes: str | Path | pd.DataFrame | None = None,
) -> dict[str, int]:
    """Filter SAM/BAM guide alignments and split valid unique vs multi-mapping.

    Rules:
    - must map to a primary contig
    - PAM must be aligned; mismatches/indels/soft clips not allowed in PAM
      except PAM's first base (`N`)
    - spacer cannot contain insertions or deletions
    - leading G soft-clipping is allowed by default
    """
    try:
        import pysam
    except ImportError as exc:
        raise ImportError("pysam is required for SAM/BAM filtering.") from exc

    primary = _resolve_allowed_contigs(primary_contigs=primary_contigs, chromsizes=chromsizes)
    input_path = Path(input_sam)
    mode = "rb" if input_path.suffix == ".bam" else "r"

    valid_by_guide: dict[str, list] = {}
    invalid_rows: list[tuple[str, str, str]] = []
    discarded_rows: list[tuple[str, str, int, int, int, str, str, int, int, str, str]] = []
    unmapped_rows: list[tuple[str, int]] = []
    valid_bed_rows: list[tuple[str, int, int, str, int, str, int, int, str]] = []
    guide_stats = defaultdict(
        lambda: {"n_aligned": 0, "n_valid": 0, "n_discarded": 0, "n_not_mapped": 0}
    )
    qname_to_sequence: dict[str, str] = {}

    with pysam.AlignmentFile(str(input_path), mode) as in_sam:
        header = in_sam.header
        for aln in in_sam.fetch(until_eof=True):
            qname = aln.query_name
            seq = aln.query_sequence
            if seq is not None and seq != "*":
                seq_u = seq.upper()
                qname_to_sequence[qname] = _revcomp(seq_u) if aln.is_reverse else seq_u
            recovered_seq = qname_to_sequence.get(qname)
            guide_id = recovered_seq if recovered_seq is not None else qname
            read_name = qname

            if aln.is_unmapped:
                invalid_rows.append((read_name, ".", "unmapped"))
                unmapped_rows.append((read_name, aln.flag))
                guide_stats[read_name]["n_not_mapped"] = 1
                continue

            contig = aln.reference_name
            guide_stats[read_name]["n_aligned"] += 1
            if primary is not None and contig not in primary:
                invalid_rows.append((read_name, contig, "non_primary_contig"))
                guide_stats[read_name]["n_discarded"] += 1
                discarded_rows.append(
                    (
                        read_name,
                        contig,
                        int(aln.reference_start) + 1,
                        int(aln.flag),
                        int(aln.mapping_quality),
                        "-" if aln.is_reverse else "+",
                        aln.cigarstring or "*",
                        int(aln.get_tag("NM")) if aln.has_tag("NM") else -1,
                        int(aln.get_tag("AS")) if aln.has_tag("AS") else -1,
                        aln.get_tag("MD") if aln.has_tag("MD") else "",
                        "non_primary_contig",
                    )
                )
                continue

            md_tag = aln.get_tag("MD") if aln.has_tag("MD") else ""
            cigartuples = aln.cigartuples or []
            query_len = aln.query_length or _query_length_from_cigartuples(cigartuples)
            mismatch_positions = _mismatch_positions_from_md(md_tag, cigartuples)
            is_valid, reason = _evaluate_alignment_layout(
                cigartuples,
                query_len,
                mismatch_positions,
                pam=pam,
                is_reverse_strand=aln.is_reverse,
                allow_leading_g_softclip=allow_leading_g_softclip,
            )
            if is_valid:
                valid_by_guide.setdefault(read_name, []).append(aln)
                guide_stats[read_name]["n_valid"] += 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
                protospacer_bounds = _protospacer_query_bounds(
                    aln,
                    query_len=query_len,
                    pam_len=len(pam),
                    allow_leading_g_softclip=allow_leading_g_softclip,
                )
                if protospacer_bounds is not None:
                    q_start, q_end = protospacer_bounds
                    bed_span = _protospacer_bed_span(aln, q_start=q_start, q_end=q_end)
                    seq_start = 0
                    if _leading_g_softclipped(
                        aln, allow_leading_g_softclip=allow_leading_g_softclip
                    ):
                        seq_start = 1
                    seq_end = len(guide_id) - len(pam)
                    protospacer_seq = guide_id[seq_start:seq_end] if seq_start < seq_end else None
                else:
                    bed_span = None
                    protospacer_seq = None

                if bed_span is not None and protospacer_seq is not None:
                    valid_bed_rows.append(
                        (
                            contig,
                            bed_span[0],
                            bed_span[1],
                            protospacer_seq,
                            int(aln.mapping_quality),
                            "-" if aln.is_reverse else "+",
                            nm_tag,
                            as_tag,
                            alias_by_guide_id.get(guide_id, read_name)
                            if alias_by_guide_id is not None
                            else read_name,
                        )
                    )
            else:
                if reason in {"pam_not_fully_aligned", "query_too_short_for_pam", "softclip_not_allowed"}:
                    normalized_reason = "discarded_tail_unaligned"
                elif reason == "pam_gg_mismatch":
                    normalized_reason = "discarded_tail_mismatch"
                elif reason in {"spacer_insertion", "spacer_deletion"}:
                    normalized_reason = "discarded_protospacer_indel"
                else:
                    normalized_reason = reason
                invalid_rows.append((read_name, contig, reason))
                guide_stats[read_name]["n_discarded"] += 1
                discarded_rows.append(
                    (
                        read_name,
                        contig,
                        int(aln.reference_start) + 1,
                        int(aln.flag),
                        int(aln.mapping_quality),
                        "-" if aln.is_reverse else "+",
                        aln.cigarstring or "*",
                        int(aln.get_tag("NM")) if aln.has_tag("NM") else -1,
                        int(aln.get_tag("AS")) if aln.has_tag("AS") else -1,
                        aln.get_tag("MD") if aln.has_tag("MD") else "",
                        normalized_reason,
                    )
                )

    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 = base_output_dir / "invalid_alignments.tsv"
    else:
        invalid_path = Path(output_invalid_tsv)
    invalid_path.parent.mkdir(parents=True, exist_ok=True)
    with invalid_path.open("w", encoding="utf-8") as handle:
        handle.write("guide_id\tcontig\treason\n")
        for row in invalid_rows:
            handle.write("\t".join(row) + "\n")

    if output_valid_bed == "auto":
        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)
    with valid_bed_path.open("w", encoding="utf-8") as handle:
        for row in valid_bed_rows:
            handle.write("\t".join(map(str, row)) + "\n")

    if output_discarded_tsv == "auto":
        discarded_path = base_output_dir / "discarded_alignments.tsv"
    else:
        discarded_path = Path(output_discarded_tsv)
    discarded_path.parent.mkdir(parents=True, exist_ok=True)
    with discarded_path.open("w", encoding="utf-8") as handle:
        handle.write(
            "read_name\tchromosome\tpos1\tflag\tmapq\tstrand\tcigar\tNM\tAS\tMD\treason\n"
        )
        for row in discarded_rows:
            handle.write("\t".join(map(str, row)) + "\n")

    if output_unmapped_tsv == "auto":
        unmapped_path = base_output_dir / "unmapped.tsv"
    else:
        unmapped_path = Path(output_unmapped_tsv)
    unmapped_path.parent.mkdir(parents=True, exist_ok=True)
    with unmapped_path.open("w", encoding="utf-8") as handle:
        handle.write("read_name\tflag\n")
        for row in unmapped_rows:
            handle.write("\t".join(map(str, row)) + "\n")

    if output_guide_log_tsv == "auto":
        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)
    with guide_log_path.open("w", encoding="utf-8") as handle:
        handle.write("guide_id\tn_aligned\tn_valid\tn_discarded\tn_not_mapped\n")
        for guide_id in sorted(guide_stats):
            s = guide_stats[guide_id]
            handle.write(
                f"{guide_id}\t{s['n_aligned']}\t{s['n_valid']}\t{s['n_discarded']}\t{s['n_not_mapped']}\n"
            )

    n_guides_unique_valid = 0
    n_guides_multi_valid = 0
    n_guides_aligned_none_valid = 0
    n_guides_unmapped = 0
    n_guides_one_valid_plus_invalid = 0

    for s in guide_stats.values():
        if s["n_not_mapped"] == 1 and s["n_aligned"] == 0:
            n_guides_unmapped += 1
            continue
        if s["n_valid"] == 0 and s["n_aligned"] > 0:
            n_guides_aligned_none_valid += 1
            continue
        if s["n_valid"] == 1:
            n_guides_unique_valid += 1
            if s["n_discarded"] > 0:
                n_guides_one_valid_plus_invalid += 1
            continue
        if s["n_valid"] > 1:
            n_guides_multi_valid += 1

    if output_summary_tsv == "auto":
        summary_path = base_output_dir / "alignment_summary.tsv"
    else:
        summary_path = Path(output_summary_tsv)
    summary_path.parent.mkdir(parents=True, exist_ok=True)
    with summary_path.open("w", encoding="utf-8") as handle:
        handle.write("metric\tcount\n")
        handle.write(f"guides_unique_valid\t{n_guides_unique_valid}\n")
        handle.write(f"guides_multi_valid\t{n_guides_multi_valid}\n")
        handle.write(f"guides_aligned_none_valid\t{n_guides_aligned_none_valid}\n")
        handle.write(f"guides_unmapped\t{n_guides_unmapped}\n")
        handle.write(f"guides_one_valid_plus_invalid\t{n_guides_one_valid_plus_invalid}\n")

    n_valid_guides = len(valid_by_guide)
    n_unique_guides = sum(1 for v in valid_by_guide.values() if len(v) == 1)
    n_multi_guides = sum(1 for v in valid_by_guide.values() if len(v) > 1)
    n_valid_alignments = sum(len(v) for v in valid_by_guide.values())

    return {
        "valid_guides": n_valid_guides,
        "unique_guides": n_unique_guides,
        "multi_guides": n_multi_guides,
        "valid_alignments": n_valid_alignments,
        "invalid_alignments": len(invalid_rows),
        "guides_unique_valid": n_guides_unique_valid,
        "guides_multi_valid": n_guides_multi_valid,
        "guides_aligned_none_valid": n_guides_aligned_none_valid,
        "guides_unmapped": n_guides_unmapped,
        "guides_one_valid_plus_invalid": n_guides_one_valid_plus_invalid,
    }