crispr_analysis_utils.igvf_export¶

crispr_analysis_utils.igvf_export ¶

Write the guide design in the IGVF per-guide metadata format.

The three steps before this one between them know everything the format asks: alignment gives each guide its site and its PAM, annotation gives the genes, and element mapping gives the element and the guide's role.

What a guide targets is decided from the genome, not from the design's own label for the element. The exception is a validated enhancer, since that an element was shown to regulate a gene is experimental knowledge no annotation can supply.

COLUMNS `module-attribute` ¶

COLUMNS = {'guide_id': 'Unique name of the guide.', 'spacer': "The targeting sequence of the gRNA, read 5' to 3', in DNA bases.", 'oligo': "The sequence the library ordered, without the PAM: the spacer plus any base 5' of it, such as a G added for U6 transcription.", 'n_oligos': 'How many library entries share this spacer. Above one, they are one feature in the counts and the reads cannot be told apart, so the counts are a sum.', 'targeting': 'TRUE when the guide targets a locus relevant to the experiment. FALSE for a safe-targeting or non-targeting control.', 'type': "The guide's role: targeting, safe-targeting, non-targeting or positive control.", 'guide_chr': "Chromosome of the guide's target sequence.", 'guide_start': "Start of the guide's target sequence, 0-based.", 'guide_end': 'End of it, half-open, and without the PAM.', 'strand': 'Strand of the genomic sequence being targeted.', 'pam': "The protospacer adjacent motif as the genome holds it, 5' to 3' on the targeted strand, immediately 3' of the spacer.", 'element_location': 'genic when the element overlaps a protein-coding gene, intergenic when it does not. Independent of genomic_element, which tests promoter and enhancer before it ever looks at the gene body.', 'genomic_element': 'What the targeted element is: promoter, enhancer, distal element, exon or intron. Decided from the genome.', 'intended_target_name': "Names the target: the gene's ENSEMBL id for a promoter, the exon's ENSEMBL id for an exon, and the coordinates otherwise.", 'intended_target_chr': 'Chromosome of the target element.', 'intended_target_start': 'Start of the target element, 0-based.', 'intended_target_end': 'End of it, half-open.', 'putative_target_genes': "Every gene in ENSEMBL associated with the element, comma-separated. A measured pair can name several, as a bidirectional promoter does; an inferred one names the nearest TSS's gene alone.", 'reporter': 'Sequence of the reporter element. Unused here.', 'imperfect': 'Sequence of the matching perfect guide. Unused here.'}

The columns of the IGVF per-guide metadata table, in order.

ELEMENTS `module-attribute` ¶

ELEMENTS = (PROMOTER, ENHANCER, DISTAL, EXON, INTRON)

The genomic_element values this writes. The format allows more.

TYPES `module-attribute` ¶

TYPES = (TARGETING, SAFE, NON, POSITIVE)

The type values this writes.

PROMOTER_WINDOW `module-attribute` ¶

PROMOTER_WINDOW = 2000

Within this of a TSS, an element is a promoter.

The format defines distal element as lying more than 2 kb from any annotated TSS, so the same threshold decides the other side.

VALIDATED `module-attribute` ¶

VALIDATED = 'e2g_validated'

The one element category taken from the design rather than the genome.

POSITIVE_CATEGORIES `module-attribute` ¶

POSITIVE_CATEGORIES = ('promoter', VALIDATED)

The design's element kinds meant as positive controls.

ExportError ¶

Bases: RuntimeError

The inputs cannot produce a valid table.

Source code in src/crispr_analysis_utils/igvf_export.py

class ExportError(RuntimeError):
    """The inputs cannot produce a valid table."""

RunResult `dataclass` ¶

What run produced.

Attributes:

outdir (Path) –

The output folder.
table (DataFrame) –

The IGVF table, one row per guide.
elements (DataFrame) –

How each element was classified.
metadata (dict) –

The record written under run_metadata/.

Source code in src/crispr_analysis_utils/igvf_export.py

