課程筆記施工中。刚操。。再芋。菊霜。

1.SCENIC簡介

目前單細(xì)胞轉(zhuǎn)錄組領(lǐng)域用的比較多的細(xì)胞聚類方法大多是直接從基因表達(dá)矩陣推斷，但是對于多樣本合并分析济赎，很多情況下會(huì)出現(xiàn)難以解決的批次效應(yīng)鉴逞，例如：有些癌癥多樣本的聚類結(jié)果大多每個(gè)樣本單獨(dú)分成一群;對于發(fā)育樣本，發(fā)育前期和后期細(xì)胞類型可能存在較大差異司训，某些樣本特異的細(xì)胞群构捡，難以判斷是批次效應(yīng)產(chǎn)生的還是真正的生物學(xué)效應(yīng)。

細(xì)胞異質(zhì)性的決定于細(xì)胞轉(zhuǎn)錄狀態(tài)的差異豁遭，轉(zhuǎn)錄狀態(tài)又是由基因調(diào)控網(wǎng)絡(luò)(GRN, gene regulatory network)決定并維持穩(wěn)定的叭喜。因此分析單細(xì)胞的基因調(diào)控網(wǎng)絡(luò)有助于深入挖掘細(xì)胞異質(zhì)性背后的生物學(xué)意義。

2017年發(fā)表在Nature Methods雜志上的Single-cell regulatory network inference and clustering SCENIC算法蓖谢，利用scRNA-seq數(shù)據(jù)捂蕴，同時(shí)進(jìn)行基因調(diào)控網(wǎng)絡(luò)重建和細(xì)胞狀態(tài)鑒定，應(yīng)用于腫瘤和小鼠大腦單細(xì)胞圖譜數(shù)據(jù)闪幽，提出并證明了順式調(diào)控網(wǎng)絡(luò)分析能夠用于指導(dǎo)轉(zhuǎn)錄因子和細(xì)胞狀態(tài)的鑒定啥辨。SCENIC通過使用生物學(xué)驅(qū)動(dòng)的features自動(dòng)清除腫瘤樣本特異性等批次效應(yīng)。

1.1 SCENIC工作流程

SCENIC工作流程主要有四步：

用GENIE3（隨機(jī)森林) 或GRNBoost (Gradient Boosting) 推斷轉(zhuǎn)錄因子與候選靶基因之間的共表達(dá)模塊盯腌。每個(gè)模塊包含一個(gè)轉(zhuǎn)錄因子及其靶基因溉知，純粹基于共表達(dá)。
使用RcisTarget分析每個(gè)共表達(dá)模塊中的基因腕够，以鑒定enriched motifs级乍；僅保留TF motif富集的模塊和targets，每個(gè)TF及其潛在的直接targets gene被稱作一個(gè)調(diào)節(jié)子（regulon）
使用AUCell評估每個(gè)細(xì)胞中每個(gè)regulon的活性帚湘，AUCell分?jǐn)?shù)用于生成Regulon活性矩陣玫荣，通過為每個(gè)regulon設(shè)置AUC閾值，可以將該矩陣進(jìn)行二值化（0|1大诸，on|off）捅厂，這將確定Regulon在哪些細(xì)胞中處于“打開”狀態(tài)。
細(xì)胞聚類：基于regulons的活性鑒定穩(wěn)定的細(xì)胞狀態(tài)并對結(jié)果進(jìn)行探索资柔。

2.SCENIC分析前的準(zhǔn)備

2.1 SCENIC安裝

詳細(xì)教程：https://rawcdn.githack.com/aertslab/SCENIC/701cc7cc4ac762b91479b3bd2eaf5ad5661dd8c2/inst/doc/SCENIC_Setup.html
安裝一些可能用到的依賴包

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::version()
# If your bioconductor version is previous to 3.9, see the section bellow

## Required
BiocManager::install(c("AUCell", "RcisTarget"))
BiocManager::install(c("GENIE3")) # Optional. Can be replaced by GRNBoost

## Optional (but highly recommended):
# To score the network on cells (i.e. run AUCell):
BiocManager::install(c("zoo", "mixtools", "rbokeh"))
# For various visualizations and perform t-SNEs:
BiocManager::install(c("DT", "NMF", "pheatmap", "R2HTML", "Rtsne"))
# To support paralell execution (not available in Windows):
BiocManager::install(c("doMC", "doRNG"))
# To export/visualize in http://scope.aertslab.org
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)

安裝SCENIC包

if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCENIC") 
packageVersion("SCENIC")

2.2 相關(guān)數(shù)據(jù)庫的下載

