空间转录组的细胞通讯分析

从单细胞空间转录组 (ST) 数据中，推断由分泌蛋白介导的细胞间通讯。示例数据为CosMx平台的肝癌ST样本。

创建SpaCET对象

创建对象需准备三类输入数据：

count矩阵：基因名 (行) × 细胞ID (列)。
空间坐标：细胞ID (行) × 坐标 (列)，第1、2 列分别为X、Y (单位mm)。
meta信息 (可选)：细胞ID (行) × 注释 (列)。

library(SecAct)
library(SpaCET)

# https://hpc.nih.gov/~Jiang_Lab/SecAct_Package/LIHC_CosMx_data.rda
load("LIHC_CosMx_data.rda")

counts[1:5,1:5]
5 x 5 sparse Matrix of class "dgCMatrix"
      c_2_100_733 c_2_101_240 c_2_101_339 c_2_101_452 c_2_102_179
AATK            .           .           .           1           2
ABL1            .           2           .           .           1
ABL2            .           .           1           .           .
ACACB           1           1           2           .           .
ACE             .           .           1           .           .

metaData[1:5,]
                  cellType         niche
c_2_100_733     Tumor_core tumor subtype
c_2_101_240 Tumor_boundary     interface
c_2_101_339     Tumor_core         tumor
c_2_101_452     Tumor_core tumor subtype
c_2_102_179     Tumor_core         tumor

spotCoordinates[1:5,]
            x_slide_mm y_slide_mm
c_2_100_733    6.70816    9.03340
c_2_101_240    7.50024    8.86156
c_2_101_339    7.60560    8.73604
c_2_101_452    7.62984    9.17260
c_2_102_179    7.86824    8.99368

SpaCET_obj <- create.SpaCET.object(counts, spotCoordinates, metaData = metaData, platform = "CosMx")
SpaCET_obj <- SpaCET.quality.control(SpaCET_obj, min.genes = 50)
SpaCET.visualize.spatialFeature(SpaCET_obj, spatialType = "QualityControl", spatialFeatures = c("UMI", "Gene"), pointSize = 0.1)

使用对象里面的细胞类型注释信息，查看各细胞类型的空间分布：

library(ggplo2)

my_cols <- c('B'='#C88888','Erythrocyte'='#fe666d','T.alpha.beta'='#B95FBB','T.gamma.delta'='#3288bd','NK'='#bb8761','Hepatocyte'='#63636d','Cholangiocyte'='#de77ae',
'Endothelial'='#D4D915','Fibroblast'='#66c2a5','Macrophage'='#ff9a36','Tumor_core'='#A4DFF2','Tumor_boundary'='blue','Other'='#cccccc')

SpaCET.visualize.spatialFeature(SpaCET_obj, spatialType = "metaData", spatialFeatures = "cellType", colors= my_cols, pointSize = 0.1) + guides(color=guide_legend(override.aes = list(size=3)))

推断分泌蛋白活性

运行 SecAct.activity.inference.ST 为每个spot推断分泌蛋白的活性 (内置数据里面的1000多个分泌蛋白)。结果保存在SpaCET_obj@results$SecAct_output$SecretedProteinActivity，包含以下四项：1) beta：回归系数；2) se：标准误；3) zscore：beta/se；4) pvalue：双侧置换检验p值。

SpaCET_obj <- SecAct.activity.inference.ST(inputProfile = SpaCET_obj, scale.factor = 1000, sigFilter = T)
SpaCET_obj@results$SecAct_output$SecretedProteinActivity$zscore[1:5, 1:3]
       c_2_100_733 c_2_101_240 c_2_101_339
ACE    -1.85461220   -4.425002   -4.400708
ADIPOQ  0.41980845   -3.286004   -4.155328
ADM2    0.06245601    4.793775    5.301861
ALCAM   0.71253379    3.315182    3.481940
ANGPT1 -0.44917882   -3.697340   -4.605047

细胞通讯推断

