使用pyscenic做转录因子分析

前面我们已经是完整的展示了使用R语言的scenic做转录因子分析的流程，但是在公开课演示100个细胞走这个流程，发现居然是耗时20min，实在是不能忍。哪怕是再不喜欢python，也得学一学pyscenic了！

毕竟在文章《A scalable SCENIC workflow for single-cell gene regulatory network analysis》，有这个时间消耗对比，不服不行！

时间消耗对比

一般来说pyscenic是在Linux相关服务器上面操作，所以我也是在自己的ubuntu操作系统演示：

安装自己的conda，每个用户独立操作

安装方法代码如下：

# 首先下载文件，20M/S的话需要几秒钟即可
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件，记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh 
#  安装成功后需要更新系统环境变量文件
source ~/.bashrc


安装好conda后需要设置镜像。

conda config --add channels r 
conda config --add channels conda-forge 
conda config --add channels bioconda
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes 


使用conda安装pyscenic的一些依赖

极端情况下，建议每个软件都可以独立给一个环境，虽然是这样会导致很多冗余的软件，不过，确实很方便。

# 需要一些依赖,尤其是这个python 版本
conda create -n pyscenic python=3.7 
conda activate pyscenic

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

pip install pyscenic 


有时候会遇到网络报错，换个时间点就好了

pip._vendor.urllib3.exceptions.ReadTimeoutError: HTTPSConnectionPool(host='files.pythonhosted.org', port=443): Read timed out.


然后，有意思的地方就来了，我测试了五台服务器，发现并不是所有的服务器都可以成功安装它！后面我们会讲解一个docker的解决方案，假如你确实只有的conda无法安装，也可以等明天的教程哈！

提取单细胞表达量矩阵csv忘记并且导入Linux服务器

首先我们对文章《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》的单细胞转录组数据进行降维聚类分群，然后提取fibo这个子亚群，然后再随机挑取1000个fibo细胞，这样的表达量矩阵进行后续分析。

在seurat里面将矩阵筛选，然后输出成csv，再用python读入，然后打包成 loom

#注意矩阵一定要转置，不然会报错
write.csv(t(as.matrix(fibo@assays$RNA@counts)),file = "fibo_1000.csv")


自己制作一个python 脚本，可以命名为 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("fibo_1000.csv");
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 ，这个命令就会读取当前文件夹的  fibo_1000.csv 文件，进行一些转为，保存为 sample.loom 文件，供后续流程。（这个时候的sample.loom 文件是后续一切分析的开始）

再次强调，使用conda安装的这个pyscenic，它依赖于一系列的python模块，就会在这里报错！比如我遇到的就是pandas的报错：

ImportError: cannot import name 'DtypeArg' from 'pandas._typing' (/home/x10/miniconda3/envs/pyscenic/lib/python3.7/site-packages/pandas/_typing.py)


不过，我使用了docker镜像的pyscenic，所以绕过了这个报错！

需要自己准备好如下所示的文件：

 54M 7月  18 11:18 fibo_1000.csv
1.1G 7月  18 11:30 hg19-tss-centered-10kb-7species.mc9nr.feather
12K 7月  18 11:29 hs_hgnc_tfs.txt
99M 7月  18 11:29 motifs-v9-nr.hgnc-m0.001-o0.0.tbl 


pyscenic 的3个步骤之 grn

代码如下所示，复制粘贴后运行即可：

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


转录因子的基因列表文件下载 ：pySCENIC/resources at master · aertslab/pySCENIC · GitHub

pyscenic 的3个步骤之 cistarget

代码如下所示，复制粘贴后运行即可：

同样的需要下载数据库文件，https://resources.aertslab.org/cistarget/

pyscenic ctx \
adj.sample.tsv \
hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather \
--annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl \
--expression_mtx_fname sample1.loom \
--mode "dask_multiprocessing" \
--output reg.csv \
--num_workers 3 \
--mask_dropouts


pyscenic 的3个步骤之 AUCell

这个步骤超级快，代码如下所示，复制粘贴后运行即可：

pyscenic aucell \
sample.loom \
reg.csv \
--output sample_SCENIC.loom \
--num_workers 3


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

首先再看看文章《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》，其降维聚类分群分群，并且拿出来fibo细胞亚群走这个转录因子分析，得到如下所示：

原文fibo细胞亚群走这个转录因子分析

