Signature scoring

Score a gene programme per cell against a background of expression-matched control genes. cell_cycle_scoring is add_module_score run twice, on the S and G2/M lists, plus the discrete phase call.

The control genes are drawn at random from expression bins, so these scores carry an RNG. Against R Seurat the per-cell phase call is 96.6 % concordant and the continuous scores correlate at Pearson ≥ 0.998 — the residual is the control draw, and nothing else. Cell-cycle vignette.

add_module_score

add_module_score(seurat, features: Union[Sequence[str], Sequence[Sequence[str]], dict], pool: Optional[list[str]] = None, nbin: int = 24, ctrl: int = 100, name: str = 'Cluster', assay: Optional[str] = None, layer: str = 'data', seed: int = 1, search: bool = False) -> 'object'

Score one or more gene programs per cell.

Mirrors R's AddModuleScore(): each program's score is the mean expression of its genes minus the mean expression of a control set drawn from the same average-expression bins (so highly/lowly expressed genes are controlled for).

Parameters:

• features (Union[Sequence[str], Sequence[Sequence[str]], dict]) –

  a single gene list, a list of gene lists, or a name->list dict.

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

  genes to sample controls from (default: all features).

• nbin (int, default: 24 ) –

  number of average-expression bins (Seurat default 24).

• ctrl (int, default: 100 ) –

  control genes sampled per program gene (Seurat default 100).

• name (str, default: 'Cluster' ) –

  metadata column prefix; programs become {name}1, {name}2 … (or the dict keys when features is a dict).

• seed (int, default: 1 ) –

  RNG seed for control-gene sampling.

• search (bool, default: False ) –

  if True, resolve program genes not found verbatim by a case/punctuation-insensitive match (local UpdateSymbolList).

Returns:

• ``seurat``, with one metadata column added per program. –

Source code in truecell/module_score.py

def add_module_score(
    seurat,
    features: Union[Sequence[str], Sequence[Sequence[str]], dict],
    pool: Optional[list[str]] = None,
    nbin: int = 24,
    ctrl: int = 100,
    name: str = "Cluster",
    assay: Optional[str] = None,
    layer: str = "data",
    seed: int = 1,
    search: bool = False,
) -> "object":
    """Score one or more gene programs per cell.

    Mirrors R's AddModuleScore(): each program's score is the mean expression of
    its genes minus the mean expression of a control set drawn from the same
    average-expression bins (so highly/lowly expressed genes are controlled for).

    Parameters
    ----------
    features : a single gene list, a list of gene lists, or a name->list dict.
    pool     : genes to sample controls from (default: all features).
    nbin     : number of average-expression bins (Seurat default 24).
    ctrl     : control genes sampled per program gene (Seurat default 100).
    name     : metadata column prefix; programs become ``{name}1``, ``{name}2`` …
               (or the dict keys when ``features`` is a dict).
    seed     : RNG seed for control-gene sampling.
    search   : if True, resolve program genes not found verbatim by a
               case/punctuation-insensitive match (local ``UpdateSymbolList``).

    Returns
    -------
    ``seurat``, with one metadata column added per program.
    """
    rng = np.random.default_rng(seed)

    # Normalise `features` into (labels, list-of-lists).
    if isinstance(features, dict):
        labels = list(features.keys())
        programs = [list(v) for v in features.values()]
    elif len(features) > 0 and isinstance(features[0], (list, tuple, set)):
        programs = [list(p) for p in features]
        labels = [f"{name}{i + 1}" for i in range(len(programs))]
    else:
        programs = [list(features)]
        labels = [f"{name}1"]

    mat, feat_names = _assay_data(seurat, assay, layer)
    feat_idx = {f: i for i, f in enumerate(feat_names)}
    pool = list(pool) if pool is not None else list(feat_names)
    pool = [g for g in pool if g in feat_idx]

    # Average expression per pooled gene, then equal-frequency bins.
    pool_rows = [feat_idx[g] for g in pool]
    if sp.issparse(mat) or is_lazy(mat):
        data_avg = np.asarray(mat[pool_rows, :].mean(axis=1)).flatten()
    else:
        data_avg = np.asarray(mat)[pool_rows, :].mean(axis=1)
    # Tiny jitter breaks ties so qcut can form `nbin` equal-frequency bins.
    jitter = rng.standard_normal(len(data_avg)) / 1e30
    bins = pd.qcut(data_avg + jitter, q=min(nbin, len(pool)), labels=False, duplicates="drop")
    gene_to_bin = {g: int(b) for g, b in zip(pool, bins)}
    bin_to_genes: dict[int, list[str]] = {}
    for gene_name, bin_idx in gene_to_bin.items():
        bin_to_genes.setdefault(bin_idx, []).append(gene_name)

    n_cells = mat.shape[1]
    for label, genes in zip(labels, programs):
        used = _resolve_symbols(genes, feat_names, search)
        if not used:
            seurat.meta_data[label] = np.zeros(n_cells)
            continue

        # Control gene set: per program-gene, sample `ctrl` from its bin.
        # A dict, not a set: a mean depends on the order its terms are added,
        # and Python randomises str hashing per process, so iterating a set of
        # gene names gave a control score that differed in its last bits from
        # one run to the next. First-seen order is also R's — `AddModuleScore`
        # applies `unique()` to the sampled names and indexes the matrix with
        # the result.
        ctrl_genes: dict[str, None] = {}
        for g in used:
            b = gene_to_bin.get(g)
            if b is None:
                continue
            candidates = bin_to_genes.get(b, [])
            if not candidates:
                continue
            size = min(ctrl, len(candidates))
            picked = rng.choice(candidates, size=size, replace=False)
            ctrl_genes.update(dict.fromkeys(picked.tolist()))

        feat_scores = _mean_over_rows(mat, [feat_idx[g] for g in used])
        if ctrl_genes:
            ctrl_scores = _mean_over_rows(mat, [feat_idx[g] for g in ctrl_genes])
        else:
            ctrl_scores = np.zeros(n_cells)
        seurat.meta_data[label] = feat_scores - ctrl_scores

    return seurat