在计算出分泌蛋白活性后，SecAct 还能进一步从整个组织切片中推断出的这些信号活动来估算出共识模式。此模块主要包含以下两个过程：

筛选显著分泌蛋白：计算蛋白信号活性与其邻居细胞 RNA 表达总和的Spearman 相关性，BH校正；阈值 r > 0.05 且 FDR < 0.01。
检验分泌蛋白是否介导两个细胞类型间的细胞通讯：首先，我们根据空间邻近性 (距离 < 20 μm) 确定了所有相邻的细胞对，其中一对细胞分别来自细胞类型 1 和细胞类型 2 。在这些相邻细胞对中，将能够进行通讯的细胞对定义为：细胞类型1表达信号分子s的RNA (count>0)，且细胞类型2对信号分子s表现出信号活性 (activity score > 0)。随后，我们计算了细胞类型1→s→细胞类型2的通讯评分，该评分被定义为通讯细胞对数量与总相邻细胞对数量的比值。p值表示随机比值大于或等于实际比值的概率，通过1000次随机化计算得出并由Benjamini–Hochberg校正。在每次随机化中，对每种细胞类型的细胞ID进行随机置换。若ratio >0.2 & FDR<0.01时，则细胞类型1→s→细胞类型2的通讯具有统计学意义。

SpaCET_obj <- SecAct.CCC.scST(SpaCET_obj, cellType_meta = "cellType", scale.factor = 1000, radius = 0.02, ratio_cutoff = 0.2, padj_cutoff = 0.01)
head(SpaCET_obj@results$SecAct_output$SecretedProteinCCC, 3)
                                        sender secretedProtein       receiver sender_count receiver_count neighboringCellPairs communicatingCellPairs     ratio          pv      pv.adj
Tumor_core_APOA1_Tumor_boundary     Tumor_core           APOA1 Tumor_boundary       318078          71718                47724                  17244 0.3613276 0.000999001 0.001801119
Tumor_boundary_APP_Tumor_core   Tumor_boundary             APP     Tumor_core        71718         318078                47724                  16814 0.3523175 0.000999001 0.001801119
Tumor_core_CCL15_Tumor_boundary     Tumor_core           CCL15 Tumor_boundary       318078          71718                47724                  10048 0.2105440 0.000999001 0.001801119

可视化细胞通讯

细胞通讯的可视化方式有两种：热图与圈图。热图里面的数字表示从发送者到接收者的显著分泌蛋白数量。

SecAct.CCC.heatmap(SpaCET_obj, row.sorted=T, column.sorted=T, colors_cellType=my_cols)

SecAct.CCC.circle(SpaCET_obj, colors_cellType=my_cols)

可以从对象里面选择感兴趣的细胞通讯并进行可视化。用SecAct.CCC.dot可视化时需设定发送者、分泌蛋白和接收者参数：

cellTypes <- c("Tumor_boundary", "Fibroblast", "Macrophage", "Endothelial")
secretedProtein <- c("BGN","COL1A1","COL1A2","DCN","IGFBP5","LGALS1","LGALS9","LYZ","LUM","MGP","SPP1","THBS1","THBS2")
SecAct.CCC.dot(SpaCET_obj, colors_cellType = my_cols, sender = cellTypes, secretedProtein = secretedProtein, receiver = cellTypes)

可视化 signaling velocity

用SecAct.signaling.velocity.scST可以展示感兴趣的细胞通讯的signaling velocity情况。作者居然在函数内部直接引用了变量my_cols，所以可以通过该变量定义一个命名向量来设置颜色：

my_cols <- c(Fibroblast='green',Tumor_boundary='blue',Other='grey')
SecAct.signaling.velocity.scST(SpaCET_obj, sender = "Fibroblast", secretedProtein = "THBS2", receiver = "Tumor_boundary", cellType_meta = "cellType") + guides(color=guide_legend(override.aes = list(size=3)))

空间转录组的细胞通讯分析