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

2022-07-09 本文已影响0人青青青山

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

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工作流程主要有四步：

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

w001：以每个TF为核心保留weight>0.001的基因形成共表达模块；
w005：以每个TF为核心保留weight>0.005的基因形成共表达模块；
top50：以每个TF为核心保留weight值top50的基因形成共表达模块；
top5perTarget：每个基因保留weight值top5的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
top10perTarget：每个基因保留weight值top10的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
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()是一个多种功能打包的函数，它包含的子功能有：

计算regulon在每个细胞中AUC值，最后得到一个以regulon为行细胞为列的矩阵。结果目录：int/3.4_regulonAUC.Rds
计算每个regulon的AUC阈值，细胞中regulonAUC值>阈值，代表此regulon在细胞中处于激活状态，否则代表沉默状态。结果目录：int/3.5_AUCellThresholds.Rds
使用regulonAUC矩阵对细胞进行降维聚类
用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