trajectory inference

For more details on scFates, please see the documentation

Environment setup¶

scFates requires R (for the mgcv package used internally) and rpy2 to bridge Python and R. The easiest way to get everything working is with conda.

Local (conda): Create the environment from the provided scFates_environment.yml:

conda env create -f scFates_environment.yml
conda activate scFates

Pip only: Use scFates_requirements.txt:

pip install -r scFates_requirements.txt

Note: The pip-only route requires that you have R installed separately with the mgcv package available, and that R_HOME is set correctly.

Google Colab: Run the cell below — it will detect Colab and install the required packages automatically.

In [48]:

import sys
if 'google.colab' in sys.modules:
    print("Colab detected — installing required packages...")
    !pip install -q scanpy==1.10.4 anndata==0.11.3 scFates==1.0.9 rpy2==3.5.11 leidenalg==0.10.2 igraph==0.11.8
else:
    print("Local environment detected — skipping install.")

import sys if 'google.colab' in sys.modules: print("Colab detected — installing required packages...") !pip install -q scanpy==1.10.4 anndata==0.11.3 scFates==1.0.9 rpy2==3.5.11 leidenalg==0.10.2 igraph==0.11.8 else: print("Local environment detected — skipping install.")

Local environment detected — skipping install.

Point interpreter to R, import packages, fetch data.

In [49]:

import os, sys
os.environ['R_HOME'] = sys.exec_prefix+"/lib/R/"

import os, sys os.environ['R_HOME'] = sys.exec_prefix+"/lib/R/"

In [50]:

print(os.environ['R_HOME'])

print(os.environ['R_HOME'])

/opt/homebrew/Caskroom/miniforge/base/envs/scFates/lib/R/

In [51]:

import scanpy as sc
import scFates as scf
import numpy as np
import pandas as pd

import scanpy as sc import scFates as scf import numpy as np import pandas as pd

In [52]:

!wget https://compstemcellbio.s3.amazonaws.com/2024/TI/adata_Stam_250_022224.h5ad

!wget https://compstemcellbio.s3.amazonaws.com/2024/TI/adata_Stam_250_022224.h5ad

--2026-03-11 23:13:15--  https://compstemcellbio.s3.amazonaws.com/2024/TI/adata_Stam_250_022224.h5ad
54.231.136.65, 16.15.185.77, 52.217.225.97, ...pstemcellbio.s3.amazonaws.com)... 
Connecting to compstemcellbio.s3.amazonaws.com (compstemcellbio.s3.amazonaws.com)|54.231.136.65|:443... connected.
200 OKequest sent, awaiting response... 
Length: 30948532 (30M) [binary/octet-stream]
Saving to: ‘adata_Stam_250_022224.h5ad.1’

adata_Stam_250_0222 100%[===================>]  29.51M  8.98MB/s    in 3.3s    

2026-03-11 23:13:19 (8.84 MB/s) - ‘adata_Stam_250_022224.h5ad.1’ saved [30948532/30948532]

In [53]:

adata = sc.read_h5ad("adata_Stam_250_022224.h5ad")

adata = sc.read_h5ad("adata_Stam_250_022224.h5ad")

In [54]:

adata

adata

Out[54]:

AnnData object with n_obs × n_vars = 2750 × 36601
    obs: 'day', 'sample', 'experiment', 'n_counts', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand'

In [55]:

adata.obs['sample'].value_counts()

adata.obs['sample'].value_counts()

Out[55]:

sample
cd34.day0    250
ery.day2     250
ery.day4     250
ery.day6     250
ery.day8     250
ery.day11    250
mk.day2      250
mk.day4      250
mk.day6      250
mk.day8      250
mk.day11     250
Name: count, dtype: int64

Note that this partcular dataset is already normalized (your HW data will not be). Now find HVG

In [56]:

sc.pp.highly_variable_genes(adata , min_mean=0.0125, max_mean=6, min_disp=0.25)

sc.pp.highly_variable_genes(adata , min_mean=0.0125, max_mean=6, min_disp=0.25)

Perform PCA, plot embedding, and color with marker genes for CD34 HSPCs, Erythrocytes, and Megakaryocytes (Human)

In [57]:

