空間信息在空間轉(zhuǎn)錄組中的運(yùn)用
Giotto|| 空間表達(dá)數(shù)據(jù)分析工具箱
SPOTlight || 用NMF解卷積空間表達(dá)數(shù)據(jù)
Seurat 新版教程：分析空間轉(zhuǎn)錄組數(shù)據(jù)
單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析|| scanpy教程：空間轉(zhuǎn)錄組數(shù)據(jù)分析
10X Visium:空間轉(zhuǎn)錄組樣本制備到數(shù)據(jù)分析
定量免疫浸潤在單細(xì)胞研究中的應(yīng)用

stLearn: integrating spatial location, tissue morphology and gene expression to find cell types, cell-cell interactions and spatial trajectories within undissociated tissues

空間轉(zhuǎn)錄組學(xué)是一種新興的技術(shù)向挖，它將空間信息和組織形態(tài)以及表達(dá)量融合成空間表達(dá)數(shù)據(jù)。整合這三種類型的數(shù)據(jù)極大地豐富了人們的想象力剿涮，寄希破譯細(xì)胞類型在其原生背景下的新狀態(tài)谋右。在這里我們向大家演示stLearn：一個(gè)綜合分析以上三種數(shù)據(jù)類型的python庫驻仅，stLearn首先像分析單細(xì)胞轉(zhuǎn)錄組一樣識(shí)別細(xì)胞類型，不同的是stLearn可以在空間中重建組織內(nèi)的細(xì)胞類型演變（Spatial trajectory inference ），并識(shí)別具有細(xì)胞間相互作用（Spatial cell-cell interaction）的區(qū)域丁稀。

stLearn的創(chuàng)新之處在于：

• SME normalisation ： (Spatial Morphological gene Expression normalisation应民，空間基因表達(dá)均一化)是一種組織內(nèi)的均一化策略话原，它利用鄰域信息(空間位置)和形態(tài)距離對基因表達(dá)數(shù)據(jù)進(jìn)行均一化。試圖把空間信息應(yīng)用到數(shù)據(jù)分析的全流程诲锹。

• SMEclust：stLearn使用SME歸一化數(shù)據(jù)繁仁，進(jìn)行無監(jiān)督聚類，將相似的點(diǎn)聚到聚類中归园，并根據(jù)組織中聚類的空間信息找到子類。與其他方法相比，SMEclust具有更高的空間聚類性能搀崭。

  
  

• spatial trajectory inference ： 這部分是此演示文本的重點(diǎn)粗卜，我們將在下面描述之。

• spatial cell-cell interactions： 為了研究細(xì)胞-細(xì)胞相互作用(CCI)桥爽，stLearn將細(xì)胞類型多樣性和L-R共表達(dá)結(jié)合成一個(gè)統(tǒng)計(jì)量朱灿，可用于自動(dòng)掃描整個(gè)組織切片。stLearn分別計(jì)算了L-R表達(dá)和細(xì)胞類型多樣性聚谁，最后將這兩種分析合并母剥，以發(fā)現(xiàn)對細(xì)胞-細(xì)胞相互作用有重要作用的L-R對。

下面是stLearn分析流程的框架：

切入今天我們主要討論的主題： 空間軌跡推斷。

我們曾經(jīng)表明：單細(xì)胞的一切分析 加前綴Spatial 都是一個(gè)新的分析點(diǎn)环疼。trajectory 當(dāng)然不會(huì)例外习霹。

