接上一篇，pySCENIC分析完之后就可以进行可视化了，个人认为最重要的有三个图：rss点图，rank图，转录因子表达水平热图，3个图结合着看。
1.loom文件读入R，提取数据
library(SCopeLoomR)
library(SCENIC)
library(AUCell)
loom <- open_loom("/yourpath/SCENIC.loom")
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")
regulonAUC <- get_regulons_AUC(loom,column.attr.name = 'RegulonsAUC')
2.可视化点图：寻找cluster特异性转录因子
#提取细胞metadata信息
cellinfo <-sc@meta.data[,c("cluster_name","group","orig.ident","nFeature_RNA","nCount_RNA")]
colnames(cellinfo)=c('celltype', 'group','orig.ident','nGene' ,'nUMI')
#计算细胞特异性TF
cellTypes <- as.data.frame(subset(cellinfo,select = 'celltype'))
selectedResolution <- "celltype"
cellAnnotation = cellTypes[colnames(regulonAUC),
selectedResolution]
cellAnnotation = na.omit(cellAnnotation)
rss <- calcRSS(AUC = getAUC(regulonAUC),
cellAnnotation = cellAnnotation)
rss = na.omit(rss)
rssPlot <- plotRSS(rss,
zThreshold = 3,#可调整
cluster_columns = FALSE,
order_rows = TRUE,
thr=0.1,
varName = "cellType",
col.low = '#330066',
col.mid = '#66CC66',
col.high = '#FFCC33')
rssPlot$rowOrder
plotly::ggplotly(rssPlot$plot)
3.可视化rank图，可以自己改一下plotRSS_oneSet参数让图更好看一点
#2.rank图
cowplot::plot_grid(plotRSS_oneSet(rss2,
setName = table(sc@active.ident)[4]%>%names(),n=3),
NULL,NULL,nrow = 2,byrow = T)
- regulon表达水平热图，转录因子表达水平同理，得到的图类似只是输入不一样
library(ggheatmap)
library(reshape2)
library(RColorBrewer)
#regulon表达水平热图
tfs <- c("Nr2f1(+)","Cebpd(+)","Hnf4g(+)","Twist1(+)","Twist2(+)","Prrx2(+)",
"Lef1(+)","Foxl2(+)","Foxp1(+)","Hey2(+)","Sox6(+)","Msx1(+)")
rss_data <- rssPlot$plot$data[which(rssPlot$plot$data$Topic %in% tfs),]
rownames(rss_data) <- rss_data[,1]
rss_data <- rss_data[,-1]
colnames(rss_data)
col_ann <- data.frame(group= c(rep("Acta2+ SMC",1),
rep("Notch3+ FB",1),
rep("Col8a1+ FB",1),
rep("Kcnma1+ SMC",1),
rep("Kdr+ EC",1),
rep("Ednrb+ EC",1),
rep("Pecam1+ EC",1),
rep("Pi16+ FB",1),
rep("Eln+ FB",1),
rep("Ly6c1+ EC",1),
rep("Pdgfra+ FB",1),
rep("Klf4+ EC",1),
rep("Ednra+ SMC",1),
rep("Angpt1+ FB",1)))
rownames(col_ann) <- colnames(rss_data)
groupcol <- colorRampPalette(brewer.pal(14,'Set3'))(14)
names(groupcol) <- c("Acta2+ SMC","Notch3+ FB","Col8a1+ FB","Kcnma1+ SMC","Kdr+ EC",
"Ednrb+ EC","Pecam1+ EC","Pi16+ FB","Eln+ FB","Ly6c1+ EC",
"Pdgfra+ FB","Klf4+ EC","Ednra+ SMC","Angpt1+ FB")
col <- list(group=groupcol)
text_columns <- sample(colnames(rss_data),0)
ggheatmap(rss_data,color=colorRampPalette(c('#1A5592','white',"#B83D3D"))(100),
cluster_rows = T,cluster_cols = F,scale = "row",
annotation_cols = col_ann,
annotation_color = col,
legendName="Relative value",
text_show_cols = text_columns)
#转录因子表达水平热图
top3tfgene <- c("Nr2f1","Cebpd","Hnf4g","Twist1","Twist2","Prrx2",
"Lef1","Foxl2","Foxp1","Hey2","Sox6","Msx1")
top3gene_cell_exp <- AverageExpression(sc,
assays = 'RNA',
features = top3tfgene,
group.by = 'celltype',
slot = 'data')
top3gene_cell_exp <- as.data.frame(top3gene_cell_exp$RNA)
top3marker_exp <- t(scale(t(top3gene_cell_exp),scale = T,center = T))
ggheatmap(top3marker_exp,color=colorRampPalette(c('#1A5592','white',"#B83D3D"))(100),
cluster_rows = T,cluster_cols = F,scale = "row",
annotation_cols = col_ann,
annotation_color = col,
legendName="Relative value",
text_show_cols = text_columns)
