pyscenic的转录因子分析结果展示之5种可视化

pyscenic的转录因子分析结果展示之5种可视化 by 生信技能树

转录因子分析作为单细胞的3大高级分析，大家应该是不再陌生，我们也多次介绍过：

SCENIC转录因子分析结果的解读
单细胞转录因子分析之SCENIC流程（R语言版本，代码运行超级慢）

但是在R里面跑这个，超级耗时，所以有使用pyscenic做转录因子分析和没想到自己会放弃conda（docker镜像的pyscenic做单细胞转录因子分析），大家可以按需取用。

运行 pyscenic的转录因子分析

我们仍然是以大家熟知的pbmc3k数据集为例。大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释，存储为Rdata文件。

# 0.安装R包 ----
# devtools::install_github('caleblareau/BuenColors')
# utils::install.packages(pkgs = "ggstatsplot")
# InstallData("pbmc3k") 

library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")  
data("pbmc3k")  
sce <- pbmc3k.final   
library(Seurat)
table(Idents(sce))
p1=DimPlot(sce,label = T)

因为这个数据集已经进行了不同单细胞亚群的标记，所以我们很容易拿到单细胞表达量矩阵和细胞对应的表型信息。首先需要在R里面的把seurat对象输出成为csv文件：

write.csv(t(as.matrix(sce@assays$RNA@counts)),
          file = "pbmc_3k.all.csv")

然后在Linux环境里面写一个 Python脚本（ csv2loom.py ）把 csv格式的表达量矩阵转为 .loom 文件

import os, sys
os.getcwd()
os.listdir(os.getcwd()) 

import loompy as lp;
import numpy as np;
import scanpy as sc;
x=sc.read_csv("pbmc_3k.csv");
row_attrs = {"Gene": np.array(x.var_names),};
col_attrs = {"CellID": np.array(x.obs_names)};
lp.create("pbmc_3k.loom",x.X.transpose(),row_attrs,col_attrs);

上面的脚本写了后，就可以运行 Python脚本（ csv2loom.py ）把 csv格式的表达量矩阵转为 .loom 文件：

conda activate pyscenic
python csv2loom.py

最后运行 run_pyscenic.sh 的脚本，命令是：

 nohup bash run_pyscenic.sh &

而 run_pyscenic.sh 的脚本，内容如下所示：

# 不同物种的数据库不一样，这里是人类是human 
dir=/home/bakdata/x10/jmzeng/pyscenic
tfs=$dir/TF/TFs_list/hs_hgnc_tfs.txt
feather=$dir/hg19-tss-centered-10kb-7species.mc9nr.feather
tbl=$dir/TF/TFs_annotation_motif/human_TFs/motifs-v9-nr.hgnc-m0.001-o0.0.tbl 
# 一定要保证上面的数据库文件完整无误哦 
input_loom=pbmc_3k.loom
ls $tfs  $feather  $tbl  

pyscenic grn \
--num_workers 20 \
--output adj.sample.tsv \
--method grnboost2 \
$input_loom  $tfs 


pyscenic ctx \
adj.sample.tsv $feather \
--annotations_fname $tbl \
--expression_mtx_fname $input_loom  \
--mode "dask_multiprocessing" \
--output reg.csv \
--num_workers 20  \
--mask_dropouts


pyscenic aucell \
$input_loom \
reg.csv \
--output out_SCENIC.loom \
--num_workers 20

最重要的文件如下所示：

10M  3 12 09:15 out_SCENIC.loom
6.7M  3 12 09:15 pbmc_3k.loom
13M  3 12 09:15 reg.csv

在R里面读取out_SCENIC.loom进行可视化

首先需要在GitHub安装SCopeLoomR和SCENIC这两个包，然后加载它们：

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)

然后提取 out_SCENIC.loom 信息


inputDir='./outputs/'
scenicLoomPath=file.path(inputDir,'out_SCENIC.loom')
library(SCENIC)
loom <- open_loom(scenicLoomPath) 

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)

然后载入前面的seurat对象，我们这里仅仅是最基础的示例数据，所以直接使用 SeuratData 包即可：

####### step2 : 加载seurat对象  #######
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_check
p <- 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$celltype
DimPlot(sce,reduction = "umap",label=T,group.by = "sub_celltype" )
ggsave('umap-by-sub_celltype.pdf')

Idents(sce) <- sce$sub_celltype


sce <- 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')

这里的代码仍然是简单的检验一下自己的降维聚类分群是否合理，方便后续合并分析。

单细胞数据分析里面最基础的就是降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释，这个大家基本上问题不大了，使用seurat标准流程即可，不过它默认出图并不好看，详见以前我们做的投票：可视化单细胞亚群的标记基因的5个方法，下面的5个基础函数相信大家都是已经烂熟于心了：

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(pbmc, features = c("MS4A1", "CD79A"))
RidgePlot(pbmc, features = c("MS4A1", "CD79A"), ncol = 1)
DotPlot(pbmc, features = unique(features)) + RotatedAxis()
DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)

现在我们就可以把pyscenic的转录因子分析结果去跟我们的降维聚类分群结合起来进行5种可视化展示。

合并的代码如下所示：


sub_regulonAUC <- regulonAUC[,match(colnames(sce),colnames(regulonAUC))]
dim(sub_regulonAUC)
sce 
#确认是否一致
identical(colnames(sub_regulonAUC), colnames(sce))

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)
head(cellTypes)
head(cellClusters)
sub_regulonAUC[1:4,1:4] 
save(sub_regulonAUC,cellTypes,cellClusters,sce,
     file = 'for_rss_and_visual.Rdata')

这个时候，我们的pbmc3k数据集里面的两千多个细胞都有表达量矩阵，也有转录因子活性打分信息。

我们知道b细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，接下来就可以对这两个进行简单的可视化啦。

首先我们需要把这两个转录因子活性信息添加到降维聚类分群后的的seurat对象里面。

regulonsToPlot = c('TCF4(+)','NR2C1(+)')
regulonsToPlot
sce@meta.data = cbind(sce@meta.data ,t(assay(   sub_regulonAUC[regulonsToPlot,])))
Idents(sce) <- sce$sub_celltype
table(Idents(sce) )

第一个可视化是DotPlot：

DotPlot(sce, features = unique(regulonsToPlot)) + RotatedAxis()

可以看到b细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)，确实是在b细胞比较独特的高，不过它们两个同时也被DC或者血小板共享：

第一个可视化是DotPlot

第二个可视化是山峦图：

RidgePlot(sce, features =  regulonsToPlot , ncol = 1)

可以看到效果其实没有前面的DotPlot好：

第二个可视化是山峦图

另外的图，效果还不如这两个，代码如下所示；

VlnPlot(sce, features =  regulonsToPlot ) 
FeaturePlot(sce, features =  regulonsToPlot )

这里就不浪费版面进行展示了。

不过，到现在为止，仅仅是帮助大家认识了使用pyscenic做转录因子分析后的结果而已，后续进行各个亚群特异性的转录因子判定，以及更多的可视化才是重点，敬请期待哦。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的，有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你

ixxmu / mp_duty