Skip to content

Differential expression tests

  • Chi-squared tests from pseudobulk de-chisq: quick and dirty \(\chi^2\) tests on pseudobulk matrices for

    • within-dataset between-factor gene enrichment
    • between-dataset within-factor gene enrichment

  • Conditional DE after pixel/molecule level inference

    • multi-conditional-de-pixel: Multi-sample cell type-specific DE using pixel-level annotations (allowing GeoJSON region selection in each sample).
    • conditional-de-region-pixel: Cell type-specific DE between two GeoJSON-defined regions using one annotation file.

  • Conditional DE for spot/cell level data

    multi-conditional-de-pois: Multi-sample (pairwise), cell type-specific DE using Poisson regression on spot-level data.

Naive chi-squared tests from pseudobulk

de-chisq runs \(\chi^2\) tests on a pseudobulk matrix for a quick and dirty check on gene enrichment at factor level. Provide one input for a 1-vs-rest test per factor, or two inputs for a pairwise comparison of each factor between datasets.

  • 1-vs-rest per factor: 

    punkst de-chisq --input pseudobulk.tsv --out de.tsv --min-fc 1.5 --max-pval 1e-3
    By default, it also reports 1-vs-neighbors results that highlight more specific genes enriched in this factor compared to k (default 3) other most similar factors.

  • Pairwise comparison between two pseudobulk matrices: 

    punkst de-chisq --input sampleA.pseudobulk.tsv sampleB.pseudobulk.tsv \
  --out de.pair.tsv --min-fc 1.5 --max-pval 1e-4

Input Format

The input is a TSV file where the first row must be a header including an arbitrary label for the feature name column (e.g. "Gene") followed by K factor / cell type names. Then each row starts with the gene name followed by counts of that gene for each factor.

Such matrices will be generated by punkst topic-model and pixel-decode.

For pairwise input, the two matrices should have matching factor columns. Feature names are intersected between the two files.

Required Parameters

--input - Input TSV file. Provide one file for 1-vs-rest tests, or two files for pairwise comparison. --out - Output TSV file.

Optional Parameters

--min-count-per-feature - Minimum total count for a feature to be considered. Default: 100.

--max-pval - Max p-value for output. Default: 1e-3.

--min-fc - Minimum fold change. Default: 1.5.

--min-count - Minimum observed count for a (feature, factor) pair to be tested. Default: 10.

--pseudocount - Pseudocount added to each cell for the Chi-squared test. Default: 0.5.

--threads - Number of threads. Default: 1. (normally unnecessary)

--neighbor-k - Number of nearest neighbor columns to aggregate as background (single input only). Default: 3. Must be smaller than K-1. When set, writes an additional output file with 1-vs-neighbors results.

--confusion - (Pairwise mode only; experimental heuristics) Two confusion matrices to deconvolve the two input pseudobulk matrices before comparison, meant to reduce the bias due to the difference of spatial environment or contaminations (from neighboring factors) between datasets.

Output

Single input output columns:

Feature, Factor, Chi2, FoldChange, log10pval

When --neighbor-k > 0, an additional file is written with the same columns and the suffix .1vsNeighbors.tsv (if --out ends with .tsv, that extension is stripped first).

Two-input output columns:

Feature, Factor, Chi2, FoldChange, log10pval, Count1, Count2

Rows are sorted by factor (ascending) and Chi2 (descending) within each factor.

Conditional DE after pixel level inference

The pixel-level conditional DE commands test feature enrichment after pixels/molecules have been assigned to factors by pixel-decode or packaged into annotated PMTiles. (we use "pixel-level" to refer to both pixel-mode and single-molecule output from pixel-decode)

  • multi-conditional-de-pixel compares groups of datasets or dataset-region pairs.
  • conditional-de-region-pixel compares two GeoJSON-defined regions in one dataset.

Pixel-mode annotation input is joined with transcript tiles. Single-molecule annotation input and annotated point PMTiles already contain feature, count, factor, and probability information, so transcript tiles are not needed.

The test

Records are aggregated into grid units of size --grid-size, stratified by factor assignment, and tested with a Binomial model:

\(X^{(k)}_{im} \sim\) Binom \((N^{(k)}_{im}, \pi^{(k)}_{im})\) for cell type \(k\), bin \(i\) and gene \(m\), where both \(X\) and \(N\) are soft-aggregated counts from pixel level cell type assignments.

