10X单细胞（10X空间转录组）细胞通讯分析联合大全

阅读量：

作者，追风少年i

周四了，一周最好的时光，不知道大家感觉怎么样了？2022年上半年的最后一天，时光总是匆匆，不会为自己停留片刻。白哥有一句台词，之前做的再好也回不去了，项目交出去的那一刻，便不再属于自己，那如果新的项目做不成，就该考虑换个工作了~~~~你再有辉煌也都是过去的事情了，要是因为过去的事出不来的话，我们就会越来越平庸。

半年了，我们要做总结，希望大家都能在自己擅长的领域有所成就~~~

七种方法的整合

复制代码

 show_methods()

    
 #>  [1] "connectome"      "logfc"           "natmi"           "sca"            
    
 #>  [5] "cellphonedb"     "cytotalk"        "call_squidpy"    "call_cellchat"  
    
 #>  [9] "call_connectome" "call_sca"        "call_italk"      "call_natmi"

其中重点涉及

考虑空间位置的通讯分析手段---CellphoneDB（V3.0）

10X单细胞（10X空间转录组）通讯分析之CytoTalk（从头构建信号转导网络）

数据库包括

复制代码

 # Resource currently included in OmniPathR (and hence `liana`) include:

    
 show_resources()
    
 #>  [1] "Default"          "Consensus"        "Baccin2019"       "CellCall"        
    
 #>  [5] "CellChatDB"       "Cellinker"        "CellPhoneDB"      "CellTalkDB"      
    
 #>  [9] "connectomeDB2020" "EMBRACE"          "Guide2Pharma"     "HPMR"            
    
 #> [13] "ICELLNET"         "iTALK"            "Kirouac2010"      "LRdb"            
    
 #> [17] "Ramilowski2015"   "OmniPath"

当然，我们要进一步扩充，首先来看LIANA与NicheNet的联合分析

LIANA（LIgand-receptor AAnalysis frAmework）是一个框架，能够使用不同的资源和方法对来自单细胞转录组学的配体-受体相互作用进行优先排序。它允许用户系统地生成关于来自给定细胞类型的哪些配体与另一种细胞类型的受体结合的假设。与 LIANA 相比，NicheNet 旨在深化将配体与一组转录靶标连接起来的细胞内机制，从而广泛使用来自多个来源的先验知识。 LIANA 和 NicheNet 并不相互排斥，但在某些情况下可能相当互补，因为它们旨在探索细胞间和细胞内通信的不同方面。

其中Nichenet的参考链接在NicheNet，值得大家深入学习。

复制代码

 library(tidyverse)

    
 library(liana)
    
 library(nichenetr)
    
 library(Seurat)
    
 library(ggrepel)
    
 library(cowplot)
    
 options(timeout=300)

加载数据

复制代码

 # single-cell expression matrix described in Puram et al. 2017

    
 hnscc_expression <- readRDS(url("https://zenodo.org/record/3260758/files/hnscc_expression.rds"))
    
 # model weights
    
 ligand_target_matrix <- readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))

感兴趣的细胞进行下游分析

复制代码

 expression <- hnscc_expression$expression

    
 sample_info <- hnscc_expression$sample_info
    
 colnames(sample_info) <- make.names(colnames(sample_info))
    
  
    
 # filter samples based on vignette's information and add cell type
    
 tumors_remove <-  c("HN10", "HN", "HN12", "HN13", "HN24", "HN7", "HN8", "HN23")
    
 sample_info <- sample_info %>%
    
   subset( !(tumor %in% tumors_remove) & Lymph.node == 0) %>%
    
   # fix some cell type identity names
    
   mutate(cell_type = ifelse(classified..as.cancer.cell == 1, "Tumor", non.cancer.cell.type)) %>%
    
   subset(cell_type %in% c("Tumor", "CAF"))
    
  
    
 # cell ID as rownames
    
 rownames(sample_info) <- sample_info$cell
    
  
    
 # subset expression to selected cells
    
 expression <- expression[sample_info$cell, ]
    
  
    
 # gene set of interest
    
 geneset_oi <- read_tsv(url("https://zenodo.org/record/3260758/files/pemt_signature.txt"), col_types = cols(), col_names = "gene") %>%
    
   pull(gene) %>%
    
   .[. %in% rownames(ligand_target_matrix)]
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-17/hSAcWFGuJNmyn9bgHIdD8Mw2z5aV.png)

Run LIANA

复制代码

 # create seurat object

    
 seurat_object <- Seurat::CreateAssayObject(counts = expm1(t(expression))) %>%
    
   Seurat::CreateSeuratObject(., meta.data = sample_info) %>%
    
   Seurat::NormalizeData()
    
  
    
 # set cell identity to cell type
    
 Idents(seurat_object) <- seurat_object@meta.data$cell_type
    
  
    
 liana_results <- liana_wrap(seurat_object) %>%
    
   liana_aggregate()
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-17/94JIof7piqHDxBuMaE1Q02lWNKnG.png)

