ixxmu
commented
2 days ago 单细胞全流程之转录因子活性预测：decoupleR by 爱喝可乐的天天
最近nature在10月23日上线了两篇关于心衰中免疫细胞和成纤维细胞间通讯的文章：
两篇都与经典的炎症因子IL1β相关
第一篇主要是发现了表达成纤维细胞激活蛋白（FAP）和periostin（POSTN）的成纤维细胞亚群，这些细胞与心肌纤维化有关，通过配体-受体分析和空间转录组学预测，CCR2+巨噬细胞和成纤维细胞之间的IL-1β信号传导在空间上定义的生态位中驱动FAP/POSTN成纤维细胞的出现，从而在心脏纤维化中发挥作用。
第二篇文章主要发现了在浸润性Cx3cr1+巨噬细胞中条件性敲除转录共激活因子Brd4可以改善小鼠的心力衰竭并显著减少成纤维细胞的激活，通过单细胞染色质可及性分析和BRD4占据位点的分析，研究人员发现了一个靠近IL-1β的大增强子，以及一系列精确的压力依赖性调控元件，这些元件控制IL-1β的表达，分泌的IL-1β激活了成纤维细胞中依赖于RELA（也称为p65）的增强子，导致人类心脏成纤维细胞中促纤维化反应。
有时间我们再来仔细看看这两篇文章的具体内容。
今天主要想分享的是第一篇文章中用到的DoRothEA转录因子活性预测分析，我也是第一次看到。翻看git之后发现DoRothEA已经被decoupleR取代，我们直接来看下decoupleR的使用方法，这是原文的图：
安装
BiocManager::install("decoupleR")
BiocManager::install("OmnipathR")#获取转录因子数据必须安装
decoupleR还可以做通路活性推断PROGENy，这个分析我们之前已经分享过了，感兴趣的朋友可查看单细胞全流程之Progeny通路活性推断
分析流程
非常简单，获取TF的背景数据，一步运行即可
library(dplyr)
library(Seurat)
library(tibble)
library(pheatmap)
library(tidyr)
library(viper)
library(decoupleR)
library(ggplot2)
library(patchwork)
library("OmnipathR")
sce <- readRDS("sc_AC_anno.rds")
Idents(sce) <- "celltype"
net <- get_collectri(organism='human', split_complexes=FALSE)
mat <- as.matrix(sce@assays$RNA@data)
plan("multisession",workers = 4)#为了提高分析速度，这里设置多个核，大家可以自行修改
acts <- run_ulm(mat, net,minsize = 5,
.source='source', .target='target',.mor='mor')
接着就是把预测结果添加到seurat对象中
sce[['tfsulm']] <- acts %>%
pivot_wider(id_cols = 'source',
names_from = 'condition',
values_from = 'score') %>%
column_to_rownames('source') %>%
Seurat::CreateAssayObject(.)
DefaultAssay(object = sce) <- "tfsulm"
sce <- ScaleData(sce)
sce@assays$tfsulm@data <- sce@assays$tfsulm@scale.data
可以用featureplot函数可视化某一个转录因子的预测结果
p1 <- DimPlot(sce, reduction = "umap",
label = TRUE, pt.size = 0.5) +
NoLegend() + ggtitle('Cell types')
p1
p2 <- (FeaturePlot(sce, features = c("ATF3")) &
scale_colour_gradient2(low = 'blue', mid = 'white', high = 'red')) +
ggtitle('ATF3 activity')
p2
DefaultAssay(object = sce) <- "RNA"
p3 <- FeaturePlot(sce, features = c("ATF3")) +
ggtitle('ATF3 expression')
p3
p1 + p2 + p3
ggsave(filename = "ATF3.pdf",height = 4,width = 12)
从结果我们可以发现在各类巨噬细胞中尽管表达量相差不大但是活性却差异明显。
查看整体的分析结果，以排名前25的转录因子为例
n_tfs <- 25
#提取活性数据并修改数据格式
df <- t(as.matrix(sce@assays$tfsulm@data)) %>%
as.data.frame() %>%
mutate(cluster = Idents(sce)) %>%
pivot_longer(cols = -cluster, #指定除cluster列以外的都需要转换
names_to = "source", #列名放到source列中
values_to = "score") %>% #值放到score列中
group_by(cluster, source) %>% #将数据按cluster和source两列进行分组
summarise(mean = mean(score))
#获得不同簇中活性差异较大的TFs
tfs <- df %>%
group_by(source) %>%
summarise(std = sd(mean)) %>%
arrange(-abs(std)) %>%
head(n_tfs) %>%
pull(source)
#修改数据格式
top_acts_mat <- df %>%
filter(source %in% tfs) %>%
pivot_wider(id_cols = 'cluster',
names_from = 'source',
values_from = 'mean') %>%
column_to_rownames('cluster') %>%
as.matrix()
#绘图
palette_length = 100
my_color = colorRampPalette(c("Darkblue","white","red"))(palette_length)
my_breaks <- c(seq(-3, 0, length.out=ceiling(palette_length/2) + 1),
seq(0.05, 3, length.out=floor(palette_length/2)))
#显示
pheatmap(top_acts_mat,
border_color = NA,
color=my_color,
breaks = my_breaks)
ggsave(filename = "top25.pdf",height = 4,width = 8)
到此就结束了，可以尝试一下哈~
https://mp.weixin.qq.com/s/402Ov0GwaJh6gm2KrrRDpw