我们的复现分成2步，首先对进行降维聚类分群，代码如下所示：

rm(list = ls())  
library(tidyverse)
library(Seurat)
library(data.table)
ct=fread('fibo_1000.csv',data.table = F)
ct[1:4,1:4]
rownames(ct)=ct[,1]
ct=ct[,-1]
ct=t(ct)
fibo_sce <- CreateSeuratObject(counts = ct, 
                           project = "fibo_1000" )
fibo_sce
#QC
fibo_sce[["percent.mt"]] <- PercentageFeatureSet(fibo_sce, pattern = "^MT-")
# Visualize QC metrics as a violin plot
VlnPlot(fibo_sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
fibo_sce <- subset(fibo_sce, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 25)
#Normalization
fibo_sce <- NormalizeData(fibo_sce, normalization.method = "LogNormalize", scale.factor = 10000)
#高变基因选择
fibo_sce <- FindVariableFeatures(fibo_sce, selection.method = "vst", nfeatures = 2000)
#标准化
all.genes <- rownames(fibo_sce)
fibo_sce <- ScaleData(fibo_sce, features = all.genes)
#去除MT，重新进行标准化
fibo_sce <- ScaleData(fibo_sce, vars.to.regress = "percent.mt")
#PCA
fibo_sce <- RunPCA(fibo_sce, features = VariableFeatures(object = fibo_sce), verbose = FALSE)
#聚类
fibo_sce <- FindNeighbors(fibo_sce, dims = 1:20)
fibo_sce <- FindClusters(fibo_sce, resolution = 0.5)
#可视化
fibo_sce <- RunUMAP(fibo_sce, dims = 1:20)
fibo_sce<- RunTSNE(fibo_sce, dims = 1:20)
## after we run UMAP and TSNE, there are more entries in the reduction slot
str(fibo_sce@reductions)
DimPlot(fibo_sce, reduction = "umap", label = TRUE)
#查找marker基因
fibo_sce.markers <- FindAllMarkers(fibo_sce, only.pos = TRUE, min.pct = 0.25, 
                               logfc.threshold = 0.25)
saveRDS(fibo_sce, "./fibo_1000.rds")
# 挑选合适的基因，进行可视化


出图如下：

我们的降维聚类分群

文章采取了RGS5和PDGFRA来区分iCAF和mCAF，我们姑且先按照作者的来！我们前面的降维聚类分群操作采取resolution = 0.5的时候确实得到了太多亚群，不过还是尊重作者，进行二分类！

rm(list = ls())  
library(SCENIC)
packageVersion("SCENIC")  
library(SCopeLoomR)
scenicLoomPath='sample_SCENIC.loom'
loom <- open_loom(scenicLoomPath)
# Read information from loom file: 
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonAUC <- get_regulons_AUC(loom, column.attr.name="RegulonsAUC") 
sce=readRDS("./fibo_1000.rds")
sce
library(pheatmap) 
n=t(scale(t( getAUC(regulonAUC[,] )))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]
dim(n) 
ac=data.frame(group= as.character( Idents(sce)))
n[1:4,1:4]
n=n[,colnames(n) %in% colnames(sce)]
rownames(ac)=colnames(n) 
cg=read.table('choose_tf.txt')[,1]
cg
cg_n=n[rownames(n) %in% cg,]
pheatmap(cg_n,show_colnames =F,show_rownames = T,
         annotation_col=ac)
table(ac$group)
# 尊重作者，进行二分类！
ac$group=ifelse(ac$group %in% c(2:5,7,9),'mCAF','iCAF')
pheatmap(cg_n,show_colnames =F,show_rownames = T,
         annotation_col=ac)
pheatmap(cg_n,show_colnames =F,show_rownames = T,
         annotation_col=ac,
         filename = 'heatmap_choose_regulon.png')
dev.off()

出图如下：

复现

可以看到，作者的10个转录因子，在我们这个里面仅仅是剩下了 7 个，不知道是不是因为前面随机抽样1000个细胞导致的。虽然有转录因子的缺失，但是转录组因子的规律并没有变化，在iCAF和mCAF这个亚群特异性激活的转录因子保持原文的样子。

号外: 我们提供单细胞数据分析服务哦

• 明码标价之10X技术单细胞（2.5万每个）（标准100G测序数据）
• 明码标价之10X转录组原始测序数据的cellranger流程
• 明码标价之单细胞转录组的质控降维聚类分群和生物学注释