cosg 1.2.0

COSG v1.2.0 - analytic significance for the COSG specificity score

on Python PyPI

COSG now reports significance alongside specificity. Marker tables carry pvals, pvals_adj, zscores and neg_log10_pvals when you ask for them; the default is off and the output is byte-for-byte unchanged when you do not.

cosg.cosg(adata, groupby='CellTypes',
          calculate_pvalues=True,        # default False
          pvalue_method='spa',           # 'normal' | 'exact' | 'permutation'
          pvalue_fdr_method='fdr_bh')    # or 'fdr_by'

adata.uns['cosg']['pvals']             # raw p
adata.uns['cosg']['pvals_adj']         # Benjamini-Hochberg within each group
adata.uns['cosg']['neg_log10_pvals']   # -log10 p, no floor
adata.uns['cosg']['zscores']           # standardised effect

This closes issue #6.

What is being tested

For gene g and a group c of n_c cells, COSG's raw cosine is T / (||x|| * sqrt(n_c)) where T is the gene's sum over the group. Under the exchangeability null — this gene's expression is independent of the grouping — ||x|| and n_c are fixed, so testing the cosine is testing T, and T is a simple random sample sum without replacement from a fixed population. That has exact moments, so no permutations are sampled and there is no p-value floor.

The p-value calibrates the raw cosine, so it is identical across mu by design: mu says how harshly to discount a gene that is also high elsewhere, which is a preference about ranking genes that are genuinely associated. The null contains no mu, and chance has one answer per gene per group.

Why not the normal approximation

Because it fails exactly where markers are decided. Its accuracy is governed by the expected number of expressing cells inside the group, not by the group size or the gene's global expressing count. Measured against an exact enumeration of the permutation null (20,000 cells, 80,000 gene × group entries under random labels):

A sparse, specific gene is not a corner case; it is the interesting kind, and BH at FDR 0.01 over 20,000 genes engages exactly those thresholds.

How it is computed

Three regimes, chosen internally from the data:

• the normal approximation where the tail cannot change a decision;
• a conditional (double) saddlepoint once the p is small — a saddlepoint has relative error, which is what a tail needs, where an Edgeworth correction has absolute error and still misses by an order of magnitude;
• exact enumeration for genes expressed in very few cells.

The saddlepoint reads a gene only through sums over its value multiset, so each gene is summarised once as (value, multiplicity) atoms on a fixed log grid — a few tens of atoms whether the gene is expressed in 60 cells or 10,000. The compressed tail agrees with the uncompressed one to 0.1% or better on normalised data and exactly on counts.

On integer data (raw or count-split) the sum lives on a lattice, and the continuous tail formula approximates a step function from the wrong side. The tail is evaluated at t - 1/2, which brings it from 0.81–0.93 to 0.94–1.00 times an exact enumeration. Integrality is detected from the data; you do not declare it.

neg_log10_pvals is computed in log space. It matters: on a clean cell type a third of the top 200 markers reach the float64 floor (~1e-308), where every pvals entry reads as the same number. On a test fixture where pvals is exactly zero, neg_log10_pvals reads 751 and 1804 for two markers the p-value cannot separate at all.

Cost

Human MTG, 30,192 genes, 24 subclasses, each arm in its own process:

Peak memory is the same as the moments pass any p-value needs — the atom summary costs nothing measurable. Both the in-memory and the streaming (cytome) path compute the same summary and agree gene by gene.

One caveat worth reading before you use these

If the group labels came from clustering the same expression matrix, the p-values are anti-conservative — the labels were chosen to separate the data they are now tested against (post-selection inference; Gao, Bien & Witten 2022). This is a property of the design, not of the approximation: sampling permutations instead would not fix it, and it applies to every marker test computed this way, not to COSG specifically.

Valid uses are labels carrying independent information — curated annotation, another modality, a reference mapping — or count splitting (Neufeld et al. 2022): derive the labels on one split and pass the other to layer=.

Also in this release

pvalue_method='permutation' samples the null and is the validation oracle, not a shipping default; its resolution is bounded by 1/(n_permutations+1). batch_key stratifies the null (labels are permuted within batch) on the in-memory path; the streaming path raises rather than silently computing an unstratified null.

This release also carries the accumulated 1.1.2 and 1.1.3 work.

Documentation

Five worked tutorials, each executed end to end on real data:

• Marker genes and their significance — the columns above, IQR normalisation for comparing scores across cell types, the dendrogram variants, and the double-dipping caveat demonstrated rather than asserted.
• COSG on a cytome — the same markers streamed from disk, and the four output_format shapes.
• COSG across batches — batch_key measured: 86% of top-20 markers hold, and the ones that move name your most dissociation-sensitive cell types.
• COSG on the GPU — device='gpu' across five matrix sizes.
• COSG on spatial data — organ markers on a whole embryo section, plotted back into tissue space.

Installation

pip install --upgrade cosg

Citation

Dai M., Pei X., Wang X.-J. Accurate and fast cell marker gene identification with COSG. Briefings in Bioinformatics (2022) 23, bbab579.