什么是tangram？

tangram是一種映射單細(xì)胞表達(dá)譜數(shù)據(jù)到空間表達(dá)譜的方法，它收集同一解剖區(qū)域或者組織類型的單細(xì)胞數(shù)據(jù)和空間數(shù)據(jù)。通過整合，該方法在scRNAseq表達(dá)譜匹配的基礎(chǔ)上創(chuàng)建一組新的空間轉(zhuǎn)錄組數(shù)據(jù)，并可投射scRNAseq的注釋(e.g.細(xì)胞類型、模塊）到不同的空間區(qū)域。

使用Tangram可以做什么？

Tangram最常用的功能是在空間中分辨細(xì)胞類型，并對空間轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行校正：因為相較于空間轉(zhuǎn)錄組(Visium 或 Slide-seq)，scRNA-seq數(shù)據(jù)較少出現(xiàn)dropout. Tangram可以生成一組包含更多基因的空間轉(zhuǎn)錄譜，因此，我們可以在空間上觀察轉(zhuǎn)錄模塊(program)的表達(dá)，以此來發(fā)現(xiàn)配受體互作和細(xì)胞通訊連接機(jī)制。在細(xì)胞分割允許的情況下，Tangram還可以用來對空間數(shù)據(jù)去卷積。如果碰到多模態(tài)的單細(xì)胞數(shù)據(jù)，Tangram還可以在空間上可視化其他模態(tài)數(shù)據(jù)，比較染色質(zhì)開放性(chromatin accessibility).

Tangram和其他空間去卷積/映射方法的不同？

驗證性。其他方法通過已知的細(xì)胞類型分布模式或者頻率來驗證映射的合理性，做的都是sanity checks(合理性檢查），但如用于探索性研究將不再有用。而tangram的好處是映射結(jié)果可以通過holdout genes(預(yù)留出的基因/test transcriptome）來進(jìn)行檢驗。

當(dāng)scRNAseq和spatial data 來自不同的樣本，可以用tangram嗎?

可行的。發(fā)明者對映射方法進(jìn)行了很巧妙的調(diào)整，可以對某一細(xì)胞類型的average cells進(jìn)行映射，而不是只映射單個細(xì)胞，這種平均的方式可以消除不同樣本帶來的生物信號變化。記住，這需要再scRNA-seq前進(jìn)行注釋來犧牲單細(xì)胞的高分辨度，并對函數(shù)傳入mode = cluster.

不多說了，上代碼

import squidpy as sq
import numpy as np
import pandas as pd
from anndata import AnnData
import pathlib
import matplotlib.pyplot as plt
import matplotlib as mpl
import skimage
import seaborn as sns
import tangram as tg

sc.logging.print_header()
print(f"squidpy=={sq.__version__}")

%load_ext autoreload
%autoreload 2
%matplotlib inline

加載數(shù)據(jù)

• Squidy加載小鼠大腦皮層的單細(xì)胞數(shù)據(jù)為adata_sc，空間轉(zhuǎn)錄組數(shù)據(jù)為adata_st.
• 只剪取包含大腦皮層的clusters, 單細(xì)胞數(shù)據(jù)的預(yù)處理步驟參考:cite:tasic2018shared

adata_st = sq.datasets.visium_fluo_adata_crop()
adata_st = adata_st[
    adata_st.obs.cluster.isin([f"Cortex_{i}" for i in np.arange(1, 5)])
].copy()
img = sq.datasets.visium_fluo_image_crop()
adata_sc = sq.datasets.sc_mouse_cortex()

可視化空間和單細(xì)胞數(shù)據(jù)

adata_st.obs

fig, axs = plt.subplots(1, 2, figsize=(20, 5))
sc.pl.spatial(
    adata_st, color="cluster", alpha=0.7, frameon=False, show=False, ax=axs[0]
)
sc.pl.umap(
    adata_sc, color="cell_subclass", size=10, frameon=False, show=False, ax=axs[1]
)
plt.tight_layout()

圖片.png

使用Tangram進(jìn)行空間映射前要先確定一組訓(xùn)練基因，數(shù)目在100-1000個之間，并具有高的信號質(zhì)量。有時候可以通過改變不同的候選訓(xùn)練集印記來觀察結(jié)果的變化。
在這里，我們使用了1401個標(biāo)記基因作為訓(xùn)練基因集(標(biāo)記基因選取為每群細(xì)胞特異表達(dá)的前100個基因）

sc.tl.rank_genes_groups(adata_sc, groupby="cell_subclass", use_raw=False)
markers_df = pd.DataFrame(adata_sc.uns["rank_genes_groups"]["names"]).iloc[0:100, :]
markers = list(np.unique(markers_df.melt().value.values))
len(markers)

接著我們使用下面一組函數(shù)傳入訓(xùn)練基因集, 函數(shù)對基因進(jìn)行篩選，如果改基因在兩個數(shù)據(jù)集中全部是0值，或者不同時出現(xiàn)在兩個數(shù)據(jù)集，該基因?qū)⒈蝗コ?/p>

tg.pp_adatas(adata_sc, adata_st, genes=markers)

兩個數(shù)據(jù)集從1401個基因中篩選得到1280個。

Find alignment
為了對scRNA-seq表達(dá)譜進(jìn)行最佳的空間匹配，使用map_cells_to_space函數(shù)：-num_epochs確定映射迭代次數(shù)；隨著不斷地迭代，當(dāng)評分進(jìn)入平臺期程序隨之停止，評分基于mapped cells vs spatial data之間訓(xùn)練基因集表達(dá)的相似性。

• 默認(rèn)的mapping mode 是mode='cells'，這推薦用于GPU環(huán)境。
• 也可用mode='clusters'，它將同一cluster的細(xì)胞進(jìn)行平均，使用cluster_label傳入標(biāo)簽，速度更快，并適用于空間數(shù)據(jù)和單細(xì)胞數(shù)據(jù)來自不同樣本。
  -使用uniform參數(shù)如果spatial voxels是單細(xì)胞分辨率(例如MERFISH）
  -rna_count_based默認(rèn)spot中細(xì)胞密度與RNA分子數(shù)目是呈線性相關(guān)的。

ad_map = tg.map_cells_to_space(adata_sc, adata_st,
    mode="cells",
#     mode="clusters",
#     cluster_label='cell_subclass',  # .obs field w cell types
    density_prior='rna_count_based',
    num_epochs=500,
    # device="cuda:0",
    device='cpu',
)

映射后的ad_map儲存成AnnData格式，具有如下的結(jié)構(gòu)：

• The cell-by-spot matrix X contains the probability of cell i to be in spot j.
• The obs dataframe contains the metadata of the single cells.
• The var dataframe contains the metadata of the spatial data.
• The uns dictionary contains a dataframe with various information about the training genes (saved as train_genes_df).

細(xì)胞類型映射
接著，使用project_cell_annotation將細(xì)胞類型注釋映射到空間中，并用plot_cell_annotation可視化，perc來設(shè)定colormap的范圍以去除離群值。

tg.project_cell_annotations(ad_map, adata_st, annotation="cell_subclass")
annotation_list = list(pd.unique(adata_sc.obs['cell_subclass']))
tg.plot_cell_annotation_sc(adata_st, annotation_list,perc=0.02)

圖片.png

接著檢驗mapping是否成功，在每個基因水平，如果想查看每個基因的simlarity scores，使用ad_map.uns['train_genes_df']查看

tg.plot_training_scores(ad_map, bins=20, alpha=.5)

圖片.png

New spatial data via aligned single cells

如果mapping mode 是mode='cells'，我們可以生成一個新的包含經(jīng)匹配單細(xì)胞的數(shù)據(jù)，該函數(shù)project_genes接受ad_map以及相應(yīng)的單細(xì)胞矩陣adata_sc作為輸入。得到的ad_ge是voxel-by-gene AnnData，與adata_st近似，但是包含mapped single cell單細(xì)胞的表達(dá)矩陣而不是Visium，接下來的分析都圍繞ad_ge。

ad_ge = tg.project_genes(adata_map=ad_map, adata_sc=adata_sc)
ad_ge

• 這里的ad_ge: AnnData object with n_obs × n_vars = 324 × 36826
• 而ad_map: AnnData object with n_obs × n_vars = 21697 × 324
  注意兩者區(qū)別：
  ad_map的n_obs表示單細(xì)胞數(shù)據(jù)中的不同細(xì)胞，顯示其在spot中的概率。
  ad_ge的n_obs表示不同spot，n_vars代表不同基因，表示單細(xì)胞mapped之后新的空間矩陣。

挑選幾個評分低的tranning gene

genes = ['rragb', 'trim17', 'eno1b']
ad_map.uns['train_genes_df'].loc[genes]

使用plot_genes函數(shù)傳入兩個voxel-by-gene AnnData，分別是adata_measured實際檢測的空間矩陣adata_st，以及上一步得到的預(yù)測的空間矩陣adata_predicted: ad_ge。

tg.plot_genes_sc(genes, adata_measured=adata_st, adata_predicted=ad_ge, perc=0.02)

圖片.png

以上結(jié)果可以解釋為什么這些基因的類似評分較低。因為有些基因在不同檢測技術(shù)下具有不同的程度的稀疏度(sparity)，通常scRNA-seq的稀疏性相較于空間轉(zhuǎn)錄組更高。由于空間轉(zhuǎn)錄組技術(shù)存在dropouts效應(yīng)，由于某些基因信號缺陷，Tangram不能很好得在空間上對這些基因進(jìn)行匹配。但只要這些tranning基因有較高的測量質(zhì)量，我們可以相信映射的結(jié)果并使用tangram預(yù)測其在空間位置上的表達(dá)。
我們也可以用這種插值方法推斷一些沒有通過空間技術(shù)檢測到的基因的表達(dá)，前提是這些基因通過單細(xì)胞測序檢測到。

genes=['loc102633833', 'gm5700', 'gm8292']
tg.plot_genes_sc(genes, adata_measured=adata_st, adata_predicted=ad_ge, perc=0.02)

圖片.png

也可以通過compare_spatial_geneexp對所有基因similarity scores進(jìn)行計算，需要傳入不同的兩個空間矩陣ad_ge和adata_st，Training genes會在is_training列用布爾值標(biāo)記。如果我們也傳入single cell AnnData，函數(shù)將返回sparsity_sc (single cell data sparsity) 和sparsity_diff (spatial data sparsity - single cell data sparsity) 兩列值，用于plot_test_scores作圖。

df_all_genes = tg.compare_spatial_geneexp(ad_ge, adata_st, adata_sc)
df_all_genes

test_genes的預(yù)測情況可以可視化：

# sns.scatterplot(data=df_all_genes, x='score', y='sparsity_sp', hue='is_training', alpha=.5);  # for legacy
tg.plot_auc(df_all_genes);

圖片.png

這張圖是tangram中用于檢驗驗證基因效果最重要的圖。圖中每個點代表一個基因，x軸表示基因的相似性評分，y軸表示基因在空轉(zhuǎn)數(shù)據(jù)中的稀疏性。不出所料，高稀疏性的基因相似性評分較低，證明tangram在test genes中有較佳的預(yù)測效果。我們可以用曲線下面積去評估預(yù)測效果。

以下這些基因稀疏性高，但能很好得預(yù)測出來。

genes=['tfap2b', 'zic4']
tg.plot_genes_sc(genes, adata_measured=adata_st, adata_predicted=ad_ge, perc=0.02)

圖片.png

某些非稀疏基因呈現(xiàn)出特異的表達(dá)模式，tangram可以更加凸顯這一模式。

genes = ['cd34', 'rasal1']
tg.plot_genes_sc(genes, adata_measured=adata_st, adata_predicted=ad_ge, perc=0.02)

圖片.png

有些基因沒有深入闡明功能，很難通過空轉(zhuǎn)技術(shù)檢測到，而tangram可以提供預(yù)測。

genes = ['gm33027', 'gm5431']
tg.plot_genes_sc(genes[:5], adata_measured=adata_st, adata_predicted=ad_ge, perc=0.02)

圖片.png