Spatial trajectory inference analysis tutorial

In this tutorial, we are using both spatial information and gene expression profile to perform spatial trajectory inference to explore the progression of Ductal Carcinoma in situ (DCIS) - Invasive Ductal Cancer (IDC)

Source: https://support.10xgenomics.com/spatial-gene-expression/datasets

1. Preparation

We are trying to keep it focus on spatial trajectory inference then every step from reading data, processing to clustering, we will give the code here to easier for user to use.

Read and preprocess data

import stlearn as st
st.settings.set_figure_params(dpi=120)
# Reading data
data = st.Read10X(path="/home/d.pham/10X/BCBA/")
# Save raw_count
data.layers["raw_count"] = data.X
# Preprocessing
st.pp.filter_genes(data,min_cells=3)
st.pp.normalize_total(data)
st.pp.log1p(data)
# Keep raw data
data.raw = data
st.pp.scale(data)

Variable names are not unique. To make them unique, call `.var_names_make_unique`.
Variable names are not unique. To make them unique, call `.var_names_make_unique`.

Normalization step is finished in adata.X
Log transformation step is finished in adata.X
Scale step is finished in adata.X

Clustering data

# Run PCA
st.em.run_pca(data,n_comps=50,random_state=0)
# Tiling image
st.pp.tiling(data,out_path="tiling",crop_size = 40)
# Using Deep Learning to extract feature
st.pp.extract_feature(data)
# Apply stSME spatial-PCA option
st.spatial.morphology.adjust(data,use_data="X_pca",radius=50,method="mean")
st.pp.neighbors(data,n_neighbors=25,use_rep='X_pca_morphology',random_state=0)
st.tl.clustering.louvain(data,random_state=0)

Tiling image:   3%|▎          [ time left: 00:03 ]

PCA is done! Generated in adata.obsm['X_pca'], adata.uns['pca'] and adata.varm['PCs']

Tiling image: 100%|██████████ [ time left: 00:00 ]
Extract feature: 100%|██████████ [ time left: 00:00 ]
Adjusting data:   0%|           [ time left: 00:24 ]

The morphology feature is added to adata.obsm['X_morphology']!

Adjusting data: 100%|██████████ [ time left: 00:00 ]

The data adjusted by morphology is added to adata.obsm['X_pca_morphology']
Created k-Nearest-Neighbor graph in adata.uns['neighbors']
Applying Louvain clustering ...
Louvain clustering is done! The labels are stored in adata.obs['louvain']

st.pl.cluster_plot(data,use_label="louvain",image_alpha=1,size=7)

st.add.annotation(data,label_list=['Fatty tissue,immune/lymphoid 1 MALAT1+',
                                   'Invasive cancer,fibrous tissue 1 CXCL14+',
                                   'Invasive cancer,fibrous tissue 2 CRISP3+',
                                   'Invasive cancer,fibrous tissue, fatty tissue',
                                   'Fatty tissue,immune/lymphoid 2 IGKC+',
                                   'Fibrous tissue',
                                   'Invasive cancer,fibrous tissue (DCIS)',
                                   'Fatty tissue, Fibrous tissue',
                                   'Invasive cancer,immune/lymphoid (IDC)' ,
                                   'Invasive cancer,fatty tissue 3 MUC5B+',
                                   'Fatty tissue'],
                 use_label="louvain")
st.pl.cluster_plot(data,use_label="louvain_anno",image_alpha=1,size=7)

The annotation is added to adata.obs['louvain_anno']

2. Spatial trajectory inference

Choosing root

3733 is the index of the spot that we chose as root. It in the DCIS cluster (6). We recommend the root spot should be in the end/begin of a cluster in UMAP space. You can find min/max point of a cluster in UMAP as root.

data.uns["iroot"] = st.spatial.trajectory.set_root(data,use_label="louvain",cluster=6,use_raw=True)
st.spatial.trajectory.pseudotime(data,eps=50,use_rep="X_pca",use_label="louvain")

