Skip to content

crispr_analysis_utils.guide_annotation.db

crispr_analysis_utils.guide_annotation.db

The annotation store: a GENCODE GTF reduced to what annotating needs.

Reading a full GENCODE takes most of a minute, so a run reads this instead. The store is parquet, which pyranges1 loads in about a second.

FEATURES module-attribute 

FEATURES = ('gene', 'transcript', 'exon', 'CDS', 'UTR')

The feature rows the store keeps. Everything else is dropped.

GENE_TYPES module-attribute 

GENE_TYPES = ('protein_coding',)

The gene types a guide is annotated against.

Every call of genic, exon, intron and promoter follows these. Genes of other types are kept as one row each, so a protospacer overlapping one can be reported without that overlap deciding where the guide sits.

COLUMNS module-attribute 

COLUMNS = ('Chromosome', 'Start', 'End', 'Strand', 'Feature', 'ensg', 'gene_symbol', 'enst', 'exon_id', 'exon_number', 'is_reference', 'gene_type')

The store's columns. The first four are what pyranges needs.

StoreError

Bases: RuntimeError

The store is missing, unreadable, or holds no usable annotation.

Source code in src/crispr_analysis_utils/guide_annotation/db.py 
48
49
class StoreError(RuntimeError):
    """The store is missing, unreadable, or holds no usable annotation."""

Store dataclass

A loaded annotation store.

Attributes:

  • features (DataFrame) –

    Every kept row: genes, transcripts, exons, CDS and UTR.

  • provenance (dict) –

    What the store was built from, as build_store recorded it.

Source code in src/crispr_analysis_utils/guide_annotation/db.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
@dataclass
class Store:
    """A loaded annotation store.

    Attributes
    ----------
    features
        Every kept row: genes, transcripts, exons, CDS and UTR.
    provenance
        What the store was built from, as `build_store` recorded it.
    """

    features: pd.DataFrame
    provenance: dict

    @property
    def gene_types(self) -> tuple[str, ...]:
        """The gene types this store's annotation is called against."""
        return tuple(self.provenance.get("gene_types") or GENE_TYPES)

    def rows(
        self, feature: str, *, reference_only: bool = False, primary: bool = True
    ) -> pd.DataFrame:
        """The rows of one feature.

        Parameters
        ----------
        feature
            One of `FEATURES`.
        reference_only
            Keep only the reference transcript's rows.
        primary
            Keep only genes of `gene_types`. On by default, so a store that
            also carries other gene types cannot widen what counts as genic.
        """
        rows = self.features[self.features["Feature"] == feature]
        if primary:
            rows = rows[rows["gene_type"].isin(self.gene_types)]
        if reference_only:
            rows = rows[rows["is_reference"]]
        return rows

    def other_genes(self) -> pd.DataFrame:
        """Gene rows whose type is not one the annotation is called against."""
        genes = self.features[self.features["Feature"] == "gene"]
        return genes[~genes["gene_type"].isin(self.gene_types)]

    @property
    def contigs(self) -> set[str]:
        """Every contig the store holds."""
        return set(self.features["Chromosome"].unique())

gene_types property 

gene_types: tuple[str, ...]

The gene types this store's annotation is called against.

contigs property 

contigs: set[str]

Every contig the store holds.

rows 

rows(feature: str, *, reference_only: bool = False, primary: bool = True) -> pd.DataFrame

The rows of one feature.

Parameters:

  • feature (str) –

    One of FEATURES.

  • reference_only (bool, default: False ) –

    Keep only the reference transcript's rows.

  • primary (bool, default: True ) –

    Keep only genes of gene_types. On by default, so a store that also carries other gene types cannot widen what counts as genic.

Source code in src/crispr_analysis_utils/guide_annotation/db.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def rows(
    self, feature: str, *, reference_only: bool = False, primary: bool = True
) -> pd.DataFrame:
    """The rows of one feature.

    Parameters
    ----------
    feature
        One of `FEATURES`.
    reference_only
        Keep only the reference transcript's rows.
    primary
        Keep only genes of `gene_types`. On by default, so a store that
        also carries other gene types cannot widen what counts as genic.
    """
    rows = self.features[self.features["Feature"] == feature]
    if primary:
        rows = rows[rows["gene_type"].isin(self.gene_types)]
    if reference_only:
        rows = rows[rows["is_reference"]]
    return rows

other_genes 

other_genes() -> pd.DataFrame

Gene rows whose type is not one the annotation is called against.

Source code in src/crispr_analysis_utils/guide_annotation/db.py 
94
95
96
97
def other_genes(self) -> pd.DataFrame:
    """Gene rows whose type is not one the annotation is called against."""
    genes = self.features[self.features["Feature"] == "gene"]
    return genes[~genes["gene_type"].isin(self.gene_types)]

build_store 

build_store(gtf: str | PathLike[str], out: str | PathLike[str], *, gene_types: tuple[str, ...] = GENE_TYPES) -> Store

Reduce a GENCODE GTF to an annotation store.

Parameters:

  • gtf (str | PathLike[str]) –

    A GENCODE GTF, gzipped or not. It must hold transcript and exon rows: a genes-only GTF cannot say exon, intron or transcript.

  • out (str | PathLike[str]) –

    The parquet file to write.

  • gene_types (tuple[str, ...], default: GENE_TYPES ) –

    The gene types a guide is annotated against: genic, exon, intron and promoter all follow these, and every feature row is kept for them. Genes of every other type are kept as one row each, so an overlap can be reported without it deciding where the guide sits.

Returns:

  • Store –

    The store just written.