鏈接：https://resources.aertslab.org/cistarget/
感覺這個(gè)數(shù)據(jù)庫國內(nèi)連接有些問題焙贷，很難下載下來，還可以嘗試 wget等下載方式

##1, For human:
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather",
"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather")
# mc9nr: Motif collection version 9: 24k motifs

##2, For mouse:
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-7species.mc9nr.feather",
"https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-7species.mc9nr.feather")
# mc9nr: Motif collection version 9: 24k motifs

##3, For fly:
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/drosophila_melanogaster/dm6/flybase_r6.02/mc8nr/gene_based/dm6-5kb-upstream-full-tx-11species.mc8nr.feather")
# mc8nr: Motif collection version 8: 20k motifs

##4, download
dir.create("cisTarget_databases");   #創(chuàng)建一個(gè)文件夾保存數(shù)據(jù)庫
setwd("cisTarget_databases")
#如果3個(gè)參考數(shù)據(jù)庫都想下載贿堰，每次設(shè)置變量dbFiles后辙芍，都要運(yùn)行以下代碼
for(featherURL in dbFiles)
{
  download.file(featherURL, destfile=basename(featherURL)) # saved in current dir
}

2.3 相關(guān)數(shù)據(jù)的準(zhǔn)備

library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
#加載Seurat標(biāo)準(zhǔn)流程跑完的數(shù)據(jù)
load(file = '../section-01-cluster/basic.sce.pbmc.Rdata')

sce=pbmc 
table( Idents(sce ))
table(sce@meta.data$seurat_clusters) 
table(sce@meta.data$orig.ident) 

#挑選感興趣的亞群
levels(Idents(sce))
sce = sce[, Idents(sce) %in% c( "FCGR3A+ Mono", "CD14+ Mono"  )] 
levels(Idents(sce))

#尋找兩個(gè)亞群的差異基因
markers_df <- FindMarkers(object = sce, ident.1 = 'FCGR3A+ Mono', ident.2 = 'CD14+ Mono',#logfc.threshold = 0,min.pct = 0.25)
head(markers_df)
#篩選差異基因
cg_markers_df=markers_df[abs(markers_df$avg_log2FC) >1,]
dim(cg_markers_df)
cg_markers_df=cg_markers_df[order(cg_markers_df$avg_log2FC),]

#由于后續(xù)計(jì)算量很大，計(jì)算資源不足可以抽樣
sce=subset(sce, downsample = 10)
sce

3.運(yùn)行SCENIC

3.1 準(zhǔn)備counts矩陣和細(xì)胞的meta信息


table(Idents(sce))
phe=sce@meta.data   
mat=sce@assays$RNA@counts

mat[1:4,1:4]
exprMat =as.matrix(mat) 
dim(exprMat)
exprMat[1:4,1:4] 
head(phe)

cellInfo <-  phe[,c('seurat_clusters','nCount_RNA' ,'nFeature_RNA' )]
colnames(cellInfo)=c('CellType', 'nGene' ,'nUMI')
head(cellInfo)
table(cellInfo$CellType)
cellInfo$CellType=Idents(sce)
table(cellInfo$CellType)

3.2 run-SCENIC

db='~/cisTarget_databases'
list.files(db)
# 保證cisTarget_databases 文件夾下面有下載好 的文件
#創(chuàng)建scenicOptions對象 存儲(chǔ) SCENIC 設(shè)置
scenicOptions <- initializeScenic(org="hgnc", 
                                  dbDir=db , nCores=4) 
#nCores 線程數(shù)羹与；dbDIR 存放數(shù)據(jù)庫的目錄沸手。

saveRDS(scenicOptions, file="int/scenicOptions.Rds") 
saveRDS(cellInfo, file="int/cellInfo.Rds")

相關(guān)性回歸分析

#過濾基因
genesKept <- geneFiltering(exprMat, scenicOptions)
length(genesKept)

exprMat_filtered <- exprMat[genesKept, ]
exprMat_filtered[1:4,1:4]
dim(exprMat_filtered)

#計(jì)算輸入表達(dá)式矩陣的 spearman 相關(guān)性并以 SCENIC 格式保存結(jié)果
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered_log <- log2(exprMat_filtered+1) 

# 最耗費(fèi)時(shí)間的就是這個(gè)步驟,推斷轉(zhuǎn)錄因子與候選靶基因之間的共表達(dá)模塊
runGenie3(exprMat_filtered_log, scenicOptions)
#runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)需要注意nParts參數(shù)外遇，它的作用是把表達(dá)矩陣分成n份分開計(jì)算，目的是防止數(shù)據(jù)量大時(shí)內(nèi)存不夠契吉。以上代碼運(yùn)行后，int目錄下有不少中間結(jié)果產(chǎn)生诡渴，簡要解釋一下：1.2_corrMat.Rds：基因之間的相關(guān)性矩陣;1.3_GENIE3_weightMatrix_part_1.Rds等：GENIE3的中間結(jié)果;1.4_GENIE3_linkList.Rds：GENIE3最終結(jié)果捐晶，是把“1.3_”開頭的文件合并在一起。

推斷共表達(dá)模塊
上一步計(jì)算了轉(zhuǎn)錄因子與每一個(gè)基因的相關(guān)性妄辩，接下來需要過濾低相關(guān)性的組合形成共表達(dá)基因集（模塊）惑灵。作者嘗試了多種策略（標(biāo)準(zhǔn)）過濾低相關(guān)性TF-Target，研究發(fā)現(xiàn)沒有一種最佳策略眼耀，因此他們的建議是6種過濾標(biāo)準(zhǔn)都用英支。這6種方法分別是：

w001：以每個(gè)TF為核心保留weight>0.001的基因形成共表達(dá)模塊；
w005：以每個(gè)TF為核心保留weight>0.005的基因形成共表達(dá)模塊哮伟；
top50：以每個(gè)TF為核心保留weight值top50的基因形成共表達(dá)模塊干花；
top5perTarget：每個(gè)基因保留weight值top5的TF得到精簡的TF-Target關(guān)聯(lián)表，然后把基因分配給TF構(gòu)建共表達(dá)模塊楞黄；
top10perTarget：每個(gè)基因保留weight值top10的TF得到精簡的TF-Target關(guān)聯(lián)表池凄，然后把基因分配給TF構(gòu)建共表達(dá)模塊；
top50perTarget：每個(gè)基因保留weight值top50的TF得到精簡的TF-Target關(guān)聯(lián)表鬼廓，然后把基因分配給TF構(gòu)建共表達(dá)模塊肿仑；

### Build and score the GRN
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)