為了發(fā)現(xiàn)組織過程——例如，回答哪個(gè)癌細(xì)胞或克隆最先出現(xiàn)炫隶，或者癌癥是如何進(jìn)化的——stLearn提供了一種稱為偽時(shí)空(pseudo-space-time, PST)軌跡分析的算法淋叶。PST是scRNAseq數(shù)據(jù)分析中常用的偽時(shí)間(pseudotime)概念的一個(gè)擴(kuò)展，它被設(shè)計(jì)用來檢測生物過程伪阶，這些生物過程可以從跨組織轉(zhuǎn)錄狀態(tài)的梯度變化來推斷煞檩。在PST中，進(jìn)化軌跡是基于組織內(nèi)細(xì)胞的空間環(huán)境和轉(zhuǎn)錄組圖譜重建的栅贴。在SMEclust分析之后斟湃，我們實(shí)現(xiàn)了PST算法來尋找組織局部層次(即單個(gè)亞群的子群之間的關(guān)系)和全局層次(即亞群之間的關(guān)系)，分別對應(yīng)亞群內(nèi)部和亞群之間的關(guān)系檐薯。

stLearn首先使用基于全組織SME 均一化基因表達(dá)數(shù)據(jù)的PAGA軌跡分析凝赛，用于發(fā)現(xiàn)亞群內(nèi)的聯(lián)系。然后坛缕，利用穩(wěn)健的軌跡推斷方法——擴(kuò)散偽時(shí)間(diffusion pseudotime 墓猎， DPT)方法計(jì)算偽時(shí)間。DPT方法使用類似隨機(jī)游走的方法來測量細(xì)胞到細(xì)胞的轉(zhuǎn)移赚楚。對于軌跡的方向毙沾，用戶可以根據(jù)組織中正在研究的生物過程來定義根節(jié)點(diǎn)。然后結(jié)合基因表達(dá)值和物理距離計(jì)算偽時(shí)空距離( pseudo-space-time distance宠页，PSTD)左胞。在此基礎(chǔ)上，構(gòu)建有向圖举户，像monocle一樣尋找分支的方法類似：利用有向最小生成樹算法對圖進(jìn)行優(yōu)化罩句，找出圖中連接節(jié)點(diǎn)最短有根樹和分支×舱可見门烂，空間軌跡推斷也是一種排序分析，只是構(gòu)建距離矩陣的方法不同兄淫，這里的距離用到了空間信息屯远。

工具安裝： https://stlearn.readthedocs.io/en/latest/installation.html
數(shù)據(jù)下載：https://support.10xgenomics.com/spatial-gene-expression/datasets/1.0.0/V1_Adult_Mouse_Brain?

安裝分析工具和下載示例數(shù)據(jù)之后我們開始做演示：

import stlearn as st
from pathlib import Path
import scanpy as sc
import matplotlib.pyplot as plt
from matplotlib import cm as cm 
import numpy as  np
import matplotlib as mpl

