单细胞分析工具--pySCENIC转录因子分析

SCENIC全称Single-Cell rEgulatory Network Inference and Clustering，即转录因子分析，用于构建GRN，gene regulatory network基因调控网络。2017年于Nat Methods首次发表R包版本(https://doi:10.1038/nmeth.4463.)，后又于2020年NATURE PROTOCOLS发表python版本(https://doi.org/10.1038/s41596-020-0336-2)，分析速度得到大幅度提升，而分析方法基本相同。

这篇笔记首先简单了解下分析步骤与原理，然后参考文章提供代码进行实操

1、TF转录因子

转录因子(Transcription factor)是指能够识别、结合在某基因上游特异核苷酸序列上的蛋白质;
这类蛋白通过介导RNA聚合酶与DNA模板的结合，从而调控下游靶基因的转录；也可以和其它转录因子形成转录因子复合体来影响基因的转录。
转录因子特异识别、结合的DNA序列称为转录因子结合位点(Transcription factor binding site,TFBS)，长度在5~20bp范围。由于同一个TF可以调控多个，而其具体结合到每个靶基因的TFBS不完全相同，但具有一定的保守性。
基于此，使用motif来概括TFBS的特征序列，具体是使用PWM(Postion Weight Matrix)表示。如下例图，即有的位点是某一种特定碱基，而有的位点是多种可能。

同一个基因上游可能存在多种motif，受到多个TF调控；而同一个TF也可能识别多种motif，从而形成复杂的基因表达转录调控网络。

2、SCENIC分析步骤

2.1 单细胞数据预处理

输入：单细胞测序数据；输出：单细胞表达矩阵

这一步主要就是读入单细胞数据，过滤低质量的细胞、低表达的基因；然后将单细胞count表达矩阵(行名是细胞，列名是基因)存储为csv格式（或者loom）;

文章建议使用Scanpy预处理，因为更适合大的单细胞数据集。但个人还是更习惯于Seurat的流程。

2.2 网络推断Network inference

输入：单细胞表达矩阵，转录因子列表；输出：TF-gene.csv文件。第1列是转录因子，第2列是共表达的基因，第3列是共表达强度(weight)

针对转录因子列表里的每一个转录因子，探索在单细胞表达矩阵中哪些基因与之存在共表达关系。

通过基于树(tree-based)的回归模型，根据转录因子的表达去预测基因的表达情况。具体算法可选择GRNBoost2或者GENIE3，前者的计算效率会更高。

2.3 得到TF-regulon

输入：上一步得到的TF-gene.csv，motif-gene注释数据，motif-TF注释数据，单细胞表达数据
输出：motif富集结果或者regulon数据

分为如下几个子步骤

(1) Module generation

首先根据第二步得到的TF-gene.csv文件，计算中TF -module(one Translate factor and its target genes)。scenic提供下面三种思路计算得到。

Based on a percentile score; the default is to select sets based on the 75th or 90th percentile of a factor’s targets.

Based on the top N targets for each factor, with a default of 50 targets.

Based on the top N regulators for a target; the default settings select sets of 5, 10, and 50 regulators.

按照如上默认值设置的话，有 2 + 1 + 3 = 6种思路的module结果。所以在最终得到的结果里，可能同一个TF有多个module。

根据调控基因与转录因子的相关性方向可分为 activating module(ρ ≥ +0.03)与 repressing module(ρ ≤ −0.03)。由于负相关模块(TF介导对靶基因的转录抑制)的数量一般较少，默认后续分析只保留正相关的activating module进行后续分析。

最后可限制module所包含的最少基因数，默认为20。

(2) Motif enrichment

在提供的motif-gene注释文件中，提供了每个motif的基因排名，即哪些基因的上游区域含有最类似该motif的区域。

基于此计算出每个module中所富集的motif有哪些。然后根据motif-TF注释数据，判断module 所富集到的motif是否为该module TF的motif。再然后，通过module的motif预测模块中真正的靶基因，去除其中的假阳性或者间接靶点。

最后通过合并同一个TF的所有保留基因，作为一个Regulon(one TF and its target gene)

2.4 Regulon活性评分

基于上一步得到的所有Regulon(可认为是Gene Set)，结合单细胞表达矩阵去评价每一个细胞对于每一个Regulon的表达程度评分。

不同细胞对于同一Regulon的表达活性是可以比较的，但同一细胞的不同Regulon的表达与否是不好比较的。

文章提供一种分析思路，将Regulon的活性得分二值化，即是否激活(On/Off)。这主要通过对那些具有 bimodal distribution双峰分布曲线的Regulon进行混合高斯分布拟合，确定合适的分割阈值，从而将细胞划分为两种状态。

原先参考代码是pySCENIC文章里提供的示范代码；如今版本以及数据库均已更新，现主要参考官方手册：

https://pyscenic.readthedocs.io/en/latest/index.html

3、pySCENIC代码实操

原文的预处理步骤是使用Scanpy进行单细胞预处理的，我还是比较熟悉Seurat流程，因此对文章给的代码有所删改。

3.1 准备conda环境

1
2
3
4


conda create -y -n pyscenic python=3.10
conda activate pyscenic

pip install pyscenic

3.2 下载参考数据

（1）motif–gene注释数据(genome ranking database)

https://resources.aertslab.org/cistarget/databases/ 提供了不同物种以及不同基因组版本的注释数据，可根据需要选择；

1
2


wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather

（2）motif–TF注释数据

https://resources.aertslab.org/cistarget/motif2tf/ 提供了不同物种的注释数据，可根据需要选择；

1

wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl

（3）转录因子列表

https://resources.aertslab.org/cistarget/tf_lists/ 提供了不同物种的注释数据，可根据需要选择；

1

wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt

3.3 示例单细胞表达矩阵

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


##下载单细胞示例数据
# wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_v3/pbmc_10k_v3_filtered_feature_bc_matrix.tar.gz
# tar xvf pbmc_10k_v3_filtered_feature_bc_matrix.tar.gz

library(Seurat)
sce = Read10X(data.dir = "filtered_feature_bc_matrix")
sce = CreateSeuratObject(sce)

## 省略质控相关步骤

sce_count = GetAssayData(sce[["RNA"]], slot = "counts")
sce_count[1:4,1:4] #行名基因，列名细胞
write.csv(sce_count, file="sce_count.csv")

Note：

上面输出、保存的表达矩阵行名是基因，列名是细胞；

而SCENIC分析需要的表达矩阵格式为：行名是细胞，列名是基因；

因此可在之后pyscenic分析时添加参数 –transpose，进行转置。

3.4 scenic分析(命令行)

分别对应上述的2.2~2.4简介

1
2
3
4
5
6
7


## expression matrix
f_ex_matrix_csv=/path/to/sce_count.csv

## referrence databases
f_db_names=$(ls /path/to/*feather | tr "\n" " ")
f_motif_path=/path/to/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl
f_tfs=/path/to/hs_hgnc_tfs.txt

(1) pyscenic grn

1
2
3
4
5
6
7
8
9


pyscenic grn \
--output adj.tsv \
--method grnboost2 \
--transpose \
${f_ex_matrix_csv} \
${f_tfs} \
--num_workers 2

head adj.tsv

(2) pyscenic ctx

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


pyscenic ctx adj.tsv \
${f_db_names} \
--annotations_fname ${f_motif_path} \
--expression_mtx_fname ${f_ex_matrix_csv} \
--transpose \
--output sce_regulon.gmt \
--mask_dropouts \
--mode "dask_multiprocessing" \
--num_workers 20


less -SN sce_regulon.gmt
#每一行是一个regulon：第一列是转录因子，第二列是说明性文字，第三列往后是target gene

(3) pyscenic aucell

1
2
3
4
5
6
7
8


pyscenic aucell \
${f_ex_matrix_csv} \
sce_regulon.gmt \
--transpose \
--output sce_regulon_AUC.csv \
--num_workers 20

less -SN regulon_AUC.csv

将AUCell结果二值化：on/off。需要进入python环境

1
2
3
4
5
6
7
8


from pyscenic import binarization
import pandas as pd
auc_mtx = pd.read_csv("sce_regulon_AUC.csv", index_col="Regulon")
# time consuming
auc_binary = binarization.binarize(auc_mtx, num_workers=20)
# 结果返回一个tuple，包含两个元素：第一个是二值化后的表达矩阵；第二个是每个转录因子二值化的阈值
auc_binary[0]
auc_binary[1]

Note :

（1）输入单细胞表达矩阵也支持loom格式；输出结果中也支持不同类型的数据格式，具体可参看命令帮助文档。

（2）官方手册中也提供了多种不同使用方法，例如在python中的分析流程；docker环境中的使用方法。

1、TF转录因子#

2、SCENIC分析步骤#

2.1 单细胞数据预处理#

2.2 网络推断Network inference#

2.3 得到TF-regulon#

(1) Module generation#

(2) Motif enrichment#

2.4 Regulon活性评分#

3、pySCENIC代码实操#

3.1 准备conda环境#

3.2 下载参考数据#

3.3 示例单细胞表达矩阵#

3.4 scenic分析(命令行)#

(1) pyscenic grn#

(2) pyscenic ctx#

(3) pyscenic aucell#