适用背景
之前写了一篇文章介绍四步完成单细胞数据调控网络流程分析，但生成的结果没有进行可视化，所以本篇文章将介绍可视化的内容。
本篇可视化主要分成三个部分：
- 1、常规热图
- 2、基于Seurat对象作图
- 3、网络图
SCENIC/pySCENIC结果具体展示什么呢？
据我初步了解，主要展示regulons的活性值，所谓regulons其实应该就是TF，另外可以展示这些regulons调控的基因网络图。
可视化举例
做热图主要分为两种，一种是把细胞分组求regulons的活性平均值/中位数，另一种是展示所有细胞regulons的活性平均值。
-
1、细胞分组求regulons的活性平均值/中位数
图1
-
2、所有细胞的活性值
图2
-
3、基于Seurat的气泡图
图3
-
4、基于Seurat的峰图
图4
5、基于Seurat的小提琴图
图5
6、基于Seurat的点图
图6
7、regulon网络调控图
图7
数据准备
以经典的pbmc3k数据集为例，获取代码如下：
library(Seurat)
library(SeuratData)
data("pbmc3k")
sce <- pbmc3k.final
saveRDS(sce,'pbmc3k.rds')
这个数据集已做好降维聚类和注释，注释列为seurat_annotations：
> head(sce)
orig.ident nCount_RNA nFeature_RNA seurat_annotations percent.mt
AAACATACAACCAC pbmc3k 2419 779 Memory CD4 T 3.0177759
AAACATTGAGCTAC pbmc3k 4903 1352 B 3.7935958
AAACATTGATCAGC pbmc3k 3147 1129 Memory CD4 T 0.8897363
AAACCGTGTATGCG pbmc3k 980 521 NK 1.2244898
AAACGCACTGGTAC pbmc3k 2163 781 Memory CD4 T 1.6643551
AAACGCTGACCAGT pbmc3k 2175 782 CD8 T 3.8160920
AAACGCTGGTTCTT pbmc3k 2260 790 CD8 T 3.0973451
AAACGCTGTAGCCA pbmc3k 1275 532 Naive CD4 T 1.1764706
AAACGCTGTTTCTG pbmc3k 1103 550 FCGR3A+ Mono 2.9011786
AAACTTGAAAAACG pbmc3k 3914 1112 B 2.6315789
RNA_snn_res.0.5 seurat_clusters group1 sub_celltype
AAACATACAACCAC 1 1 Memory CD4 T Memory CD4 T
AAACATTGAGCTAC 3 3 B B
AAACATTGATCAGC 1 1 Memory CD4 T Memory CD4 T
AAACCGTGTATGCG 6 6 NK NK
AAACGCACTGGTAC 1 1 Memory CD4 T Memory CD4 T
AAACGCTGACCAGT 4 4 CD8 T CD8 T
AAACGCTGGTTCTT 4 4 CD8 T CD8 T
AAACGCTGTAGCCA 0 0 Naive CD4 T Naive CD4 T
AAACGCTGTTTCTG 5 5 FCGR3A+ Mono FCGR3A+ Mono
AAACTTGAAAAACG 3 3 B B
然后按照之前的文章跑一下流程：
Rscript get_count_from_seurat.R -i pbmc3k.rds -d seurat_annotations -s 200 -l out
python create_loom_file_from_scanpy.py -i out.csv
sh pyscenic_from_loom.sh -i out.loom -n 10
Rscript calcRSS_by_scenic.R -l aucell.loom -m metadata_subset.xls -c seurat_annotations
生成文件有：
aucell.loom
ctx.csv
grn.tsv
metadata_subset.xls
out.csv
out.loom
regulon_RSS.Rdata
seurat_annotations_rss.rds
subset.rds
加载需要的包：
library(Seurat)
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
library(dplyr)
library(KernSmooth)
library(RColorBrewer)
library(plotly)
library(BiocParallel)
library(pheatmap)
library(cowplot)
library(ggpubr)
library(ggsci)
library(ggplot2)
library(tidygraph)
library(ggraph)
library(stringr)
colpalettes<-unique(c(pal_npg("nrc")(10),pal_aaas("default")(10),pal_nejm("default")(8),pal_lancet("lanonc")(9),
pal_jama("default")(7),pal_jco("default")(10),pal_ucscgb("default")(26),pal_d3("category10")(10),
pal_locuszoom("default")(7),pal_igv("default")(51),
pal_uchicago("default")(9),pal_startrek("uniform")(7),
pal_tron("legacy")(7),pal_futurama("planetexpress")(12),pal_rickandmorty("schwifty")(12),
pal_simpsons("springfield")(16),pal_gsea("default")(12)))
len <- 100
incolor<-c(brewer.pal(8, "Dark2"),brewer.pal(12, "Paired"),brewer.pal(8, "Set2"),brewer.pal(9, "Set1"),colpalettes,rainbow(len))
主函数plot_pyscenic
plot_pyscenic函数参数介绍：
- inloom='aucell.loom'，pyscenic得到的aucell.loom文件，可更改。
- incolor=incolor，调色板，可自定义
- inrss='seurat_annotations_rss.rds'，计算细胞类型特异性regulons得到的_rss.rds后缀文件
- inrds='subset.rds'，用于分析pyscenic的rds，这里分组随机取了200个细胞
- infun='median'，展示的热图用中位数还是平均值，默认用中位数，平均值改成mean。
- ct.col='seurat_annotations'，分组的列，这里用seurat_annotations
- inregulons=NULL，可以指定感兴趣的regulons，如果没有指定，默认使用前5的
- ingrn='grn.tsv'
- ntop1=5，热图展示的活性值top regulon，默认为前5.
- ntop2=50，网络图展示regulons调控的权重top基因，默认为前50。
plot_pyscenic <- function(inloom='aucell.loom',incolor=incolor,inrss='seurat_annotations_rss.rds',inrds='subset.rds',infun='median', ct.col='seurat_annotations',inregulons=NULL,ingrn='grn.tsv',ntop1=5,ntop2=50){
###load data
loom <- open_loom(inloom)
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')
regulonAucThresholds <- get_regulon_thresholds(loom)
embeddings <- get_embeddings(loom)
close_loom(loom)
rss <- readRDS(inrss)
sce <- readRDS(inrds)
##calculate fc
df = do.call(rbind,
lapply(1:ncol(rss), function(i){
dat= data.frame(
regulon = rownames(rss),
cluster = colnames(rss)[i],
sd.1 = rss[,i],
sd.2 = apply(rss[,-i], 1, get(infun))
)
}))
df$fc = df$sd.1 - df$sd.2
#select top regulon
ntopg <- df %>% group_by(cluster) %>% top_n(ntop, fc)
anrow = data.frame( path = ntopg$cluster)
#plot rss by cluster
lcolor <- incolor[1:length(unique(ntopg$cluster))]
names(lcolor) <- unique(anrow$path)
annotation_colors <- list(path=lcolor)
pn1 = rss[ntopg$regulon,]
pn2 = rss[unique(ntopg$regulon),]
rownames(pn1) <- make.unique(rownames(pn1))
rownames(anrow) <- rownames(pn1)
scale='row'
hei <- ceiling(length(ntopg$regulon)*0.16)
pdf(paste0('regulon_activity_in_',ct.col,'.pdf'),10,hei)
pheatmap(pn1,annotation_row = anrow,scale=scale,annotation_colors=annotation_colors,show_rownames = T)
pheatmap(pn2,scale=scale,show_rownames = T)
dev.off()
hei <- ceiling(length(rownames(rss))*0.16)
pdf(paste0('all_regulons_activity_in_',ct.col,'.pdf'),10,hei)
pheatmap(rss,annotation_row = anrow,scale=scale,annotation_colors=annotation_colors,show_rownames = T)
dev.off()
#plot rss by all cells
if (is.null(inregulons)){
inregulons <- unique(ntopg$regulon)
}
pn3=as.matrix(regulonAUC@assays@data$AUC)
regulon <- inregulons
pn3 <- pn3[regulon,]
sce$group1=sce@meta.data[,ct.col]
ancol = data.frame(sce@meta.data[,c('group1')])
colnames(ancol) <- c('group1')
rownames(ancol) <- colnames(sce)
lcolor <- incolor[1:length(unique(ntopg$cluster))]
names(lcolor) <- unique(ntopg$cluster)
annotation_colors <- list(group1 =lcolor)
df1 <- ancol
df1$cell <- rownames(df1)
df1 <- df1[order(df1$group1),]
pn3 <- pn3[,rownames(df1)]
torange=c(-2,2)
pn3 <- scales::rescale(pn3,to=torange)
scale='none'
hei <- ceiling(length(unique(regulon))*0.16)
pdf(paste0('regulon_activity_in_allcells.pdf'),10,hei)
pheatmap(pn3,annotation_col = ancol,scale=scale,annotation_colors=annotation_colors,show_rownames = T,show_colnames = F,cluster_cols=F)
#pheatmap(pn3,scale=scale,show_rownames = T, show_colnames = F,cluster_cols=F)
dev.off()
regulonsToPlot = inregulons
sce$sub_celltype <- sce@meta.data[,ct.col]
sub_regulonAUC <- regulonAUC[,match(colnames(sce),colnames(regulonAUC))]
cellClusters <- data.frame(row.names = colnames(sce),
seurat_clusters = as.character(sce$seurat_clusters))
cellTypes <- data.frame(row.names = colnames(sce),
celltype = sce$sub_celltype)
sce@meta.data = cbind(sce@meta.data ,t(sub_regulonAUC@assays@data@listData$AUC[regulonsToPlot,]))
Idents(sce) <- sce$sub_celltype
nlen <- length(regulonsToPlot)
hei <- ceiling(nlen)*0.2
blu<-colorRampPalette(brewer.pal(9,"Blues"))(100)
pdf('regulons_activity_in_dotplot.pdf',hei,10)
print(DotPlot(sce, features = unique(regulonsToPlot)) + coord_flip()+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))+
scale_color_gradientn(colours = blu)
)
dev.off()
hei=ceiling(nlen/4)*4
pdf('regulons_activity_in_umap.pdf',16,hei)
print(RidgePlot(sce, features = regulonsToPlot , ncol = 4))
print(VlnPlot(sce, features = regulonsToPlot,pt.size = 0 ))
print(FeaturePlot(sce, features = regulonsToPlot))
dev.off()
grn <- read.table(ingrn,sep='\t',header=T,stringsAsFactors=F)
inregulons1=gsub('[(+)]','',inregulons)
c1 <- which(grn$TF %in% inregulons1)
grn <- grn[c1,]
#edge1 <- data.frame()
#node1 <- data.frame()
pdf(paste0(ntop,'_regulon_netplot.pdf'),10,10)
for (tf in unique(grn$TF)) {
tmp <- subset(grn,TF==tf)
if (dim(tmp)[1] > ntop) {
tmp <- tmp[order(tmp$importance,decreasing=T),]
tmp <- tmp[1:ntop,]
}
node2 <- data.frame(tmp$target)
node2$node.size=1.5
node2$node.colour <- 'black'
colnames(node2) <- c('node','node.size','node.colour')
df1 <- data.frame(node=tf,node.size=2,node.colour='#FFDA00')
node2 <- rbind(df1,node2)
edge2 <- tmp
colnames(edge2) <- c('from','to','edge.width')
edge2$edge.colour <- "#1B9E77"
torange=c(0.1,1)
edge2$edge.width <- scales::rescale(edge2$edge.width,to=torange)
graph_data <- tidygraph::tbl_graph(nodes = node2, edges = edge2, directed = T)
p1 <- ggraph(graph = graph_data, layout = "stress", circular = TRUE) + geom_edge_arc(aes(edge_colour = edge.colour, edge_width = edge.width)) +
scale_edge_width_continuous(range = c(1,0.2)) +geom_node_point(aes(colour = node.colour, size = node.size))+ theme_void() +
geom_node_label(aes(label = node,colour = node.colour),size = 3.5, repel = TRUE)
p1 <- p1 + scale_color_manual(values=c('#FFDA00','black'))+scale_edge_color_manual(values=c("#1B9E77"))
print(p1)
}
dev.off()
}
使用方法及结果
使用示例：
plot_pyscenic(inloom='aucell.loom',incolor=incolor,inrss='seurat_annotations_rss.rds',inrds='subset.rds',infun='median', ct.col='seurat_annotations',inregulons=NULL,ingrn='grn.tsv',ntop1=5,ntop2=50)
生成文件如下：
分组热图文件 regulon_activity_in_seurat_annotations.pdf
所有细胞热图文件 regulon_activity_in_allcells.pdf
基于Seurat的气泡图文件 regulons_activity_in_dotplot.pdf
基于Seurat的峰图、小提琴图和散点图文件 regulons_activity_in_umap.pdf
网络图文件 50_regulon_netplot.pdf
小结与讨论
这是我目前想到的可视化类型，欢迎大家评论区补充。
*以上作图参考了pySCENIC的转录因子分析及数据可视化的系列文章和使用pyscenic做转录因子分析。
网友评论