单细胞测序单细胞测序单细胞转录组

单细胞转录组高级分析二:转录调控网络分析

2021-04-28  本文已影响0人  Seurat_Satija

SCENIC简介

组织内细胞异质性的基础是细胞转录状态的差异,转录状态的特异性又是由转录因子主导的基因调控网络(GRNs)决定并维持稳定的。因此分析单细胞的GRNs有助于深入挖掘细胞异质性背后的生物学意义,并为疾病的诊断、治疗以及发育分化的研究提供有价值的线索。然而单细胞转录组数据具有背景噪音高、基因检出率低和表达矩阵稀疏性的特点,给传统统计学和生物信息学方法推断高质量的GRNs带来了挑战。
Single-cell regulatory network inference and clustering (SCENIC)是一种专为单细胞数据开发的GRNs算法,它的创新之处在于引入了转录因子motif序列验证统计学方法推断的基因共表达网络,从而识别高可靠性的由转录因子主导的GRNs。SCENIC相关的文章2017年首先发表于nature methods,2020年又将流程整理后发表于nature protocls。需要深入了解分析原理和流程的朋友可以参考这两篇文章:

SCENIC : single-cell regulatory network inference and clustering

A scalable SCENIC workflow for single-cell gene regulatory network analysis

SCENIC分析流程
官方介绍的主要分析有四步:

  1. GENIE3/GRNBoost:基于共表达情况鉴定每个TF的潜在靶点;
  2. RcisTarget:基于DNA-motif 分析选择潜在的直接结合靶点;
  3. AUCell:分析每个细胞的regulons活性;
  4. 细胞聚类:基于regulons的活性鉴定稳定的细胞状态并对结果进行探索

分析流程可以细化为以下步骤:

图片

SCENIC安装

详细的安装教程参阅官网:

https://rawcdn.githack.com/aertslab/SCENIC/701cc7cc4ac762b91479b3bd2eaf5ad5661dd8c2/inst/doc/SCENIC_Setup.html

安装代码如下

在安装SCENIC之前,请按照以下代码安装一些依赖包:


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)

检查一下核心依赖包的版本,并确保版本符合以下要求:

AUCell >=1.4.1 (minimum 1.2.4);

RcisTarget>=1.2.0 (minimum 1.0.2);

GENIE3>=1.4.0 (minimum 1.2.1).

packageVersion("AUCell")
packageVersion("RcisTarget")
packageVersion("GENIE3")

安装SCENIC包

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

还需要下载一些数据库

可以在此网址用浏览器下载:https://resources.aertslab.org/cistarget/

也可以使用以下代码在R环境中直接下载:

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

SCENIC分析准备

SCENIC分析需要输入单细胞表达矩阵和细胞meta信息,正式分析之前还要设置运行线程数,参考数据库目录、版本,细胞meta信息存放目录等配置信息。


library(Seurat)
library(tidyverse)
library(patchwork)
library(SCENIC)
rm(list=ls())

##==分析准备==##
dir.create("SCENIC")
dir.create("SCENIC/int")
scRNA <- readRDS("scRNA.rds")
setwd("~/project/10xDemo2/SCENIC") 
##准备细胞meta信息
cellInfo <- data.frame(scRNA@meta.data)
colnames(cellInfo)[which(colnames(cellInfo)=="orig.ident")] <- "sample"
colnames(cellInfo)[which(colnames(cellInfo)=="seurat_clusters")] <- "cluster"
colnames(cellInfo)[which(colnames(cellInfo)=="celltype_Monaco")] <- "celltype"
cellInfo <- cellInfo[,c("sample","cluster","celltype")]
saveRDS(cellInfo, file="int/cellInfo.Rds")
##准备表达矩阵
#为了节省计算资源,随机抽取1000个细胞的数据子集
subcell <- sample(colnames(scRNA),1000)
scRNAsub <- scRNA[,subcell]
saveRDS(scRNAsub, "scRNAsub.rds")
exprMat <- as.matrix(scRNAsub@assays$RNA@counts)
##设置分析环境
mydbDIR <- "./cisTarget"
mydbs <- c("hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather",
           "hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather")
names(mydbs) <- c("500bp", "10kb")
scenicOptions <- initializeScenic(org="hgnc", 
                                  nCores=8,
                                  dbDir=mydbDIR, 
                                  dbs = mydbs,
                                  datasetTitle = "HNSCC")
saveRDS(scenicOptions, "int/scenicOptions.rds")

scenicOptions <- initializeScenic()这一步要注意自己的情况设置参数:

共表达网络计算

相关性回归分析

SCENIC正式分析的第一步是计算转录因子与每个基因的相关性,此步骤消耗的计算资源非常大,作者建议两个策略处理大数据集:

为了减少初学者的学习负担,本文采用第一种策略。

##==转录调控网络推断==##
##基因过滤
#过滤标准是基因表达量之和>细胞数*3%,且在1%的细胞中表达
genesKept <- geneFiltering(exprMat, scenicOptions, 
              minCountsPerGene = 3 * 0.01 * ncol(exprMat), 
              minSamples = ncol(exprMat) * 0.01)
exprMat_filtered <- exprMat[genesKept, ]
##计算相关性矩阵
runCorrelation(exprMat_filtered, scenicOptions)
##TF-Targets相关性回归分析
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)
#这一步消耗的计算资源非常大,个人电脑需要几个小时的运行时间

runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)需要注意nParts参数,它的作用是把表达矩阵分成n份分开计算,目的是防止数据量大时内存不够。以上代码运行后,int目录下有不少中间结果产生,简要解释一下:

查看1.4_GENIE3_linkList.Rds文件,head(1.4_GENIE3_linkList)

> TF Target    weight
2003093  XBP1   MZB1 0.1616139
2003094  XBP1  DERL3 0.1489923
2003095  XBP1  PRDX4 0.1393732
1003616 RPS4X  RPL34 0.1285798
1254129 RPS4X   RPS6 0.1265229
1       RPS4X  RPL13 0.1256756</pre>

TF是转录因子名称,Target是潜在靶基因的名字,weight是TF与Target之间的相关性权重。

推断共表达模块

上一步计算了转录因子与每一个基因的相关性,接下来需要过滤低相关性的组合形成共表达基因集(模块)。作者尝试了多种策略(标准)过滤低相关性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构建共表达模块;

##推断共表达模块
runSCENIC_1_coexNetwork2modules(scenicOptions)

主要运行结果是int目录下的1.6_tfModules_asDF.Rds

> Target    TF method corr
1   BCDIN3D ABCF2   w001    1
2    STRADB ABCF2   w001    1
3 HIST1H2BG ABCF2   w001    1
4      GID4 ABCF2   w001    1
5    LPCAT1 ABCF2   w001   -1
6      ZHX2 ABCF2   w001    1</pre>

method是上面提到的6种方法,corr是runCorrelation(exprMat_filtered, scenicOptions)命令得到的,1代表激活,-1代表抑制,0代表中性,SCENIC只会采用corr值为1的数据用于后续分析。

Motif验证共表达模块

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

  1. 对每个共表达模块进行motif富集分析,保留显著富集的motif;此项分析依赖gene-motif评分(排行)数据库,其行为基因、列为motif、值为排名,就是我们下载的cisTarget数据库。

  2. 使用数据库对motif进行TF注释,注释结果分高、低可信度 。数据库直接注释和同源基因推断的TF是高可信结果,使用motif序列相似性注释的TF是低可信结果。

  3. 用保留的motif对共表达模块内的基因进行打分(同样依据cisTarget数据库),识别显著高分的基因(理解为motif离这些基因的TSS很近);

  4. 删除共表达模块内与motif评分不高的基因,剩下的基因集作者称为调控单元(regulon)


##推断转录调控网络(regulon)
runSCENIC_2_createRegulons(scenicOptions)
#以上代码可增加参数coexMethod=c("w001", "w005", "top50", "top5perTarget", "top10perTarget", "top50perTarget"))
#默认6种方法的共表达网络都计算,可以少选几种方法以减少计算量

内存不够的电脑在这一步运行时会报错,重要的结果保存在output文件夹:

Step2_MotifEnrichment.tsv

图片

这是各个共表达模块显著富集motif的注释信息,表头官方解释如下:

Step2_MotifEnrichment_preview.html

图片

显著富集的motif注释信息,表头含义同上个文件。

Step2_regulonTargetsInfo.tsv

图片

最终的regulon文件,将第一个文件的数据整合在了一起,表头含义如下:

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

Regulon活性评分与可视化

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

##==regulon活性评分与可视化==##

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活性tSNE图更有价值,举例说明如下:

图片

左图展示的是某个regulon在所有细胞中AUC值分布情况,红色虚线是自动分配的阈值分界线;左中图蓝色点是regulonAUC值大于阈值的细胞,灰色点是regulonAUC值小于阈值的细胞;右中图是regulonAUC原始值的tSNE图;右图是regulon中转录因子的表达量tSNE图。

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")

二进制转换及衍生分析

runSCENIC_4_aucell_binarize(scenicOptions, exprMat=exprMat_all)

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

SCENIC结果可视化

runSCENIC_3_scoreCells()和runSCENIC_4_aucell_binarize()运行之后会生成一些可视化的heatmap图与tSNE图,但是他们既不容易与seurat分析的结果联系起来,又不容易调整图形参数和分析内容。我们可以调用SCENIC的分析结果,使用seurat和pheatmap进行可视化。

图片

热图图例显示不全,尺寸不可调;中图和右图是runSCENIC_3与runSCENIC_4得到的tSNE图,与seurat的tSNE图很难联系起来。

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
scRNAauc <- AddMetaData(scRNA, AUCmatrix)
scRNAauc@assays$integrated <- NULL
saveRDS(scRNAauc,'scRNAauc.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
scRNAbin <- AddMetaData(scRNA, BINmatrix)
scRNAbin@assays$integrated <- NULL
saveRDS(scRNAbin, 'scRNAbin.rds')

##利用Seurat可视化AUC
dir.create('scenic_seurat')
#FeaturePlot
p1 = FeaturePlot(scRNAauc, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p2 = FeaturePlot(scRNAbin, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p3 = DimPlot(scRNA, 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
图片

用山脊图和小提琴图展示CEBPB调控网络的AUC值

pheatmap可视化SCENIC结果****

#RidgePlot&VlnPlot
p1 = RidgePlot(scRNAauc, features = "CEBPB_extended_2290g", group.by="celltype_Monaco") + 
               theme(legend.position='none')
p2 = VlnPlot(scRNAauc, 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)
图片

有兴趣的朋友还可以做个附加题:把CEBPB主导的基因调控网络做GO和KEGG富集分析。

转自:单细胞转录组高级分析二:转录调控网络分析

上一篇下一篇

猜你喜欢

热点阅读