接上一篇肋乍，pySCENIC分析完之后就可以進(jìn)行可視化了，個(gè)人認(rèn)為最重要的有三個(gè)圖：rss點(diǎn)圖镊逝，rank圖，轉(zhuǎn)錄因子表達(dá)水平熱圖废士，3個(gè)圖結(jié)合著看溪厘。
1.loom文件讀入R骡显，提取數(shù)據(jù)
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.可視化點(diǎn)圖：尋找cluster特異性轉(zhuǎn)錄因子
#提取細(xì)胞metadata信息
cellinfo <-sc@meta.data[,c("cluster_name","group","orig.ident","nFeature_RNA","nCount_RNA")]
colnames(cellinfo)=c('celltype', 'group','orig.ident','nGene' ,'nUMI')
#計(jì)算細(xì)胞特異性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,#可調(diào)整
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參數(shù)讓圖更好看一點(diǎn)
#2.rank圖
cowplot::plot_grid(plotRSS_oneSet(rss2,
setName = table(sc@active.ident)[4]%>%names(),n=3),
NULL,NULL,nrow = 2,byrow = T)
- regulon表達(dá)水平熱圖，轉(zhuǎn)錄因子表達(dá)水平同理聋亡，得到的圖類似只是輸入不一樣
library(ggheatmap)
library(reshape2)
library(RColorBrewer)
#regulon表達(dá)水平熱圖
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)
#轉(zhuǎn)錄因子表達(dá)水平熱圖
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)