Motif驗(yàn)證共表達(dá)模塊

經(jīng)過上述分析每個(gè)轉(zhuǎn)錄因子都找到了強(qiáng)相關(guān)的靶基因，很多基因調(diào)控網(wǎng)絡(luò)分析到此就結(jié)束了碎税。SCENIC的創(chuàng)新之處是對此結(jié)果提出了質(zhì)疑尤慰，并通過以下步驟修剪共表達(dá)模塊形成有生物學(xué)意義的調(diào)控單元（regulons）：

對每個(gè)共表達(dá)模塊進(jìn)行motif富集分析，保留顯著富集的motif雷蹂；此項(xiàng)分析依賴gene-motif評分（排行）數(shù)據(jù)庫伟端，其行為基因、列為motif萎河、值為排名荔泳，就是我們下載的cisTarget數(shù)據(jù)庫。
使用數(shù)據(jù)庫對motif進(jìn)行TF注釋虐杯，注釋結(jié)果分高玛歌、低可信度。數(shù)據(jù)庫直接注釋和同源基因推斷的TF是高可信結(jié)果擎椰，使用motif序列相似性注釋的TF是低可信結(jié)果支子。
用保留的motif對共表達(dá)模塊內(nèi)的基因進(jìn)行打分（同樣依據(jù)cisTarget數(shù)據(jù)庫），識(shí)別顯著高分的基因（理解為motif離這些基因的TSS很近）达舒；
刪除共表達(dá)模塊內(nèi)與motif評分不高的基因值朋，剩下的基因集作者稱為調(diào)控單元（regulon）叹侄。

# 斷轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)（regulon）
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions,coexMethod=c("top5perTarget")) # Toy run settings

結(jié)果保存在output文件夾：Step2_regulonTargetsInfo.tsv,表頭含義如下：

TF：轉(zhuǎn)錄因子名稱
gene：TF靶基因名稱
nMotif：靶基因在數(shù)據(jù)庫的motif數(shù)量
bestMotif：最顯著富集的motif名稱
NES：標(biāo)準(zhǔn)富集分?jǐn)?shù)，分值越高越顯著
highConfAnnot：是不是高可信注釋
Genie3Weight：TF與靶基因的相關(guān)性權(quán)重