Raises:

  • StoreError –

    The GTF holds no gene of those types, or no exons.

Source code in src/crispr_analysis_utils/guide_annotation/db.py 
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def build_store(
    gtf: str | os.PathLike[str],
    out: str | os.PathLike[str],
    *,
    gene_types: tuple[str, ...] = GENE_TYPES,
) -> Store:
    """Reduce a GENCODE GTF to an annotation store.

    Parameters
    ----------
    gtf
        A GENCODE GTF, gzipped or not. It must hold `transcript` and `exon`
        rows: a genes-only GTF cannot say exon, intron or transcript.
    out
        The parquet file to write.
    gene_types
        The gene types a guide is annotated against: genic, exon, intron and
        promoter all follow these, and every feature row is kept for them.
        Genes of every other type are kept as one row each, so an overlap can
        be reported without it deciding where the guide sits.

    Returns
    -------
    Store
        The store just written.

    Raises
    ------
    StoreError
        The GTF holds no gene of those types, or no exons.
    """
    import pyranges1 as pr

    gtf = Path(gtf)
    if not gtf.is_file():
        raise FileNotFoundError(f"The GTF {gtf} does not exist.")

    # duplicate_attr is not optional: a GENCODE row carries several `tag`
    # keys, and without it the reader keeps one and drops most of MANE_Select.
    frame = pr.read_gtf(gtf, duplicate_attr=True)
    # Genes of the annotated types keep every feature row; genes of any other
    # type keep one row each, enough to report an overlap without that overlap
    # deciding where a guide sits.
    primary = frame["gene_type"].isin(gene_types)
    frame = frame.loc[
        (primary & frame["Feature"].isin(FEATURES))
        | (~primary & (frame["Feature"] == "gene"))
    ]
    kept_primary = frame["gene_type"].isin(gene_types)
    if not kept_primary.any():
        raise StoreError(
            f"{gtf} holds no {' or '.join(gene_types)} gene. Is it the right GTF?"
        )
    if not ((frame["Feature"] == "exon") & kept_primary).any():
        raise StoreError(
            f"{gtf} holds no exon rows, so it cannot say exon, intron or"
            " transcript. A genes-only GTF is not enough; use the full"
            " GENCODE annotation."
        )

    tags = frame["tag"].fillna("")
    frame = frame.assign(
        ensg=_unversioned(frame["gene_id"]),
        enst=_unversioned(frame["transcript_id"]),
        exon_id=_optional(frame, "exon_id"),
        gene_symbol=frame["gene_name"],
        mane_select=tags.str.contains("MANE_Select"),
        ensembl_canonical=tags.str.contains("Ensembl_canonical"),
    )

    reference = _reference_transcripts(frame)
    frame = frame.assign(
        is_reference=frame["ensg"].map(reference).eq(frame["enst"])
        & frame["enst"].notna()
    )

    provenance = {
        "gtf": str(gtf),
        "gene_types": list(gene_types),
        "other_gene_types": sorted(
            frame["gene_type"][~frame["gene_type"].isin(gene_types)].unique().tolist()
        ),
        "cau_version": __version__,
        "counts": {
            feature: int((frame["Feature"] == feature).sum()) for feature in FEATURES
        },
        "genes": int(
            frame["ensg"][
                (frame["Feature"] == "gene") & frame["gene_type"].isin(gene_types)
            ].nunique()
        ),
        "other_genes": int(
            frame["ensg"][
                (frame["Feature"] == "gene") & ~frame["gene_type"].isin(gene_types)
            ].nunique()
        ),
        "reference_transcripts": int(len(reference)),
        "mane_select": int(
            frame["mane_select"][frame["Feature"] == "transcript"].sum()
        ),
    }

    out = Path(out)
    out.parent.mkdir(parents=True, exist_ok=True)
    # A plain DataFrame, so a built store and a loaded one are the same type.
    table = pd.DataFrame(frame[list(COLUMNS)]).reset_index(drop=True)
    _write_parquet(table, out, provenance)
    return Store(features=table, provenance=provenance)

load_store 

load_store(path: str | PathLike[str]) -> Store

Load a store written by build_store.

Parameters:

  • path (str | PathLike[str]) –

    The parquet file.

Returns:

  • Store –

    Its rows and its provenance.

Raises:

  • StoreError –

    The file is missing, or was not written by build_store.

Source code in src/crispr_analysis_utils/guide_annotation/db.py 
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
def load_store(path: str | os.PathLike[str]) -> Store:
    """Load a store written by `build_store`.

    Parameters
    ----------
    path
        The parquet file.

    Returns
    -------
    Store
        Its rows and its provenance.

    Raises
    ------
    StoreError
        The file is missing, or was not written by `build_store`.
    """
    import pyarrow.parquet as pq

    path = Path(path)
    if not path.is_file():
        raise StoreError(
            f"The annotation store {path} does not exist. Build one with"
            " `cau guide-annotation build-db`."
        )
    table = pq.read_table(path)
    missing = set(COLUMNS) - set(table.column_names)
    if missing:
        # gene_type arrived with the other-gene-type reporting, so a store
        # built before it is simply old rather than foreign.
        rebuild = (
            "Rebuild it with `cau guide-annotation build-db`."
            if missing == {"gene_type"}
            else "Rebuild it."
        )
        raise StoreError(
            f"{path} is missing {', '.join(sorted(missing))}, so it was not"
            f" written by this version of the tool. {rebuild}"
        )
    raw = (table.schema.metadata or {}).get(_METADATA_KEY)
    return Store(
        features=table.to_pandas(),
        provenance=json.loads(raw) if raw else {},
    )