Skip to content

Demultiplexing and pooled screens

Two workflows that both start from a second assay of oligo counts and end in a per-cell call written back into meta_data.

Cell hashing assigns pooled samples from hashtag counts. hto_demux is Seurat's HTODemux — CLR, cluster into k = n_hashtags + 1, per-hashtag negative-binomial background threshold, then singlet / doublet / negative. multiseq_demux is the MULTI-seq alternative, a Gaussian-KDE quantile threshold per barcode. On the cross-species ground truth they are 99.81 % call-concordant with R.

Mixscape separates real CRISPR knockouts from escapers. calc_perturb_sig subtracts each cell's nearest non-targeting controls; run_mixscape then fits the two-component mixture per guide; mixscape_lda builds the supervised map on which each guide population separates.

The CLR margin defaults are deliberate

hto_demux and multiseq_demux normalize across features (margin 1), not across cells. That is what HTODemux does, and it is not the same default as the general-purpose CLR path. Changing it to 2 to "make them consistent" would break agreement with Seurat.

Cell hashing

hto_demux 

hto_demux(seurat, assay: str = 'HTO', positive_quantile: float = 0.99, init: Optional[int] = None, nstarts: int = 10, kfunc: str = 'clara', nsamples: int = 100, normalize: bool = True, margin: int = 1, seed: int = 42, verbose: bool = False)