sc.tl.pca(adata, mask_var='highly_variable')
sc.pl.pca(adata, color=["CD34", "HBA1", "PF4"], s=20)

sc.tl.pca(adata, mask_var='highly_variable') sc.pl.pca(adata, color=["CD34", "HBA1", "PF4"], s=20)

In [58]:

sc.pl.pca(adata, color='sample')

sc.pl.pca(adata, color='sample')

cf.tl.tree() fits a principal tree to this embedding.

See documentation for parameters

In [59]:

adata

adata

Out[59]:

AnnData object with n_obs × n_vars = 2750 × 36601
    obs: 'day', 'sample', 'experiment', 'n_counts', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'hvg', 'pca', 'sample_colors'
    obsm: 'X_pca'
    varm: 'PCs'

In [60]:

adata.obsm['X_pca'].shape

adata.obsm['X_pca'].shape

Out[60]:

(2750, 50)

In [61]:

scf.tl.tree(adata,method="ppt",Nodes=10,use_rep="X_pca")

scf.tl.tree(adata,method="ppt",Nodes=10,use_rep="X_pca")

inferring a principal tree --> parameters used 
    10 principal points, sigma = 0.1, lambda = 1, metric = euclidean
    fitting:   0%|                                       | 0/50 [00:00<?, ?it/s]
    converged
    finished (0:00:00) --> added 
    .uns['ppt'], dictionnary containing inferred tree.
    .obsm['X_R'] soft assignment of cells to principal points.
    .uns['graph']['B'] adjacency matrix of the principal points.
    .uns['graph']['F'] coordinates of principal points in representation space.

Plot the tree on top of the pca embedding with scf.pl.graph()

In [62]:

scf.pl.graph(adata, basis='pca')

scf.pl.graph(adata, basis='pca')

Why does the trajectory bend back (top left)?

In [63]:

sc.pl.pca(adata, color=["sample"],projection="3d", components = ['1,2,3', '1,3,2'], s=20)

sc.pl.pca(adata, color=["sample"],projection="3d", components = ['1,2,3', '1,3,2'], s=20)

We should restrict tree finding to 1st 3 PCs

In [64]:

scf.tl.tree(adata, method="ppt", Nodes=10, use_rep="X_pca", ndims_rep = 3)

scf.tl.tree(adata, method="ppt", Nodes=10, use_rep="X_pca", ndims_rep = 3)

inferring a principal tree --> parameters used 
    10 principal points, sigma = 0.1, lambda = 1, metric = euclidean
    fitting:   0%|                                       | 0/50 [00:00<?, ?it/s]
    converged
    finished (0:00:00) --> added 
    .uns['ppt'], dictionnary containing inferred tree.
    .obsm['X_R'] soft assignment of cells to principal points.
    .uns['graph']['B'] adjacency matrix of the principal points.
    .uns['graph']['F'] coordinates of principal points in representation space.

In [65]:

scf.pl.graph(adata, basis='pca')

scf.pl.graph(adata, basis='pca')

Explore other parameters:

• nodes
• ppt_lambda
• ppt_sigma

In [66]:

scf.tl.tree(adata,method="ppt",Nodes=50,use_rep="X_pca", ndims_rep = 3)
scf.pl.graph(adata, basis='pca')

scf.tl.tree(adata,method="ppt",Nodes=50,use_rep="X_pca", ndims_rep = 3) scf.pl.graph(adata, basis='pca')

inferring a principal tree --> parameters used 
    50 principal points, sigma = 0.1, lambda = 1, metric = euclidean
    fitting:  40%|███████████▌                 | 20/50 [00:00<00:00, 115.97it/s]
    converged
    finished (0:00:00) --> added 
    .uns['ppt'], dictionnary containing inferred tree.
    .obsm['X_R'] soft assignment of cells to principal points.
    .uns['graph']['B'] adjacency matrix of the principal points.
    .uns['graph']['F'] coordinates of principal points in representation space.

In [67]:

scf.tl.tree(adata,method="ppt",Nodes=100,use_rep="X_pca", ppt_lambda=100,ppt_sigma=0.75,ppt_nsteps=100, ndims_rep = 3)
scf.pl.graph(adata, basis='pca')

