crispr_analysis_utils.guide_annotation.annotate¶

crispr_analysis_utils.guide_annotation.annotate ¶

Place each aligned site in a gene: genic or intergenic, exon or intron.

Overlaps are decided on the protospacer, so a guide reaching into an exon counts as exonic. Distances and the position along a CDS are measured from the target site, which is where the break falls.

Every interval operation ignores strand: a guide on one strand targets a gene on either.

DELIMITER `module-attribute` ¶

DELIMITER = ';'

Separates the ids in a list-valued column.

A semicolon rather than a comma, because a comma already separates ids inside the IGVF format's own putative_target_genes.

FEATURES `module-attribute` ¶

FEATURES = (EXON, INTRON)

Where a site falls, in the order qc.tsv counts them.

SITE_COLUMNS `module-attribute` ¶

SITE_COLUMNS = ('chr', 'start', 'end', 'strand', 'target_site')

What annotate_sites needs of every site.

DISTANCE_FROM `module-attribute` ¶

DISTANCE_FROM = (MIDPOINT, TARGET_SITE)

Where a distance is measured from.

The midpoint of the protospacer by default. A dCas9 occupies the whole protospacer rather than acting at one base, so its centre is the more honest reference for how far a guide sits from a TSS. target_site measures from the boundary 3 bases 5' of the PAM instead, which is where a nuclease cuts.

AnnotationError ¶

Bases: RuntimeError

The sites and the annotation cannot be compared at all.

Source code in src/crispr_analysis_utils/guide_annotation/annotate.py

class AnnotationError(RuntimeError):
    """The sites and the annotation cannot be compared at all."""

annotate_sites ¶

annotate_sites(sites: DataFrame, store: Store, promoters: DataFrame, *, promoter_window: int = 1000, distance_from: str = MIDPOINT) -> pd.DataFrame

Annotate every site against the genes of the store.

Parameters:

sites (DataFrame) –

One row per site, with SITE_COLUMNS at least. Any other column is carried through to the output untouched.
store (Store) –

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

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

How far either side of a TSS counts as its promoter.

Returns:

DataFrame –

One row per site per overlapping gene; a site in no gene keeps one row, intergenic, with the gene columns empty. The input's own columns come first.

Raises:

AnnotationError –

No site contig appears in the store, so nothing could ever overlap.

Source code in src/crispr_analysis_utils/guide_annotation/annotate.py

def annotate_sites(
    sites: pd.DataFrame,
    store: Store,
    promoters: pd.DataFrame,
    *,
    promoter_window: int = 1000,
    distance_from: str = MIDPOINT,
) -> pd.DataFrame:
    """Annotate every site against the genes of the store.

    Parameters
    ----------
    sites
        One row per site, with `SITE_COLUMNS` at least. Any other column is
        carried through to the output untouched.
    store
        The annotation store, as `load_store` returns it.
    promoters
        Transcription start sites, as `read_promoters` returns them.
    promoter_window
        How far either side of a TSS counts as its promoter.

    Returns
    -------
    pandas.DataFrame
        One row per site per overlapping gene; a site in no gene keeps one
        row, `intergenic`, with the gene columns empty. The input's own
        columns come first.

    Raises
    ------
    AnnotationError
        No site contig appears in the store, so nothing could ever overlap.
    """
    import pyranges1 as pr

    missing = set(SITE_COLUMNS) - set(sites.columns)
    if missing:
        raise AnnotationError(f"The sites are missing {', '.join(sorted(missing))}.")
    if promoter_window < 0:
        raise ValueError(
            f"promoter_window must not be negative, got {promoter_window}."
        )
    if distance_from not in DISTANCE_FROM:
        raise ValueError(
            f"distance_from must be one of {', '.join(DISTANCE_FROM)}, got"
            f" {distance_from!r}."
        )

    sites = sites.reset_index(drop=True)
    _check_contigs(sites, store)

    spacers = pr.PyRanges(
        {
            "Chromosome": sites["chr"].to_numpy(),
            "Start": sites["start"].to_numpy(),
            "End": sites["end"].to_numpy(),
            "Strand": sites["strand"].to_numpy(),
            "_row": np.arange(len(sites)),
        }
    )
    # The CDS position asks where along the coding sequence Cas9 acts, so it
    # stays on the target site whatever distances are measured from.
    targets = pr.PyRanges(
        {
            "Chromosome": sites["chr"].to_numpy(),
            "Start": sites["target_site"].to_numpy(),
            "End": sites["target_site"].to_numpy() + 1,
            "Strand": sites["strand"].to_numpy(),
            "_row": np.arange(len(sites)),
        }
    )

    # Where distances are measured from. The midpoint rounds down, so a
    # protospacer of even length takes the base left of centre.
    anchor_at = (
        (sites["start"].to_numpy() + sites["end"].to_numpy()) // 2
        if distance_from == MIDPOINT
        else sites["target_site"].to_numpy()
    )
    anchors = pr.PyRanges(
        {
            "Chromosome": sites["chr"].to_numpy(),
            "Start": anchor_at,
            "End": anchor_at + 1,
            "Strand": sites["strand"].to_numpy(),
            "_row": np.arange(len(sites)),
        }
    )

    pairs = _gene_pairs(spacers, store, pr)
    pairs = pairs.merge(
        _exon_calls(spacers, store, pr), on=["_row", "ensg"], how="left"
    )
    pairs = pairs.merge(
        _cds_position(targets, store, pr), on=["_row", "ensg"], how="left"
    )

    # Two axes, not one. Whether a gene is there at all, and where inside it:
    # a guide can be genic and in a promoter and in an exon at once, which a
    # single column could not say.
    pairs["location"] = np.where(pairs["ensg"].isna(), INTERGENIC, GENIC)
    pairs["feature"] = np.where(
        pairs["ensg"].isna(),
        "",
        np.where(pairs["enst"].notna(), EXON, INTRON),
    )
    pairs["region"] = np.where(
        pairs["feature"] == EXON,
        np.where(pairs["cds_position"].notna(), CDS, UTR),
        "",
    )
    pairs["exonic_any_transcript"] = pairs["exonic_any_transcript"].fillna(False)
    # How far into the gene the target site falls, from the gene's own 5'
    # end. Near zero is the annotated start; a large one is deeper in, which
    # is what distinguishes a canonical promoter from an alternative one.
    five_prime = np.where(
        pairs["gene_strand"] == "+", pairs["gene_start"], pairs["gene_end"] - 1
    )
    toward = anchor_at[pairs["_row"].to_numpy()] - five_prime
    pairs["distance_to_gene_start"] = np.where(
        pairs["gene_strand"].isna(),
        np.nan,
        np.where(pairs["gene_strand"] == "+", toward, -toward),
    )

    pairs["sense"] = np.where(
        pairs["gene_strand"].isna(),
        "",
        np.where(pairs["strand_site"] == pairs["gene_strand"], SENSE, ANTISENSE),
    )

    tss = _tss_columns(anchors, spacers, promoters, promoter_window, pr)
    others = _other_gene_overlaps(spacers, store, pr, len(sites))
    frame = sites.join(
        pairs.set_index("_row").drop(columns=["strand_site", "gene_start", "gene_end"])
    )
    frame = frame.merge(tss, left_index=True, right_index=True, how="left")
    frame = frame.merge(others, left_index=True, right_index=True, how="left")
    frame["gc_percent"] = _gc_percent(frame)
    frame["max_polyt_run"] = _max_polyt(frame)
    return frame.reset_index(drop=True)

crispr_analysis_utils.guide_annotation.annotate¶