Resolve gene dynamics by lineage

A pseudotime places every cell on one axis. When a process branches, that axis pools the branches: a gene that rises toward one fate is averaged with cells that went another way, and its trend looks weak. A fate map adds, for each cell, a probability of reaching each candidate endpoint. Here we use those probabilities as weights, so that each lineage counts every cell by how strongly it connects to that lineage’s endpoint, the approach of CellRank.

We continue the developing-pancreas example of Explore candidate cell fates with four endocrine fates, Alpha, Beta, Delta, and Epsilon, and ask three questions:

1. Which genes change along each lineage, and which does a single pseudotime miss?

2. Which genes go with commitment to each fate where the lineages part?

3. How do expression programs differ between the lineages?

Open the prepared analysis¶

Map four candidate fates¶

As in Order cell states with pseudotime, the Ductal cells share a total source weight of -1 and the endocrine populations share +1, which orients the graph. We then freeze the published cell types as sink labels and map the four endocrine fates, so every Alpha, Beta, Delta, and Epsilon cell becomes an endpoint of its own lineage.

Summarize the probabilities within each published cell type:

Endpoint cells have probability one for their own fate. Earlier cells spread their probability across the fates and lean toward Beta, the largest endocrine population, while Delta, the smallest, receives little probability upstream. These probabilities describe this graph and these endpoints, not the observed future of any cell.

Follow genes along each lineage¶

trajectory.markers correlates each gene with pseudotime over all scored cells. trajectory.lineage_markers correlates each gene with pseudotime once per lineage, weighting every cell by its probability of that lineage’s fate. Cells of the shared trunk count in every lineage, and cells drawn to another fate count little.

Each lineage’s p-values use its effective number of cells, Kish’s $(\sum w)^2 / \sum w^2$ of its weights $w$. Delta amounts to a few hundred equally weighted cells and the other lineages to well over a thousand, so its statistics are less certain. P-values also treat cells as independent: use them to rank genes, not as evidence across samples.

Compare these correlations with the pooled pseudotime for two known markers of each fate:

The pooled pseudotime barely registers Sst and even ranks Hhex as falling, because most cells late in pseudotime are Alpha and Beta cells that do not express them. Weighted by the Delta fate, both rise. Ghrl and Mboat4 behave the same way for Epsilon. A missing value marks a gene left untested in a lineage, because its expressing cells carry less than 10 cells’ worth of that lineage’s weight.

Genes shared by all endocrine cells rise along every lineage, so a lineage’s strongest correlations are often shared ones. To find genes specific to one lineage, compare each lineage with the others:

Each column now names genes that rise toward one fate more than toward the others, such as Gcg toward Alpha, Pdx1, Nkx6-1, and Slc2a2 toward Beta, Hhex toward Delta, and Ghrl toward Epsilon. These are candidates to inspect, not evidence that a gene drives its lineage.

Find genes that go with commitment¶

Pseudotime correlations follow genes along a lineage. Lineage drivers ask a different question: among cells that have not committed, which genes rise with the probability of one fate? trajectory.drivers correlates each gene with each fate probability and tests it with Fisher’s z-transform. It leaves out the endpoint cells by default, because their probabilities are fixed at zero or one and would mostly recover the markers that defined the endpoints.

The lineages part in the Ngn3-high progenitors and pre-endocrine cells, so we restrict the search to them.

Known regulators appear without being asked for. Hhex, a transcription factor that delta cells need, leads the Delta drivers with Sst and Pyy; Mafb and the secretory-granule genes Chgb and Pcsk2 lead Beta; and Arx, a determinant of alpha cells, is among the strongest Alpha drivers. The Epsilon drivers mix ribosomal genes with Cdkn1a, which marks cells leaving the cell cycle, a reminder that drivers can follow cell state as much as fate. The intervals are narrow because hundreds of cells contribute, but all of them come from one graph and one sample: a driver is a candidate regulator of commitment, not a demonstrated cause.

Compare expression programs across lineages¶

trajectory.lineage_aggregation summarizes genes along every lineage. Along pseudotime, each lineage’s fate weight is cut into 50 bins of equal weight, so each lineage runs from the shared start to its endpoint, and neighbouring bins are smoothed together. Each gene is scaled once across all lineages, which keeps its lineages comparable, and genes with similar profiles across the lineages form modules. We profile the highly variable genes of the saved analysis; the full gene set adds thousands of flat or sparse profiles.

Read each panel from the shared start on the left to the endpoint on the right. Modules that fade on every panel are progenitor programs; modules that rise on one panel only are fate programs. The strip above each panel shows the bins’ mean pseudotime: bins near an endpoint span little pseudotime, because pseudotime levels off inside the endpoints. The Delta panel is noisier than the others because its lineage carries the least fate weight.

Module numbers are labels, not stages. For individual genes, plot the profiles of every lineage together. The x-axis is each lineage’s progress, the share of its fate weight already passed, so the endpoints are as wide as the trunk:

Neurog3 peaks in the shared trunk of every lineage before the lineages part, while the hormone genes rise toward one endpoint each. Arx rises toward Alpha and Epsilon, and Pdx1 toward Beta. The Neurog3 peak sits later on the Epsilon lineage because progress counts fate weight: the Epsilon endpoint is small, so progenitors carry most of that lineage’s weight and fill more of its axis. Compare the order of events within a lineage, not positions across lineages. Profiles are scaled per gene across lineages, so the curves show where a gene is relatively high, not its absolute abundance. Pass axis="pseudotime" to plot against pseudotime instead.

Common mistakes and limitations¶

• Reading shared programs as fate-specific: the trunk counts in every lineage, so a lineage’s strongest correlations include pan-endocrine genes. Compare lineages, as above, or use drivers.

• Treating fate probabilities as observed fates: they depend on the graph, the pseudotime, and the chosen endpoints, and a missing endpoint redistributes its probability among the others. Compare plausible endpoint sets as in Trajectory validation checklist.

• Including endpoint cells in drivers: their probabilities are fixed by the model, which makes their driver scores circular.

• Over-reading p-values: they assume independent cells, and the effective number of cells accounts only for unequal weights.

Continue to Draw lineage curves and check candidate sinks to draw the lineages over a UMAP and to check candidate endpoints with a cluster tree.