scf.tl.tree(adata,method="ppt",Nodes=100,use_rep="X_pca", ppt_lambda=100,ppt_sigma=0.75,ppt_nsteps=100, ndims_rep = 3) scf.pl.graph(adata, basis='pca')

inferring a principal tree --> parameters used 
    100 principal points, sigma = 0.75, lambda = 100, metric = euclidean
    fitting:  30%|████████▋                    | 30/100 [00:00<00:00, 99.91it/s]
    converged
    finished (0:00:00) --> added 
    .uns['ppt'], dictionnary containing inferred tree.
    .obsm['X_R'] soft assignment of cells to principal points.
    .uns['graph']['B'] adjacency matrix of the principal points.
    .uns['graph']['F'] coordinates of principal points in representation space.

Set a root node with scf.tl.root()

Then compute Pseudotime with scf.tl.pseudotime()

In [68]:

scf.tl.root(adata, 96)

scf.tl.root(adata, 96)

node 96 selected as a root --> added
    .uns['graph']['root'] selected root.
    .uns['graph']['pp_info'] for each PP, its distance vs root and segment assignment.
    .uns['graph']['pp_seg'] segments network information.

tl.pseudotime also assigns cells to:

obs. "seg", "milestones"

plot PCA colored by ^^

In [69]:

scf.tl.pseudotime(adata,seed=42)

scf.tl.pseudotime(adata,seed=42)

projecting cells onto the principal graph
    finished (0:00:00) --> added
    .obs['edge'] assigned edge.
    .obs['t'] pseudotime value.
    .obs['seg'] segment of the tree assigned.
    .obs['milestones'] milestone assigned.
    .uns['pseudotime_list'] list of cell projection from all mappings.

In [70]:

sc.pl.pca(adata, color=["seg", "milestones"])

sc.pl.pca(adata, color=["seg", "milestones"])

In [71]:

sc.pl.pca(adata, color=["seg", "milestones"],projection="3d", components = ['1,2,3', '1,3,2'], s=20)

sc.pl.pca(adata, color=["seg", "milestones"],projection="3d", components = ['1,2,3', '1,3,2'], s=20)

You can overlay trajectory, PT, on the embedding with scf.pl.trajectory()

In [72]:

scf.pl.trajectory(adata, basis='pca')

scf.pl.trajectory(adata, basis='pca')

And you can use scany plotting to show pseudotime for each cell

In [73]:

sc.pl.pca(adata, color='t')
sc.pl.pca(adata, color=["t"],projection="3d", components = ['1,2,3', '1,3,2'], s=20)

sc.pl.pca(adata, color='t') sc.pl.pca(adata, color=["t"],projection="3d", components = ['1,2,3', '1,3,2'], s=20)

cluster cells into a real tree with scf.tl.dendrogram() and plot with scf.pl.dendrogram()

In [74]:

scf.tl.dendrogram(adata)

scf.tl.dendrogram(adata)

Generating dendrogram of tree
    segment : 100%|███████████████████████████████| 3/3 [00:00<00:00,  3.22it/s]
    finished (0:00:00) --> added 
    .obsm['X_dendro'], new embedding generated.
    .uns['dendro_segments'] tree segments used for plotting.

In [75]:

scf.pl.dendrogram(adata, color=['HBA1'])

scf.pl.dendrogram(adata, color=['HBA1'])

Test for genes whose expresion is associated with PT and branching

First, subset the data so that this is quick

In [76]:

ad2 = adata.copy()
ad2 = ad2[:,ad2.var['highly_variable']]
# to speed up execution of test_association, sub-sample the data as below
# ad2 = ad2[:,np.random.choice(ad2.n_vars, 500, replace=False)]

# you can also specify the number of jobs to run in parallel if you system supports this with the n_jobs parameter

ad2 = adata.copy() ad2 = ad2[:,ad2.var['highly_variable']] # to speed up execution of test_association, sub-sample the data as below # ad2 = ad2[:,np.random.choice(ad2.n_vars, 500, replace=False)] # you can also specify the number of jobs to run in parallel if you system supports this with the n_jobs parameter

Now perform the test

In [77]:

scf.tl.test_association(ad2, fdr_cut=0.25, A_cut=.4, n_jobs = 20)

