参考生信技能树：
pyscenic的转录因子分析结果展示之5种可视化、
pyscenic的转录因子分析结果展示之各个单细胞亚群特异性激活转录因子
上一篇见：
pySCENIC的转录因子分析及数据可视化（一）
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)
library(scRNAseq)

load('for_rss_and_visual.Rdata')
head(cellTypes) 
sub_regulonAUC[1:4,1:2] 
dim(sub_regulonAUC)
#[1]  220 2638

值得一提的是这个pbmc3k数据集的两千多个细胞，其实就220个转录因子结果（曾老师教程里面是208个）。

1. TF活性均值

看看不同单细胞亚群的转录因子活性平均值

# Split the cells by cluster:
selectedResolution <- "celltype" # select resolution
cellsPerGroup <- split(rownames(cellTypes), 
                       cellTypes[,selectedResolution]) 

# 去除extened regulons
sub_regulonAUC <- sub_regulonAUC[onlyNonDuplicatedExtended(rownames(sub_regulonAUC)),] 
dim(sub_regulonAUC)
#[1]  220 2638 #似乎没啥区别

# Calculate average expression:
regulonActivity_byGroup <- sapply(cellsPerGroup,
                                  function(cells) 
                                    rowMeans(getAUC(sub_regulonAUC)[,cells]))

# Scale expression. 
# Scale函数是对列进行归一化，所以要把regulonActivity_byGroup转置成细胞为行，基因为列
# 参考：https://www.jianshu.com/p/115d07af3029
regulonActivity_byGroup_Scaled <- t(scale(t(regulonActivity_byGroup),
                                          center = T, scale=T)) 
# 同一个regulon在不同cluster的scale处理
dim(regulonActivity_byGroup_Scaled)
#[1] 220   9
regulonActivity_byGroup_Scaled=regulonActivity_byGroup_Scaled[]
regulonActivity_byGroup_Scaled=na.omit(regulonActivity_byGroup_Scaled)

2. 热图查看TF分布

pheatmap(regulonActivity_byGroup_Scaled)

这里碰到一个怪事，一开始用pheatmap:pheatmap(regulonActivity_byGroup_Scaled) ，这样写报错：Error in pheatmap:pheatmap(regulonActivity_byGroup_Scaled) : NA/NaN argument。

image.png

可以看到，确实每个单细胞亚群都是有 自己的特异性的激活的转录因子。

3. rss查看特异TF

不过，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)

image.png

PAX5(+)和REL(+)的确在B细胞里面高表达。

4. 其他查看TF方式

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)

image.png

也可以按照sd计算df的sd.2

image.png

或者按照mean计算df的sd.2

image.png

至此， pySCENIC的转录因子分析及数据可视化教程复现结束，感谢曾老师和生信技能树。