@dataclass
class RunResult:
    """What `run` produced.

    Attributes
    ----------
    outdir
        The output folder.
    table
        The IGVF table, one row per guide.
    elements
        How each element was classified.
    metadata
        The record written under ``run_metadata/``.
    """

    outdir: Path
    table: pd.DataFrame
    elements: pd.DataFrame
    metadata: dict = field(default_factory=dict)

classify_elements ¶

classify_elements(elements: DataFrame, store, promoters: DataFrame, *, promoter_window: int = PROMOTER_WINDOW, genes: Mapping[str, list[str]] | None = None) -> pd.DataFrame

Decide what each element is, and what it targets.

Parameters:

elements (DataFrame) –

The elements, as element_mapping.read_elements returns them.
store –

The annotation store, as guide_annotation.load_store returns it.
promoters (DataFrame) –

Transcription start sites, as guide_annotation.read_promoters returns them.
promoter_window (int, default: PROMOTER_WINDOW ) –

Within this many bases of a TSS, an element is a promoter.
genes (Mapping[str, list[str]] | None, default: None ) –

Element id to the ENSEMBL ids known to be regulated by it, for the elements where that was measured rather than inferred.

Returns:

DataFrame –

One row per element: genomic_element, intended_target_name, the target's coordinates and putative_target_genes.

Notes

Everything but VALIDATED is decided from the genome. The design's own label for an element is not evidence about what the element is.

Source code in src/crispr_analysis_utils/igvf_export.py

def classify_elements(
    elements: pd.DataFrame,
    store,
    promoters: pd.DataFrame,
    *,
    promoter_window: int = PROMOTER_WINDOW,
    genes: Mapping[str, list[str]] | None = None,
) -> pd.DataFrame:
    """Decide what each element is, and what it targets.

    Parameters
    ----------
    elements
        The elements, as `element_mapping.read_elements` returns them.
    store
        The annotation store, as `guide_annotation.load_store` returns it.
    promoters
        Transcription start sites, as `guide_annotation.read_promoters`
        returns them.
    promoter_window
        Within this many bases of a TSS, an element is a promoter.
    genes
        Element id to the ENSEMBL ids known to be regulated by it, for the
        elements where that was measured rather than inferred.

    Returns
    -------
    pandas.DataFrame
        One row per element: ``genomic_element``, ``intended_target_name``,
        the target's coordinates and ``putative_target_genes``.

    Notes
    -----
    Everything but `VALIDATED` is decided from the genome. The design's own
    label for an element is not evidence about what the element is.
    """
    import pyranges1 as pr

    if elements.empty:
        raise ExportError("There are no elements to classify.")

    frame = elements.drop_duplicates("element_id").reset_index(drop=True)
    ranges = pr.PyRanges(frame[["Chromosome", "Start", "End", "element_id"]])

    in_gene = _overlapping(ranges, store.rows("gene"), pr)
    in_exon = _overlapping(ranges, store.rows("exon"), pr)
    near, nearest_gene = _nearest_tss(ranges, promoters, pr)

    frame["genomic_element"] = [
        _element_kind(row, in_gene, in_exon, near, promoter_window)
        for row in frame.itertuples()
    ]
    frame["nearest_ensg"] = frame["element_id"].map(nearest_gene)
    # Whether the element itself sits in a gene. genomic_element cannot say:
    # promoter and enhancer win on precedence before the gene body is tested,
    # so this is computed from the same overlap and kept beside it.
    frame["element_location"] = np.where(
        frame["element_id"].isin(in_gene), GENIC, INTERGENIC
    )

    exon_ids = _exon_ids(ranges, store, pr)
    introns = _intron_spans(ranges, store, pr)

    known = {} if genes is None else dict(genes)
    frame["putative_target_genes"] = [
        ",".join(_genes_for(row, known)) for row in frame.itertuples()
    ]
    frame["intended_target_name"] = [
        _target_name(row, known, exon_ids, introns) for row in frame.itertuples()
    ]
    frame["intended_target_chr"] = frame["Chromosome"]
    frame["intended_target_start"] = frame["Start"]
    frame["intended_target_end"] = frame["End"]
    return frame.drop(columns=["nearest_ensg"])

