空间通过细胞空间排布分析银屑病的严重程度

作者，Evil Genius

今天我们复习，文章在Spatial transcriptomics stratifies psoriatic disease severity by emergent cellular ecosystems，2023.6.2发表于Science Immunology,IF 30.63，样本类型：健康的皮肤和银屑病。其中在研究的过程就是借助了其他人的单细胞数据辅助研究自己的空间转录组数据，思路值得借鉴。

这个文献之前分享过在如何利用数据库的单细胞数据分析自有的空间转录组，包括在空转系列课上也重点讲过，其中我们需要学习的地方

1、空间整合分析如何判断使用CCA还是harmony？ 2、空间整合合理的结果是如何判断的？ 3、MIA的用法 4、如何借助公共的单细胞数据辅助分析自己的空间样本？也就是MIA在数据不匹配的情况下的实际运用 5、SpaceFold + BayesPrism的联合用法

我们来回顾一下

• 文章运用空间转录组学(ST)分析了25例健康、活动性病变和临床未受损伤的皮肤活检，并与公开的单细胞转录组学数据相结合，揭示了健康和炎症皮肤之间免疫微孔的显著差异。
• 轻度和严重形式的PsO（人类牛皮癣）具有不同的分子特征，严重的PsO可能会深刻地改变远端未受影响的皮肤部位的细胞和代谢组成。

背景知识

• 人体皮肤构成了对周围环境的物理和免疫屏障。
• 牛皮癣(PsO)是最普遍的免疫介导的炎症性疾病之一，影响人口的2%至3%。PsO的临床特征是皮肤红斑斑块，这些斑块要么局限于一小块皮肤区域，要么覆盖身体的大部分，并且在厚度(硬结)、鳞屑(脱屑)和发红(红斑)方面变化
• 空间结构数据和健康和病理细胞相互作用发生的独特环境或“微生态”是具有挑战性的。考虑到细胞组织和相互作用组在协调健康和疾病中的重要性，在组织尺度上未知地绘制这些特征可能会揭示银屑病严重程度和/或其系统性表现的新特征。

结果1、健康人皮肤的空间图谱

After Seurat anchor–based integration, we implemented graphbased clustering of the healthy skin spot transcriptome. Dimensionality reduction, performed using uniform manifold approximation and projection (UMAP), revealed 12 distinct clusters. Cluster annotations based on differentially expressed (DE) marker genes had notable concordance with histopathologic tissue features. Clusters were grouped by the skin’s three primary histologically defined layers (epidermis, dermis, and hypodermis) and captured discrete structures such as the hair follicle Harmony-based integration and graph-based clustering analysis yielded 11 distinct clusters and similarly delineated the distinct skin layers. Serial analysis of tissue sections yielded highly reproducible results, and each of the 12 clusters was represented in all healthy samples when using both anchor-based and Harmony-based integration analyses. We thus performed subsequent validations and analyses of healthy skin using anchor-based integration. Across all samples, epidermal and hair follicle clusters 3, 7, and 9 had the highest number of UMIs, followed by cluster 6 (upper follicle and perifollicular dermis) and dermis 2, highlighting regions that were the most cell dense and/or transcriptionally active

• 利用公开可用的scRNA-seq数据获得ST spot内的细胞分辨率并验证数据分析。
• 为了进一步证实细胞类型的空间分布，使用了多模态交叉分析(MIA)，这是一种计算方法，通过超几何测试，根据每个ST和scRNA-seq cluster之间上调DE/标记基因的重叠程度，推断组织区域中特定细胞类型的富集。
• 一旦数据存在严重的不合理现象，需要人工矫正，单细胞数据如果存在细胞类型缺失的情况下也要矫正。

结果2、滤泡周围免疫生态位监测健康皮肤

• 表征健康皮肤的免疫图谱，使用的方法是MIA，以确定ST cluster中免疫细胞的富集程度。
• DC淋巴细胞介导的对健康皮肤毛囊和血管周围空间的免疫监视。

结果3、银屑病病变中表皮-真皮相互作用失调

• 再次比较了Harmony和Seurat基于锚点的整合方法，然后对健康和银屑病病变和非病变皮肤spot转录组进行基于图的聚类和UMAP可视化，以生成区域特异性聚类。根据DE基因和组织学特征对Harmony和anchor集成umap中的cluster进行了注释。对两个umap的交叉检验表明，Harmony整合的cluster具有更大的结构特异性，因此使用Harmony样本整合执行下游分析。根据表皮、真皮层、毛囊皮脂腺和脂肪区域对空间spot cluster进行分组，并追踪病变和非病变皮肤的关键组织病理学特征。
• 比较不同样本组之间cluster的富集差异，基因特征，通讯，代谢等差异。

