參考生信技能樹：
pyscenic的轉(zhuǎn)錄因子分析結(jié)果展示之5種可視化、
pyscenic的轉(zhuǎn)錄因子分析結(jié)果展示之各個(gè)單細(xì)胞亞群特異性激活轉(zhuǎn)錄因子
上一篇見：
pySCENIC的轉(zhuǎn)錄因子分析及數(shù)據(jù)可視化（一）
pySCENIC的轉(zhuǎn)錄因子分析及數(shù)據(jù)可視化（二）

本次內(nèi)容主要各個(gè)單細(xì)胞亞群特異性激活轉(zhuǎn)錄因子及可視化。首先再次回顧一下pyscenic的轉(zhuǎn)錄因子分析結(jié)果

######  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

值得一提的是這個(gè)pbmc3k數(shù)據(jù)集的兩千多個(gè)細(xì)胞，其實(shí)就220個(gè)轉(zhuǎn)錄因子結(jié)果（曾老師教程里面是208個(gè)）。

1. TF活性均值

看看不同單細(xì)胞亞群的轉(zhuǎn)錄因子活性平均值

# 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 #似乎沒啥區(qū)別

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

# Scale expression. 
# Scale函數(shù)是對(duì)列進(jìn)行歸一化，所以要把regulonActivity_byGroup轉(zhuǎn)置成細(xì)胞為行，基因?yàn)榱?# 參考：http://www.itdecent.cn/p/115d07af3029
regulonActivity_byGroup_Scaled <- t(scale(t(regulonActivity_byGroup),
                                          center = T, scale=T)) 
# 同一個(gè)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)

這里碰到一個(gè)怪事，一開始用pheatmap:pheatmap(regulonActivity_byGroup_Scaled) ，這樣寫報(bào)錯(cuò)：Error in pheatmap:pheatmap(regulonActivity_byGroup_Scaled) : NA/NaN argument。

image.png

可以看到，確實(shí)每個(gè)單細(xì)胞亞群都是有 自己的特異性的激活的轉(zhuǎn)錄因子。

3. rss查看特異TF

不過(guò)，SCENIC包自己提供了一個(gè) calcRSS函數(shù)，幫助我們來(lái)挑選各個(gè)單細(xì)胞亞群特異性的轉(zhuǎn)錄因子，全稱是：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
運(yùn)行超級(jí)簡(jiǎn)單。

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細(xì)胞里面高表達(dá)。

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計(jì)算df的sd.2

image.png

或者按照mean計(jì)算df的sd.2

image.png

至此， pySCENIC的轉(zhuǎn)錄因子分析及數(shù)據(jù)可視化教程復(fù)現(xiàn)結(jié)束，感謝曾老師和生信技能樹。