build_table ¶

build_table(alignment: DataFrame, mapping: DataFrame, elements: DataFrame, *, positive_categories: tuple[str, ...] = POSITIVE_CATEGORIES) -> pd.DataFrame

Assemble the IGVF table, one row per guide.

Parameters:

alignment (DataFrame) –

guide_genome_annotation.tsv.
mapping (DataFrame) –

guide_element_mapping.tsv.
elements (DataFrame) –

What classify_elements returned.
positive_categories (tuple[str, ...], default: POSITIVE_CATEGORIES ) –

The design's element categories meant as positive controls. A guide is only called one when the genome agrees the element is a promoter, or the element is a validated enhancer.

Returns:

DataFrame –

The columns of COLUMNS, one row per guide, in library order.

Raises:

ExportError –

A guide has more than one spacer, which the format forbids.

Source code in src/crispr_analysis_utils/igvf_export.py

def build_table(
    alignment: pd.DataFrame,
    mapping: pd.DataFrame,
    elements: pd.DataFrame,
    *,
    positive_categories: tuple[str, ...] = POSITIVE_CATEGORIES,
) -> pd.DataFrame:
    """Assemble the IGVF table, one row per guide.

    Parameters
    ----------
    alignment
        ``guide_genome_annotation.tsv``.
    mapping
        ``guide_element_mapping.tsv``.
    elements
        What `classify_elements` returned.
    positive_categories
        The design's element categories meant as positive controls. A guide
        is only called one when the genome agrees the element is a promoter,
        or the element is a validated enhancer.

    Returns
    -------
    pandas.DataFrame
        The columns of `COLUMNS`, one row per guide, in library order.

    Raises
    ------
    ExportError
        A guide has more than one spacer, which the format forbids.
    """
    guides = alignment.drop_duplicates("guide_id").set_index("guide_id")
    if guides["alias"].duplicated().any():
        raise ExportError("Two guides share an alias, so guide_id would not be unique.")

    rows = _chosen_rows(mapping, elements)
    missing = guides.index.difference(rows.index)
    if len(missing):
        raise ExportError(
            f"{len(missing)} guides are in the alignment but not in the element"
            f" mapping, starting with {list(missing[:3])}. The two tables"
            " describe different libraries."
        )
    chosen = rows["element_id"]
    info = elements.set_index("element_id")

    out = pd.DataFrame(index=guides.index)
    out["guide_id"] = guides["alias"]
    out["spacer"] = guides["spacer"]
    # Carried so the submitted table says which guides share a feature: their
    # counts are a sum and are not comparable with a guide ordered once.
    for column in ("oligo", "n_oligos"):
        # An alignment table written before these existed has neither, and
        # the column is then written empty rather than the table refused.
        out[column] = guides[column] if column in guides.columns else pd.NA

    element = chosen.reindex(out.index)
    kind = element.map(info["genomic_element"])
    category = (
        element.map(info.set_index(info.index)["element_category"])
        if "element_category" in info
        else pd.Series(dtype=object)
    )
    klass = rows["target_class"].reindex(out.index)

    positive = category.isin(positive_categories) & (
        (category == VALIDATED) | (kind == PROMOTER)
    )
    out["targeting"] = (klass == "targeting").map({True: "TRUE", False: "FALSE"})
    out["type"] = np.select(
        [klass == "safe_targeting", klass == "non_targeting", positive],
        [SAFE, NON, POSITIVE],
        default=TARGETING,
    )

    mapped = klass != "non_targeting"
    for column, source in (
        ("guide_chr", "chr"),
        ("guide_start", "start"),
        ("guide_end", "end"),
        ("strand", "strand"),
    ):
        out[column] = rows[source].reindex(out.index).where(mapped)
    out["pam"] = (
        guides["mapped_sequence"]
        .str[-3:]
        .where(mapped & guides["mapped_sequence"].notna())
    )

    for column, source in (
        ("element_location", "element_location"),
        ("genomic_element", "genomic_element"),
        ("intended_target_name", "intended_target_name"),
        ("intended_target_chr", "intended_target_chr"),
        ("intended_target_start", "intended_target_start"),
        ("intended_target_end", "intended_target_end"),
        ("putative_target_genes", "putative_target_genes"),
    ):
        out[column] = element.map(info[source])
    out["reporter"] = ""
    out["imperfect"] = ""

    for column in (
        "guide_start",
        "guide_end",
        "intended_target_start",
        "intended_target_end",
    ):
        out[column] = out[column].astype("Float64").astype("Int64")
    return out.reset_index(drop=True)[list(COLUMNS)]

