单细胞转录组之转录调控网络分析

课程笔记施工中。。。。。

1.SCENIC简介

目前单细胞转录组领域用的比较多的细胞聚类方法大多是直接从基因表达矩阵推断，但是对于多样本合并分析，很多情况下会出现难以解决的批次效应，例如：有些癌症多样本的聚类结果大多每个样本单独分成一群;对于发育样本，发育前期和后期细胞类型可能存在较大差异，某些样本特异的细胞群，难以判断是批次效应产生的还是真正的生物学效应。

细胞异质性的决定于细胞转录状态的差异，转录状态又是由基因调控网络(GRN, gene regulatory network)决定并维持稳定的。因此分析单细胞的基因调控网络有助于深入挖掘细胞异质性背后的生物学意义。

2017年发表在Nature Methods杂志上的Single-cell regulatory network inference and clustering SCENIC算法，利用scRNA-seq数据，同时进行基因调控网络重建和细胞状态鉴定，应用于肿瘤和小鼠大脑单细胞图谱数据，提出并证明了顺式调控网络分析能够用于指导转录因子和细胞状态的鉴定。SCENIC通过使用生物学驱动的features自动清除肿瘤样本特异性等批次效应。

1.1 SCENIC工作流程

SCENIC工作流程主要有四步：

  1. 用GENIE3（随机森林) 或GRNBoost (Gradient Boosting) 推断转录因子与候选靶基因之间的共表达模块。每个模块包含一个转录因子及其靶基因，纯粹基于共表达。
  2. 使用RcisTarget分析每个共表达模块中的基因，以鉴定enriched motifs；仅保留TF motif富集的模块和targets，每个TF及其潜在的直接targets gene被称作一个调节子（regulon）
  3. 使用AUCell评估每个细胞中每个regulon的活性，AUCell分数用于生成Regulon活性矩阵，通过为每个regulon设置AUC阈值，可以将该矩阵进行二值化（0|1，on|off），这将确定Regulon在哪些细胞中处于“打开”状态。
  4. 细胞聚类：基于regulons的活性鉴定稳定的细胞状态并对结果进行探索。

2.SCENIC分析前的准备

2.1 SCENIC安装

详细教程：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 相关数据库的下载

链接：https://resources.aertslab.org/cistarget/
感觉这个数据库国内连接有些问题，很难下载下来，还可以尝试 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");   #创建一个文件夹保存数据库
setwd("cisTarget_databases")
#如果3个参考数据库都想下载，每次设置变量dbFiles后，都要运行以下代码
for(featherURL in dbFiles)
{
  download.file(featherURL, destfile=basename(featherURL)) # saved in current dir
}

2.3 相关数据的准备

library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
#加载Seurat标准流程跑完的数据
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))

#寻找两个亚群的差异基因
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),]

#由于后续计算量很大，计算资源不足可以抽样
sce=subset(sce, downsample = 10)
sce

3.运行SCENIC

3.1 准备counts矩阵和细胞的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 文件夹下面有下载好 的文件
#创建scenicOptions对象 存储 SCENIC 设置
scenicOptions <- initializeScenic(org="hgnc", 
                                  dbDir=db , nCores=4) 
#nCores 线程数；dbDIR 存放数据库的目录。

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

相关性回归分析

#过滤基因
genesKept <- geneFiltering(exprMat, scenicOptions)
length(genesKept)

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

#计算输入表达式矩阵的 spearman 相关性并以 SCENIC 格式保存结果
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered_log <- log2(exprMat_filtered+1) 

# 最耗费时间的就是这个步骤,推断转录因子与候选靶基因之间的共表达模块
runGenie3(exprMat_filtered_log, scenicOptions)
#runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)需要注意nParts参数，它的作用是把表达矩阵分成n份分开计算，目的是防止数据量大时内存不够。以上代码运行后，int目录下有不少中间结果产生，简要解释一下：1.2_corrMat.Rds：基因之间的相关性矩阵;1.3_GENIE3_weightMatrix_part_1.Rds等：GENIE3的中间结果;1.4_GENIE3_linkList.Rds：GENIE3最终结果，是把“1.3_”开头的文件合并在一起。

推断共表达模块
上一步计算了转录因子与每一个基因的相关性，接下来需要过滤低相关性的组合形成共表达基因集（模块）。作者尝试了多种策略（标准）过滤低相关性TF-Target，研究发现没有一种最佳策略，因此他们的建议是6种过滤标准都用。这6种方法分别是：

  1. w001：以每个TF为核心保留weight>0.001的基因形成共表达模块；
  2. w005：以每个TF为核心保留weight>0.005的基因形成共表达模块；
  3. top50：以每个TF为核心保留weight值top50的基因形成共表达模块；
  4. top5perTarget：每个基因保留weight值top5的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
  5. top10perTarget：每个基因保留weight值top10的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
  6. top50perTarget：每个基因保留weight值top50的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