• NicheNet分析空间的通讯差异。

结果4、PsO病变和非病变皮肤中免疫环境和成纤维细胞的动态重组

• 使用MIA，确认了关键的表皮和真皮结构细胞与组织学组织结构一致
• 使用MIA绘制了免疫细胞的空间分布。包括银屑病病变性和非病变性ST簇中的组织常驻记忆T细胞(Trm)
• 作为MIA的补充策略，使用反卷积、积分和投影方法BayesPrism和SpaceFold，沿着健康、病变和非病变皮肤的空间维度绘制数据集中的结构细胞和免疫细胞。
• 这种组织尺度的计算制图进一步强调了非病变皮肤和健康皮肤之间的差异，表明即使是离临床明显疾病远端的皮肤也表现出细胞和分子的改变。

结果5、无监督聚类揭示了银屑病严重程度的分层

• 利用PASI评分测量的皮肤病严重程度的先验知识和伴有的全身、皮外疾病[活动性皮肤银屑病病变和关节炎
• 分析在第一个主成分上观察到转录组变异的梯度，从病变的中度到重度皮肤，到病变的轻度皮肤，到非病变的中度到重度皮肤，最后到与健康皮肤重叠的非病变轻度皮肤。

结果6、轻度和中度至重度疾病的区别特征

• 对按轻度或中度至重度疾病分组的病变和非病变样本进行BayesPrism和SpaceFold投影。每种细胞类型的定位沿着真皮-表皮组织轴(0到1假尺度)进行投影，并使用黑素细胞作为基底层的细胞标记，我们用虚线表示分离表皮和真皮的基膜。与合并的病变皮肤分析类似，一维空间折叠投影显示病变皮肤的表皮增生跨越疾病严重程度的谱。

• 不同的分子特征分层轻度和重度PsO，临床未受损伤的皮肤在严重疾病中发生改变

总结

研究强调了ST在健康皮肤中同时绘制基因表达图谱和定义不同的、生物学相关的细胞“邻域”的能力。在这样做的过程中，发现了ST在捕获成熟脂肪细胞(成熟脂肪细胞和皮脂细胞)的基因表达方面的明显优势。分析还表明，免疫细胞在健康皮肤的滤泡周围和血管周围的niche中富集，而在炎症皮肤中，这种病灶定位明显改变。PsO病变在免疫细胞空间分布上有明显的变化，从健康皮肤中主要分布在毛囊周围，到银屑病皮肤的表面到真皮中部。因此，除了活动性PsO病变中免疫群体的组成变化外，它们的位置可能为了解它们的调节和功能提供了急需的见解。例如，分析发现稳态效应T细胞(TH/TC)和调节性T细胞(Tregs)与apc聚集在血管周围的滤泡周围区域。卵泡是共生微生物的家园，是病原体和过敏原进入的门户，并丰富表达免疫调节因子。小鼠研究表明，稳态免疫监视依赖于居住在卵泡中的共生体。了解健康和患病人体皮肤中特定的共生细菌种群以及它们如何影响淋巴细胞功能是积极研究的途径。 长期以来，皮肤免疫耐受的丧失被认为是导致PsO发生和延续的原因。因此，除了炎症细胞的涌入外，卵泡生态位周围的稳态免疫监视和耐受性机制的紊乱也可能是病变形成的关键诱因。与这一观点相一致的是，全球转录组学分析显示，与严重疾病相比，轻度PsO病变富含免疫调节特征。在临床上，已知头皮上的银屑病病变(人体毛囊最密集的皮肤)更难以治疗，并且更有可能累及全身。此外，在头皮和非头皮PsO中都观察到毛囊脱落、缩小和皮脂腺萎缩。卵泡上皮周围效应/调节性T细胞的战略性定位和比例导致(或防止)皮肤炎症的机制值得进一步研究。 将ST应用于活跃的银屑病病变和来自未受损伤的非病变皮肤的活检，能够在组织尺度上发现细胞组织和细胞间通讯的扭曲。通过这样做，说明了这个平台的价值，将公开可用的皮肤单细胞转录图谱放入其原生组织环境中。研究结果确定了两个表皮和一个真皮空间簇在病变皮肤横截面上显著富集(在某些情况下独特表达)。这些cluster表达了参与炎症(例如TH17激活)和表皮增殖和分化途径的基因，这些途径也在表皮修复过程中出现。长期以来，PsO被认为是一种伤口过度愈合的病理。在所有样本中都有pso富集cluster，因此可能代表疾病的“核心特征”，无论其严重程度或是否存在合并症。先前的工作报道了IL-17A在轻度疾病参与者样本中的表达与严重疾病相比的矛盾增加。尽管如此，IL-17阻断在临床上对轻度和严重疾病都有效，验证了既定的概念，即驱动整个疾病严重程度谱的机制沿着17型免疫轴聚集。 研究受到Visium ST技术目前可用分辨率的限制，该技术以50 μm的增量捕获基因表达。然而，我们整合了scRNA-seq数据集，以获得更高的细胞分辨率和空间背景。随着这项技术的发展，具有更高密度和更高分辨率的空间条形码珠的平台将提供更多关于健康和患病状态下细胞微环境的粒度。另一个关键的考虑因素是测序的深度和基因的辍学率。因此，我们的分析也可能没有捕获到罕见的细胞类型或非常低表达的基因。因此，我们避免了“缺乏证据”的争论，只关注基因的积极富集来得出结论。SpaceFold对组织结构做了两个假设:组织应该具有典型的结构，每个点的细胞类型分数表明这种细胞类型的物理坐标，因为它与起源点有关。因此，异质组织区域内斑点的坐标可能与物理细胞类型位置不对应，因此不能直接解释。然而，不同细胞类型分布区域之间的相对坐标仍然可以与使用组织细胞标志(例如，lc和黑素细胞区分基底表皮层和基底膜)进行比较来得出结论。最后，与病变样本相比，一些空间簇(0、1、2、3、8、15和16)在非病变皮肤中趋向于富集;然而，在这些cluster中观察到的样本间差异以及分析的低样本量使得这些差异在统计上不显著。

Collectively, our findings underscore the value of spatially resolved gene profiling in understanding emergent cellular ecosystems underlying health and disease. Thus, in addition to unearthing previously undescribed disease features, our data will serve as an invaluable resource for the research community.

方法细节

Anchor-based integration

After quality control and sample filtering, ST samples were normalized independently using variance-stabilizing transformation by running the SCTransform function in the Seurat package with default parameters as described in the vignette, setting parameter variable.features.n to 3000, i.e., the top 3000 most variable genes ranked by residual variance were used for downstream integration and cluster analysis. The top 20 PCs with a resolution of 0.6 were used to identify spatial clusters by running the FindNeighbours and FindClusters functions. The clusters were manually annotated on the basis of cluster-specific marker genes identified using the FindAllMarkers function with default parameters, including statistical test set to Wilcoxon rank-sum test, minimum cell percentage set to 0.25, and log2FC cutoff set to 0.25. Heatmap plots were generated using DoHeatmap, showing the top eight marker genes per cluster

Harmony-based integration

Both groups of data were first merged into a single Seurat object using the merge function and were normalized using the SCTransform function using default parameters, setting parameter variable.features.n to 3000, i.e., using the top 3000 most variable genes ranked by residual variance. PCA was performed on the normalized data as needed before running Harmony batch correction. Harmony batch correction was performed by running the RunHarmony Function from the Harmony package using default parameters, with assay set to SCT counts and group variable set to sample ID column in the metadata. The downstream analysis of batch-corrected data was performed using Seurat, and the top 40 dimensions from Harmony embeddings were used for UMAP-based dimensionality reduction and graph-based clustering, with the resolution set to 0.35. The clusters were then manually annotated using cluster-specific marker genes after running the FindAllMarkers function using default parameters, including a Wilcoxon rank-sum statistical test, minimum cell percentage set to 25%, and log2FC cutoff set to 0.25. Heatmaps were generated using DoHeatmap, showing the top 10 marker genes per cluster

MIA方法

MIA uses a hypergeometric test to assess overlap in marker genes between each ST cluster and scRNA-seq cluster . The top 300 up-regulated marker genes with an adjusted P value less than or equal to 0.05 and log fold change greater than 0.25 were used to calculate MIA enrichment scores. MIA uses a hypergeometric test to assess overlap in marker genes between each ST cluster and scRNA-seq cluster. The top 300 up-regulated marker genes with an adjusted P value less than or equal to 0.05 and log fold change greater than 0.25 were used to calculate MIA enrichment scores.

最后，附上MIA的分析代码,略有更新

## MULTI-MODAL INTERSECTION ANALYSIS

MIA_ENRICH <- function(stlist,sclist,total){
  overlap <- length(intersect(stlist$gene,sclist$gene))
  C <- length(sclist$gene)
  D <- length(stlist$gene)
  
  # ENRICHMENT CALCULATION
  e <- -log10(phyper(overlap,C, total - C, D, lower.tail = FALSE))
  return(e)
}
MIA_DEPLETE <- function(stlist,sclist){
  overlap <- length(intersect(stlist$gene,sclist$gene))
  C <- length(sclist$gene)
  D <- length(stlist$gene)
  
  # ENRICHMENT CALCULATION
  d <- -log10(1- (phyper(overlap,C, total - C, D, lower.tail = FALSE)))
  return(d)
}
MIA <- function(total_genes,single_cell.markers,spatial.markers)
{
  #D.SCORES <- c()
  # Perform this operation for every cell type
  #single_cell.markers <- FindAllMarkers(single_cell,assay = assay_use,logfc.threshold = 0.25,return.thresh = p_val_adj < 0.1)
  #spatial.markers <- FindAllMarkers(spatial_data,assay = assay_use,logfc.threshold = 0.25,return.thresh = p_val_adj < 0.1)
  cell.types <- single_cell.markers %>% dplyr::select(cluster) %>% unique() %>% as.list()
  spatial.regions <- spatial.markers %>% dplyr::select(cluster) %>% unique() %>% as.list()
  E.SCORES <- data.frame(spatial.regions)
  for(i in cell.types){
    for (x in i){
      e_list <- c()
    #list.append(e_list,i)
      for(y in spatial.regions){
        for(z in y){
          single_cell <- single_cell.markers %>% filter(cluster==x)
          spatial_data <- spatial.markers %>% filter(cluster==z)
          e <- MIA_ENRICH(single_cell,spatial_data,total = total_genes)
          #d <- MIA_DEPLETE(single_cell,spatial_data,total = total_genes)
          e_list <- c(e_list,e)
        }
        #D.SCORES <- append(D.SCORES,d)
      }
      E.SCORES[paste(x)] <- e_list
    #E.SCORES <- append(E.SCORES,e)
    }
  }
  #e.data <- data.frame("GA"=E.SCORES[1],"ER"=E.SCORES[2],"C0L17A1+"=E.SCORES[3])
  #d.data <- data.frame("GA"=D.SCORES[1],"ER"=D.SCORES[2],"C0L17A1+"=D.SCORES[3])
  #res <- list(E.SCORES,e.data)
  return(E.SCORES)
}

MIA_ENRICH_bullk <- function(stlist,bulklist,total){
  overlap <- length(intersect(stlist$gene,bulklist))
  C <- length(bulklist)
  D <- length(stlist$gene)
  
  # ENRICHMENT CALCULATION
  e <- -log10(phyper(overlap,C, total - C, D, lower.tail = FALSE))
  return(e)
}
MIA_bulk <- function(total_genes,markers,spatial.markers,name)
{
  spatial.regions <- spatial.markers %>% dplyr::select(cluster) %>% unique() %>% as.list()
  E.SCORES <- data.frame(spatial.regions)
  e_list <- c()
      #list.append(e_list,i)
    for(y in spatial.regions){
      for(z in y){
        spatial_data <- spatial.markers %>% filter(cluster==z)
        e <- MIA_ENRICH_bullk(spatial_data,markers,total = total_genes)
        e_list <- c(e_list,e)
      }
    }
   E.SCORES[paste(name)] <- e_list
  return(E.SCORES)
}

n_cells = nrow(combo_tbl)
for (ac in integrated_clusters) {
  ac_cells = combo_tbl %>% filter(cluster_integrated == ac) %>% pull(cell)
  n_ac = length(ac_cells)
  for (bc in independent_clusters) {
    bc_cells = combo_tbl %>% filter(cluster_independent == bc) %>% pull(cell)
    n_bc = length(bc_cells)
    n_common = length(intersect(ac_cells, bc_cells))
    n_either = length(union(ac_cells, bc_cells))
    
    # calculate jaccard index
    jaccard_index = n_common / n_either
    # message(glue("{ac} : {bc} : {length(ac_cells)} : {length(bc_cells)} : {jaccard_index}"))
    jaccard_mat[ac, bc] = round(jaccard_index, 3)
    
    # perform hypergeometric test
    ph = 1.0
    if (n_common > 0) {
      ph = phyper(n_common - 1, n_ac, n_cells - n_ac, n_bc, lower.tail = FALSE, log.p = FALSE)
    }
    # round(-log10(ph), 1)
    ph_mat[ac, bc] = round(ph, 3)
    # perform fisher's exact test
    fe = fisher.test(
      matrix(c(n_common,
               n_ac - n_common,
               n_bc - n_common,
               n_cells - n_ac - n_bc + n_common),
             2, 2),
      alternative="greater"
    )$p.value
    fe_mat[ac, bc] = round(ph, 3)
    
  }
}

生活很好，有你更好