write_table ¶

write_table(path: str | PathLike[str], table: DataFrame) -> None

Write the IGVF table, with empty cells where a field does not apply.

Source code in src/crispr_analysis_utils/igvf_export.py

def write_table(path: str | os.PathLike[str], table: pd.DataFrame) -> None:
    """Write the IGVF table, with empty cells where a field does not apply."""
    table.reindex(columns=list(COLUMNS)).to_csv(path, sep="\t", index=False, na_rep="")

write_columns_yaml ¶

write_columns_yaml(path: str | PathLike[str]) -> None

Write guide_metadata_columns.yaml, generated from COLUMNS so it cannot drift.

Source code in src/crispr_analysis_utils/igvf_export.py

def write_columns_yaml(path: str | os.PathLike[str]) -> None:
    """Write guide_metadata_columns.yaml, generated from `COLUMNS` so it cannot drift."""
    lines = ["# What each column of this folder's output holds.", "", f"{TABLE}:"]
    lines += [f'  {name}: "{" ".join(text.split())}"' for name, text in COLUMNS.items()]
    Path(path).write_text("\n".join(lines) + "\n", encoding="utf-8")

counts ¶

counts(table: DataFrame, elements: DataFrame) -> dict

What the table holds, by type and by element class.

Source code in src/crispr_analysis_utils/igvf_export.py

def counts(table: pd.DataFrame, elements: pd.DataFrame) -> dict:
    """What the table holds, by type and by element class."""
    return {
        "guides": int(len(table)),
        "type": {name: int((table["type"] == name).sum()) for name in TYPES},
        "genomic_element": {
            name: int((table["genomic_element"] == name).sum()) for name in ELEMENTS
        },
        "elements": int(len(elements)),
        "with_a_gene": int((elements["putative_target_genes"] != "").sum()),
        "missing_pam": int(((table["type"] != NON) & table["pam"].isna()).sum()),
    }

run ¶

run(alignment: str | PathLike[str], mapping: str | PathLike[str], elements: str | PathLike[str], store: str | PathLike[str], outdir: str | PathLike[str], *, promoters: str | PathLike[str] | None = None, genes: str | PathLike[str] | None = None, promoter_window: int = PROMOTER_WINDOW) -> RunResult

Write the guide design in the IGVF per-guide metadata format.

Parameters:

alignment (str | PathLike[str]) –

guide_genome_annotation.tsv.
mapping (str | PathLike[str]) –

guide_element_mapping.tsv.
elements (str | PathLike[str]) –

The elements BED the mapping used.
store (str | PathLike[str]) –

The annotation store, which supplies the genes, exons and exon ids.
outdir (str | PathLike[str]) –

Where the table goes. Nothing is written outside it.
promoters (str | PathLike[str] | None, default: None ) –

A TSS BED. Without one, the store's reference transcript starts.
genes (str | PathLike[str] | None, default: None ) –

A TSV of measured element-to-gene pairs, with an element id column and an ENSEMBL id column. Used for the elements whose gene was measured rather than inferred.
promoter_window (int, default: PROMOTER_WINDOW ) –

Within this many bases of a TSS, an element is a promoter.

