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

最新推荐文章于 2024-07-31 19:38:21 发布

追风少年ii

最新推荐文章于 2024-07-31 19:38:21 发布

阅读量380

点赞数 5

文章标签：数据库空间转录组数据分析单细胞

本文链接：https://blog.csdn.net/weixin_53637133/article/details/138837164

版权

作者，追风少年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)]

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()

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))

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))