[Python知识库][scRNA-seq]单细胞转录因子分析——SCENIC实操示例

之前的文章我们写到了用SCENIC对单细胞RNA测序数据进行转录因子分析的原理：[scRNA-seq]单细胞转录因子分析——SCENIC算法简析。在这篇推送中，我们就用数据实例来展示如何得到SCENIC分析结果。

SCENIC有R的版本，通过SCENIC实现。在http://htmlpreview.github.io/?https://github.com/aertslab/SCENIC/blob/master/inst/doc/SCENIC_Setup.html、http://htmlpreview.github.io/?https://github.com/aertslab/SCENIC/blob/master/inst/doc/SCENIC_Running.html网址中有详细的讲解。

但是SCENIC整个算法的计算量都比较大，如果使用R来获得结果，需要非常长的时间。一个10x的样本就需要计算两天左右，当数据量增多时，使用R版本就显得不太现实。因此，SCENIC的作者还开发了对应的python版本，名为pySCENIC。

使用pySCENIC分析，我们可以通过运行python代码，调用函数获得分析结果，也可以直接使用命令行、输入不同的参数得到结果。因为使用命令行更为直接了当，在这，我们就来看看怎么用命令行获得SCENIC的分析结果。

安装及文件准备

pySCENIC的安装较为简单，只需要在系统输入一下命令行，即可安装。详细信息可见https://pyscenic.readthedocs.io/en/latest/installation.html

pip install pyscenic

SCENIC还需要一些参考的数据，包括：

1. 转录起始子的信息。可以在https://resources.aertslab.org/cistarget/找到，在下面的例子中，我用了人类的两个数据hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather和hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather

2. 转录因子的信息。https://github.com/aertslab/pySCENIC/tree/master/resources，我们在下面的示例中使用的是hs_hgnc_tfs.txt

3. 包含motif信息的数据库。人类的在https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl，小鼠的在https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.mgi-m0.001-o0.0.tbl

正如我们之前介绍的，SCENIC的完整分析流程由三部分组成，grn、cistarget和AUCell，想要得到SCENIC的结果，其实也就是分别运行三行命令，得到这三部分的分析结果。下一步分析流程的输入为上一步分析流程的输出。而grn的输入为一个包含细胞基因表达矩阵的loom文件。想要获得这个loom文件，我们可以先将细胞基因表达矩阵存成一个csv文件，再在python中，将其转化成loom文件。

将细胞基因表达矩阵存成csv文件的代码如下，其关键点在于要对基因表达矩阵进行一个转置。因为我已经在Seurat里进行过初步的分析，就直接将需要用SCENIC分析的细胞存成名为相应的Seurat object，再保存它的基因表达矩阵，这一步在R中完成。

library(Seurat)
library(tidyverse)
load(<Seurat Object位置>)
write.csv(t(as.matrix(<Seurat Object>@assays$RNA@counts)),file = <csv文件位置>)

紧接着，就是将csv文件保存为loom文件，使其能符合grn的输入格式。这一步通过python完成。

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(<csv文件位置>);
row_attrs = {"Gene": np.array(x.var_names),};
col_attrs = {"CellID": np.array(x.obs_names)};
lp.create(<loom文件>,x.X.transpose(),row_attrs,col_attrs);

上述代码中<>中的内容（包括<>）需要替换成自己的数据的相应信息，下同。在得到作为输入的loom文件，以及参考数据hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather、hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather、hs_hgnc_tfs.txt、motifs-v9-nr.hgnc-m0.001-o0.0.tbl之后，我们就完成了前期的文件准备工作，可以开始SCENIC的分析。

SCENIC分析

pySCENIC的使用还是比较简单，就是通过三行命令依次跑grn、cistarget和AUCell。

pyscenic grn --num_workers 20 --output <样本名>.adj.tsv --method grnboost2 <loom文件> hs_hgnc_tfs.txt
pyscenic ctx <样本名>.adj.tsv hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather --annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl --expression_mtx_fname <loom文件> --output <样本名>.reg.csv --num_workers 20
pyscenic aucell <loom文件> <样本名>.reg.csv --output <样本名>_SCENIC.loom --num_workers 20

其中：

--num_workers: 多少个线程并行计算

--output: 输出文件名字

--method: grn的算法，可以为genie3或grnboost2a

取决于数据的大小和使用多少计算资源，完成上述三步需要的时间可能在一到几个小时不等。使用上面的参数分析一个10x的数据，除了步骤一花了4个小时，步骤二三都只需要十几分钟。而步骤三的结果文件<样本名>_SCENIC.loom就包含了SCENIC的分析结果。

在实操过程中，小L在跑步骤二时曾遇到下面的报错：
ValueError: "/home/linxy29/data/reference/pySCENIC/hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather" is a cisTarget Feather database in Feather v1 format, which is not supported anymore. Convert them with "convert_cistarget_databases_v1_to_v2.py" (https://github.com/aertslab/create_cisTarget_databases/) to Feather v2 format.

即使是用convert_cistarget_databases_v1_to_v2.py处理之后，也还是会遇到文件格式有问题的报错。这个报错可能和操作系统或者其它东西有关，因为小L安装pySCENIC时，是重新创建了一个新的conda环境，应该不存在版本问题。更换服务器后，使用相同的步骤和数据，是能够得到最终的结果而不会出现这个报错。

在得到结果文件之后，我们就可以利用它画出我们需要的图。以下是一个详细的教程：

http://htmlpreview.github.io/?https://github.com/aertslab/SCENIC/blob/master/Tutorials_JupyterNotebooks/SCENIC_tutorial_2-ExploringOutput.html

最后，祝大家吃好喝好睡好，科研快乐～

Ref:

[1]?https://cloud.tencent.com/developer/article/1854003

[2]?https://scenic.aertslab.org/tutorials/

[3]?Aibar, Sara, et al. "SCENIC: single-cell regulatory network inference and clustering." Nature methods 14.11 (2017): 1083-1086.