运行SCENIC参考文献
Cytoscape视频教程

本教包括：

一、pyscenic对单细胞进行转录因子的分析；
二、筛选不同细胞类型差异表达转录因子；
三、通过Cytoscape对不用细胞类型差异转录因子进行可视化。

1、pyscenic对单细胞进行转录因子的分析

需要一些依赖

conda create -n pyscenic python=3.7 
conda activate pyscenic 

mamba install -y numpy scanpy
mamba install -y -c anaconda cytoolz

pip install pyscenic 

数据库配置
转录因子分析还需要自行下载最新更新的数据库（https://resources.aertslab.org/cistarget/databases/）：数据库储存在/data/index_genome/cisTarget_databases/

mkdir  /data/index_genome/cisTarget_databases/
cd /data/index_genome/cisTarget_databases/
cat >download.bash
##1, For human:
wget -c https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather
wget -c https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather

##2, For mouse:
wget -c https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather
wget -c https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather

wget https://github.com/aertslab/pySCENIC/blob/master/resources/hs_hgnc_tfs.txt
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
## 下载
bash download.bash 

Step1.提取单细胞表达量矩阵

library(SCENIC)
packageVersion("SCENIC")
setwd("D:/D/temp/PySCENIC/hexiaoying")
sce <- readRDS("D:/D/metropolypus/FastMNN/Mmerge.rds")

###根据colname名字区分BD单细胞和10X单细胞测序，提取10X单细胞测序数据
colnames(sce)
sce$tech <- sapply(colnames(x = sce), function(x) unlist(strsplit(x, "_"))[1]) 
sce$IF=sapply(sce$tech, function(x) unlist(nchar(x)>10))
table(sce$IF)

Idents(sce)<-sce$IF
ten=subset(sce,idents = c("TRUE"))
table(ten$sample)

#由于pySCENIC运行比较慢，随机提取2万个细胞进行分析
subcell <- sample(colnames(ten),20000)
End.combinedsub <- ten[,subcell]

write.csv(t(as.matrix(End.combinedsub@assays$RNA@counts)),file = "for.scenic.data.csv")
saveRDS(End.combinedsub, "SCENIC.sub.rds")

Step2.pySCENIC常规运行（Python+Linux环境）

cd /mnt/d/D/temp/PySCENIC/hexiaoying

cat >change.py

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("for.scenic.data.csv",first_column_names=True);
row_attrs = {"Gene": np.array(x.var_names),};
col_attrs = {"CellID": np.array(x.obs_names)};
lp.create("sample.loom",x.X.transpose(),row_attrs,col_attrs);

python change.py

##这一步可能会出现的问题
##这一步pyscenic aucell 会报错AttributeError: 'Series' object has no attribute 'iteritems'只要把"count in values.value_counts().iteritems():" to "count in values.value_counts().items():"   解决办法参考https://github.com/dusty-nv/jetson-inference/issues/1640

cat >scenic.sh

#!/bin/sh
##运行pySCENIC
# 不同物种的数据库不一样，这里是人类是human 
dir=/mnt/d/D/linux/cisTarget_databases/ #改成自己的目录
tfs=$dir/hs_hgnc_tfs.txt
feather=$dir/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather
tbl=$dir/motifs-v9-nr.hgnc-m0.001-o0.0.tbl 
# 一定要保证上面的数据库文件完整无误哦 
input_loom=/mnt/d/D/temp/PySCENIC/hexiaoying/sample.loom
ls $tfs  $feather  $tbl  

#2.1 grn
pyscenic grn \
--num_workers 20 \
--output adj.sample.tsv \
--method grnboost2 \
sample.loom \
$tfs #转录因子文件，1839个基因的名字列表

#2.2 cistarget
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

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

bash scenic.sh

输出结果

image.png

Step3.sample_SCENIC.loom导入R里面进行可视化

##可视化
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)
library(scRNAseq)
library(patchwork)
library(ggplot2) 
library(stringr)
library(circlize)

#### 1.提取 out_SCENIC.loom 信息
pyScenicLoomFile <- file.path(getwd(), "out_SCENIC.loom")
loom <- open_loom(pyScenicLoomFile) 

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)


#### 2.加载SeuratData
seurat.data <- readRDS("D:/D/temp/PySCENIC/hexiaoying/SCENIC.sub.rds") 
Idents(seurat.data)<-seurat.data$celltype
DimPlot(seurat.data,label = T)