### 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验证共表达模块

经过上述分析每个转录因子都找到了强相关的靶基因，很多基因调控网络分析到此就结束了。SCENIC的创新之处是对此结果提出了质疑，并通过以下步骤修剪共表达模块形成有生物学意义的调控单元（regulons）：

对每个共表达模块进行motif富集分析，保留显著富集的motif；此项分析依赖gene-motif评分（排行）数据库，其行为基因、列为motif、值为排名，就是我们下载的cisTarget数据库。
使用数据库对motif进行TF注释，注释结果分高、低可信度 。数据库直接注释和同源基因推断的TF是高可信结果，使用motif序列相似性注释的TF是低可信结果。
用保留的motif对共表达模块内的基因进行打分（同样依据cisTarget数据库），识别显著高分的基因（理解为motif离这些基因的TSS很近）；
删除共表达模块内与motif评分不高的基因，剩下的基因集作者称为调控单元（regulon）。

# 断转录调控网络（regulon）
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions,coexMethod=c("top5perTarget")) # Toy run settings

结果保存在output文件夹：Step2_regulonTargetsInfo.tsv,表头含义如下：

  • TF：转录因子名称
  • gene：TF靶基因名称
  • nMotif：靶基因在数据库的motif数量
  • bestMotif：最显著富集的motif名称
  • NES：标准富集分数，分值越高越显著
  • highConfAnnot：是不是高可信注释
  • Genie3Weight：TF与靶基因的相关性权重

请特别注意这个文件，后续分析找到了有价值的regulon，需要回到这个文件找对应的转录因子和靶基因 。SCENIC中regulon的名称有两种，一种是TF名称+extended+靶基因数目，另一种是TF名称+靶基因数目。第一种是转录因子与所有靶基因组成的基因调控网络，第二种是转录因子与高可信靶基因（即highConfAnnot=TRUE的基因）组成的基因调控网络。

Regulon活性评分与可视化
每个Regulon就是一个转录因子及其直接调控靶基因的基因集，SCENIC接下来的工作就是对每个regulon在各个细胞中的活性评分。评分的基础是基因的表达值，分数越高代表基因集的激活程度越高。我们推断regulons虽然只用了1000个随机抽取的细胞，但是regulon评分的时候可以把所有细胞导进来计算。

library(doParallel)
scenicOptions <- initializeScenic(org="hgnc", dbDir=db , nCores=4) 
## 如果报错，尝试多线程重新设置成为 1 个线程
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log )

runSCENIC_3_scoreCells()是一个多种功能打包的函数，它包含的子功能有：

  1. 计算regulon在每个细胞中AUC值，最后得到一个以regulon为行细胞为列的矩阵。结果目录：int/3.4_regulonAUC.Rds
  2. 计算每个regulon的AUC阈值，细胞中regulonAUC值>阈值，代表此regulon在细胞中处于激活状态，否则代表沉默状态。结果目录：int/3.5_AUCellThresholds.Rds
  3. 使用regulonAUC矩阵对细胞进行降维聚类
  4. 用heatmap图展示regulonAUC矩阵，用t-SNE图分别展示每个regulon的活性分布情况。结果目录：output/Step3_开头的系列文件

Regulon活性二进制转换与可视化

对于细胞类型清晰的数据集而言，构建regulons并对其活性打分足够后续分析。但是，在很多情况下将评分转换为二进制的“开/关”，则既方便解释，又能最大化体现细胞类型的差异。将特定的regulon转换为“0/1”有利于探索和解释关键转录因子。将所有的regulons转换为“0/1”后创建二进制的活性矩阵，则可以用于细胞聚类，对消除技术偏倚特别有用。AUCell会自动计算可能的阈值进行“0/1”转换，作者建议在转换之前手工检查并调整这些阈值。

查看调整阈值的代码(可选 )

#使用shiny互动调整阈值
aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_all)
savedSelections <- shiny::runApp(aucellApp)
#保存调整后的阈值
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")

二进制转换及衍生分析

scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)

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

runSCENIC_4_aucell_binarize()将regulonAUC矩阵转换为二进制矩阵后，会重新降维聚类，输出的结果与runSCENIC_3_scoreCells()类似。

3.Seurat可视化SCENIC结果

把SCENIC结果中最重要的regulonAUC矩阵导入Seurat，这样得到的可视化结果更容易与我们之前的分析联系起来。

##导入原始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')

##导入二进制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)

决定单核细胞命运的regulon：转录因子CEBPB主导的调控网络

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

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二进制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已更新#「生信技能树」单细胞公开课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

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 211,948评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,371评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 157,490评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,521评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,627评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,842评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,997评论 3赞 408
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,741评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,203评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,534评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,673评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,339评论 4赞 330
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,955评论 3赞 313
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,770评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,000评论 1赞 266
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,394评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,562评论 2赞 349