Returns:

RunResult –

The table, the element classification and the record.

Source code in src/crispr_analysis_utils/igvf_export.py

def run(
    alignment: str | os.PathLike[str],
    mapping: str | os.PathLike[str],
    elements: str | os.PathLike[str],
    store: str | os.PathLike[str],
    outdir: str | os.PathLike[str],
    *,
    promoters: str | os.PathLike[str] | None = None,
    genes: str | os.PathLike[str] | None = None,
    promoter_window: int = PROMOTER_WINDOW,
) -> RunResult:
    """Write the guide design in the IGVF per-guide metadata format.

    Parameters
    ----------
    alignment
        ``guide_genome_annotation.tsv``.
    mapping
        ``guide_element_mapping.tsv``.
    elements
        The elements BED the mapping used.
    store
        The annotation store, which supplies the genes, exons and exon ids.
    outdir
        Where the table goes. Nothing is written outside it.
    promoters
        A TSS BED. Without one, the store's reference transcript starts.
    genes
        A TSV of measured element-to-gene pairs, with an element id column
        and an ENSEMBL id column. Used for the elements whose gene was
        measured rather than inferred.
    promoter_window
        Within this many bases of a TSS, an element is a promoter.

    Returns
    -------
    RunResult
        The table, the element classification and the record.
    """
    from .element_mapping import read_elements
    from .guide_annotation import load_store, promoters_from_store, read_promoters

    outdir = Path(outdir)
    for name in (RESULTS_DIR, METADATA_DIR, "logs"):
        (outdir / name).mkdir(parents=True, exist_ok=True)
    record = outdir / METADATA_DIR / f"{STEP}_{_stamp()}.json"

    state: dict = {
        "step": STEP,
        "complete": False,
        "error": None,
        "started": _now(),
        "cau_version": __version__,
        "inputs": {
            "alignment": str(alignment),
            "mapping": str(mapping),
            "elements": str(elements),
            "store": str(store),
            "genes": None if genes is None else str(genes),
        },
        "parameters": {"promoter_window": promoter_window},
        "counts": {},
    }
    handler = logging.FileHandler(outdir / "logs" / LOG_FILE, encoding="utf-8")
    handler.setFormatter(logging.Formatter(_LOG_FORMAT))
    package = logging.getLogger(__name__.rsplit(".", 1)[0])
    package.addHandler(handler)
    if package.level == logging.NOTSET:
        package.setLevel(logging.INFO)
    propagate = package.propagate
    package.propagate = False
    try:
        loaded = load_store(store)
        tss = read_promoters(promoters) if promoters else promoters_from_store(loaded)
        regions = read_elements(elements)
        align = pd.read_csv(alignment, sep="\t", dtype={"chr": str}, low_memory=False)
        maps = pd.read_csv(mapping, sep="\t", dtype={"chr": str}, low_memory=False)
        state["inputs"]["alignment_md5"] = file_md5(alignment)
        state["inputs"]["elements_md5"] = file_md5(elements)

        info = classify_elements(
            regions,
            loaded,
            tss,
            promoter_window=promoter_window,
            genes=_read_genes(genes),
        )
        table = build_table(align, maps, info)

        write_table(outdir / RESULTS_DIR / TABLE, table)
        write_columns_yaml(outdir / RESULTS_DIR / "guide_metadata_columns.yaml")
        state["counts"] = counts(table, info)
        state["complete"] = True
        state["finished"] = _now()
        write_json(record, state)
        logger.info("IGVF export complete: %s", state["counts"]["type"])
        return RunResult(outdir=outdir, table=table, elements=info, metadata=state)
    except Exception as error:
        state["error"] = f"{type(error).__name__}: {error}"
        state["finished"] = _now()
        write_json(record, state)
        logger.error("IGVF export failed: %s", state["error"])
        raise
    finally:
        package.removeHandler(handler)
        package.propagate = propagate
        handler.close()

crispr_analysis_utils.igvf_export¶