前面的筆記：pyscenic的轉(zhuǎn)錄因子分析結(jié)果展示之5種可視化 帶領(lǐng)大家回顧了一下 單細(xì)胞轉(zhuǎn)錄因子分析之SCENIC流程  ，并且重新認(rèn)識了  使用pyscenic做轉(zhuǎn)錄因子分析 后的結(jié)果。

我們根據(jù)pbmc3k數(shù)據(jù)集里面的b細(xì)胞有兩個(gè)非常出名的轉(zhuǎn)錄因子，TCF4(+) 以及NR2C1(+)，進(jìn)行了可視化。其實(shí)這兩個(gè)轉(zhuǎn)錄因子并不是先驗(yàn)知識，是我們根據(jù)這個(gè)分析結(jié)果進(jìn)行各個(gè)單細(xì)胞亞群特異性激活轉(zhuǎn)錄因子統(tǒng)計(jì)得到的。

讓我們來看看這個(gè)統(tǒng)計(jì)分析如何做，以及如何更好的可視化。如果你確實(shí)計(jì)算資源有限制，其實(shí)看看 各個(gè)單細(xì)胞亞群特異性的轉(zhuǎn)錄因子熱圖，也是很容易理解，并不一定要   使用pyscenic做轉(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)

load(   file = 'for_rss_and_visual.Rdata')
head(cellTypes) 
sub_regulonAUC[1:4,1:2] 

可以看到每個(gè)單細(xì)胞都標(biāo)記好了生物學(xué)名字，而且配套了一個(gè)轉(zhuǎn)錄因子活性矩陣，自己去看 前面的   使用pyscenic做轉(zhuǎn)錄因子分析  教程， 拿到  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

值得一提的是這個(gè)pbmc3k數(shù)據(jù)集的兩千多個(gè)細(xì)胞，其實(shí)就兩百多個(gè)轉(zhuǎn)錄因子結(jié)果哦。

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

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)) 
# 同一個(gè)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)

如下所示的簡單熱圖 ：

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

不過，SCENIC包自己提供了一個(gè) calcRSS函數(shù)，幫助我們來挑選各個(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)行超級簡單。

rss <- calcRSS(AUC=getAUC(sub_regulonAUC), 
   cellAnnotation=cellTypes[colnames(sub_regulonAUC), 
 selectedResolution]) 
rss=na.omit(rss) 
rssPlot <- plotRSS(rss)
plotly::ggplotly(rssPlot$plot)

如下所示，可以看到我們前面的筆記：pyscenic的轉(zhuǎn)錄因子分析結(jié)果展示之5種可視化 ，提到的 pbmc3k數(shù)據(jù)集里面的b細(xì)胞有兩個(gè)非常出名的轉(zhuǎn)錄因子，TCF4(+) 以及NR2C1(+)，確實(shí)是在b細(xì)胞里面名列前茅哦：

大家可以把每個(gè)單細(xì)胞亞群各自的獨(dú)特激活的轉(zhuǎn)錄因子按照我們前面的筆記：pyscenic的轉(zhuǎn)錄因子分析結(jié)果展示之5種可視化 ，拿去安裝進(jìn)行可視化哦。

比如血小板的PRDM1這個(gè)轉(zhuǎn)錄因子，你會發(fā)現(xiàn)，不同的可視化方法只有結(jié)合起來你才能對數(shù)據(jù)有一個(gè)更直觀的感受。

其它可視化（基本都是展現(xiàn)不同細(xì)胞亞群特異性轉(zhuǎn)錄因子）

我們也多次介紹過：

• SCENIC轉(zhuǎn)錄因子分析結(jié)果的解讀
• 單細(xì)胞轉(zhuǎn)錄因子分析之SCENIC流程  （R語言版本，代碼運(yùn)行超級慢）

經(jīng)過了前面的步驟，我們順利挑選到了各個(gè)單細(xì)胞亞群特異性的轉(zhuǎn)錄因子列表，就可以自己去提取它們在每個(gè)細(xì)胞的活性值，自己繪制熱圖或者其它圖表，大同小異的。

下面簡單的分享一下我自己的解決方案，

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包自己提供了一個(gè) calcRSS函數(shù)跟我的方法去找到的各個(gè)單細(xì)胞亞群的特異性的轉(zhuǎn)錄因子其實(shí)大同小異哦

確實(shí)每個(gè)單細(xì)胞亞群都是有 自己的特異性的激活的轉(zhuǎn)錄因子，仍然是推薦看看 各個(gè)單細(xì)胞亞群特異性的轉(zhuǎn)錄因子熱圖。其實(shí)很久很久以前的教程  使用pyscenic做轉(zhuǎn)錄因子分析 ，也是 同樣的可視化思想。