#可视化
sub_regulonAUC <- regulonAUC[,match(colnames(seurat.data),colnames(regulonAUC))]
dim(sub_regulonAUC)
seurat.data
#确认是否一致
identical(colnames(sub_regulonAUC), colnames(seurat.data))

cellClusters <- data.frame(row.names = colnames(seurat.data), 
                           seurat_clusters = as.character(seurat.data$celltype))
cellTypes <- data.frame(row.names = colnames(seurat.data), 
                        celltype = seurat.data$celltype)
head(cellTypes)
head(cellClusters)
sub_regulonAUC[1:4,1:4] 

#保存一下
save(sub_regulonAUC,cellTypes,cellClusters,seurat.data,
     file = 'for_rss_and_visual.Rdata')


##
#B细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)
regulonsToPlot = c('TAL1(+)','GATA1(+)')
regulonsToPlot %in% row.names(sub_regulonAUC)
seurat.data@meta.data = cbind(seurat.data@meta.data ,t(assay(sub_regulonAUC[regulonsToPlot,])))

# Vis
p1 = DotPlot(seurat.data, features = unique(regulonsToPlot)) + RotatedAxis()
p1
p2 = RidgePlot(seurat.data, features = regulonsToPlot , ncol = 2)
p2
p3 = VlnPlot(seurat.data, features = regulonsToPlot,pt.size = 0)
p3
p4 = FeaturePlot(seurat.data,features = regulonsToPlot)
p4

wrap_plots(p1,p2,p3,p4)
#可选
#scenic_res = assay(sub_regulonAUC) %>% as.matrix()
#seurat.data[["scenic"]] <- SeuratObject::CreateAssayObject(counts = scenic_res)
#seurat.data <- SeuratObject::SetAssayData(seurat.data, slot = "scale.data",
#                                  new.data = scenic_res, assay = "scenic")



###B细胞有两个非常出名的转录因子，TCF4(+) 以及NR2C1(+)
regulonsToPlot = c('TCF4(+)','NR2C1(+)')
regulonsToPlot %in% row.names(sub_regulonAUC)
seurat.data@meta.data = cbind(seurat.data@meta.data ,t(assay(sub_regulonAUC[regulonsToPlot,])))

# Vis
p1 = DotPlot(seurat.data, features = unique(regulonsToPlot)) + RotatedAxis()
p2 = RidgePlot(seurat.data, features = regulonsToPlot , ncol = 2) 
p3 = VlnPlot(seurat.data, features = regulonsToPlot,pt.size = 0)
p4 = FeaturePlot(seurat.data,features = regulonsToPlot)

wrap_plots(p1,p2,p3,p4)
#可选
#scenic_res = assay(sub_regulonAUC) %>% as.matrix()
#seurat.data[["scenic"]] <- SeuratObject::CreateAssayObject(counts = scenic_res)
#seurat.data <- SeuratObject::SetAssayData(seurat.data, slot = "scale.data",
#                                  new.data = scenic_res, assay = "scenic")




##4.1 TF活性均值
### 4.1. TF活性均值
# 看看不同单细胞亚群的转录因子活性平均值
# Split the cells by cluster:
selectedResolution <- "celltype" # select resolution
cellsPerGroup <- split(rownames(cellTypes), 
                       cellTypes[,selectedResolution])

# 去除extened regulons
sub_regulonAUC <- sub_regulonAUC[onlyNonDuplicatedExtended(rownames(sub_regulonAUC)),] 
dim(sub_regulonAUC)

# Calculate average expression:
regulonActivity_byGroup <- sapply(cellsPerGroup,
                                  function(cells) 
                                    rowMeans(getAUC(sub_regulonAUC)[,cells]))

# Scale expression. 
# Scale函数是对列进行归一化，所以要把regulonActivity_byGroup转置成细胞为行，基因为列
# 参考：//www.greatytc.com/p/115d07af3029
regulonActivity_byGroup_Scaled <- t(scale(t(regulonActivity_byGroup),
                                          center = T, scale=T)) 
# 同一个regulon在不同cluster的scale处理
dim(regulonActivity_byGroup_Scaled)
#[1] 209   9
regulonActivity_byGroup_Scaled=na.omit(regulonActivity_byGroup_Scaled)