mpl.rcParams["font.sans-serif"] = ["SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
st.settings.set_figure_params(dpi=120)
# Reading data
data = st.Read10X(path="D:\\spatial\\V1_Adult_Mouse_Brain_filtered_feature_bc_matrix\\")


data
AnnData object with n_obs × n_vars = 2698 × 31053
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'spatial'
    obsm: 'spatial'

# 可以查看數(shù)據(jù)
data.obsm['spatial'] # 略
data.var['genome'] # 略 

可見 stlearn 用的是和scanpy一樣的AnnData數(shù)據(jù)結(jié)構(gòu)，其實(shí)stlearn的很多函數(shù)也是調(diào)用scanpy的捕虽，所以scanpy的所有分析在這里也是完全做的慨丐，如：

sc.pl.highest_expr_genes(data, n_top=20)

基本上與單細(xì)胞預(yù)處理一樣的過程，質(zhì)控泄私，均一化房揭，標(biāo)準(zhǔn)化過程备闲。

# 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)

data.layers["normal_count"] = data.X
# Keep raw data
data.raw = data
data
st.pp.scale(data)
# 再存一份scale的數(shù)據(jù)
data.layers["scale_count"] = data.X

# 也存了三分?jǐn)?shù)據(jù)：
print(data.layers['raw_count'])
print(data.layers["normal_count"])
print(data.layers['scale_count'])

降維聚類：

st.em.run_pca(data,n_comps=50,random_state=0)
st.pp.neighbors(data,n_neighbors=25,use_rep='X_pca',random_state=0)
st.tl.clustering.louvain(data,random_state=0)
sc.tl.umap(data)

data
AnnData object with n_obs × n_vars = 2698 × 18908
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow', 'louvain'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'mean', 'std'
    uns: 'spatial', 'log1p', 'pca', 'neighbors', 'louvain'
    obsm: 'spatial', 'filtered_counts', 'normalized_total', 'X_pca'
    varm: 'PCs'
    layers: 'raw_count', 'normal_count', 'scale_count'
    obsp: 'distances', 'connectivities'

有了聚類信息可以看看分群情況：

labels = data.obs['louvain'].value_counts().index
students = data.obs['louvain'].value_counts()
colors =  ["#8da0cd","#fc8d62","#66c2a5","#ff7f00"]
explode = [0.1,0.1,0.1,0.1]

plt.pie(students,
        labels= labels,
        startangle=45,
        autopct="%3.1f%%",
        shadow= True,
        pctdistance= 0.7,
        labeldistance= 1.2 )
plt.title("cluster")
plt.show()

sc.pl.umap(data, color='louvain', palette=sc.pl.palettes.default_20)

st.pl.cluster_plot(data,use_label="louvain",tissue_alpha=1,spot_size=5,show_legend=True)

只畫幾個(gè)亞群：

st.pl.cluster_plot(data,use_label="louvain",tissue_alpha=1,spot_size=5,list_cluster=[1,2,3],show_legend=True)

繪制基因表達(dá)量：

st.pl.gene_plot(data,genes=["Cnp"],use_label="louvain",use_raw_count=True)

再特定亞群的表達(dá)情況：

st.pl.gene_plot(data,genes=["Cnp"],use_label="louvain",use_raw_count=True,list_clusters=[6,7])

基本的數(shù)據(jù)探索之后，我們開始做空間軌跡推斷捅暴。注意恬砂，這里和單細(xì)胞轉(zhuǎn)錄組一樣，根節(jié)點(diǎn)的選擇也是很重要的蓬痒。

import numpy as np
data.uns["iroot"] = np.flatnonzero(data.obs["louvain"]  == str(3))[50]
st.spatial.trajectory.pseudotime(data,eps=50,use_rep="X_pca")
data

AnnData object with n_obs × n_vars = 2698 × 18908
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow', 'louvain', 'sub_cluster_labels', 'dpt_pseudotime'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'mean', 'std'
    uns: 'spatial', 'log1p', 'pca', 'neighbors', 'louvain', 'umap', 'louvain_colors', 'tmp_color', 'iroot', 'paga', 'louvain_sizes', 'diffmap_evals', 'threshold_spots', 'split_node', 'global_graph', 'centroid_dict'
    obsm: 'spatial', 'filtered_counts', 'normalized_total', 'X_pca', 'X_umap', 'X_diffmap'
    varm: 'PCs'
    layers: 'raw_count', 'normal_count', 'scale_count'
    obsp: 'distances', 'connectivities'

沒有空間信息的軌跡（擬時(shí)序）泻骤。

st.pl.non_spatial_plot(data,use_label="louvain")

st.pl.trajectory.pseudotime_plot(data,list_cluster="all",show_graph=True,node_alpha=1,tissue_alpha=1,edge_alpha=0.1,node_size=5)

st.spatial.trajectory.pseudotimespace_local函數(shù)可以構(gòu)造局部軌跡, 該算法計(jì)算子類之間的時(shí)空距離∥嗌荩可以理解為亞群內(nèi)部的軌跡（異質(zhì)性）狱掂。

st.spatial.trajectory.pseudotimespace_local(data,use_label="louvain",cluster=3)
st.pl.subcluster_plot(data,use_label="louvain",cluster=3,tissue_alpha=1)

st.pl.trajectory.local_plot(data,use_cluster=3,branch_alpha=0.2,reverse=True)

可以看出指定的亞群內(nèi)各部分間的轉(zhuǎn)移軌跡，并給出對應(yīng)的score值和方向亲轨。

欲構(gòu)建全局的軌跡可以用：

st.spatial.trajectory.pseudotimespace_global(data,use_label="louvain",list_cluster=[0,3,4,5])

Screening:   1%|           [ time left: 00:15 ]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.47
Start to construct the trajectory: 3 -> 0 -> 5 -> 4

st.pl.cluster_plot(data,use_label="louvain",show_trajectory=True,list_cluster=[1,3,4,5],show_subcluster=False)

st.pl.trajectory.tree_plot(data)

有了在空間中的軌跡（形成不同的分支）趋惨，我們肯定想知道哪些基因決定了這種軌跡，stLearn可以檢測分支間過渡的marker gene惦蚊。像差異表達(dá)分析的結(jié)果一樣會(huì)返回一個(gè)gene list希柿，這就可以走向基因調(diào)控，代謝通路之中了养筒。

st.spatial.trajectory.detect_transition_markers_clades(data,clade=2,use_raw_count=True)
st.pl.trajectory.transition_markers_plot(data,top_genes=30,trajectory="clade_2")

data
AnnData object with n_obs × n_vars = 2698 × 18908
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow', 'louvain', 'sub_cluster_labels', 'dpt_pseudotime'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'mean', 'std'
    uns: 'spatial', 'log1p', 'pca', 'neighbors', 'louvain', 'umap', 'louvain_colors', 'tmp_color', 'iroot', 'paga', 'louvain_sizes', 'diffmap_evals', 'threshold_spots', 'split_node', 'global_graph', 'centroid_dict', 'draw_graph', 'nonabs_dpt_distance_matrix', 'ST_distance_matrix', 'screening_result_local', 'PTS_graph', 'screening_result_global', 'clade_2'
    obsm: 'spatial', 'filtered_counts', 'normalized_total', 'X_pca', 'X_umap', 'X_diffmap', 'X_draw_graph_fa'
    varm: 'PCs'
    layers: 'raw_count', 'normal_count', 'scale_count'
    obsp: 'distances', 'connectivities'

data.uns['clade_2']

          gene     score       p-value
6842    Nap1l5  0.640180  1.287045e-60
6440   Tmem130  0.634474  2.987392e-59
2797     Rtl8a  0.632571  8.402284e-59
2659    Pcsk1n  0.630333  2.808432e-58
5676     Uchl1  0.629674  4.000150e-58
       ...       ...           ...
3480      Gmps  0.400526  3.161635e-21
9385    Slc1a6  0.400444  3.226337e-21
10329     Lsm6  0.400368  3.286873e-21
5349    Steap2  0.400231  3.399766e-21
13802     Gfap -0.546523  2.295220e-41

思考題：

考察下圖，有了空間信息并檢測出亞群邊界后端姚，這些邊界附近的細(xì)胞有沒有什么特殊的特征（處于邊界本身就是一種特征）晕粪，請嘗試提出一些統(tǒng)計(jì)量來描述這些特征及其擬說明的生物學(xué)問題。

stLearn: integrating spatial location, tissue morphology and gene expression to find cell types, cell-cell interactions and spatial trajectories within undissociated tissues
trajectories: Classes and Methods for TrajectoryData
PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells
Tempora: cell trajectory inference using time-series single-cell RNA sequencing data
A comparison of single-cell trajectory inference methods
Inference of multiple trajectories in single cell RNA-seq data from RNA velocity
Untangling biological factors influencing trajectory inference from single cell data

https://en.wikipedia.org/wiki/Trajectory_inference
https://www.youtube.com/watch?v=Fbd08bIv4yk
https://github.com/mdozmorov/scRNA-seq_notes
https://stlearn.readthedocs.io/en/latest/Pseudo-space-time-tutorial.html