logit \((\pi^{(k)}_{im}) = a^{(k)}_m + y_i b^{(k)}_m\), where \(y_i\) is a binary indicator for the dataset (0/1). The null hypothesis is \(b^{(k)}_m = 0\).

The main p-values are computed from robust sandwich estimators of the standard errors of the estimated effect sizes. Optional permutation p-values can also be added with --perm.

Note: the intended use is when the cell type model is from external sources (e.g. apply pixel-decode with a reference-based cell type pseudobulk matrix). The interpretation is tricky if the cell types are learned from the same datasets. Either way, this is only a data exploration tool and we do not claim that it is a rigorous statistical test.

Shared Inputs And Options

Input choices:

  • Pixel-mode annotation: provide annotation files with --anno or --anno-data/--anno-index, plus transcript tiles with --pts.
  • Single-molecule annotation: provide annotation files with --anno or --anno-data/--anno-index; --pts is not needed.
  • PMTiles: provide annotated point PMTiles with --pmtiles; --pts is not used.

PMTiles input may be local files, http(s):// URLs, or unsigned/public s3://bucket/key paths when remote I/O is enabled. MLT and MVT point tiles are supported; non-point layers are ignored. For pyramid PMTiles, the default is the archive max zoom; use --pmtiles-zoom to select a specific level.

Common parameters:

--K - Number of factors in the selected annotation model. Optional for indexed annotation files whose headers store K, and optional for PMTiles that carry factor metadata. Required for older or external PMTiles without usable factor metadata. When provided, it overrides PMTiles metadata.

--features - Feature list to test. The first whitespace-delimited field on each non-comment line is used. This is preferred for PMTiles input because it avoids scanning the archive to discover features. If omitted for local deploy-cartoscope transcript PMTiles with --pmtiles-feature-field gene, the command falls back to sibling genes_bin_counts.json; otherwise it scans PMTiles rows. Optional for single-molecule annotation input with embedded feature names.

--grid-size - Grid size used to aggregate records into units.

--out - Output prefix.

Common PMTiles options:

--pmtiles-feature-field, --pmtiles-count-field - PMTiles point property names for feature and count. Defaults: feature, ct. For deploy-cartoscope transcript PMTiles, use --pmtiles-feature-field gene.

--pmtiles-factor-prefix - PMTiles K/P column prefix selecting one factor annotation group, for example celltype-pixel for columns celltype-pixel_K1, celltype-pixel_P1, etc. Optional only when the PMTiles schema contains exactly one K/P group.

Other common options:

--binary - Indicates annotation data files are binary (.bin).

--icol-x, --icol-y, --icol-feature, --icol-val - Column indices for X/Y coordinates, feature name, and count in transcript files (0-based). Defaults are 0, 1, 2, and 3; ignored for single-molecule annotation and PMTiles input.

--region-scale - Integer coordinate scale used during GeoJSON preprocessing. Default: 10.

--pseudo-rel - Relative pseudo count fraction used in the Binomial model. Default: 0.05.

--min-count-per-feature - Minimum total count for a feature to be considered. Default: 100.

--min-count - Minimum observed factor-specific count for a unit to be included. Default: 10.

--min-prob - Minimum annotation probability for a record to contribute to a factor. Default: 0.01.

--max-pval - Max p-value for output. Default: 1.

--max-pval-deconv - If at least two factors reach this p-value threshold, run deconvolution. Default: 0.05.

--min-or - Minimum odds ratio for output. Default: 1.

--perm - Number of permutations for beta calibration. Default: 0.

--min-or-perm - Minimum odds ratio for permutation testing. Default: 1.2.

--seed - Random seed for permutation.

--threads - Number of threads. Default: 1.

--aux-suff - Optional suffix inserted before .nobs.tsv and .sums.tsv.

--debug - Enable debug logging.

multi-conditional-de-pixel

multi-conditional-de-pixel compares groups of datasets. Use a contrast design file for custom comparisons or dataset-region pairs; otherwise the command builds pairwise contrasts from --labels.

Example usage with pixel-mode annotation: 