#热图查看每个单细胞亚群都是有 自己的特异性的激活的转录因子
p6<-Heatmap(
  regulonActivity_byGroup_Scaled,
  name                         = "z-score",
  col                          = colorRamp2(seq(from=-2,to=2,length=11),rev(brewer.pal(11, "Spectral"))),
  show_row_names               = TRUE,
  show_column_names            = TRUE,
  row_names_gp                 = gpar(fontsize = 6),
  clustering_method_rows = "ward.D2",
  clustering_method_columns = "ward.D2",
  row_title_rot                = 0,
  cluster_rows                 = TRUE,
  cluster_row_slices           = FALSE,
  cluster_columns              = FALSE)

pdf("TF.activity1.heatmap.pdf",width = 14,height = 27);p6;dev.off()




#4.2 rss查看特异TF
### 4.2. rss查看特异TF
rss <- calcRSS(AUC=getAUC(sub_regulonAUC), 
               cellAnnotation=cellTypes[colnames(sub_regulonAUC), selectedResolution]) 
rss=na.omit(rss) 
rssPlot <- plotRSS(rss)
plotly::ggplotly(rssPlot$plot)
ggsave("TF_activity.pdf",width = 14,height = 27) 

#4.3其他方式查看TF
rss=regulonActivity_byGroup_Scaled
head(rss)
df = do.call(rbind,
             lapply(1:ncol(rss), function(i){
               dat= data.frame(
                 path  = rownames(rss),
                 cluster = colnames(rss)[i],
                 sd.1 = rss[,i],
                 sd.2 = apply(rss[,-i], 1, median)  
               )
             }))
df$fc = df$sd.1 - df$sd.2
top5 <- df %>% group_by(cluster) %>% top_n(5, fc)
rowcn = data.frame(path = top5$cluster) 
n = rss[top5$path,] 
#rownames(rowcn) = rownames(n)
p7<-pheatmap(n,
         annotation_row = rowcn,
         show_rownames = T)
pdf("TF.activity2.heatmap.pdf",width = 14,height = 14);p7;dev.off()

##建立可以用于转录因子调控基因Cytoscape网络图的表格
#melt的用法https://www.cnblogs.com/liujiaxin2018/p/16673493.html
library(reshape2)
edges <- melt(regulons,id="regulons")
colnames(edges) <- c("to","from")



##提取参与调节的差异转录因子进行分析
##############------------------------------------pySCENIC
edges.DEG<-edges[edges$to %in% DEG.A$gene,]
levels(factor(edges.tf$from))

DEG.A.tf<-DEG.A[DEG.A$gene %in% str_trim(str_extract(edges$from, '[^(]*')),]
##str_trim(str_extract(edges$from, '[^(]*')) 删除ATF3(+)后面的(+)

###写入上调和下调标记
DEG.A.tf$change <- ifelse((DEG.A.tf$avg_log2FC>0),'up',
                     ifelse((DEG.A.tf$avg_log2FC<0),'down','nochange'))

##提取细胞对应的基因
DEG.tf<-dplyr::select(DEG.A.tf,c("cluster","gene"))
write.csv(DEG.tf,file = "DEG.tf.csv",row.names = F)

##提取基因和上下调标记
label.tf<-dplyr::select(DEG.A.tf,c("gene","change"))
write.csv(label.tf,file = "label.tf.csv",row.names = F)

##提取基因FC
fc.tf<-dplyr::select(DEG.A.tf,c("gene","avg_log2FC"))
write.csv(fc.tf,file = "fc.tf.csv",row.names = F)

##统计上下调转录因子在细胞类型中的个数
celltype.tf.chang<-dcast(celltype.tf,cluster~change)
write.csv(celltype.tf.chang,file = "celltype.tf.change.csv",row.names = F)

Step4 导入Cytoscape进行可视化（按住control键可以对结点进行选择）

1、输入文件

image.png

通过file-import-Network from File选择文件

image.png

2、其他参数全部是通过table导入

输入文件

image.png

通过file-import-Table from File选择文件

image.png

image.png

3、给细胞类型设置分组

image.png

分组的样式

image.png

4、根据第一步的分组，对选中的结点设置颜色

image.png

image.png

5、根据avg_log2FC定义node的颜色

这是avg_log2FC的表格

image.png

image.png

6、输入细胞上下调转录因子

image.png

细胞转录差异可视化

image.png

结果

image.png

7、最后导出结果

image.png

操作技巧

过滤的方式选择结点：比如直接选择之前设置的type结点

image.png

选中的结点如果叠加在一起，通过这样设置可以把结点分散开来

image.png