单细胞 | pySCENIC·转录因子分析（二）

接上一篇，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)

4. 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)