pyscenic的转录因子分析结果展示之各个单细胞亚群特异性激活转录因子
前面的笔记： pyscenic的转录因子分析结果展示之5种可视化 带领大家回顾了一下 单细胞转录因子分析之SCENIC流程 ，并且重新认识了 使用pyscenic做转录因子分析 后的结果。
我们根据pbmc3k数据集里面的b细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，进行了可视化。其实这两个转录因子并不是先验知识，是我们根据这个分析结果进行各个单细胞亚群特异性激活转录因子统计得到的。
让我们来看看这个统计分析如何做，以及如何更好的可视化。如果你确实计算资源有限制，其实看看 各个单细胞亚群特异性的转录因子热图 ，也是很容易理解，并不一定要 使用pyscenic做转录因子分析 哦。
首先再次回顾一下pyscenic的转录因子分析结果
代码如下所示：
###### step0 加载 各种R包 #####
rm(list=ls())
library(Seurat)
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
library(dplyr)
library(KernSmooth)
library(RColorBrewer)
library(plotly)
library(BiocParallel)
library(grid)
library(ComplexHeatmap)
library(data.table)
load( file = 'for_rss_and_visual.Rdata')
head(cellTypes)
sub_regulonAUC[1:4,1:2]
可以看到每个单细胞都标记好了生物学名字，而且配套了一个转录因子活性矩阵，自己去看 前面的 使用pyscenic做转录因子分析 教程， 拿到 for_rss_and_visual.Rdata 的文件哦 ：
> head(cellTypes)
celltype
AAACATACAACCAC-1 Naive CD4 T
AAACATTGAGCTAC-1 B
AAACATTGATCAGC-1 Memory CD4 T
AAACCGTGCTTCCG-1 CD14+ Mono
AAACCGTGTATGCG-1 NK
AAACGCACTGGTAC-1 Memory CD4 T
> sub_regulonAUC[1:4,1:2]
AUC for 4 regulons (rows) and 2 cells (columns).
Top-left corner of the AUC matrix:
cells
regulons AAACATACAACCAC-1 AAACATTGAGCTAC-1
ARNTL(+) 0.08163265 0.05424823
ASCL2(+) 0.21191691 0.28589650
ATF1(+) 0.04203110 0.04387755
ATF2(+) 0.24141561 0.19748137
> dim(sub_regulonAUC)
[1] 208 2638
值得一提的是这个pbmc3k数据集的两千多个细胞，其实就两百多个转录因子结果哦。
看看不同单细胞亚群的转录因子活性平均值
selectedResolution <- "celltype" # select resolution
# Split the cells by cluster:
cellsPerGroup <- split(rownames(cellTypes),
cellTypes[,selectedResolution])
sub_regulonAUC <- sub_regulonAUC[onlyNonDuplicatedExtended(rownames(sub_regulonAUC)),] # 去除extened regulons
dim(sub_regulonAUC)
# Calculate average expression:
regulonActivity_byGroup <- sapply(cellsPerGroup,
function(cells)
rowMeans(getAUC(sub_regulonAUC)[,cells]))
# Scale expression:
regulonActivity_byGroup_Scaled <- t(scale(t(regulonActivity_byGroup),
center = T, scale=T))
# 同一个regulon在不同cluster的scale处理
dim(regulonActivity_byGroup_Scaled)
regulonActivity_byGroup_Scaled=regulonActivity_byGroup_Scaled[]
regulonActivity_byGroup_Scaled=na.omit(regulonActivity_byGroup_Scaled)
pheatmap:pheatmap(regulonActivity_byGroup_Scaled)
如下所示的简单热图 ：
简单热图
可以看到，确实每个单细胞亚群都是有 自己的特异性的激活的转录因子，仍然是推荐看看 各个单细胞亚群特异性的转录因子热图 。
不过，SCENIC包自己提供了一个 calcRSS函数，帮助我们来挑选各个单细胞亚群特异性的转录因子，全称是：Calculates the regulon specificity score
参考文章：The RSS was first used by Suo et al. in: Revealing the Critical Regulators of Cell Identity in the Mouse Cell Atlas. Cell Reports (2018). doi: 10.1016/j.celrep.2018.10.045
运行超级简单。
rss <- calcRSS(AUC=getAUC(sub_regulonAUC),
cellAnnotation=cellTypes[colnames(sub_regulonAUC),
selectedResolution])
rss=na.omit(rss)
rssPlot <- plotRSS(rss)
plotly::ggplotly(rssPlot$plot)
如下所示，可以看到我们前面的笔记： pyscenic的转录因子分析结果展示之5种可视化 ，提到的 pbmc3k数据集里面的b细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，确实是在b细胞里面名列前茅哦：
挑选各个单细胞亚群特异性的转录因子
大家可以把每个单细胞亚群各自的独特激活的转录因子按照我们前面的笔记： pyscenic的转录因子分析结果展示之5种可视化 ，拿去安装进行可视化哦。
比如血小板的PRDM1这个转录因子，你会发现，不同的可视化方法只有结合起来你才能对数据有一个更直观的感受。
其它可视化（基本都是展现不同细胞亚群特异性转录因子）
我们也多次介绍过：
- SCENIC转录因子分析结果的解读
- 单细胞转录因子分析之SCENIC流程 （R语言版本，代码运行超级慢）
经过了前面的步骤，我们顺利挑选到了各个单细胞亚群特异性的转录因子列表，就可以自己去提取它们在每个细胞的活性值，自己绘制热图或者其它图表，大同小异的。
下面简单的分享一下我自己的解决方案，
library(dplyr)
rss=regulonActivity_byGroup_Scaled
head(rss)
library(dplyr)
df = do.call(rbind,
lapply(1:ncol(rss), function(i){
dat= data.frame(
path = rownames(rss),
cluster = colnames(rss)[i],
sd.1 = rss[,i],
sd.2 = apply(rss[,-i], 1, median)
df$fc = df$sd.1 - df$sd.2
top5 <- df %>% group_by(cluster) %>% top_n(5, fc)
rowcn = data.frame(path = top5$cluster)
n = rss[top5$path,]
#rownames(rowcn) = rownames(n)