punkst multi-conditional-de-pixel \
  --anno sampleA/pixel sampleB/pixel --binary \
  --pts sampleA/transcripts.tiled sampleB/transcripts.tiled \
  --labels sampleA sampleB \
  --K 12 --features features.tsv \
  --grid-size 20 --min-count 10 --perm 5000 \
  --icol-x 0 --icol-y 1 --icol-feature 2 --icol-val 3 \
  --out de/pixel_de --threads 8 --seed 1

Example usage with explicit contrast file: 

punkst multi-conditional-de-pixel \
  --contrast contrast.tsv \
  --K 12 --features features.tsv \
  --grid-size 20 --min-count 10 --perm 5000 \
  --icol-x 0 --icol-y 1 --icol-feature 2 --icol-val 3 \
  --out de/pixel_de --threads 8 --seed 1

Example usage with annotated PMTiles: 

punkst multi-conditional-de-pixel \
  --pmtiles sampleA/genes_all.pmtiles sampleB/genes_all.pmtiles \
  --labels sampleA sampleB \
  --pmtiles-feature-field gene \
  --pmtiles-factor-prefix celltype-pixel \
  --K 25 --features features.tsv \
  --grid-size 20 --min-count 10 \
  --out de/pmtiles_de --threads 8

contrast.tsv format (tab-delimited): 

anno_prefix pts_prefix  B_vs_A  C_vs_A
sampleA/pixel   sampleA/transcripts.tiled   -1  -1
sampleB/pixel   sampleB/transcripts.tiled   1   0
sampleC/pixel   sampleC/transcripts.tiled   0   1
Each contrast column uses -1 (group 0), 1 (group 1), or 0 (exclude). Contrast names come from the header.

Region-aware contrast.tsv format: 

anno    pts region  label   confusion   B_region_vs_A_region
sampleA/pixel   sampleA/transcripts.tiled   regions/A.geojson   A_region        -1
sampleB/pixel   sampleB/transcripts.tiled   regions/B.geojson   B_region        1
sampleC/pixel   sampleC/transcripts.tiled   regions/C.geojson   C_region        0
When region is present, each row is treated as one pseudo-sample defined by the annotation dataset, transcript dataset, and GeoJSON region. The anno column may contain either an annotation prefix or an explicit .bin/.tsv annotation file. Single-molecule annotation rows can leave pts empty because the annotation file already contains feature and factor information; pixel-mode rows still require pts. The same annotation factor index must refer to the same factor across all rows. The optional label column names rows in auxiliary outputs, and the optional confusion column supplies a per-row K x K matrix. Empty confusion cells are allowed and trigger on-the-fly estimation from the region-filtered pixels.

PMTiles contrast rows use a pmtiles column instead of anno/pts.

Required Parameters

Provide either --anno, --pts (or --anno-data/--anno-index, --pts), --pmtiles, or a --contrast design file.

--anno - Prefixes of pixel annotation files (from punkst pixel-decode). For each prefix, the tool expects <prefix>.tsv (or <prefix>.bin if --binary) and <prefix>.index. Cannot be combined with --contrast.

--anno-data, --anno-index - Alternative to --anno. Provide explicit data and index files (same count for each).

--pts - Prefixes of transcript data files (from punkst pts2tiles). For each prefix, the tool expects <prefix>.tsv and <prefix>.index. Cannot be combined with --contrast.

--pmtiles - Annotated point PMTiles files. Cannot be combined with --anno/--anno-data/--anno-index outside a contrast file.

--contrast - Contrast design TSV with columns: either anno/anno_prefix plus pts/pts_prefix, or pmtiles; optional region, optional label, optional confusion; and one or more contrast columns (values -1/0/1). When provided, it replaces --anno/--anno-data/--anno-index/--pts/--pmtiles, --confusion, and --labels.

For pixel-mode annotation rows, annotation and transcript tiles must use the same tile size for each dataset pair (this is guaranteed if the pixel level decoding results are generated from the corresponding transcript tiles by punkst pixel-decode).

Optional Parameters

--labels - Labels for datasets used in pairwise output; defaults to 0..N-1. Ignored when --contrast is supplied.

--confusion - Optional per-dataset confusion matrices (K x K TSV). When omitted, the command estimates one confusion matrix from each annotation dataset on the fly.

Output Files

Given --out PREFIX, the tool generates one main output file per contrast, plus auxiliary files.

  • For each contrast, PREFIX.CONTRAST.tsv is generated, where CONTRAST is from the --contrast file header or labelA_vs_labelB for auto-pairwise.

  • Main output columns:

    • Slice: Cell type index (0-based).
    • Feature: Feature name.
    • Beta: Observed log odds ratio, i.e. logit(Pi1) - logit(Pi0).
    • log10p: -log10(p-value) for Beta.
    • Pi0, Pi1: Observed feature proportions in group 0 and group 1.
    • TotalCount: Total feature count in this cell type and contrast.
    • Beta_deconv: Deconvolved log odds ratio adjusted for cell type mixing. This is experimental and is only non-zero when deconvolution is triggered and succeeds.
    • FC_deconv: Deconvolved fold change computed from the deconvolved latent proportions.
    • log10p_deconv: -log10(p-value) for Beta_deconv.
    • p_perm: Included only when --perm > 0.

  • Auxiliary files:

    • PREFIX.nobs.tsv: The number of units and total counts per factor and group.
    • PREFIX.sums.tsv: Nonzero total feature counts per factor and group for units passing filters.
    • PREFIX.CONTRAST.da.tsv: Differential abundance chi-square test per cell type for each contrast. Count0 and Count1 are factor-specific abundance totals for the negative/reference and positive/comparison groups; FC is the positive/comparison group abundance fraction divided by the negative/reference group abundance fraction using group totals from PREFIX.nobs.tsv.

conditional-de-region-pixel

conditional-de-region-pixel compares two groups defined by --region-neg and --region-pos in one dataset. For pixel-mode input, region membership is defined by transcript coordinates. For single-molecule and PMTiles input, region membership is defined by molecule coordinates in the annotated input.

Example usage with pixel-mode annotation: 

punkst conditional-de-region-pixel \
  --anno sample/pixel --binary \
  --pts sample/transcripts.tiled \
  --region-neg regions/control.geojson \
  --region-pos regions/treatment.geojson \
  --region-label-neg control \
  --region-label-pos treatment \
  --features features.tsv \
  --grid-size 20 --min-count 10 --perm 2000 \
  --icol-x 0 --icol-y 1 --icol-feature 2 --icol-val 3 \
  --out de/region_de --threads 8

Example usage with annotated PMTiles: 

punkst conditional-de-region-pixel \
  --pmtiles sample/genes_all.pmtiles \
  --pmtiles-feature-field gene \
  --pmtiles-factor-prefix celltype-pixel \
  --K 25 --features features.tsv \
  --region-neg regions/control.geojson \
  --region-pos regions/treatment.geojson \
  --grid-size 20 --min-count 10 \
  --out de/region_pmtiles_de --threads 8

Required Parameters

Provide either --anno, both --anno-data and --anno-index, or --pmtiles.

--anno - Prefix of the annotation files. The tool expects <prefix>.tsv (or <prefix>.bin if --binary) and <prefix>.index.

--anno-data, --anno-index - Alternative to --anno. Provide the explicit annotation data and index files.

--pmtiles - Annotated point PMTiles file. MLT and MVT point tiles are supported; non-point layers are ignored. For pyramid PMTiles, the default is the archive max zoom.

--pts - Prefix of the transcript tiled dataset. Required only for pixel-mode annotation input. The tool expects <prefix>.tsv and <prefix>.index.

--region-neg, --region-pos - GeoJSON files defining the negative/reference and positive/comparison regions. See GeoJSON Region Input for the accepted file format.

For pixel-mode input, the annotation and transcript tiled inputs must use the same tile size.

Optional Parameters

--region-label-neg, --region-label-pos - Labels used in the contrast name and auxiliary outputs. Defaults: 0, 1.

Output Files

Given --out PREFIX, the command generates:

  • PREFIX.REGIONNEG_vs_REGIONPOS.tsv, where the names come from --region-label-neg and --region-label-pos.
  • PREFIX.REGIONNEG_vs_REGIONPOS.da.tsv
  • PREFIX.nobs.tsv
  • PREFIX.sums.tsv

The main output columns are the same as multi-conditional-de-pixel:

  • Slice
  • Feature
  • Beta
  • log10p
  • Pi0
  • Pi1
  • TotalCount
  • Beta_deconv
  • FC_deconv
  • log10p_deconv
  • p_perm when --perm > 0

If the two regions overlap, transcripts in the overlap contribute to both groups.

Conditional DE for spot/cell level data

multi-conditional-de-pois performs cell-type-specific differential expression on pre-defined spatial units (e.g. hexagons from tiles2hex, other binned data, single cells, or low resolution data like Visium). It leverages a pre-trained topic model (e.g., LDA from punkst topic-model) to estimate cell type mixture proportions within each unit, then for each gene fits a Poisson regression model to estimate differential expression between specified groups of datasets.

For gene \(m\) and unit \(i\), the model is: \(Y_{im} \sim \text{Pois}(\lambda_{im})\) with rate \(\lambda_{im}\) defined by a mixture across cell types: $\(\lambda_{im} = \sum_k \theta_{ik} g^{-1}( \eta_{mk}^0 + x_i b_{mk} )\)$ where \(\theta_{ik}\) is the proportion of cell type \(k\) in unit \(i\), \(x_i \in \{\pm 1\}\) indicates the group that unit \(i\) belongs to, and the null hypothesis is \(b_{mk}=0\). \(g\) is a link function with non-negative domain. The log1p link \(g^{-1}(\eta)=c_i (\exp(\eta) - 1)\) (where \(c_i=n_i/L\)) is much more efficient than the canonical log link.

Example usage with explicit contrast file: 

punkst multi-conditional-de-pois \
  --contrast contrast.tsv \
  --model model.tsv \
  --link log1p \
  --out de/unit_de --threads 8

contrast.tsv format (tab-delimited): 

in-data in-meta B_vs_A
sampleA/hex.tsv sampleA/hex.json    -1
sampleB/hex.tsv sampleB/hex.json    1
Each contrast column uses -1 (group 0), 1 (group 1), or 0 (exclude). Contrast names come from the header. The in-data column points to unit-level count data (e.g. from tiles2hex) and in-meta points to the corresponding metadata file.

Required Parameters

Provide either --in-data and --in-meta (for automatic pairwise contrasts) or a --contrast design file.

--model - Path to the trained topic model file (e.g. LDA from topic-model).

--out - Output prefix.

--contrast - Contrast design TSV with columns: in-data, in-meta, and one or more contrast columns (values -1/0/1). When provided, it replaces --in-data, --in-meta and --labels.

Optional Parameters

--in-data - Input data files (.tsv from tiles2hex).

--in-meta - Metadata files (.json from tiles2hex).

--labels - Labels for datasets used in pairwise output; defaults to 0..N-1. Ignored when --contrast is supplied.

--link - Regression link function: log or log1p. Default: log1p.

--min-count - Minimum total count per unit to be included. Default: 50.

--min-count-per-feature - Minimum total count for a feature to be tested. Default: 100.

--min-units-per-feature - Minimum number of units with a non-zero count for a feature to be tested. Default: K * 10 where K is number of cell types.

--max-pval - Max p-value for output. Default: 1.

--se-method - Method for calculating standard errors. 1: Fisher (model-based), 2: Sandwich (robust), 3: both. Default: 1.

--size-factor - Size factor \(L\) for the Poisson regression model. Default: 10000.0.

--threads - Number of threads. Default: 1.

Output Files

Given --out PREFIX, the tool generates one main output file per contrast, plus an auxiliary file.

  1. For each contrast, PREFIX.CONTRAST.tsv is generated, where CONTRAST is from the --contrast file header or labelA_vs_labelB for auto-pairwise. Columns:

    • Feature: Feature name.
    • Factor: Cell type name.
    • Beta: log fold change for the feature in this cell type.
    • SE: Standard error of Beta. If se-method=3, an additional SE0 column with the Fisher-based SE is included.
    • log10p: -log10(p-value) for Beta.
    • Dist2bd: Distance of the beta estimate to the optimization boundary. Small values may indicate unstable estimates.
    • n: Number of units included in the test for this feature.
    • Count1, Count2: Total counts of the feature in group 1 and group 2.

  2. An auxiliary file PREFIX.CONTRAST.eta0.tsv is also generated, containing the fitted baseline (log) expression rates for each gene and cell type.