转录因子分析作为单细胞的3大高级分析，大家应该是不再陌生，我们也多次介绍过：
但是在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))
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)
#确认是否一致
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)
