pySCENIC的转录因子分析及数据可视化（二）

参考生信技能树：

上一篇见：

下一篇见：

1. 依赖包安装

## SCENIC需要一些依赖包，先安装好BiocManager::install(c("AUCell", "RcisTarget"))BiocManager::install(c("GENIE3"))BiocManager::install(c("zoo", "mixtools", "rbokeh"))BiocManager::install(c("DT", "NMF", "pheatmap", "R2HTML", "Rtsne"))BiocManager::install(c("doMC", "doRNG"))devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)devtools::install_github("aertslab/SCENIC")

网上很多说SCENIC很难安装，还好我这抄作业成功了。

library(SCENIC)packageVersion("SCENIC")#[1] '1.2.4’

2. 提取 out_SCENIC.loom 信息

##可视化rm(list=ls())library(Seurat)library(SCopeLoomR)library(AUCell)library(SCENIC)library(dplyr)library(KernSmooth)library(RColorBrewer)library(plotly)library(BiocParallel)library(grid)library(ComplexHeatmap)library(data.table)library(scRNAseq)## 提取 out_SCENIC.loom 信息#inputDir='./outputs/'#scenicLoomPath=file.path(inputDir,'out_SCENIC.loom')library(SCENIC)loom <- open_loom('out_SCENIC.loom') regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")regulons_incidMat[1:4,1:4] regulons <- regulonsToGeneLists(regulons_incidMat)regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')regulonAucThresholds <- get_regulon_thresholds(loom)tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))])embeddings <- get_embeddings(loom) close_loom(loom)rownames(regulonAUC)names(regulons)

3. 加载SeuratData

library(SeuratData) #加载seurat数据集 data("pbmc3k") sce <- pbmc3k.final table(sce$seurat_clusters)table(Idents(sce))sce$celltype = Idents(sce)library(ggplot2) genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A', 'CD19', 'CD79A', 'MS4A1' , 'IGHG1', 'MZB1', 'SDC1', 'CD68', 'CD163', 'CD14', 'TPSAB1' , 'TPSB2', # mast cells, 'RCVRN','FPR1' , 'ITGAM' , 'C1QA', 'C1QB', # mac 'S100A9', 'S100A8', 'MMP19',# monocyte 'FCGR3A', 'LAMP3', 'IDO1','IDO2',## DC3 'CD1E','CD1C', # DC2 'KLRB1','NCR1', # NK 'FGF7','MME', 'ACTA2', ## fibo 'DCN', 'LUM', 'GSN' , ## mouse PDAC fibo 'MKI67' , 'TOP2A', 'PECAM1', 'VWF', ## endo 'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' )library(stringr) genes_to_check=str_to_upper(genes_to_check)genes_to_checkp <- DotPlot(sce , features = unique(genes_to_check), assay='RNA' ) + coord_flip() + theme(axis.text.x=element_text(angle=45,hjust = 1))p ggsave('check_last_markers.pdf',height = 11,width = 11)DimPlot(sce,reduction = "umap",label=T ) sce$sub_celltype = sce$celltypeDimPlot(sce,reduction = "umap",label=T,group.by = "sub_celltype" )ggsave('umap-by-sub_celltype.pdf')Idents(sce) <- sce$sub_celltypesce <- FindNeighbors(sce, dims = 1:15)sce <- FindClusters(sce, resolution = 0.8)table(sce@meta.data$RNA_snn_res.0.8) DimPlot(sce,reduction = "umap",label=T ) ggsave('umap-by-sub_RNA_snn_res.0.8.pdf')

4. 四种可视化

sub_regulonAUC <- regulonAUC[,match(colnames(sce),colnames(regulonAUC))]dim(sub_regulonAUC)sce #确认是否一致identical(colnames(sub_regulonAUC), colnames(sce))#[1] TRUEcellClusters <- data.frame(row.names = colnames(sce), seurat_clusters = as.character(sce$seurat_clusters))cellTypes <- data.frame(row.names = colnames(sce), celltype = sce$sub_celltype)head(cellTypes)head(cellClusters)sub_regulonAUC[1:4,1:4] save(sub_regulonAUC,cellTypes,cellClusters,sce, file = 'for_rss_and_visual.Rdata')

B细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，接下来就可以对这两个进行简单的可视化。首先我们需要把这两个转录因子活性信息添加到降维聚类分群后的的seurat对象里面。

#尴尬的是TCF4(+)我这个数据里面没有，换了个PAX5(+)和REL(+)regulonsToPlot = c('PAX5(+)','REL(+)')regulonsToPlotsce@meta.data = cbind(sce@meta.data ,t(assay( sub_regulonAUC[regulonsToPlot,])))Idents(sce) <- sce$sub_celltypetable(Idents(sce))

DotPlot(sce, features = unique(regulonsToPlot)) + RotatedAxis()RidgePlot(sce, features = regulonsToPlot , ncol = 1)VlnPlot(sce, features = regulonsToPlot,pt.size = 0 ) FeaturePlot(sce, features = regulonsToPlot )

可以看到b细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，确实是在b细胞比较独特的高。

image.png

下一步将进行亚群特异性转录因子分析及可视化。

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。