查看原文
其他

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

The following article is from 生信会客厅 Author Kinesin

上期专题我们介绍了单细胞转录组数据的基础分析,然而那些分析只是揭开了组织异质性的面纱,还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容:多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正,欢迎添加Kinesin微信交流探讨!


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
## RequiredBiocManager::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.orgif (!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, downloaddir.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()这一步要注意自己的情况设置参数:

  • nCores=8,代表开启8个线程计算,你的参数要根据自己电脑CPU的情况设置;

  • mydbDIR <- "./cisTarget",设置的是我存放数据库的目录,你要填写自己存放数据库的目录;

  • dbs = mydbs,我在前面的变量设置中,把hg38的数据库文件赋值给了mydbs,你用hg19数据库或小鼠的数据库需要相应调整。


共表达网络计算


 相关性回归分析 


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


  • 采用GENIE3推断共表达模块时随机抽取少量细胞计算,计算regulons的活性时,所有细胞都代入运算;

  • 使用python版本的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.2_corrMat.Rds:基因之间的相关性矩阵

  • 1.3_GENIE3_weightMatrix_part_1.Rds等:GENIE3的中间结果

  • 1.4_GENIE3_linkList.Rds:GENIE3最终结果,是把“1.3_”开头的文件合并在一起。


查看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


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


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的注释信息,表头官方解释如下:

  • geneSet: Name of the gene set

  • motif: ID of the motif

  • NES: Normalized enrichment score of the motif in the gene-set

  • AUC: Area Under the Curve (used to calculate the NES)

  • TFinDB: Indicates whether the highlightedTFs are included within the high confidence annotation (two asterisks) or low confidence annotation (one asterisk).

  • TF_highConf: Transcription factors annotated to the motif according to ‘motifAnnot_highConfCat’.

  • TF_lowConf: Transcription factors annotated to the motif according to ‘motifAnnot_lowConfCat’.

  • enrichedGenes: Genes that are highly ranked for the given motif.

  • nErnGenes: Number of genes highly ranked

  • rankAtMax: Ranking at the maximum enrichment, used to determine the number of enriched genes.

Step2_MotifEnrichment_preview.html

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


Step2_regulonTargetsInfo.tsv

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

  • 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评分的时候可以把所有细胞导进来计算。

##==regulon活性评分与可视化==####regulons计算AUC值并进行下游分析exprMat_all <- as.matrix(scRNA@assays$RNA@counts)exprMat_all <- log2(exprMat_all+1)runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)

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$thresholdsscenicOptions@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$AUCAUCmatrix <- data.frame(t(AUCmatrix), check.names=F)RegulonName_AUC <- colnames(AUCmatrix)RegulonName_AUC <- gsub(' \\(','_',RegulonName_AUC)RegulonName_AUC <- gsub('\\)','',RegulonName_AUC)colnames(AUCmatrix) <- RegulonName_AUCscRNAauc <- AddMetaData(scRNA, AUCmatrix)scRNAauc@assays$integrated <- NULLsaveRDS(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_BINscRNAbin <- AddMetaData(scRNA, BINmatrix)scRNAbin@assays$integrated <- NULLsaveRDS(scRNAbin, 'scRNAbin.rds')
##利用Seurat可视化AUCdir.create('scenic_seurat')#FeaturePlotp1 = 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|p3ggsave('scenic_seurat/CEBPB_extended_2290g.png', plotc, width=14 ,height=4)


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


#RidgePlot&VlnPlotp1 = 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 + p2ggsave('scenic_seurat/Ridge-Vln_CEBPB_extended_2290g.png', plotc, width=10, height=8)


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


pheatmap可视化SCENIC结果 

library(pheatmap)cellInfo <- readRDS("int/cellInfo.Rds")celltype = subset(cellInfo,select = 'celltype')AUCmatrix <- t(AUCmatrix)BINmatrix <- t(BINmatrix)#挑选部分感兴趣的regulonsmy.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)


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


后记


写这篇教程我花了一周的业余时间,这几天我都在忙些什么呢?考虑到服务器的配置和使用门槛较高,大家很可能在个人电脑上运行教程的内容,我所有的代码都会在自己的笔记本上验证一遍。此前文章的代码都没有问题,这次运行SCENIC遇到过不去的坎了。SCENIC有几步的计算量非常大,我的电脑运行一次要10几个小时,连续尝试了几个晚上都报错。后来搞清楚根本原因是内存不够,最后用服务器才顺利跑完了流程。如果你也遇到以下警告和错误提示信息,请放弃分析或减少分析的细胞数。

Warning message:In mclapply(argsList, FUN, mc.preschedule = preschedule, mc.set.seed = set.seed, : scheduled cores 3, 6, 8 did not deliver results, all values of the jobs will be affected
Error in signif(trhAssignment, 3)


获取帮助


本教程的目的在于把常用的单细胞分析流程串起来,适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释,官网有详细的教程和权威的解释,翻译好的中文教程也有多个版本,有兴趣的可以搜索一下。如果您学习本教程有一定困难,可以点击 “阅读原文” 找到作者联系方式,获取帮助。


往期回顾

单细胞转录组基础分析一:分析环境搭建

单细胞转录组基础分析二:数据质控与标准化

单细胞转录组基础分析三:降维与聚类

单细胞转录组基础分析四:细胞类型鉴定

单细胞转录组基础分析五:细胞再聚类

单细胞转录组基础分析六:伪时间分析

单细胞转录组基础分析七:差异基因富集分析

单细胞转录组基础分析八:可视化工具总结

欢迎加入生信技能树小圈子

期待单细胞工具的大浪淘沙,洗尽铅华






如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程



看完记得顺手点个“在看”哦!


生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存