在看文獻Single-cell RNA sequencing of blood antigen-presenting cells in severe COVID-19 reveals multi-process defects in antiviral immunity的時候看到一個沒用過的轉(zhuǎn)錄因子預測方法DoRothEA滋饲。
Fig.6b, Heatmap of top 50 highly variable TF activities among the three severity groups; the z-scores of TF activities are colour-coded
DoRothEA是一種包含轉(zhuǎn)錄因子（TF）與其靶標相互作用的基因集。一個TF及其對應靶點的集合被定義為調(diào)節(jié)子（regulons）劫乱。DoRothEA regulons 收集了文獻膊畴，ChIP-seq peaks掘猿，TF結(jié)合位點基序，從基因表達推斷相互作用等不同類型的互作證據(jù)唇跨。TF和靶標之間的互作可信度根據(jù)支持的證據(jù)數(shù)量劃分為A-E五個等級稠通，A是最可信，E為可信度低买猖。
DoRothEA可以用于bulk RNAseq和scRNAseq的數(shù)據(jù)改橘。
DoRothEA regulon可以與幾種統(tǒng)計方法結(jié)合使用，從而產(chǎn)生一種功能分析工具玉控，以從基因表達數(shù)據(jù)推斷TF活性飞主。通過不考慮TF本身的基因表達，而是考慮其直接轉(zhuǎn)錄靶標的mRNA水平來計算活性高诺。
R包安裝和載入
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("dorothea")
## We load the required packages
library(dorothea)
library(dplyr)
library(Seurat)
library(tibble)
library(pheatmap)
library(tidyr)
library(viper)
導入注釋好的pbmc數(shù)據(jù)集
pbmc <- readRDS("pbmc.rds")
DimPlot(pbmc,label = T,repel = T)
計算細胞的TF活性
## We read Dorothea Regulons for Human:
dorothea_regulon_human <- get(data("dorothea_hs", package = "dorothea"))
##如果是小鼠碌识，就用
##dorothea_regulon_mouse <- get(data("dorothea_mm", 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"))
## We compute Viper Scores
pbmc <- run_viper(pbmc, regulon,
options = list(method = "scale", minsize = 4,
eset.filter = FALSE, cores = 1,
verbose = FALSE))
這一步之后，assays中除了"RNA"以外懒叛，多了一個"dorothea"丸冕。
266x2638的矩陣是266個TFx2638個細胞
隨后可以用TFx細胞的矩陣對細胞進行重新聚類，方法和用GENEx細胞的矩陣做聚類一樣薛窥。
## We compute the Nearest Neighbours to perform cluster
DefaultAssay(object = pbmc) <- "dorothea"
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc, features = rownames(pbmc), verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:10, verbose = FALSE)
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:10, umap.method = "uwot", metric = "cosine")
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25,
logfc.threshold = 0.25, verbose = FALSE)
## Assigning cell type identity to clusters
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T",
"FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
轉(zhuǎn)錄因子的聚類圖
每個細胞群的TF活性（相當于每個細胞群的bulk RNAseq）
## We transform Viper scores, scaled by seurat, into a data frame to better
## handling the results
viper_scores_df <- GetAssayData(pbmc, slot = "scale.data",
assay = "dorothea") %>%
data.frame(check.names = F) %>%
t()
## 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) #也可以使用其他的分類信息
## 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))
根據(jù)前一步計算的score胖烛，選擇在細胞群間變化最大的20個TFs進行可視化
## 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)
## 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)
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)))
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)