cell_cycle_scoring

cell_cycle_scoring(seurat, s_features: Optional[list[str]] = None, g2m_features: Optional[list[str]] = None, assay: Optional[str] = None, layer: str = 'data', set_ident: bool = False, nbin: int = 24, ctrl: int = 100, seed: int = 1) -> 'object'

Score S and G2/M phases and assign a discrete phase per cell.

Mirrors R's CellCycleScoring(): runs AddModuleScore for the S and G2/M gene sets, writes S.Score / G2M.Score to metadata, and assigns Phase: G1 when both scores are ≤ 0, otherwise whichever of S / G2M is larger. Defaults to the Tirosh 2016 human gene sets (CC_GENES).

If set_ident is True, the active identity is set to Phase.

Source code in truecell/module_score.py

def cell_cycle_scoring(
    seurat,
    s_features: Optional[list[str]] = None,
    g2m_features: Optional[list[str]] = None,
    assay: Optional[str] = None,
    layer: str = "data",
    set_ident: bool = False,
    nbin: int = 24,
    ctrl: int = 100,
    seed: int = 1,
) -> "object":
    """Score S and G2/M phases and assign a discrete phase per cell.

    Mirrors R's CellCycleScoring(): runs AddModuleScore for the S and G2/M gene
    sets, writes ``S.Score`` / ``G2M.Score`` to metadata, and assigns ``Phase``:
    ``G1`` when both scores are ≤ 0, otherwise whichever of S / G2M is larger.
    Defaults to the Tirosh 2016 human gene sets (``CC_GENES``).

    If ``set_ident`` is True, the active identity is set to ``Phase``.
    """
    s_features = s_features if s_features is not None else CC_GENES["s_genes"]
    g2m_features = g2m_features if g2m_features is not None else CC_GENES["g2m_genes"]

    add_module_score(
        seurat,
        features={"S.Score": s_features, "G2M.Score": g2m_features},
        nbin=nbin, ctrl=ctrl, assay=assay, layer=layer, seed=seed,
    )

    s = seurat.meta_data["S.Score"].values.astype(float)
    g2m = seurat.meta_data["G2M.Score"].values.astype(float)

    phase = np.empty(len(s), dtype=object)
    for i in range(len(s)):
        if s[i] <= 0 and g2m[i] <= 0:
            phase[i] = "G1"
        elif s[i] > g2m[i]:
            phase[i] = "S"
        else:
            phase[i] = "G2M"
    seurat.meta_data["Phase"] = phase

    if set_ident:
        seurat.idents = list(phase)
    return seurat

The bundled gene lists

CC_GENES module-attribute

CC_GENES = {'s_genes': ['MCM5', 'PCNA', 'TYMS', 'FEN1', 'MCM7', 'MCM4', 'RRM1', 'UNG', 'GINS2', 'MCM6', 'CDCA7', 'DTL', 'PRIM1', 'UHRF1', 'MLF1IP', 'HELLS', 'RFC2', 'RPA2', 'NASP', 'RAD51AP1', 'GMNN', 'WDR76', 'SLBP', 'CCNE2', 'UBR7', 'POLD3', 'MSH2', 'ATAD2', 'RAD51', 'RRM2', 'CDC45', 'CDC6', 'EXO1', 'TIPIN', 'DSCC1', 'BLM', 'CASP8AP2', 'USP1', 'CLSPN', 'POLA1', 'CHAF1B', 'BRIP1', 'E2F8'], 'g2m_genes': ['HMGB2', 'CDK1', 'NUSAP1', 'UBE2C', 'BIRC5', 'TPX2', 'TOP2A', 'NDC80', 'CKS2', 'NUF2', 'CKS1B', 'MKI67', 'TMPO', 'CENPF', 'TACC3', 'FAM64A', 'SMC4', 'CCNB2', 'CKAP2L', 'CKAP2', 'AURKB', 'BUB1', 'KIF11', 'ANP32E', 'TUBB4B', 'GTSE1', 'KIF20B', 'HJURP', 'CDCA3', 'HN1', 'CDC20', 'TTK', 'CDC25C', 'KIF2C', 'RANGAP1', 'NCAPD2', 'DLGAP5', 'CDCA2', 'CDCA8', 'ECT2', 'KIF23', 'HMMR', 'AURKA', 'PSRC1', 'ANLN', 'LBR', 'CKAP5', 'CENPE', 'CTCF', 'NEK2', 'G2E3', 'GAS2L3', 'CBX5', 'CENPA']}