By default, LIANA will score the ligand-receptor interactions in all the possible directions within the two cell types of interest. This includes: Autocrine signaling (e.g. CAFs -> CAFs), CAFs -> Tumor cells and Tumor cells -> CAFs. As we are only interested in the CAFs -> Tumor cell direction, we filter the results and visualize the top 50 interactions according to the consensus/aggregate rank across methods. The aggregate rank itself can be interpreted as the significance of preferential enrichment for the interactions.

复制代码

 # filter results to cell types of interest

    
 caf_tumor_results <- liana_results %>%
    
   subset(source == "CAF" & target == "Tumor") %>%
    
   dplyr::rename(ligand=ligand.complex, receptor=receptor.complex)
    
  
    
 # filter results to top N interactions
    
 n <- 50
    
 top_n_caf_tumor <- caf_tumor_results %>%
    
   arrange(aggregate_rank) %>%
    
   slice_head(n = n) %>%
    
   mutate(id = fct_inorder(paste0(ligand, " -> ", receptor)))
    
  
    
 # visualize median rank
    
 top_n_caf_tumor %>%
    
   ggplot(aes(y = aggregate_rank, x = id)) +
    
   geom_bar(stat = "identity") +
    
   xlab("Interaction") + ylab("LIANA's aggregate rank") +
    
   theme_cowplot() +
    
   theme(axis.text.x = element_text(size = 8, angle = 60, hjust = 1, vjust = 1))
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-17/iQaK5LCA0ftHyNrsnjmI62pxPc4W.png)

Run NicheNet using LIANA’s ligands

将 LIANA 与 NicheNet 结合的关键方面是可以使用 LIANA 优先考虑的配体作为 NicheNet 的潜在配体集 。不会评估受体也在受体细胞类型中表达的所有表达配体，我们将仅探索那些被 LIANA 中包含的方法优先考虑的配体。因此，选择形成前面所示相互作用的配体。

复制代码

 # get ligands and filter to those included in NicheNet's ligand-target matrix

    
 ligands <- unique(top_n_caf_tumor$ligand)
    
 ligands <- ligands[ligands %in% colnames(ligand_target_matrix)]
    
  
    
 background_genes <- expression[sample_info$cell[sample_info$cell_type == "Tumor"], ] %>%
    
   apply(2,function(x){10*(2**x - 1)}) %>%
    
   apply(2,function(x){log2(mean(x) + 1)}) %>%
    
   .[. >= 4] %>%
    
   names()

并使用前面提到的 pEMT 基因集执行 NicheNet 预测配体活动

复制代码

 nichenet_activities <- predict_ligand_activities(

    
   geneset = geneset_oi,
    
   background_expressed_genes = background_genes,
    
   ligand_target_matrix = ligand_target_matrix, potential_ligands = ligands
    
 )

可视化

复制代码

 # prepare data for visualization

    
 vis_liana_nichenet <- top_n_caf_tumor %>%
    
   inner_join(nichenet_activities, by = c("ligand" = "test_ligand")) %>%
    
   arrange(pearson) %>%
    
   mutate(ligand = fct_inorder(ligand))
    
  
    
 # prepare NicheNet figure
    
 nichenet_scores_plot <- vis_liana_nichenet %>%
    
   group_by(ligand) %>%
    
   summarize(pearson = mean(pearson)) %>%
    
   ggplot(aes(y = ligand, x = pearson)) +
    
   geom_bar(stat = "identity") +
    
   ggtitle("NicheNet") +
    
   xlab("Pearson's score") +
    
   theme_cowplot() +
    
   theme(axis.text.y = element_blank(),
    
     axis.ticks.y = element_blank(),
    
     axis.title.y = element_blank(),
    
     axis.line.y = element_line(color = "white"),
    
     plot.title = element_text(hjust = 0.5),
    
     axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))
    
  
    
 # prepare LIANA figure
    
 liana_receptor_heatmap <- vis_liana_nichenet %>%
    
   ggplot(aes(y = ligand, x = receptor, fill = aggregate_rank)) +
    
   geom_tile() +
    
   theme_cowplot() +
    
   ggtitle("LIANA") +
    
   ylab("Ligand") + xlab("Receptor") +
    
   theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1),
    
     plot.title = element_text(hjust = 0.5),
    
     panel.grid.major = element_line(colour = "gray", linetype = 2),
    
     legend.position = "left")
    
  
    
 # combine plots
    
 plot_grid(liana_receptor_heatmap, nichenet_scores_plot,
    
       align = "h", nrow = 1, rel_widths = c(0.8,0.3))
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-17/q3vO9ELk2p5aXuSGV7rzjFCIwxDc.png)