前面的笔记：pyscenic的转录因子分析结果展示之5种可视化 带领大家回顾了一下 单细胞转录因子分析之SCENIC流程 ，并且重新认识了 使用pyscenic做转录因子分析 后的结果。
我们根据pbmc3k数据集里面的b细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，进行了可视化。其实这两个转录因子并不是先验知识，是我们根据这个分析结果进行各个单细胞亚群特异性激活转录因子统计得到的。
让我们来看看这个统计分析如何做，以及如何更好的可视化。如果你确实计算资源有限制，其实看看 各个单细胞亚群特异性的转录因子热图，也是很容易理解，并不一定要 使用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这个转录因子，你会发现，不同的可视化方法只有结合起来你才能对数据有一个更直观的感受。
我们也多次介绍过：
经过了前面的步骤，我们顺利挑选到了各个单细胞亚群特异性的转录因子列表，就可以自己去提取它们在每个细胞的活性值，自己绘制热图或者其它图表，大同小异的。
下面简单的分享一下我自己的解决方案，
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)
pheatmap(n,
annotation_row = rowcn,
show_rownames = T)
可以看到， SCENIC包自己提供了一个 calcRSS函数跟我的方法去找到的各个单细胞亚群的特异性的转录因子其实大同小异哦
确实每个单细胞亚群都是有 自己的特异性的激活的转录因子，仍然是推荐看看 各个单细胞亚群特异性的转录因子热图。其实很久很久以前的教程 使用pyscenic做转录因子分析 ，也是 同样的可视化思想。
我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你
联系客服