請?zhí)貏e注意這個(gè)文件昨登，后續(xù)分析找到了有價(jià)值的regulon趾代，需要回到這個(gè)文件找對應(yīng)的轉(zhuǎn)錄因子和靶基因。SCENIC中regulon的名稱有兩種丰辣，一種是TF名稱+extended+靶基因數(shù)目撒强，另一種是TF名稱+靶基因數(shù)目。第一種是轉(zhuǎn)錄因子與所有靶基因組成的基因調(diào)控網(wǎng)絡(luò)笙什，第二種是轉(zhuǎn)錄因子與高可信靶基因（即highConfAnnot=TRUE的基因）組成的基因調(diào)控網(wǎng)絡(luò)飘哨。

Regulon活性評分與可視化
每個(gè)Regulon就是一個(gè)轉(zhuǎn)錄因子及其直接調(diào)控靶基因的基因集，SCENIC接下來的工作就是對每個(gè)regulon在各個(gè)細(xì)胞中的活性評分琐凭。評分的基礎(chǔ)是基因的表達(dá)值芽隆，分?jǐn)?shù)越高代表基因集的激活程度越高。我們推斷regulons雖然只用了1000個(gè)隨機(jī)抽取的細(xì)胞统屈，但是regulon評分的時(shí)候可以把所有細(xì)胞導(dǎo)進(jìn)來計(jì)算胚吁。

library(doParallel)
scenicOptions <- initializeScenic(org="hgnc", dbDir=db , nCores=4) 
## 如果報(bào)錯(cuò)，嘗試多線程重新設(shè)置成為 1 個(gè)線程
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log )

runSCENIC_3_scoreCells()是一個(gè)多種功能打包的函數(shù)鸿吆，它包含的子功能有：

計(jì)算regulon在每個(gè)細(xì)胞中AUC值囤采，最后得到一個(gè)以regulon為行細(xì)胞為列的矩陣。結(jié)果目錄：int/3.4_regulonAUC.Rds
計(jì)算每個(gè)regulon的AUC閾值惩淳，細(xì)胞中regulonAUC值>閾值蕉毯，代表此regulon在細(xì)胞中處于激活狀態(tài)，否則代表沉默狀態(tài)思犁。結(jié)果目錄：int/3.5_AUCellThresholds.Rds
使用regulonAUC矩陣對細(xì)胞進(jìn)行降維聚類
用heatmap圖展示regulonAUC矩陣代虾，用t-SNE圖分別展示每個(gè)regulon的活性分布情況。結(jié)果目錄：output/Step3_開頭的系列文件

Regulon活性二進(jìn)制轉(zhuǎn)換與可視化

對于細(xì)胞類型清晰的數(shù)據(jù)集而言激蹲，構(gòu)建regulons并對其活性打分足夠后續(xù)分析棉磨。但是，在很多情況下將評分轉(zhuǎn)換為二進(jìn)制的“開/關(guān)”学辱，則既方便解釋乘瓤，又能最大化體現(xiàn)細(xì)胞類型的差異。將特定的regulon轉(zhuǎn)換為“0/1”有利于探索和解釋關(guān)鍵轉(zhuǎn)錄因子策泣。將所有的regulons轉(zhuǎn)換為“0/1”后創(chuàng)建二進(jìn)制的活性矩陣衙傀，則可以用于細(xì)胞聚類，對消除技術(shù)偏倚特別有用萨咕。AUCell會(huì)自動(dòng)計(jì)算可能的閾值進(jìn)行“0/1”轉(zhuǎn)換统抬，作者建議在轉(zhuǎn)換之前手工檢查并調(diào)整這些閾值。

查看調(diào)整閾值的代碼(可選 )

#使用shiny互動(dòng)調(diào)整閾值
aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_all)
savedSelections <- shiny::runApp(aucellApp)
#保存調(diào)整后的閾值
newThresholds <- savedSelections$thresholds
scenicOptions@fileNames$int["aucell_thresholds",1] <- "int/newThresholds.Rds"
saveRDS(newThresholds, file=getIntName(scenicOptions, "aucell_thresholds"))
saveRDS(scenicOptions, file="int/scenicOptions.Rds")

二進(jìn)制轉(zhuǎn)換及衍生分析

scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)

export2loom(scenicOptions, exprMat)
saveRDS(scenicOptions, file="int/scenicOptions.Rds")

runSCENIC_4_aucell_binarize()將regulonAUC矩陣轉(zhuǎn)換為二進(jìn)制矩陣后，會(huì)重新降維聚類聪建，輸出的結(jié)果與runSCENIC_3_scoreCells()類似钙畔。

3.Seurat可視化SCENIC結(jié)果

把SCENIC結(jié)果中最重要的regulonAUC矩陣導(dǎo)入Seurat，這樣得到的可視化結(jié)果更容易與我們之前的分析聯(lián)系起來金麸。

