關(guān)于本教程的數(shù)據(jù)說明趾访，見：http://www.reibang.com/p/2ac31d3a5560
stLearn是python語言編寫的值骇，更新速度賊快，網(wǎng)上很多筆記教程已經(jīng)運(yùn)行不通了后豫，前兩天一口氣更新了三個版本：
今天再看的時候這三個版本已經(jīng)沒啦时捌，又又又更新了一個：
本教程基于0.4.0版本，請安裝：pip install stlearn==0.4.0
分析基本原理
主要用于鑒定空間轉(zhuǎn)錄組切片上的熱點(diǎn)區(qū)域，即CCI hotspot霞赫，cell cell interaction hotspot，這個區(qū)域具有特征：These regions with high cell diversity and L-R activities are considered as hotspots in the tissue with most likely cell-cell interaction activities肥矢。
第一點(diǎn)：細(xì)胞類型多端衰；第二點(diǎn)：配體-受體共表達(dá)。因此甘改，在分析之前旅东，最好已經(jīng)對空間轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行了spot細(xì)胞類型注釋。
分析步驟：
- 1.加載已知的配體受體基因?qū)?/li>
- 2.鑒定配受體對顯著互作的spot
- 3.對于每個L-R對和每個細(xì)胞類型-細(xì)胞類型組合：count the instances where neighbours of a signficant spot for that LR pair link two given cell types
- 4.擾動細(xì)胞類型標(biāo)簽楼誓，鑒定顯著的互作(pvalue<0.05)
- 5.可視化CCI結(jié)果
一玉锌、配體-受體分析
stLearn CCI(Cell-Cell Interaction)流程的第一步是配體受體(LR)分析。
這個分析從候選配體-受體數(shù)據(jù)庫中鑒定配體-受體相互作用的顯著spot疟羹。
運(yùn)行時將嚴(yán)重依賴于可用的數(shù)據(jù)集和計(jì)算資源; 請注意主守，分析支持多線程。
1. 數(shù)據(jù)加載和預(yù)處理
數(shù)據(jù)目錄結(jié)構(gòu)：
注意：LR分析不需要log1p數(shù)據(jù)
conda activate stlearn
import stlearn as st
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Loading raw data #
data_dir = "./"
data = st.Read10X(data_dir)
data.var_names_make_unique()
st.add.image(adata=data,
imgpath=data_dir+"spatial/tissue_hires_image.png",
library_id="V1_Breast_Cancer_Block_A_Section_1", visium=True)
# Basic normalisation #
st.pp.filter_genes(data, min_cells=3)
st.pp.normalize_total(data) # NOTE: no log1p
作者提示到：對于上面的細(xì)胞注釋類型信息榄融，是使用Seurat生成的参淫。數(shù)據(jù)下載：tutorials
對于每一個spot可以沒有細(xì)胞類型打分，如果有細(xì)胞類型打分愧杯，do need to add the dominant cell type to the adata.obs slot with the same key as the cell type scores added to the adata.uns slot涎才。
# Adding the label transfer results, #
spot_mixtures = pd.read_csv(f'{data_dir}/tutorials/label_transfer_bc.txt', index_col=0, sep='\t')
labels = spot_mixtures.loc[:,'predicted.id'].values.astype(str)
spot_mixtures = spot_mixtures.drop(['predicted.id','prediction.score.max'],axis=1)
spot_mixtures.columns = [col.replace('prediction.score.', '')
for col in spot_mixtures.columns]
# Note the format
print(labels)
print(spot_mixtures)
# Check is in correct order
print('Spot mixture order correct?: ', np.all(spot_mixtures.index.values==data.obs_names.values))
# NOTE: using the same key in data.obs & data.uns
# Adding the dominant cell type labels per spot
data.obs['cell_type'] = labels
data.obs['cell_type'] = data.obs['cell_type'].astype('category')
# Adding the cell type scores
data.uns['cell_type'] = spot_mixtures
st.pl.cluster_plot(data, use_label='cell_type')
plt.savefig('./s1.cluster_plot.jpg')
plt.close()
細(xì)胞注釋后的標(biāo)簽圖：
2.運(yùn)行配體受體分析
運(yùn)行：st.tl.cci.run這一步耗時比較久，建議在集群上進(jìn)行分析力九，n_pairs教程推薦10000次
# Loading the LR databases available within stlearn (from NATMI)
lrs = st.tl.cci.load_lrs(['connectomeDB2020_lit'], species='human')
print(len(lrs))
# Running the analysis #
# min_spots：Filter out any LR pairs with no scores for less than min_spots
# distance： None defaults to spot+immediate neighbours; distance=0 for within-spot mode
# n_pairs： Number of random pairs to generate; low as example, recommend ~10,000
# n_cpus：Number of CPUs for parallel. If None, detects & use all available.
# 這一步耗時比較久耍铜，如果選擇n_pairs=10000，可能要好幾個小時跌前，可以多選擇幾個cpu
st.tl.cci.run(data, lrs, min_spots = 20, distance=None, n_pairs=100, n_cpus=15)
# A dataframe detailing the LR pairs ranked by number of significant spots.
lr_info = data.uns['lr_summary']
print('\n', lr_info)
Spot neighbour indices保存在data.obsm['spot_neighbours'] 與 data.obsm['spot_neigh_bcs']中
結(jié)果保存在：
- lr_scores 保存在 adata.obsm['lr_scores']
- p_vals 保存在 adata.obsm['p_vals']
- p_adjs 保存在 adata.obsm['p_adjs']
- -log10(p_adjs) 保存在 adata.obsm['-log10(p_adjs)']
- lr_sig_scores 保存在 adata.obsm['lr_sig_scores']
每個spot結(jié)果在adata.obsm的列與adata.uns['lr_summary']中的行名一樣且順序一致
LR結(jié)果的summary在data.uns['lr_summary']中：
3.P-value矯正
可以使用不同的方法矯正 p 值; p 值已經(jīng)通過運(yùn)行 st.tl.cci.run 進(jìn)行了矯正棕兼。不同的矯正方法差別在于correct_axis參數(shù)，可以是LR抵乓，spot或者none伴挚。
st.tl.cci.adj_pvals(data, correct_axis='spot', pval_adj_cutoff=0.05, adj_method='fdr_bh')
4.LRs排序并可視化
# Showing the rankings of the LR from a global and local perspective.
# Ranking based on number of significant hotspots.
st.pl.lr_summary(data, n_top=500)
plt.savefig('./s2.lr_summary500.jpg')
plt.close()
st.pl.lr_summary(data, n_top=50, figsize=(10,3))
plt.savefig('./s2.lr_summary50.jpg')
plt.close()
根據(jù)每個LR顯著spots數(shù)展示LR top 50, 500：
5.診斷圖
LR 分析的一個關(guān)鍵是識別顯著的hotspot時控制 LR 的表達(dá)水平和頻率，所以我們的診斷圖應(yīng)該是hotspot區(qū)的 配受體對與其表達(dá)水平和表達(dá)頻率之間不相關(guān)灾炭。
下面的診斷圖可以檢查這個情況茎芋；如果相關(guān)，需要設(shè)置更大的n_pairs擾動次數(shù)蜈出。
st.pl.lr_diagnostics(data, figsize=(10,2.5))
plt.savefig('./s2.lr_diagnostics.jpg')
plt.close()
- 左圖: LR 表達(dá)水平(LR 對中非零點(diǎn)平均基因中位數(shù)表達(dá))與 LR 排名的關(guān)系
- 右圖: LR 表達(dá)頻率(LR 對中每個基因的平均零點(diǎn)比例)與 LR 排名之間的關(guān)系
這個例子中田弥，LR 表達(dá)頻率和significant spots 數(shù)量之間的有比較弱的相關(guān)性，因此铡原，這n _ pairs 參數(shù)應(yīng)該設(shè)置得更高皱蹦，以創(chuàng)建更準(zhǔn)確的背景分布(文獻(xiàn)中中使用了n _ pairs=10,000)煤杀。
st.pl.lr_n_spots(data, n_top=50, figsize=(11, 4),max_text=100)
plt.savefig('./s2.lr_n_spots50.jpg')
plt.close()
st.pl.lr_n_spots(data, n_top=500, figsize=(11, 4), max_text=100)
plt.savefig('./s2.lr_n_spots500.jpg')
plt.close()
根據(jù)每個配受體對在多少個spot里面lr打分顯著進(jìn)行排序，展示top500沪哺、top50：
6.Biologic Plots (可選)
配受體對的功能富集分析沈自，結(jié)果保存在：adata.uns['lr_go']
r_path = "/share/nas5/zhangj/biosoft/miniconda3/envs/R4/bin/Rscript"
st.tl.cci.run_lr_go(data, r_path)
# save
st.pl.lr_go(data, lr_text_fp={'weight': 'bold', 'size': 10}, rot=15, figsize=(12,3.65),n_top=15, show=False)
plt.savefig('./s3.lr_go.jpg')
plt.close()
二、L-R可視化
LR _ result _ plot 在空間數(shù)據(jù)中繪制指定 LR 對的分析結(jié)果辜妓，需要的值存儲在：data.obsm中枯途，包括: lr_scores, p_vals, p_adjs, -log10(p_adjs), lr_sig_scores
# Just choosing one of the top from lr_summary
best_lr = data.uns['lr_summary'].index.values[0]
best_lr
stats = ['lr_scores', 'p_vals', 'p_adjs', '-log10(p_adjs)']
fig, axes = plt.subplots(ncols=len(stats), figsize=(16,6))
for i, stat in enumerate(stats):
st.pl.lr_result_plot(data, use_result=stat, use_lr=best_lr, show_color_bar=False, ax=axes[i])
axes[i].set_title(f'{best_lr} {stat}')
plt.savefig('./s4.lr_result_plot.jpg')
plt.close()
fig, axes = plt.subplots(ncols=2, figsize=(8,6))
st.pl.lr_result_plot(data, use_result='-log10(p_adjs)', use_lr=best_lr, show_color_bar=False, ax=axes[0])
st.pl.lr_result_plot(data, use_result='lr_sig_scores', use_lr=best_lr, show_color_bar=False, ax=axes[1])
axes[0].set_title(f'{best_lr} -log10(p_adjs)')
axes[1].set_title(f'{best_lr} lr_sig_scores')
plt.savefig('./s4.lr_result_plot-1.jpg')
plt.close()
這兩張圖一般是文獻(xiàn)中出現(xiàn)的比較多的圖，選擇的配受體對的共表達(dá)打分籍滴，顯著性：
打分結(jié)合顯著性后的結(jié)果：
三酪夷、LR 說明可視化
這些圖旨在幫助解釋細(xì)胞之間cross-talk的方向性
st.pl.lr_plot(data, best_lr, inner_size_prop=0.1, outer_mode='binary', pt_scale=5, use_label=None, show_image=True,sig_spots=False)
plt.savefig('./s5.lr_plot.jpg')
plt.close()
st.pl.lr_plot(data, best_lr, outer_size_prop=1, outer_mode='binary', pt_scale=20, use_label=None, show_image=True,sig_spots=True)
plt.savefig('./s5.lr_plot-1.jpg')
plt.close()
紅色為配體，綠色為受體孽惰，藍(lán)色為共表達(dá)晚岭。有助于了解配體/受體在哪里以及在多大程度上表達(dá)。理想情況是：受體位于細(xì)胞表面勋功，配體從細(xì)胞表面滲透出來坦报。
具有顯著性的spot：
看表達(dá)連續(xù)值情況：
# All spots #
st.pl.lr_plot(data, best_lr,
inner_size_prop=0.04, middle_size_prop=.07, outer_size_prop=.4,
outer_mode='continuous', pt_scale=60,
use_label=None, show_image=True,
sig_spots=False)
plt.savefig('./s5.lr_plot-2.jpg')
plt.close()
# Only significant spots #
st.pl.lr_plot(data, best_lr,
inner_size_prop=0.04, middle_size_prop=.07, outer_size_prop=.4,
outer_mode='continuous', pt_scale=60,
use_label=None, show_image=True,
sig_spots=True)
plt.savefig('./s5.lr_plot-3.jpg')
plt.close()
這僅在放大并希望同時顯示細(xì)胞信息和互作方向時有用：
st.pl.lr_plot(data, best_lr,
inner_size_prop=0.08, middle_size_prop=.3, outer_size_prop=.5,
outer_mode='binary', pt_scale=50,
show_image=True, arrow_width=10, arrow_head_width=10,
sig_spots=True, show_arrows=True)
plt.savefig('./s5.lr_plot-4.jpg')
plt.close()
這個圖使用交互式繪圖效果更好，交互式繪圖參考：https://stlearn.readthedocs.io/en/latest/tutorials/Interactive_plot.html
這里看不出來圖上的細(xì)節(jié)：
我們還可以觀察到配體或受體在哪里與占優(yōu)勢的點(diǎn)細(xì)胞類型同時表達(dá)/共表達(dá)：
st.pl.lr_plot(data, best_lr,
inner_size_prop=0.08, middle_size_prop=.3, outer_size_prop=.5,
outer_mode='binary', pt_scale=150,
use_label='cell_type', show_image=True,
sig_spots=True)
plt.savefig('./s5.lr_plot-5.jpg')
plt.close()
外點(diǎn)顯示配體(紅色)狂鞋、受體(綠色)和共表達(dá)(藍(lán)色)的表達(dá)片择。內(nèi)點(diǎn)由占參與互作的主要細(xì)胞類型著色：
四、細(xì)胞間互作分析CCIs
Calls significant celltype-celltype interactions based on cell-type data randomisation骚揍。
確定了LR 共表達(dá)的顯著spot區(qū)域后字管，現(xiàn)在可以確定顯著的互作細(xì)胞類型，st.tl.cci.run_cci參數(shù)含義如下：
- min_spots: Specifies the minimum number of spots where LR score present to include in subsequent analysis
- spot_mixtures: # If True will use the label transfer scores, 前面導(dǎo)入的細(xì)胞類型打分信不，so spots can have multiple cell types if score>cell_prop_cutoff
- cell_prop_cutoff: # Spot considered to have cell type if score>0.2
- sig_spots: # Only consider neighbourhoods of spots which had significant LR scores.
- n_perms: # Permutations of cell information to get background, recommend ~1000
# Running the counting of co-occurence of cell types and LR expression hotspots #
# Spot cell information either in data.obs or data.uns
# 這一步運(yùn)行也會比較耗時嘲叔，可以適當(dāng)增加n_perms次數(shù)，推薦使用1000
st.tl.cci.run_cci(data, 'cell_type', min_spots=3, spot_mixtures=True, cell_prop_cutoff=0.2, sig_spots=True, n_perms=100 )
Significant counts of cci_rank interactions for each LR pair stored in dictionary per_lr_cci_cell_type抽活。
五借跪、診斷圖: 檢查互作與細(xì)胞類型種類數(shù)的相關(guān)性
理論上cci_rank與cell type 種類應(yīng)該無關(guān)，如果有相關(guān)酌壕，可以增加st.tl.cci.run_cci函數(shù)中的n_perms值，增大擾動次數(shù)歇由。
st.pl.cci_check(data, 'cell_type')
plt.savefig('./s6.cci_check.jpg')
plt.close()
六卵牍、CCI 可視化
如果你更喜歡使用R來進(jìn)行繪圖，可以保存 anndata 對象中的鄰接矩陣沦泌，并使用 CellChat 包的可視化功能來實(shí)現(xiàn)可視化糊昙。
# Visualising the no. of interactions between cell types across all LR pairs #
pos_1 = st.pl.ccinet_plot(data, 'cell_type', return_pos=True)
plt.savefig('./s7.ccinet_plot.jpg')
plt.close()
# Just examining the cell type interactions between selected pairs #
lrs = data.uns['lr_summary'].index.values[0:3]
for best_lr in lrs[0:3]:
st.pl.ccinet_plot(data, 'cell_type', best_lr, min_counts=2,figsize=(10,7.5), pos=pos_1)
plt.savefig(fname=f's8.ccinet_plot_{best_lr}.jpg')
plt.close()
總互作：
圈的大小：某個細(xì)胞群參與的spot互作數(shù)谢谦；箭頭大惺臀：兩個細(xì)胞群間的總互作數(shù)：
CCI和弦圖
可視化少數(shù)細(xì)胞類型之間的互作時和弦圖非常有用
st.pl.lr_chord_plot(data, 'cell_type')
plt.savefig('./s9.lr_chord_plot.jpg')
plt.close()
lrs = data.uns['lr_summary'].index.values[0:3]
for lr in lrs:
st.pl.lr_chord_plot(data, 'cell_type', lr)
plt.savefig(fname=f'./s9.lr_chord_plot_{lr}.jpg')
plt.close()
總的和弦圖：
單個配受體：
Heatmap Visualisations
LR-CCI-Map
我們還提供了一些熱度圖的可視化萝衩，這樣你就可以同時觀察每個celltype-celltype 的多個 LR 對之間的相互作用
# This will automatically select the top interacting CCIs and their respective LRs #
st.pl.lr_cci_map(data, 'cell_type', lrs=None, min_total=100, figsize=(20,8))
plt.savefig('./s10.lr_cci_map.jpg')
plt.close()
# You can also put in your own LR pairs of interest #
lrs = data.uns['lr_summary'].index.values[0:11]
st.pl.lr_cci_map(data, 'cell_type', lrs=lrs, min_total=100, figsize=(20,12))
plt.savefig('./s10.lr_cci_map-1.jpg')
plt.close()
自動選取top LR對：
也可以用戶指定配體受體后：
CCI Maps
Heatmap visualising sender->receivers of cell type interactions
# of interactions 指的是一個點(diǎn)與接受細(xì)胞類型表達(dá)配體和源細(xì)胞類型表達(dá)受體在同一鄰域的次數(shù)。
st.pl.cci_map(data, 'cell_type')
plt.savefig('./s11.cci_map.jpg')
plt.close()
lrs = data.uns['lr_summary'].index.values[0:3]
for lr in lrs[0:3]:
st.pl.cci_map(data, 'cell_type', lr)
plt.savefig(fname='s11.cci_map_{lr}.jpg')
plt.close()
總的：
指定配受體的：
七没咙、Spatial cell type interactions
j結(jié)合LR分析與CCI分析猩谊，我們現(xiàn)在可以看到這些細(xì)胞在組織的哪個部位相互通訊。
best_lr = lrs[0]
### This will plot with simple black arrows ####
st.pl.lr_plot(data, best_lr, outer_size_prop=1, outer_mode=None,
pt_scale=40, use_label='cell_type', show_arrows=True,
show_image=True, sig_spots=False, sig_cci=True,
arrow_head_width=4,
arrow_width=1, cell_alpha=.8 )
plt.savefig('./s12.lr_plot_{best_lr}.jpg')
plt.close()
### This will colour the spot by the mean LR expression in the spots connected by arrow
st.pl.lr_plot(data, best_lr, outer_size_prop=1, outer_mode=None,
pt_scale=10, use_label='cell_type', show_arrows=True,
show_image=True, sig_spots=False, sig_cci=True,
arrow_head_width=4, arrow_width=2,arrow_cmap='YlOrRd', arrow_vmax=1.5)
plt.savefig('./s12.lr_plot_{best_lr}.jpg')
plt.close()
# 保存數(shù)據(jù)
# 設(shè)置備份地址
data.filename = './LR.h5ad'
# 查看是否備份成功
print(data.isbacked) # True
八祭刚、Visualisation Tips
空間切片中的 CCIs 具有豐富的信息牌捷，以上哪種可視化方式會有用，將取決于你想要強(qiáng)調(diào)的生物學(xué)和關(guān)鍵方面涡驮。根據(jù)我們的經(jīng)驗(yàn)暗甥，顯示感興趣的 LR 統(tǒng)計(jì)數(shù)據(jù)，然后使用細(xì)胞類型信息和箭頭繪制圖表是很有用的捉捅。
在后一種情節(jié)中撤防，最好是突出顯示感興趣的區(qū)域(如上面的箭頭) ，然后放大顯示棒口。因此寄月，你可以突出顯示預(yù)計(jì)將發(fā)生有趣 cci 的代表區(qū)域。請參閱“ Interactive stLearn”教程陌凳，使用交互式散景應(yīng)用程序輕松快速地創(chuàng)建此類可視化效果剥懒。
參考資料：