... storing 'feature_types' as categorical
... storing 'genome' as categorical

Created k-Nearest-Neighbor graph in adata.uns['neighbors']
All available trajectory paths are stored in adata.uns['available_paths'] with length < 4 nodes

/home/d.pham/UQ/ST/Working/stLearn/stlearn/spatials/trajectory/pseudotime.py:196: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  v = np.array(list(dic.values()))
/home/d.pham/UQ/ST/Working/stLearn/stlearn/spatials/trajectory/pseudotime.py:196: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  v = np.array(list(dic.values()))

Spatial trajectory inference - global level

We run the global level of pseudo-time-space (PSTS) method to reconstruct the spatial trajectory between cluster 6 (DCIS) and 8 (lesions IDC)

st.spatial.trajectory.pseudotimespace_global(data,use_label="louvain",list_clusters=["6","8"])

Screening:   2%|▏          [ time left: 00:08 ]

Screening PTS global graph...

Screening: 100%|██████████ [ time left: 00:00 ]
Calculating: 100%|██████████ [ time left: 00:00 ]

Calculate the graph dissimilarity using Laplacian matrix...
The optimized weighting is: 0.76
Start to construct the trajectory: 6 -> 8

st.pl.cluster_plot(data_tmp,use_label="louvain",show_trajectories=True,list_clusters=["6","8"],show_subcluster=True)

st.pl.trajectory.tree_plot(data)

Transition markers detection

Based on the spatial trajectory/tree plot, we can see 2 clades are started from sub-cluster 6 and 13. Then we run the function to detect the highly correlated genes with the PSTS values.

st.spatial.trajectory.detect_transition_markers_clades(data,clade=6,use_raw_count=False,cutoff_spearman=0.4)

Detecting the transition markers of clade_6...
Transition markers result is stored in adata.uns['clade_6']

st.spatial.trajectory.detect_transition_markers_clades(data,clade=13,use_raw_count=False,cutoff_spearman=0.4)

Detecting the transition markers of clade_13...
Transition markers result is stored in adata.uns['clade_13']

st.spatial.trajectory.detect_transition_markers_clades(data,clade=16,use_raw_count=False,cutoff_spearman=0.4)

Detecting the transition markers of clade_16...
Transition markers result is stored in adata.uns['clade_16']

data.uns['clade_6'] = data.uns['clade_6'][data.uns['clade_6']['gene'].map(lambda x: "RPL"  not in x)]
data.uns['clade_13'] = data.uns['clade_13'][data.uns['clade_13']['gene'].map(lambda x: "RPL" not in x)]
data.uns['clade_16'] = data.uns['clade_16'][data.uns['clade_16']['gene'].map(lambda x: "RPL" not in x)]
data.uns['clade_6'] = data.uns['clade_6'][data.uns['clade_6']['gene'].map(lambda x: "RPS" not in x)]
data.uns['clade_13'] = data.uns['clade_13'][data.uns['clade_13']['gene'].map(lambda x: "RPS" not in x)]
data.uns['clade_16'] = data.uns['clade_16'][data.uns['clade_16']['gene'].map(lambda x: "RPS" not in x)]

For the transition markers plot, genes from left side (red) are negatively correlated with the spatial trajectory and genes from right side (blue) are positively correlated with the spatial trajectory.

st.pl.trajectory.transition_markers_plot(data,top_genes=30,trajectory="clade_6")

st.pl.trajectory.transition_markers_plot(data,top_genes=30,trajectory="clade_13")

st.pl.trajectory.transition_markers_plot(data,top_genes=30,trajectory="clade_16")

st.spatial.trajectory.compare_transitions(data,trajectories=["clade_13","clade_16"])

The result of comparison between clade_13 and clade_16 stored in 'adata.uns['compare_result']'

st.pl.trajectory.DE_transition_plot(data)

We also provide a function to compare the transition markers between two clades.