空转和单细胞确实有着不同的目的与意义，更多的分析更适合单个样本逐一进行，所以以其中单一样本为例。
COMMOT推断细胞通讯信号流
import os
import scanpy as sc
import matplotlib.pyplot as plt
import commot as ct
from pathlib import Path
os.getcwd()
# 读取空转数据
adata_dct = {}
for i in Path("../").glob("rawdata_st/*"):
_s = str(i).split('_')[2]
_a = sc.read_visium(i,library_id=_s)
_a.obs.index = [_s + "_" + bc[0:len(bc)] for bc in _a.obs.index.tolist()]
_a.var_names_make_unique()
adata_dct[_s] = _a
st_adata = sc.concat(adata_dct,label="library_id",uns_merge="unique")
# 数据标准化和归一化 并保存原始数据
st_adata.raw = st_adata
sc.pp.normalize_total(st_adata, inplace=True)
sc.pp.log1p(st_adata)
st_adata_dis500 = st_adata.copy()
# 载入细胞通讯数据库并筛选信号
df_cellchat = ct.pp.ligand_receptor_database(species='human', signaling_type='Secreted Signaling',
database='CellChat')
print(df_cellchat.shape)
df_cellchat_filtered = ct.pp.filter_lr_database(df_cellchat, st_adata_dis500, min_cell_pct=0.05)
print(df_cellchat_filtered.shape)
# 选取单一样本进行推断
ad = st_adata_dis500[st_adata_dis500.obs.library_id == 'D2', :].copy()
a = ad.uns.get('spatial').get('D2')
ad.uns['spatial'] = {}
ad.uns['spatial']['D2'] = a
ct.tl.spatial_communication(ad,
database_name='cellchat', df_ligrec=df_cellchat_filtered, dis_thr=500, heteromeric=True, pathway_sum=True) # 空间距离限制为 500
# 判断信号通路的空间方向并可视化 （以TGF-β为例）
ct.tl.communication_direction(ad, database_name='cellchat', pathway_name='TGFb', k=5)
ct.pl.plot_cell_communication(ad, database_name='cellchat', pathway_name='TGFb', plot_method='grid', background_legend=True,
scale=0.000001, ndsize=8, grid_density=0.4, summary='sender', background='image', clustering='leiden', cmap='Alphabet',
normalize_v = True, normalize_v_quantile=0.995) # plot_method='stream' 流线型呈现
SPATA推断细胞轨迹
R
library(Seurat)
library(tidyverse)
library(patchwork)
library(SPATA) # SPATA2和SeuratV5兼容，之前对象都为V4建立，所以以SPATA演示
# 读取数据
spata_obj=SPATA::initiateSpataObject_10X('../liver-spatial/rawdata_st/GSM5764425_D2/',sample_name='D2') # SPATA旧版本需10x数据以outs形式读入
# 建立SPATA对象
spata_obj <- createTrajectories(object = spata_obj) # 进入Shiny交互界面
# 获取保存的轨迹名
getTrajectoryNames(object = spata_obj)
# 可视化随轨迹方向上的分群数量变化
plotTrajectory(object = spata_obj,
trajectory_name = "TGF",
color_to = "seurat_clusters",
pt_size = 1.7,
pt_clrp = "npg",
display_image = FALSE) +
theme(legend.position = "top")
SPATA结合Cell2location/RCTD
# 添加自己想要的Feature进行可视化
# 可结合cell2location的结果
# barcodes进行匹配修改
cell2location <- read.csv("../liver-spatial/st_cell2location_res.csv",row.names = 1) # 可换成RCTD结果
head(rownames(cell2location))
"H2_AAACATTTCCCGGATT-1"
head(spata_obj@fdata$barcodes)
"AAACTGCTGGCTCCAA-1_D2"
for ( i in 1:length(rownames(cell2location))){
cell2location$barcodes[i] <- str_split(rownames(cell2location)[i],"_")[[1]][2]
}
for ( i in 1:length(rownames(cell2location))){
cell2location$barcodes[i] <- paste(cell2location$barcodes[i],str_split(rownames(cell2location)[i],"_")[[1]][1],sep = "_")
}
head(cell2location$barcodes)
"AAACATTTCCCGGATT-1_H2"
# 仅读取D2样本
table(spata_obj@fdata$barcodes %in% cell2location$barcodes) # 查看数目是否匹配
TRUE
532
spata_obj@fdata <- merge(spata_obj@fdata,cell2location,by="barcodes")
# 以cell2location的细胞比例结果可视化
getFeatureNames(spata_obj)
plotTrajectoryFeatures(object = spata_obj,
trajectory_name = "TGF",
features = "Hep2",
smooth_method = "loess",
smooth_span = 0.2,
smooth_se = TRUE)
all_genes <- getGenes(spata_obj)
all_gene_sets <- getGeneSets(spata_obj)
Results <- grep('TGF', all_gene_sets, value = T)
[1] "BC_TGFB_PATHWAY"
[2] "RCTM_TGF_BETA_RECEPTOR_SIGNALING_IN_EMT_EPITHELIAL_TO_MESENCHYMAL_TRANSITION"
[3] "RCTM_DOWNREGULATION_OF_TGF_BETA_RECEPTOR_SIGNALING"
[4] "RCTM_TGF_BETA_RECEPTOR_SIGNALING_ACTIVATES_SMADS"
[5] "RCTM_SIGNALING_BY_TGF_BETA_RECEPTOR_COMPLEX"
[6] "RCTM_SIGNALING_BY_TGF_BETA_RECEPTOR_COMPLEX_IN_CANCER"
[7] "RCTM_LOSS_OF_FUNCTION_OF_TGFBR1_IN_CANCER"
[8] "RCTM_SIGNALING_BY_TGF_BETA_FAMILY_MEMBERS"
[9] "HM_TGF_BETA_SIGNALING"
plotTrajectoryGeneSets(
object = spata_obj,
trajectory_name = "TGF",
gene_sets = Results,
display_trajectory_parts = FALSE) + # results in missing vertical lines
ggplot2::theme(legend.position = "top")
SPATA拟合分析判定关键变量
# 寻找沿轨迹的有趣变量，进行拟合
atdf_gs <- assessTrajectoryTrends(object = spata_obj,
trajectory_name = "TGF",
variables = all_genes)
atdf_gene_sets <- assessTrajectoryTrends(object = spata_obj,
trajectory_name = "TGF",
variables = all_gene-sets)
atdf_gs
atdf_gene_sets
# 不同模型的可视化
plotTrajectoryFit(object = spata_obj,
trajectory_name = "TGF",
variable = "TGFA",
display_residuals = TRUE) +
ggplot2::theme(legend.position = "top")
# 报错，查看源码发现是left_join的参数设置有误
trace(SPATA::plotTrajectoryFit, edit = T) 将key改为by
# 以残差下AUC面积评估拟合模型
# 判定TGFA随着轨迹是Late peak的表达模式
