10X单细胞（10X空间转录组）转录因子活性分析之DoRothEA
发布时间
阅读量：
阅读量
hello
很高兴与大家分享今天为大家介绍一款专为转录因子活性预测而设计的工具软件DoRothEA。它广泛应用于多篇高水平的研究文献中，请大家参考[DoRothEA](https://links.jianshu.com/go?to=https://bioconductor.org/packages/release/data/experiment/vignettes/dorothea inst/doc/single_cell_vignette.html#introduction)。
先来看看介绍
首先介绍数据库DoRothEA的概念：DoRothEA是一种由转录因子及其靶标相互作用组成的基因集合。每个转录因子与其对应的靶点集合被定义为调节子（regulons）。通过整合来自文献中的数据、ChIP-seq峰值以及转录因子结合位点基准等信息，并基于基因表达结果推断出的各种交互证据来构建DoRothEA所包含的regulons集。为了量化这些交互的重要程度，我们根据支持的证据数量将TF与靶标的互作可信度划分为A-E五级等级划分依据的是支持的证据数量
TF 的活动度是基于其目标基因的 mRNA 表达水平计算得出的。 由此可见，可将其视为反映当前转录活性水平的重要指标。
看一看代码案例
安装和加载
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为例)
## Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
## Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
min.cells = 3, min.features = 200)前处理（可选，如果读取的rds已经做过处理，这一步就不需要了）
## Identification of mithocondrial genes
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
## Filtering cells following standard QC criteria.
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 &
percent.mt < 5)
## Normalizing the data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize",
scale.factor = 10000)
pbmc <- NormalizeData(pbmc)
## Identify the 2000 most highly variable genes
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
## In addition we scale the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
降维聚类（可选，Seurat的方法，通常我们前面都已经分析过了）
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = 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()
图片.png
评估细胞中的TF活性水平，在案例研究中我们采用了run_viper()包装函数，在DoRothEA的regulons上应用VIPER算法来评估TF活性。该功能能够处理多种类型的输入数据格式：如矩阵、数据框、表达式集以及Seurat对象。当输入为Seurat对象时(run_viper会返回一个具有dorothea assay的新对象)此assay记录了slot数据中相应的TF活性信息
## 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"))
## We compute Viper Scores
pbmc <- run_viper(pbmc, regulon,
options = list(method = "scale", minsize = 4,
eset.filter = FALSE, cores = 1,
verbose = FALSE))
随后我们采用 Seurat 按照类似的方式运用 TFs activity分数来进行细胞聚类分析。
## 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()
图片.png
对于每个细胞群中的 TF 活性（相当于 bulk RNAseq），我们基于先前计算出的 DoRothEA regulons 的 VIPER 得分进行评估，并通过分析这些 TF 的活动水平来区分不同的细胞群。
对于每个细胞群中的 TF 活性（相当于 bulk RNAseq），我们基于先前计算出的 DoRothEA regulons 的 VIPER 得分进行评估，并通过分析这些 TF 的活动水平来区分不同的细胞群。
## 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))
选择在细胞群间变化最大的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)
图片.png
感觉还挺好，方便，能说明一些生物学的问题
生活很好，有你更好
全部评论 (0)
还没有任何评论哟~