今天学习另外一个转录因子活性预测工具：DoRothEA。很多文章用的也是这个工具，看起来比pySCENIC好像快很多。其实，单细胞转录因子分析本质上是计算一种specific的得分，DoRothEA计算的是 Viper 得分。

======安装=====

官网地址是：

https://saezlab.github.io/dorothea/

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("dorothea")

或者用devtools安装：
devtools::install_github("saezlab/dorothea")

还需要一个依赖包：
BiocManager::install("viper")

=====例子测试====

还是用经典的pbmc的例子。

加载需要的库：
library(Seurat)
library(dorothea)
library(tidyverse)
library(viper)

把pbmc的数据load进来，分析和前面的很多例子一样

pbmc <- readRDS("pbmc.rds")
str(pbmc@meta.data)

下面，加载这个包自带的human的数据库

## We read Dorothea Regulons for Human:
dorothea_regulon_human <- get(data("dorothea_hs", package = "dorothea"))

## We obtain the regulons based on interactions with confidence level A, B and C
regulon <- dorothea_regulon_human %>% dplyr::filter(confidence %in% c("A","B","C"))

然后计算TF的活性，通过使用函数 run_viper() 在 DoRothEA 的regulons上运行 VIPER 以获得 TFs activity。 在 seurat 对象的情况下，该函数返回相同的 seurat 对象，其中包含一个名为 dorothea 的assay，其中包含slot数据中的 TFs activity。

pbmc <- run_viper(pbmc, regulon,options = list(method = "scale", minsize = 4, 
                  eset.filter = FALSE, cores = 1, verbose = FALSE))
Assays(pbmc)

image.png

可以看到，这个pbmc数据集里面的约 三千个细胞都有自己的266个转录因子的活性得分。

image.png

下面，我们试试看FindAllMarkers函数获取各个单细胞亚群特异性的转录因子。

DefaultAssay(object = pbmc) <- "dorothea"
table(Idents(pbmc))

image.png

pbmc.markers <- FindAllMarkers(object = pbmc, 
                              only.pos = TRUE, 
                              min.pct = 0.25, 
                              thresh.use = 0.25)                                                          
DT::datatable(pbmc.markers)
write.csv(pbmc.markers,file='pbmc.markers.csv')

下面，我们尝试热图展示：

library(dplyr) 
pbmc.markers$fc = pbmc.markers$pct.1 - pbmc.markers$pct.2
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(10, fc)
pbmc@assays$dorothea@data[1:4,1:4]
top10 =top10[top10$fc > 0.3,] 
pbmc <- ScaleData(pbmc)
DoHeatmap(pbmc,top10$gene,size=3,slot = 'scale.data')

从热图可以看出，其实还是有一大部分TF在不同的细胞类型中特异性表达的。

image.png

然后，我们尝试利用这个包来找下特异性比较好的TF regulons。官网有详细的文档介绍，比如如何选择TF activity per cell population，根据 previously computed VIPER scores on DoRothEA’s regulons.

首先也是获取Viper得分矩阵 ，根据不同细胞亚群进行归纳汇总 ：

viper_scores_df <- GetAssayData(pbmc, slot = "scale.data", assay = "dorothea") %>%
data.frame(check.names = F) %>% t()
viper_scores_df[1:4,1:4]

image.png

## We create a data frame containing the cells and their clusters
CellsClusters <- data.frame(cell = names(Idents(pbmc)), 
                            cell_type = as.character(Idents(pbmc)),
                            check.names = F)
head(CellsClusters)

## We create a data frame with the Viper score per cell and its clusters
viper_scores_clusters <- viper_scores_df  %>%
  data.frame() %>% 
  rownames_to_column("cell") %>%
  gather(tf, activity, -cell) %>%
  inner_join(CellsClusters)
  
  
## We summarize the Viper scores by cellpopulation
summarized_viper_scores <- viper_scores_clusters %>% 
  group_by(tf, cell_type) %>%
  summarise(avg = mean(activity),
            std = sd(activity))
            
# For visualization purposes, we select the 20 most variable TFs across clusters according to our scores.
head(summarized_viper_scores)


## We select the 20 most variable TFs. (20*9 populations = 180)
highly_variable_tfs <- summarized_viper_scores %>%
  group_by(tf) %>%
  mutate(var = var(avg))  %>%
  ungroup() %>%
  top_n(180, var) %>%
  distinct(tf)
highly_variable_tfs

## We prepare the data for the plot
summarized_viper_scores_df <- summarized_viper_scores %>%
  semi_join(highly_variable_tfs, by = "tf") %>%
  dplyr::select(-std) %>%   
  spread(tf, avg) %>%
  data.frame(row.names = 1, check.names = FALSE) 

                  

最后，整理好的数据如下所示：

summarized_viper_scores_df[1:4,1:4]

image.png

最终数据是挑选的top 20个TF，它们在不同的单细胞亚群里面的变化最大。这样就可以可视化这些细胞类型特异表达的TF regulon了。

有了数据，就很容易可视化：

palette_length = 100
my_color = colorRampPalette(c("Darkblue", "white","red"))(palette_length)
my_breaks <- c(seq(min(summarized_viper_scores_df), 0, 
                   length.out=ceiling(palette_length/2) + 1),
               seq(max(summarized_viper_scores_df)/palette_length, 
                   max(summarized_viper_scores_df), 
                   length.out=floor(palette_length/2)))
library(pheatmap)
viper_hmap <- pheatmap(t(summarized_viper_scores_df),fontsize=14, 
                       fontsize_row = 10, 
                       color=my_color, breaks = my_breaks, 
                       main = "DoRothEA (ABC)", angle_col = 45,
                       treeheight_col = 0,  border_color = NA) 

image.png