關于本教程的數(shù)據(jù)說明，見：http://www.itdecent.cn/p/2ac31d3a5560
stLearn是python語言編寫的，更新速度賊快，網(wǎng)上很多筆記教程已經(jīng)運行不通了，前兩天一口氣更新了三個版本：
今天再看的時候這三個版本已經(jīng)沒啦，又又又更新了一個：
本教程基于0.4.0版本，請安裝：pip install stlearn==0.4.0
分析基本原理
主要用于鑒定空間轉(zhuǎ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。
第一點：細胞類型多；第二點：配體-受體共表達。因此，在分析之前，最好已經(jīng)對空間轉(zhuǎn)錄組數(shù)據(jù)進行了spot細胞類型注釋。
分析步驟：
- 1.加載已知的配體受體基因?qū)?/li>
- 2.鑒定配受體對顯著互作的spot
- 3.對于每個L-R對和每個細胞類型-細胞類型組合：count the instances where neighbours of a signficant spot for that LR pair link two given cell types
- 4.擾動細胞類型標簽，鑒定顯著的互作(pvalue<0.05)
- 5.可視化CCI結(jié)果
一、配體-受體分析
stLearn CCI(Cell-Cell Interaction)流程的第一步是配體受體(LR)分析。
這個分析從候選配體-受體數(shù)據(jù)庫中鑒定配體-受體相互作用的顯著spot。
運行時將嚴重依賴于可用的數(shù)據(jù)集和計算資源; 請注意，分析支持多線程。
1. 數(shù)據(jù)加載和預處理
數(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
作者提示到：對于上面的細胞注釋類型信息，是使用Seurat生成的。數(shù)據(jù)下載：tutorials
對于每一個spot可以沒有細胞類型打分，如果有細胞類型打分，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()
細胞注釋后的標簽圖：
2.運行配體受體分析
運行：st.tl.cci.run這一步耗時比較久，建議在集群上進行分析，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)通過運行 st.tl.cci.run 進行了矯正。不同的矯正方法差別在于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 分析的一個關鍵是識別顯著的hotspot時控制 LR 的表達水平和頻率，所以我們的診斷圖應該是hotspot區(qū)的 配受體對與其表達水平和表達頻率之間不相關。
下面的診斷圖可以檢查這個情況；如果相關，需要設置更大的n_pairs擾動次數(shù)。
st.pl.lr_diagnostics(data, figsize=(10,2.5))
plt.savefig('./s2.lr_diagnostics.jpg')
plt.close()
- 左圖: LR 表達水平(LR 對中非零點平均基因中位數(shù)表達)與 LR 排名的關系
- 右圖: LR 表達頻率(LR 對中每個基因的平均零點比例)與 LR 排名之間的關系
這個例子中，LR 表達頻率和significant spots 數(shù)量之間的有比較弱的相關性，因此，這n _ pairs 參數(shù)應該設置得更高，以創(chuàng)建更準確的背景分布(文獻中中使用了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打分顯著進行排序，展示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)的比較多的圖，選擇的配受體對的共表達打分，顯著性：
打分結(jié)合顯著性后的結(jié)果：
三、LR 說明可視化
這些圖旨在幫助解釋細胞之間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()
紅色為配體，綠色為受體，藍色為共表達。有助于了解配體/受體在哪里以及在多大程度上表達。理想情況是：受體位于細胞表面，配體從細胞表面滲透出來。
具有顯著性的spot：
看表達連續(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()
這僅在放大并希望同時顯示細胞信息和互作方向時有用：
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
這里看不出來圖上的細節(jié)：
我們還可以觀察到配體或受體在哪里與占優(yōu)勢的點細胞類型同時表達/共表達：
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()
外點顯示配體(紅色)、受體(綠色)和共表達(藍色)的表達。內(nèi)點由占參與互作的主要細胞類型著色：
四、細胞間互作分析CCIs
Calls significant celltype-celltype interactions based on cell-type data randomisation。
確定了LR 共表達的顯著spot區(qū)域后，現(xiàn)在可以確定顯著的互作細胞類型，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, 前面導入的細胞類型打分，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
# 這一步運行也會比較耗時，可以適當增加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。
五、診斷圖: 檢查互作與細胞類型種類數(shù)的相關性
理論上cci_rank與cell type 種類應該無關，如果有相關，可以增加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來進行繪圖，可以保存 anndata 對象中的鄰接矩陣，并使用 CellChat 包的可視化功能來實現(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()
總互作：
圈的大小：某個細胞群參與的spot互作數(shù)；箭頭大?。簝蓚€細胞群間的總互作數(shù)：
CCI和弦圖
可視化少數(shù)細胞類型之間的互作時和弦圖非常有用
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 指的是一個點與接受細胞類型表達配體和源細胞類型表達受體在同一鄰域的次數(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)在可以看到這些細胞在組織的哪個部位相互通訊。
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ù)
# 設置備份地址
data.filename = './LR.h5ad'
# 查看是否備份成功
print(data.isbacked) # True
八、Visualisation Tips
空間切片中的 CCIs 具有豐富的信息，以上哪種可視化方式會有用，將取決于你想要強調(diào)的生物學和關鍵方面。根據(jù)我們的經(jīng)驗，顯示感興趣的 LR 統(tǒng)計數(shù)據(jù)，然后使用細胞類型信息和箭頭繪制圖表是很有用的。
在后一種情節(jié)中，最好是突出顯示感興趣的區(qū)域(如上面的箭頭) ，然后放大顯示。因此，你可以突出顯示預計將發(fā)生有趣 cci 的代表區(qū)域。請參閱“ Interactive stLearn”教程，使用交互式散景應用程序輕松快速地創(chuàng)建此類可視化效果。
參考資料：