##導(dǎo)入原始regulonAUC矩陣
AUCmatrix <- readRDS("int/3.4_regulonAUC.Rds")
AUCmatrix <- AUCmatrix@assays@data@listData$AUC
AUCmatrix <- data.frame(t(AUCmatrix), check.names=F)
RegulonName_AUC <- colnames(AUCmatrix)
RegulonName_AUC <- gsub(' \\(','_',RegulonName_AUC)
RegulonName_AUC <- gsub('\\)','',RegulonName_AUC)
colnames(AUCmatrix) <- RegulonName_AUC
sceauc <- AddMetaData(sce, AUCmatrix)
sceauc@assays$integrated <- NULL
saveRDS(sceauc,'sceauc.rds')

##導(dǎo)入二進(jìn)制regulonAUC矩陣
BINmatrix <- readRDS("int/4.1_binaryRegulonActivity.Rds")
BINmatrix <- data.frame(t(BINmatrix), check.names=F)
RegulonName_BIN <- colnames(BINmatrix)
RegulonName_BIN <- gsub(' \\(','_',RegulonName_BIN)
RegulonName_BIN <- gsub('\\)','',RegulonName_BIN)
colnames(BINmatrix) <- RegulonName_BIN
scebin <- AddMetaData(sce, BINmatrix)
scebin@assays$integrated <- NULL
saveRDS(scebin, 'scebin.rds')

##利用Seurat可視化AUC
dir.create('scenic_seurat')
#FeaturePlot
p1 = FeaturePlot(sceauc, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p2 = FeaturePlot(scebin, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p3 = DimPlot(sce, reduction = 'tsne', group.by = "celltype_Monaco", label=T)
plotc = p1|p2|p3
ggsave('scenic_seurat/CEBPB_extended_2290g.png', plotc, width=14 ,height=4)

決定單核細(xì)胞命運(yùn)的regulon：轉(zhuǎn)錄因子CEBPB主導(dǎo)的調(diào)控網(wǎng)絡(luò)

#RidgePlot&VlnPlot
p1 = RidgePlot(sceauc, features = "CEBPB_extended_2290g", group.by="celltype_Monaco") + 
               theme(legend.position='none')
p2 = VlnPlot(sceauc, features = "CEBPB_extended_2290g", pt.size = 0, group.by="celltype_Monaco") + 
             theme(legend.position='none')
plotc = p1 + p2
ggsave('scenic_seurat/Ridge-Vln_CEBPB_extended_2290g.png', plotc, width=10, height=8)

pheatmap可視化SCENIC結(jié)果

library(pheatmap)
cellInfo <- readRDS("int/cellInfo.Rds")
celltype = subset(cellInfo,select = 'celltype')
AUCmatrix <- t(AUCmatrix)
BINmatrix <- t(BINmatrix)
#挑選部分感興趣的regulons
my.regulons <- c('ETS1_2372g','ETV7_981g','IRF7_239g','XBP1_854g','ATF4_37g',
                 'KLF13_78g','ATF6_129g','CREB3L2_619g','TAGLN2_13g',
                 'STAT1_extended_1808g','CEBPB_extended_2290g','IRF5_extended_422g',
                 'SPI1_1606g','HMGA1_14g','SPIB_1866g','IRF8_348g','BCL11A_136g',
                 'EBF1_40g','MAF_45g','BATF_131g','FOXP3_55g','TBX21_388g',
                 'EOMES_extended_101g','TCF7_extended_31g','LEF1_extended_49g')
myAUCmatrix <- AUCmatrix[rownames(AUCmatrix)%in%my.regulons,]
myBINmatrix <- BINmatrix[rownames(BINmatrix)%in%my.regulons,]
#使用regulon原始AUC值繪制熱圖
pheatmap(myAUCmatrix, show_colnames=F, annotation_col=celltype,
         filename = 'scenic_seurat/myAUCmatrix_heatmap.png',
         width = 6, height = 5)
#使用regulon二進(jìn)制AUC值繪制熱圖
pheatmap(myBINmatrix, show_colnames=F, annotation_col=celltype,
         filename = 'scenic_seurat/myBINmatrix_heatmap.png',
         color = colorRampPalette(colors = c("white","black"))(100),
         width = 6, height = 5)

參考來源

#section 3已更新#「生信技能樹」單細(xì)胞公開課2021_嗶哩嗶哩_bilibili

致謝

I thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

THE END

單細(xì)胞轉(zhuǎn)錄組之轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)分析