Demultiplex pooled samples from hashtag counts (Seurat's HTODemux).

Mirrors HTODemux(object, assay = "HTO", positive.quantile = 0.99). Each hashtag's positive/negative cutoff is learned by fitting a negative binomial to the tag's background — the cluster in which it is least expressed — and thresholding at positive_quantile. Cells positive for zero / one / many hashtags are called Negative / Singlet / Doublet.

The object is mutated in place: five <assay>_* metadata columns plus hash.ID are written (matching Seurat), the active identity is set to hash.ID, and the learned cutoffs are stashed in obj.misc["hto_demux"].

Parameters:

  • seurat –

    a Truecell object carrying a hashtag assay.

  • assay (str, default: 'HTO' ) –

    the hashtag assay to demultiplex (default "HTO").

  • positive_quantile (float, default: 0.99 ) –

    quantile of each tag's fitted background at which the positive cutoff is set (Seurat default 0.99).

  • init (Optional[int], default: None ) –

    number of clusters; default n_hashtags + 1.

  • nstarts (int, default: 10 ) –

    kmeans only — restarts (n_init). Seurat uses 100; 10 is a faster, usually-equivalent default. Ignored by clara.

  • kfunc (str, default: 'clara' ) –

    "clara" (default, Seurat's k-medoids — see truecell._clara) or "kmeans". The two rarely disagree on which cluster is a tag's background, so the calls usually match either way (~1% of cells differ on synthetic panels, rising with tag count). Both scale linearly in cells; clara costs a roughly constant 4× (~1.3 s vs ~0.3 s at 100k cells), so the choice is about matching R, not speed.

  • nsamples (int, default: 100 ) –

    clara only — sub-samples to draw (Seurat's default, 100). Ignored by kmeans.

  • normalize (bool, default: True ) –

    CLR-normalize the counts internally for clustering and margins (default). Set False to use the assay's existing data layer (e.g. a prior normalize_data(method="CLR")).

  • margin (int, default: 1 ) –

    CLR margin when normalize is True — 1 (per hashtag across cells; Seurat's default, and what the hashing vignette normalizes with) or 2 (per cell across hashtags).

  • seed (int, default: 42 ) –

    kmeans only — random seed. Has no effect on clara (the default), which draws from its own generator that R cannot seed either — see truecell._clara.

  • verbose (bool, default: False ) –

    print each hashtag's learned cutoff.

Returns:

  • Truecell –

    seurat, with the classification metadata and hash.ID identity.

Source code in truecell/hto.py 
 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
103
104
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
def hto_demux(
    seurat,
    assay: str = "HTO",
    positive_quantile: float = 0.99,
    init: Optional[int] = None,
    nstarts: int = 10,
    kfunc: str = "clara",
    nsamples: int = 100,
    normalize: bool = True,
    margin: int = 1,
    seed: int = 42,
    verbose: bool = False,
):
    """Demultiplex pooled samples from hashtag counts (Seurat's ``HTODemux``).

    Mirrors ``HTODemux(object, assay = "HTO", positive.quantile = 0.99)``. Each
    hashtag's positive/negative cutoff is learned by fitting a negative binomial
    to the tag's background — the cluster in which it is least expressed — and
    thresholding at ``positive_quantile``. Cells positive for zero / one / many
    hashtags are called ``Negative`` / ``Singlet`` / ``Doublet``.

    The object is mutated in place: five ``<assay>_*`` metadata columns plus
    ``hash.ID`` are written (matching Seurat), the active identity is set to
    ``hash.ID``, and the learned cutoffs are stashed in ``obj.misc["hto_demux"]``.

    Parameters
    ----------
    seurat            : a :class:`~truecell.Truecell` object carrying a hashtag assay.
    assay             : the hashtag assay to demultiplex (default ``"HTO"``).
    positive_quantile : quantile of each tag's fitted background at which the
                        positive cutoff is set (Seurat default 0.99).
    init              : number of clusters; default ``n_hashtags + 1``.
    nstarts           : ``kmeans`` only — restarts (``n_init``). Seurat uses 100;
                        10 is a faster, usually-equivalent default. Ignored by
                        ``clara``.
    kfunc             : ``"clara"`` (default, Seurat's k-medoids — see
                        :mod:`truecell._clara`) or ``"kmeans"``. The two rarely
                        disagree on which cluster is a tag's background, so the
                        calls usually match either way (~1% of cells differ on
                        synthetic panels, rising with tag count). Both scale
                        linearly in cells; ``clara`` costs a roughly constant 4×
                        (~1.3 s vs ~0.3 s at 100k cells), so the choice is about
                        matching R, not speed.
    nsamples          : ``clara`` only — sub-samples to draw (Seurat's default,
                        100). Ignored by ``kmeans``.
    normalize         : CLR-normalize the counts internally for clustering and
                        margins (default). Set False to use the assay's existing
                        ``data`` layer (e.g. a prior ``normalize_data(method="CLR")``).
    margin            : CLR margin when ``normalize`` is True — 1 (per hashtag
                        across cells; Seurat's default, and what the hashing
                        vignette normalizes with) or 2 (per cell across hashtags).
    seed              : ``kmeans`` only — random seed. Has no effect on ``clara``
                        (the default), which draws from its own generator that R
                        cannot seed either — see :mod:`truecell._clara`.
    verbose           : print each hashtag's learned cutoff.

    Returns
    -------
    Truecell
        ``seurat``, with the classification metadata and ``hash.ID`` identity.
    """
    if kfunc not in ("kmeans", "clara"):
        raise NotImplementedError(
            f"kfunc={kfunc!r} is not supported; choose from 'kmeans' or 'clara'."
        )

    counts, data, feats, cells = _hto_matrices(seurat, assay, normalize, margin)
    n_htos, n_cells = data.shape
    if n_htos < 2:
        raise ValueError(
            f"HTODemux needs at least 2 hashtags; assay {assay!r} has {n_htos}."
        )

    labels, ncenters = _cluster_cells(data, init, nstarts, seed, kfunc, nsamples)

    # Average (de-logged) expression of each hashtag within each cluster, so the
    # least-expressing cluster can be read off as that hashtag's background.
    expd = np.expm1(data)
    avg = np.full((n_htos, ncenters), np.inf)
    for c in range(ncenters):
        mask = labels == c
        if mask.any():
            avg[:, c] = expd[:, mask].mean(axis=1)

    # Per-hashtag negative-binomial threshold on its background cluster.
    discrete = np.zeros((n_htos, n_cells), dtype=bool)
    cutoffs: dict[str, float] = {}
    for i in range(n_htos):
        neg_cluster = int(np.argmin(avg[i]))
        values_use = counts[i, labels == neg_cluster]
        cutoff = _positive_cutoff(values_use, positive_quantile)
        cutoffs[feats[i]] = cutoff
        discrete[i] = counts[i] > cutoff
        if verbose:
            print(f"Cutoff for {feats[i]}: {cutoff:g} reads")

    _write_classification(
        seurat, assay, data, discrete, feats, cells,
    )
    seurat.misc.setdefault("hto_demux", {})[assay] = {
        "cutoffs": cutoffs,
        "ncenters": ncenters,
        "positive_quantile": positive_quantile,
    }
    return seurat

multiseq_demux 

multiseq_demux(seurat, assay: str = 'HTO', quantile: float = 0.7, autothresh: bool = False, maxiter: int = 5, qrange: Optional[Sequence[float]] = None, normalize: bool = True, margin: int = 1, verbose: bool = False)

Demultiplex pooled samples from barcode counts (Seurat's MULTIseqDemux).

Mirrors MULTIseqDemux(object, assay = "HTO", quantile = 0.7). For each barcode a Gaussian-kernel-density threshold is placed a fraction quantile of the way between its background and positive modes; cells positive for zero / one / many barcodes are called Negative / Singlet / Doublet.

The object is mutated in place: MULTI_ID and MULTI_classification metadata columns are written, the active identity is set to MULTI_ID, and the learned thresholds are stashed in obj.misc["multiseq_demux"].

Parameters:

  • seurat –

    a Truecell object carrying a hashtag/barcode assay.

  • assay (str, default: 'HTO' ) –

    the barcode assay to demultiplex (default "HTO").

  • quantile (float, default: 0.7 ) –

    fraction between each barcode's background and positive modes at which its positive cutoff is placed (Seurat default 0.7). Ignored when autothresh is True.

  • autothresh (bool, default: False ) –

    sweep qrange for the quantile that maximizes the singlet rate, iteratively removing negatives and re-thresholding the remainder (up to maxiter rounds). Overrides quantile.

  • maxiter (int, default: 5 ) –

    maximum auto-threshold rounds (default 5).

  • qrange (Optional[Sequence[float]], default: None ) –

    quantiles swept when autothresh is True (default 0.1, 0.15, … 0.9).

  • normalize (bool, default: True ) –

    CLR-normalize the counts internally (default). Set False to use the assay's existing data layer (e.g. a prior normalize_data(method="CLR")).

  • margin (int, default: 1 ) –

    CLR margin when normalize is True — 1 (per barcode across cells; Seurat's default) or 2 (per cell across barcodes).

  • verbose (bool, default: False ) –

    print each auto-threshold round's chosen quantile.

Returns:

  • Truecell –

    seurat, with the MULTI_ID classification and identity.

Source code in truecell/multiseq.py 
 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
103
104
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
def multiseq_demux(
    seurat,
    assay: str = "HTO",
    quantile: float = 0.7,
    autothresh: bool = False,
    maxiter: int = 5,
    qrange: Optional[Sequence[float]] = None,
    normalize: bool = True,
    margin: int = 1,
    verbose: bool = False,
):
    """Demultiplex pooled samples from barcode counts (Seurat's ``MULTIseqDemux``).

    Mirrors ``MULTIseqDemux(object, assay = "HTO", quantile = 0.7)``. For each
    barcode a Gaussian-kernel-density threshold is placed a fraction ``quantile``
    of the way between its background and positive modes; cells positive for zero /
    one / many barcodes are called ``Negative`` / ``Singlet`` / ``Doublet``.

    The object is mutated in place: ``MULTI_ID`` and ``MULTI_classification``
    metadata columns are written, the active identity is set to ``MULTI_ID``, and
    the learned thresholds are stashed in ``obj.misc["multiseq_demux"]``.

    Parameters
    ----------
    seurat     : a :class:`~truecell.Truecell` object carrying a hashtag/barcode assay.
    assay      : the barcode assay to demultiplex (default ``"HTO"``).
    quantile   : fraction between each barcode's background and positive modes at
                 which its positive cutoff is placed (Seurat default 0.7). Ignored
                 when ``autothresh`` is True.
    autothresh : sweep ``qrange`` for the quantile that maximizes the singlet rate,
                 iteratively removing negatives and re-thresholding the remainder
                 (up to ``maxiter`` rounds). Overrides ``quantile``.
    maxiter    : maximum auto-threshold rounds (default 5).
    qrange     : quantiles swept when ``autothresh`` is True (default
                 ``0.1, 0.15, … 0.9``).
    normalize  : CLR-normalize the counts internally (default). Set False to use the
                 assay's existing ``data`` layer (e.g. a prior
                 ``normalize_data(method="CLR")``).
    margin     : CLR margin when ``normalize`` is True — 1 (per barcode across cells;
                 Seurat's default) or 2 (per cell across barcodes).
    verbose    : print each auto-threshold round's chosen quantile.

    Returns
    -------
    Truecell
        ``seurat``, with the ``MULTI_ID`` classification and identity.
    """
    _, data, feats, cells = _hto_matrices(seurat, assay, normalize, margin)
    n_bc, n_cells = data.shape
    if n_bc < 2:
        raise ValueError(
            f"MULTIseqDemux needs at least 2 barcodes; assay {assay!r} has {n_bc}."
        )

    # A local array rather than rebinding the Sequence[float]|None parameter.
    qs = (np.round(np.arange(0.1, 0.9 + 1e-9, 0.05), 4) if qrange is None
          else np.asarray(qrange, dtype=float))

    if autothresh:
        calls, thresholds, q_used = _auto_classify(
            data, feats, qs, maxiter, verbose
        )
    else:
        calls, thresholds = _classify_cells(data, feats, quantile)
        q_used = quantile

    _write_calls(seurat, calls, cells)
    seurat.misc.setdefault("multiseq_demux", {})[assay] = {
        "thresholds": thresholds,
        "quantile": float(q_used),
        "autothresh": bool(autothresh),
    }
    return seurat

Mixscape

calc_perturb_sig 

calc_perturb_sig(seurat, assay: str = 'RNA', features: Optional[Sequence[str]] = None, layer: str = 'data', labels: str = 'gene', nt_class: str = 'NT', split_by: Optional[str] = None, num_neighbors: int = 20, reduction: str = 'pca', ndims: int = 15, new_assay: str = 'PRTB')

Compute each cell's local perturbation signature (Seurat's CalcPerturbSig).

Mirrors CalcPerturbSig(object, assay, gd.class, nt.cell.class, reduction, ndims, num.neighbors, new.assay.name = "PRTB"). For every cell, the mean expression of its num_neighbors nearest non-targeting (NT) control cells — in the first ndims dimensions of reduction — is subtracted from its own expression. The residual (the deviation from the controls the cell most resembles) is stored as a new assay, ready for run_mixscape.

Parameters:

  • seurat –

    a Truecell object with a guide-assignment metadata column and a computed reduction.

  • assay (str, default: 'RNA' ) –

    source expression assay (default "RNA").

  • features (Optional[Sequence[str]], default: None ) –

    genes to include (default: the assay's variable features, or all features if none are set).

  • layer (str, default: 'data' ) –

    expression layer to difference (default "data", the log-normalized values).

  • labels (str, default: 'gene' ) –

    metadata column holding each cell's target-gene / guide class.

  • nt_class (str, default: 'NT' ) –

    the value in labels marking non-targeting control cells.

  • split_by (Optional[str], default: None ) –

    optional metadata column; neighbours are found only within the same group (e.g. replicate), so batch is not mistaken for signal.

  • num_neighbors (int, default: 20 ) –

    NT neighbours averaged per cell (Seurat default 20); capped at the number of NT cells available in the group.

  • reduction (str, default: 'pca' ) –

    reduction whose embedding defines "nearest" (default "pca").

  • ndims (int, default: 15 ) –

    leading dimensions of reduction to use (default 15).

  • new_assay (str, default: 'PRTB' ) –

    name of the perturbation-signature assay to create (default "PRTB").

Returns:

  • Truecell –

    seurat, with the perturbation-signature assay new_assay attached.

Source code in truecell/mixscape.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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
def calc_perturb_sig(
    seurat,
    assay: str = "RNA",
    features: Optional[Sequence[str]] = None,
    layer: str = "data",
    labels: str = "gene",
    nt_class: str = "NT",
    split_by: Optional[str] = None,
    num_neighbors: int = 20,
    reduction: str = "pca",
    ndims: int = 15,
    new_assay: str = "PRTB",
):
    """Compute each cell's local perturbation signature (Seurat's ``CalcPerturbSig``).

    Mirrors ``CalcPerturbSig(object, assay, gd.class, nt.cell.class, reduction,
    ndims, num.neighbors, new.assay.name = "PRTB")``. For every cell, the mean
    expression of its ``num_neighbors`` nearest non-targeting (NT) control cells
    — in the first ``ndims`` dimensions of ``reduction`` — is subtracted from its
    own expression. The residual (the deviation from the controls the cell most
    resembles) is stored as a new assay, ready for :func:`run_mixscape`.

    Parameters
    ----------
    seurat        : a :class:`~truecell.Truecell` object with a guide-assignment
                    metadata column and a computed ``reduction``.
    assay         : source expression assay (default ``"RNA"``).
    features      : genes to include (default: the assay's variable features, or
                    all features if none are set).
    layer         : expression layer to difference (default ``"data"``, the
                    log-normalized values).
    labels        : metadata column holding each cell's target-gene / guide class.
    nt_class      : the value in ``labels`` marking non-targeting control cells.
    split_by      : optional metadata column; neighbours are found only within the
                    same group (e.g. replicate), so batch is not mistaken for signal.
    num_neighbors : NT neighbours averaged per cell (Seurat default 20); capped at
                    the number of NT cells available in the group.
    reduction     : reduction whose embedding defines "nearest" (default ``"pca"``).
    ndims         : leading dimensions of ``reduction`` to use (default 15).
    new_assay     : name of the perturbation-signature assay to create (default
                    ``"PRTB"``).

    Returns
    -------
    Truecell
        ``seurat``, with the perturbation-signature assay ``new_assay`` attached.
    """
    from .assay5 import create_assay5_object

    emb = seurat.embeddings(reduction)
    if ndims is not None:
        emb = emb[:, :ndims]
    emb = np.asarray(emb, dtype=float)

    data, feats, cells = _layer_matrix(seurat.assays[assay], layer)
    if features is not None:
        keep = [f for f in features if f in set(feats)]
        idx = [feats.index(f) for f in keep]
        data = data[idx, :]
        feats = keep
    if not feats:
        raise ValueError("No features selected for the perturbation signature.")

    labels_vec = _aligned_meta(seurat, labels, cells)
    nt_mask = labels_vec == nt_class
    if not nt_mask.any():
        raise ValueError(
            f"No non-targeting cells: column {labels!r} has no value {nt_class!r}."
        )

    groups = _split_groups(seurat, split_by, cells)

    signature = np.array(data, dtype=float, copy=True)
    for gidx in groups:
        nt_local = gidx[nt_mask[gidx]]
        if nt_local.size == 0:
            # No controls in this split — fall back to the whole-dataset NT mean.
            control = data[:, nt_mask].mean(axis=1, keepdims=True)
            signature[:, gidx] = data[:, gidx] - control
            continue
        k = int(min(num_neighbors, nt_local.size))
        neigh = _nt_neighbor_means(emb, gidx, nt_local, data, k)
        signature[:, gidx] = data[:, gidx] - neigh

    prtb = create_assay5_object(
        data=sp.csc_matrix(signature),
        feature_names=list(feats),
        cell_names=list(cells),
        key=f"{new_assay.lower()}_",
    )
    prtb.variable_features = list(feats)
    seurat.assays[new_assay] = prtb
    seurat.misc.setdefault("calc_perturb_sig", {})[new_assay] = {
        "assay": assay,
        "reduction": reduction,
        "ndims": ndims,
        "num_neighbors": num_neighbors,
        "n_features": len(feats),
    }
    return seurat

run_mixscape 

run_mixscape(seurat, assay: str = 'PRTB', labels: str = 'gene', nt_class: str = 'NT', de_assay: str = 'RNA', layer: str = 'data', min_de_genes: int = 5, min_cells: int = 5, logfc_threshold: float = 0.25, min_pct: float = 0.05, pval_cutoff: float = 0.05, iter_num: int = 10, prtb_type: str = 'KO', new_class: str = 'mixscape_class', de_test: str = 'wilcox', seed: int = 0, verbose: bool = False)

Classify perturbed vs. escaping cells per guide (Seurat's RunMixscape).

Mirrors RunMixscape(object, assay = "PRTB", labels = "gene", nt.class.name = "NT", de.assay = "RNA", min.de.genes = 5, iter.num = 10, prtb.type = "KO"). Operating on the perturbation signature from calc_perturb_sig, each target gene's cells are split into knockout (KO) and non-perturbed (NP) by an iterative two-component Gaussian mixture over their projection onto the gene's perturbation vector (see the module docstring). NT cells stay NT.

The object is mutated in place: mixscape_class (also set as the active identity), mixscape_class.global, and mixscape_class_p_<type> metadata columns are written, and per-gene bookkeeping is stashed in obj.misc["mixscape"].

Parameters:

  • seurat –

    a Truecell object carrying the assay perturbation signature and a labels guide column.

  • assay (str, default: 'PRTB' ) –

    perturbation-signature assay (default "PRTB").

  • labels (str, default: 'gene' ) –

    metadata column of per-cell target-gene / guide class.

  • nt_class (str, default: 'NT' ) –

    value in labels marking non-targeting controls.

  • de_assay (str, default: 'RNA' ) –

    assay used for the gene-vs-NT differential expression (default "RNA").

  • layer (str, default: 'data' ) –

    signature layer to project (default "data").

  • min_de_genes (int, default: 5 ) –

    a gene needs at least this many DE genes to be testable; otherwise all its cells are NP (Seurat default 5).

  • min_cells (int, default: 5 ) –

    a gene needs at least this many cells; otherwise NP.

  • logfc_threshold (float, default: 0.25 ) –

    |avg_log2FC| DE cutoff passed to find_markers (0.25).

  • min_pct (float, default: 0.05 ) –

    min.pct DE cutoff passed to find_markers (0.05).

  • pval_cutoff (float, default: 0.05 ) –

    adjusted-p cutoff a DE gene must clear (0.05).

  • iter_num (int, default: 10 ) –

    maximum mixture-refinement rounds per gene (Seurat 10).

  • prtb_type (str, default: 'KO' ) –

    label for the perturbed class (default "KO"; use e.g. "KD" for a knock-down screen). Also names the posterior column mixscape_class_p_<type>.

  • new_class (str, default: 'mixscape_class' ) –

    base name for the output columns / identity (default "mixscape_class").

  • de_test (str, default: 'wilcox' ) –

    find_markers test for the gene-vs-NT DE (default "wilcox").

  • seed (int, default: 0 ) –

    random state for the Gaussian mixture (determinism).

  • verbose (bool, default: False ) –

    print each gene's DE-gene count and final KO count.

Returns:

  • Truecell –

    seurat, with the mixscape_class classification and identity.

Source code in truecell/mixscape.py 
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
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
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
def run_mixscape(
    seurat,
    assay: str = "PRTB",
    labels: str = "gene",
    nt_class: str = "NT",
    de_assay: str = "RNA",
    layer: str = "data",
    min_de_genes: int = 5,
    min_cells: int = 5,
    logfc_threshold: float = 0.25,
    min_pct: float = 0.05,
    pval_cutoff: float = 0.05,
    iter_num: int = 10,
    prtb_type: str = "KO",
    new_class: str = "mixscape_class",
    de_test: str = "wilcox",
    seed: int = 0,
    verbose: bool = False,
):
    """Classify perturbed vs. escaping cells per guide (Seurat's ``RunMixscape``).

    Mirrors ``RunMixscape(object, assay = "PRTB", labels = "gene",
    nt.class.name = "NT", de.assay = "RNA", min.de.genes = 5, iter.num = 10,
    prtb.type = "KO")``. Operating on the perturbation signature from
    :func:`calc_perturb_sig`, each target gene's cells are split into knockout
    (``KO``) and non-perturbed (``NP``) by an iterative two-component Gaussian
    mixture over their projection onto the gene's perturbation vector (see the
    module docstring). NT cells stay ``NT``.

    The object is mutated in place: ``mixscape_class`` (also set as the active
    identity), ``mixscape_class.global``, and ``mixscape_class_p_<type>`` metadata
    columns are written, and per-gene bookkeeping is stashed in
    ``obj.misc["mixscape"]``.

    Parameters
    ----------
    seurat          : a :class:`~truecell.Truecell` object carrying the ``assay``
                      perturbation signature and a ``labels`` guide column.
    assay           : perturbation-signature assay (default ``"PRTB"``).
    labels          : metadata column of per-cell target-gene / guide class.
    nt_class        : value in ``labels`` marking non-targeting controls.
    de_assay        : assay used for the gene-vs-NT differential expression
                      (default ``"RNA"``).
    layer           : signature layer to project (default ``"data"``).
    min_de_genes    : a gene needs at least this many DE genes to be testable;
                      otherwise all its cells are NP (Seurat default 5).
    min_cells       : a gene needs at least this many cells; otherwise NP.
    logfc_threshold : ``|avg_log2FC|`` DE cutoff passed to ``find_markers`` (0.25).
    min_pct         : ``min.pct`` DE cutoff passed to ``find_markers`` (0.05).
    pval_cutoff     : adjusted-p cutoff a DE gene must clear (0.05).
    iter_num        : maximum mixture-refinement rounds per gene (Seurat 10).
    prtb_type       : label for the perturbed class (default ``"KO"``; use e.g.
                      ``"KD"`` for a knock-down screen). Also names the posterior
                      column ``mixscape_class_p_<type>``.
    new_class       : base name for the output columns / identity
                      (default ``"mixscape_class"``).
    de_test         : ``find_markers`` test for the gene-vs-NT DE (default
                      ``"wilcox"``).
    seed            : random state for the Gaussian mixture (determinism).
    verbose         : print each gene's DE-gene count and final KO count.

    Returns
    -------
    Truecell
        ``seurat``, with the ``mixscape_class`` classification and identity.
    """
    sig, sig_feats, cells = _layer_matrix(seurat.assays[assay], layer)
    sig_feat_idx = {f: i for i, f in enumerate(sig_feats)}

    labels_vec = _aligned_meta(seurat, labels, cells)
    nt_idx = np.where(labels_vec == nt_class)[0]
    if nt_idx.size == 0:
        raise ValueError(
            f"No non-targeting cells: column {labels!r} has no value {nt_class!r}."
        )

    genes = sorted(
        {g for g in labels_vec.tolist() if isinstance(g, str) and g != nt_class}
    )
    if not genes:
        raise ValueError(f"No target genes to test in column {labels!r}.")

    n_cells = len(cells)
    mixscape_class = np.array(
        [str(x) for x in labels_vec], dtype=object
    )                                                # NT stays "NT"; genes filled below
    global_class = np.where(labels_vec == nt_class, "NT", "NP").astype(object)
    p_prtb = np.full(n_cells, np.nan, dtype=float)
    bookkeeping: dict[str, dict] = {}

    # find_markers reads the active identity, so drive the gene-vs-NT DE off the
    # guide labels — restored afterwards.
    saved_ident = pd.Categorical(list(seurat.idents))
    seurat.idents = list(labels_vec)
    try:
        for gene in genes:
            gene_local = np.where(labels_vec == gene)[0]
            # Counts plus a "scores" DataFrame-or-None, so not dict[str, int].
            info: dict[str, Any] = {
                "n_cells": int(gene_local.size), "n_de": 0, "n_iter": 0, "n_ko": 0,
            }

            if gene_local.size < min_cells:
                _assign(mixscape_class, global_class, gene_local, gene, "NP", prtb_type)
                info["scores"] = None
                bookkeeping[gene] = info
                if verbose:
                    print(f"[mixscape] {gene}: {gene_local.size} cells < min_cells → NP")
                continue

            de_genes = _de_genes(
                seurat, gene, nt_class, de_assay, de_test,
                logfc_threshold, min_pct, pval_cutoff, sig_feat_idx,
            )
            info["n_de"] = len(de_genes)
            if len(de_genes) < min_de_genes:
                _assign(mixscape_class, global_class, gene_local, gene, "NP", prtb_type)
                info["scores"] = None
                bookkeeping[gene] = info
                if verbose:
                    print(f"[mixscape] {gene}: {len(de_genes)} DE genes < min → NP")
                continue

            de_rows = [sig_feat_idx[g] for g in de_genes]
            ko_pos, post, n_iter, score = _mixscape_em(
                sig, de_rows, nt_idx, gene_local, iter_num, seed,
            )
            info["n_iter"] = n_iter
            info["n_ko"] = int(len(ko_pos))
            # R's gv data.frame: one `pvec` column plus the guide-label column,
            # indexed by cell — what PlotPerturbScore reads back out.
            if score is None:
                info["scores"] = None
            else:
                score_cells = [
                    cells[i] for i in np.concatenate([nt_idx, gene_local])
                ]
                info["scores"] = pd.DataFrame(
                    {
                        "pvec": score,
                        labels: [nt_class] * nt_idx.size + [gene] * gene_local.size,
                    },
                    index=score_cells,
                )

            ko_set = set(int(p) for p in ko_pos)
            for pos, cell in enumerate(gene_local):
                if pos in ko_set:
                    mixscape_class[cell] = f"{gene} {prtb_type}"
                    global_class[cell] = prtb_type
                else:
                    mixscape_class[cell] = f"{gene} NP"
                    global_class[cell] = "NP"
                p_prtb[cell] = post[pos]
            bookkeeping[gene] = info
            if verbose:
                print(
                    f"[mixscape] {gene}: {len(de_genes)} DE genes, "
                    f"{len(ko_set)}/{gene_local.size} {prtb_type} in {n_iter} iters"
                )
    finally:
        seurat.idents = saved_ident

    target = seurat.cell_names()

    def put(col, values):
        seurat.meta_data[col] = (
            pd.Series(list(values), index=cells).reindex(target).values
        )

    put(new_class, mixscape_class)
    put(f"{new_class}.global", global_class)
    put(f"{new_class}_p_{prtb_type.lower()}", p_prtb)
    seurat.idents = list(
        pd.Series(list(mixscape_class), index=cells).reindex(target).values
    )
    seurat.misc.setdefault("mixscape", {})[assay] = {
        "genes": bookkeeping,
        "nt_class": nt_class,
        "prtb_type": prtb_type,
    }
    return seurat

mixscape_lda 

mixscape_lda(seurat, labels: str = 'gene', nt_class: str = 'NT', assay: str = 'PRTB', de_assay: str = 'RNA', layer: str = 'data', npcs: int = 10, logfc_threshold: float = 0.25, min_pct: float = 0.1, pval_cutoff: float = 0.05, de_test: str = 'wilcox', reduction_name: str = 'lda', reduction_key: str = 'LDA_', scale_max: float = 10.0, seed: int = 42, verbose: bool = False)

Linear-discriminant projection that separates the guide classes (Seurat's MixscapeLDA).

Mirrors MixscapeLDA(object, pc.assay = "PRTB", labels = "gene", nt.class.name = "NT", npcs = 10), which asks a complementary question to run_mixscape: not which cells are perturbed but how do the whole guide populations differ from one another and from control. It builds a single supervised 2-D-ish map on which every guide class (and NT) forms its own cloud, the classic mixscape LDA plot. The only prerequisite is a perturbation-signature assay from calc_perturb_sig; the mixscape KO/NP calls are not used — cells are grouped by their raw guide label.

The construction (Seurat's PrepLDA → RunLDA):

  1. Per-guide feature blocks. For each target gene, its cells are tested against NT (on de_assay) to find that guide's response genes; a guide with fewer than npcs + 1 such genes is dropped (it cannot support npcs components). Restricted to those genes on the perturbation-signature assay, a PCA is fit on that guide's cells plus the NT cells, and then every cell in the dataset is projected onto that guide's npcs-dim subspace. Each surviving guide thus contributes npcs columns describing where all cells fall along its perturbation axes.
  2. One LDA over the concatenation. The per-guide blocks are stacked side by side into one cell × (guides · npcs) matrix and a single linear discriminant analysis (sklearn LinearDiscriminantAnalysis) is fit with the guide label as the class — finding the n_classes − 1 directions that best separate the guide populations (including NT). The discriminant scores are stored as a reduction (default "lda", key "LDA_"), and the per-cell class assignment and posteriors are written to metadata (lda_assignments and LDAP_<class>).

Two choices differ from a literal reading of R, both documented:

  • The per-guide subspace is read from the signature's data layer, scaled against the guide-plus-NT reference, in place of Seurat's ScaleData → RunPCA → ProjectCellEmbeddings chain. The composition is the same map: centre/scale each response gene by the reference cells' mean and SD, project through the reference PCA loadings.
  • The leave-one-out CV posterior (MASS lda(..., CV = TRUE)) that Seurat stashes in misc is not computed; only the resubstitution assignment and posterior are kept, which is all the plot and the metadata columns use.

Parameters:

  • seurat –

    a Truecell object carrying the assay perturbation signature (from calc_perturb_sig) and a labels guide column.

  • labels (str, default: 'gene' ) –

    metadata column of per-cell target-gene / guide class.

  • nt_class (str, default: 'NT' ) –

    value in labels marking non-targeting controls.

  • assay (str, default: 'PRTB' ) –

    perturbation-signature assay projected for the LDA features (default "PRTB").

  • de_assay (str, default: 'RNA' ) –

    assay for the per-guide guide-vs-NT differential expression (default "RNA").

  • layer (str, default: 'data' ) –

    signature layer to project (default "data").

  • npcs (int, default: 10 ) –

    per-guide PCA components (Seurat default 10); a guide needs at least npcs + 1 DE genes to contribute.

  • logfc_threshold (float, default: 0.25 ) –

    |avg_log2FC| DE cutoff passed to find_markers (0.25).

  • min_pct (float, default: 0.1 ) –

    min.pct DE cutoff passed to find_markers (0.1, as in Seurat's TopDEGenesMixscape).

  • pval_cutoff (float, default: 0.05 ) –

    adjusted-p cutoff a DE gene must clear (0.05).

  • de_test (str, default: 'wilcox' ) –

    find_markers test for the guide-vs-NT DE ("wilcox").

  • reduction_name (str, default: 'lda' ) –

    key under which the LDA reduction is stored (default "lda").

  • reduction_key (str, default: 'LDA_' ) –

    prefix for the discriminant dimension names ("LDA_").

  • scale_max (float, default: 10.0 ) –

    clip scaled values to ±scale_max before PCA, matching Seurat's ScaleData (default 10; None to disable).

  • seed (int, default: 42 ) –

    random state for the per-guide PCA (determinism).

  • verbose (bool, default: False ) –

    print each guide's DE-gene count and whether it contributed.

Returns:

  • Truecell –

    seurat, with the reduction_name LDA reduction and the lda_assignments / LDAP_<class> metadata columns.

Source code in truecell/mixscape.py 
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
def mixscape_lda(
    seurat,
    labels: str = "gene",
    nt_class: str = "NT",
    assay: str = "PRTB",
    de_assay: str = "RNA",
    layer: str = "data",
    npcs: int = 10,
    logfc_threshold: float = 0.25,
    min_pct: float = 0.1,
    pval_cutoff: float = 0.05,
    de_test: str = "wilcox",
    reduction_name: str = "lda",
    reduction_key: str = "LDA_",
    scale_max: float = 10.0,
    seed: int = 42,
    verbose: bool = False,
):
    """Linear-discriminant projection that separates the guide classes (Seurat's ``MixscapeLDA``).

    Mirrors ``MixscapeLDA(object, pc.assay = "PRTB", labels = "gene",
    nt.class.name = "NT", npcs = 10)``, which asks a complementary question to
    :func:`run_mixscape`: not *which cells are perturbed* but *how do the whole
    guide populations differ from one another and from control*. It builds a
    single supervised 2-D-ish map on which every guide class (and NT) forms its
    own cloud, the classic mixscape LDA plot. The only prerequisite is a
    perturbation-signature assay from :func:`calc_perturb_sig`; the mixscape KO/NP
    calls are **not** used — cells are grouped by their raw guide label.

    The construction (Seurat's ``PrepLDA`` → ``RunLDA``):

    1. **Per-guide feature blocks.** For each target gene, its cells are tested
       against NT (on ``de_assay``) to find that guide's response genes; a guide
       with fewer than ``npcs + 1`` such genes is dropped (it cannot support
       ``npcs`` components). Restricted to those genes on the perturbation-signature
       ``assay``, a PCA is fit on that guide's cells **plus** the NT cells, and then
       **every** cell in the dataset is projected onto that guide's ``npcs``-dim
       subspace. Each surviving guide thus contributes ``npcs`` columns describing
       where all cells fall along *its* perturbation axes.
    2. **One LDA over the concatenation.** The per-guide blocks are stacked side by
       side into one cell × (guides · ``npcs``) matrix and a single linear
       discriminant analysis (``sklearn`` ``LinearDiscriminantAnalysis``) is fit
       with the guide label as the class — finding the ``n_classes − 1`` directions
       that best separate the guide populations (including NT). The discriminant
       scores are stored as a reduction (default ``"lda"``, key ``"LDA_"``), and
       the per-cell class assignment and posteriors are written to metadata
       (``lda_assignments`` and ``LDAP_<class>``).

    Two choices differ from a literal reading of R, both documented:

    * **The per-guide subspace is read from the signature's ``data`` layer, scaled
      against the guide-plus-NT reference**, in place of Seurat's ``ScaleData`` →
      ``RunPCA`` → ``ProjectCellEmbeddings`` chain. The composition is the same map:
      centre/scale each response gene by the reference cells' mean and SD, project
      through the reference PCA loadings.
    * **The leave-one-out CV posterior** (MASS ``lda(..., CV = TRUE)``) that Seurat
      stashes in ``misc`` is not computed; only the resubstitution assignment and
      posterior are kept, which is all the plot and the metadata columns use.

    Parameters
    ----------
    seurat          : a :class:`~truecell.Truecell` object carrying the ``assay``
                      perturbation signature (from :func:`calc_perturb_sig`) and a
                      ``labels`` guide column.
    labels          : metadata column of per-cell target-gene / guide class.
    nt_class        : value in ``labels`` marking non-targeting controls.
    assay           : perturbation-signature assay projected for the LDA features
                      (default ``"PRTB"``).
    de_assay        : assay for the per-guide guide-vs-NT differential expression
                      (default ``"RNA"``).
    layer           : signature layer to project (default ``"data"``).
    npcs            : per-guide PCA components (Seurat default 10); a guide needs at
                      least ``npcs + 1`` DE genes to contribute.
    logfc_threshold : ``|avg_log2FC|`` DE cutoff passed to ``find_markers`` (0.25).
    min_pct         : ``min.pct`` DE cutoff passed to ``find_markers`` (0.1, as in
                      Seurat's ``TopDEGenesMixscape``).
    pval_cutoff     : adjusted-p cutoff a DE gene must clear (0.05).
    de_test         : ``find_markers`` test for the guide-vs-NT DE (``"wilcox"``).
    reduction_name  : key under which the LDA reduction is stored (default
                      ``"lda"``).
    reduction_key   : prefix for the discriminant dimension names (``"LDA_"``).
    scale_max       : clip scaled values to ``±scale_max`` before PCA, matching
                      Seurat's ``ScaleData`` (default 10; ``None`` to disable).
    seed            : random state for the per-guide PCA (determinism).
    verbose         : print each guide's DE-gene count and whether it contributed.

    Returns
    -------
    Truecell
        ``seurat``, with the ``reduction_name`` LDA reduction and the
        ``lda_assignments`` / ``LDAP_<class>`` metadata columns.
    """
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

    sig, sig_feats, cells = _layer_matrix(seurat.assays[assay], layer)
    sig_feat_idx = {f: i for i, f in enumerate(sig_feats)}

    labels_vec = _aligned_meta(seurat, labels, cells)
    nt_idx = np.where(labels_vec == nt_class)[0]
    if nt_idx.size == 0:
        raise ValueError(
            f"No non-targeting cells: column {labels!r} has no value {nt_class!r}."
        )

    genes = sorted(
        {g for g in labels_vec.tolist() if isinstance(g, str) and g != nt_class}
    )
    if not genes:
        raise ValueError(f"No target genes to test in column {labels!r}.")

    blocks: list[np.ndarray] = []
    feat_names: list[str] = []
    genes_used: list[str] = []

    # find_markers reads the active identity — drive the guide-vs-NT DE off the
    # guide labels, restored afterwards.
    saved_ident = pd.Categorical(list(seurat.idents))
    seurat.idents = list(labels_vec)
    try:
        for gene in genes:
            gene_local = np.where(labels_vec == gene)[0]
            de_genes = _de_genes(
                seurat, gene, nt_class, de_assay, de_test,
                logfc_threshold, min_pct, pval_cutoff, sig_feat_idx,
            )
            if len(de_genes) < npcs + 1:
                if verbose:
                    print(
                        f"[mixscape_lda] {gene}: {len(de_genes)} DE genes "
                        f"< npcs+1={npcs + 1} → skipped"
                    )
                continue
            de_rows = [sig_feat_idx[g] for g in de_genes]
            block = _lda_guide_block(
                sig, de_rows, nt_idx, gene_local, npcs, scale_max, seed,
            )
            blocks.append(block)
            feat_names.extend(f"{gene}_PC_{k + 1}" for k in range(block.shape[1]))
            genes_used.append(gene)
            if verbose:
                print(
                    f"[mixscape_lda] {gene}: {len(de_genes)} DE genes → "
                    f"{block.shape[1]} PCs"
                )
    finally:
        seurat.idents = saved_ident

    if not blocks:
        raise ValueError(
            f"No guide reached npcs+1={npcs + 1} DE genes; lower `npcs` or the DE "
            f"thresholds."
        )

    features = np.hstack(blocks)                     # cells × (guides · npcs)
    y = np.array([str(v) for v in labels_vec], dtype=object)
    n_comp = min(len(set(y.tolist())) - 1, features.shape[1])

    lda = LinearDiscriminantAnalysis(n_components=n_comp)
    lda.fit(features, y)
    embeddings = lda.transform(features)             # cells × n_comp
    loadings = np.asarray(lda.scalings_)[:, :n_comp]
    assignments = lda.predict(features)
    posterior = lda.predict_proba(features)          # cells × n_classes
    classes = [str(c) for c in lda.classes_]

    from .dimreduc import DimReduc

    dr = DimReduc(
        cell_embeddings=embeddings,
        cell_names=list(cells),
        feature_loadings=loadings,
        feature_names=list(feat_names),
        assay_used=assay,
        key=reduction_key,
        misc={
            "assignments": list(assignments),
            "posterior": posterior,
            "classes": classes,
            "genes_used": genes_used,
            "npcs": npcs,
        },
    )
    seurat.reductions[reduction_name] = dr

    target = seurat.cell_names()

    def put(col, values):
        seurat.meta_data[col] = (
            pd.Series(list(values), index=cells).reindex(target).values
        )

    put("lda_assignments", assignments)
    for j, cls in enumerate(classes):
        put(f"LDAP_{cls}", posterior[:, j])
    return seurat