scf.tl.test_association(ad2, fdr_cut=0.25, A_cut=.4, n_jobs = 20)

test features for association with the trajectory
    single mapping : 100%|██████████████████| 4639/4639 [00:47<00:00, 97.59it/s]
    found 699 significant features (0:00:47) --> added
    .var['p_val'] values from statistical test.
    .var['fdr'] corrected values from multiple testing.
    .var['st'] proportion of mapping in which feature is significant.
    .var['A'] amplitue of change of tested feature.
    .var['signi'] feature is significantly changing along pseudotime.
    .uns['stat_assoc_list'] list of fitted features on the graph for all mappings.

/opt/homebrew/Caskroom/miniforge/base/envs/scFates/lib/python3.12/site-packages/scFates/tools/test_association.py:413: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.
  adata.var[c] = stat_assoc[c]

In [78]:

# ad2.var['signi'].value_counts()

# ad2.var['signi'].value_counts()

In [79]:

ad3 = ad2.copy()
scf.tl.fit(ad3, n_jobs=20) # fits genes to trajectory
ad3.var

ad3 = ad2.copy() scf.tl.fit(ad3, n_jobs=20) # fits genes to trajectory ad3.var

fit features associated with the trajectory
    single mapping : 100%|████████████████████| 699/699 [00:14<00:00, 48.70it/s]
    finished (adata subsetted to keep only fitted features!) (0:00:14) --> added
    .layers['fitted'], fitted features on the trajectory for all mappings.
    .raw, unfiltered data.

Out[79]:

	Accession	Chromosome	End	Start	Strand	highly_variable	means	dispersions	dispersions_norm	p_val	A	fdr	st	signi
HP1BP3	ENSG00000127483	1	20787323	20742679	-	True	0.769912	1.865717	0.448787	1.181477e-16	0.499296	5.480873e-13	1	True
EIF4G3	ENSG00000075151	1	21176888	20806292	-	True	0.450619	1.914048	0.730141	1.092165e-93	0.827920	5.066555e-90	1	True
HNRNPR	ENSG00000125944	1	23344336	23303771	-	True	0.582793	1.787895	0.273731	3.026882e-40	0.529806	1.404171e-36	1	True
ID3	ENSG00000117318	1	23559501	23557926	-	True	0.331609	2.872770	4.198693	4.870580e-212	1.142615	2.259462e-208	1	True
RHCE	ENSG00000188672	1	25430192	25362249	-	True	0.428462	1.796913	0.306358	0.000000e+00	1.337406	0.000000e+00	1	True
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
PGRMC1	ENSG00000101856	X	119244466	119236245	+	True	0.555356	1.794529	0.297733	1.834102e-19	0.440285	8.508400e-16	1	True
RAB33A	ENSG00000134594	X	130184870	130171962	+	True	0.470637	1.811213	0.358095	1.758523e-75	0.713661	8.157787e-72	1	True
PHF6	ENSG00000156531	X	134428791	134373253	+	True	0.425871	1.936639	0.811870	8.659728e-20	0.405178	4.017248e-16	1	True
FHL1	ENSG00000022267	X	136211359	136146702	+	True	0.591447	2.977816	4.578740	2.898772e-143	1.173347	1.344740e-139	1	True
ATP6AP1	ENSG00000071553	X	154436516	154428645	+	True	0.499192	1.795571	0.301502	1.386872e-26	0.513673	6.433700e-23	1	True

699 rows × 14 columns

In [80]:

ad3.var['st'].value_counts()

ad3.var['st'].value_counts()

Out[80]:

st
1    699
Name: count, dtype: int64

In [81]:

scf.pl.single_trend(ad3,"HBA1",basis="pca")

scf.pl.single_trend(ad3,"HBA1",basis="pca")

In [82]:

scf.pl.single_trend(ad3,"ID3",basis="pca")

scf.pl.single_trend(ad3,"ID3",basis="pca")

In [83]:

scf.pl.single_trend(ad3,"PF4",basis="pca")

scf.pl.single_trend(ad3,"PF4",basis="pca")

In [84]:

scf.pl.single_trend(ad3,"PF4",basis="dendro")

scf.pl.single_trend(ad3,"PF4",basis="dendro")

In [85]:

scf.pl.single_trend(ad3,"PF4",basis="dendro")
scf.pl.single_trend(ad3,"HBA1",basis="dendro")

scf.pl.single_trend(ad3,"PF4",basis="dendro") scf.pl.single_trend(ad3,"HBA1",basis="dendro")

In [86]:

scf.tl.test_fork(ad3,root_milestone="96",milestones=["66","92"],n_jobs=20,rescale=True)

scf.tl.test_fork(ad3,root_milestone="96",milestones=["66","92"],n_jobs=20,rescale=True)

testing fork
    single mapping
    Differential expression: 100%|████████████| 699/699 [00:13<00:00, 53.42it/s]
    test for upregulation for each leave vs root
    upreg 66: 100%|██████████████████████████| 141/141 [00:00<00:00, 949.07it/s]
    upreg 92: 100%|█████████████████████████| 558/558 [00:00<00:00, 1709.37it/s]
    finished (0:00:13) --> added 
    .uns['96->66<>92']['fork'], DataFrame with fork test results.

In [87]:

scf.tl.branch_specific(ad3,root_milestone="96",milestones=["66","92"],effect=0.25)

scf.tl.branch_specific(ad3,root_milestone="96",milestones=["66","92"],effect=0.25)

    branch specific features: 92: 276, 66: 48
    finished --> updated 
    .uns['96->66<>92']['fork'], DataFrame updated with additionnal 'branch' column.

In [88]:

g1=scf.pl.trends(ad3,
                 root_milestone="96",
                 milestones=["66", "92"],
                 branch="66",
                 plot_emb=False,ordering="max",return_genes=True, n_features=12)

g1=scf.pl.trends(ad3, root_milestone="96", milestones=["66", "92"], branch="66", plot_emb=False,ordering="max",return_genes=True, n_features=12)

In [89]:

g2=scf.pl.trends(ad3,
                 root_milestone="96",
                 milestones=["66", "92"],
                 branch="92",
                 plot_emb=False,ordering="max",return_genes=True, n_features=12)

g2=scf.pl.trends(ad3, root_milestone="96", milestones=["66", "92"], branch="92", plot_emb=False,ordering="max",return_genes=True, n_features=12)

In [90]:

gg=g1.tolist()+g2.tolist()

gg=g1.tolist()+g2.tolist()

In [91]:

scf.pl.trends(ad3,gg,figsize=(4,4),annot="milestones",n_features=12,
        plot_emb=False,ordering=None,return_genes=True)

scf.pl.trends(ad3,gg,figsize=(4,4),annot="milestones",n_features=12, plot_emb=False,ordering=None,return_genes=True)

Out[91]:

Index(['APOE', 'GSTM3', 'EZR', 'HMGB2', 'APOC1', 'NMU', 'MT1X', 'S100A6',
       'FAM178B', 'MT1G',
       ...
       'SPARC', 'TNFAIP8', 'MTSS1', 'FRMD4B', 'GLA', 'NT5C3A', 'TUBB1',
       'ENDOD1', 'SQSTM1', 'ITM2B'],
      dtype='object', length=324)

In [92]:

scf.pl.matrix(ad3,gg,norm="minmax",cmap="RdBu_r",colorbar=False)

scf.pl.matrix(ad3,gg,norm="minmax",cmap="RdBu_r",colorbar=False)

In [93]:

scf.tl.cluster(ad3,n_neighbors=10,metric="correlation")

scf.tl.cluster(ad3,n_neighbors=10,metric="correlation")

Clustering features using fitted layer
    finished (0:00:00) --> added 
    .var['clusters'] identified modules.

In [94]:

ad3.var.clusters.unique()

ad3.var.clusters.unique()

Out[94]:

['5', '1', '4', '8', '0', '3', '6', '7', '9', '2']
Categories (10, object): ['0', '1', '2', '3', ..., '6', '7', '8', '9']

In [95]:

for c in ad3.var["clusters"].unique():
    scf.pl.trends(ad3,features=ad3.var_names[ad3.var.clusters==c],basis="pca")

for c in ad3.var["clusters"].unique(): scf.pl.trends(ad3,features=ad3.var_names[ad3.var.clusters==c],basis="pca")

In [ ]: