GSET 🌉

GSET (Gene Signals from Embedding Trajectories) is a Python package that selects a compact set of genes from time-resolved single-cell data. The selected genes are meant as the gene context of a mechanistic gene regulatory network (GRN) inference method such as CardamomOT, whose cost grows quickly with the number of genes.

The method follows every gene along two embedding trajectories. A co-expression embedding (a graph neural network on the co-expression graph of each timepoint) shows the changes of the co-expression neighbourhood of a gene, measured by a KL divergence. An embedding of the expression distribution (a DeepSets encoder) shows the changes of the activity of a gene, measured by an optimal transport (OT) residual, which is the displacement of a gene after the global movement of all genes is removed. The genes are grouped into programs, and lagged comparisons of the two signals between consecutive programs retrieve the bridge genes that link them.

Main functionalities:

1. Gene programs from the temporal co-expression embeddings;
2. Bridge genes between consecutive programs, from the KL signal without lag and the OT signal with lag 1;
3. A compact gene set of a given size: the bridge genes and the core genes of the programs.

Installation

pip install gset

Python 3.11 or newer. The package depends on numpy, scipy, pandas, scikit-learn, torch, torch_geometric and pot. anndata is needed to read an AnnData object.

Basic usage

import anndata as ad
from gset import select_genes

adata = ad.read_h5ad("my_data.h5ad")    # raw counts in X, the time of every cell in adata.obs["time"]
genes = adata.var_names[select_genes(adata, budget=30, n_programs=3)]

select_genes returns the indices of the selected genes, sorted. By default a third of the budget goes to the core genes of the programs and two thirds to the bridge genes, and the programs are found by clustering, with their number given.

About the data

The data are single-cell snapshots of the same genes at several timepoints. The cells collected at the same time form a snapshot, and successive snapshots are made of different cells, because the measurement destroys the cells. The cells of different timepoints are therefore not matched. The input is a table of raw counts (cells x genes) with the time of every cell, or a list of count matrices, one per timepoint in the order of the time. With several replicates of the same experiment, give a list with one such list per replicate.

Settings

Programs are chosen in one of three ways: n_programs (clustering), query_genes, or labels and order.

# the numbers of cores and bridges (the budget is their sum), and of the two channels of the bridges
genes = select_genes(adata, n_cores=10, n_bridges=20, n_kl=14, n_ot=6, n_programs=3)

# programs built around one query gene each, in the order in which the programs become active
queries = [adata.var_names.get_loc("Lox"), adata.var_names.get_loc("Birc5")]
genes = select_genes(adata, budget=30, query_genes=queries, program_size=10)

# programs of your own: the program of every gene (0, 1, 2, ...) and the order of the programs
genes = select_genes(adata, budget=30, labels=labels, order=[0, 1, 2])

# the data as a list of count matrices (cells x genes), one per timepoint
genes = select_genes([X_t0, X_t1, X_t2], budget=30, n_programs=3)

All the settings are described in help(select_genes).

Encoders

The co-expression embedding has two encoders. The default one ("simulation") takes the correlation profile of a gene as its features and is meant for panels of tens to a few hundred genes. The other one ("real_data") takes the leading singular vectors of the correlation matrix as features and adds self-loops to the graph, for thousands of genes:

genes = select_genes(adata, budget=30, n_programs=3, coexpr_encoder="real_data")

The expression embedding (DeepSets) is the same in both cases.

Gene Program Explorer

The genes related to a query gene, at the static level and with the lagged signals, can also be explored interactively in the Gene Program Explorer: https://rmsexplorer.org

Reproducibility

The training of the neural networks gives slightly different numbers from one run to the next when it uses several threads. To get the same genes every time, call make_reproducible() once before select_genes:

from gset import make_reproducible, select_genes

make_reproducible()