前面我们用自带例子和人类的pbmc数据测试了SCENIC的R版本的安装和使用。但是，通常情况下runGenie3这一步，稍微大一点点的数据集，运行时间都很长。比如pbmc这个数据，还不算太大，这一步都运行了近15个小时 。
所以，我们今天测试一下pySCENIC。pySCENIC是通过python实现的一个快速的SCENIC pipeline，pySCENIC支持多线程运行，如果你的服务器条件满足的话，可以大大的节省我们的运算时间。
========利用conda安装========
注：pySCENIC需要python 3.6以上版本解释器。
conda create -y -n pyscenic python=3.7
conda activate pyscenic
conda install -y numpy
conda install -y -c bioconda cytoolz
pip install pyscenic
========下载数据库所需要的文件======
数据库根据调节特征（即转录因子）对用户感兴趣的物种的整个基因组 进行 排名。排名数据库通常以feather格式存储，可以从cisTargetDBs下载。和前面下载的一样。（需要注意的是：需注意版本问题，一个是 pySCENIC < 0.12.0 和ctxcore < 0.2.0软件版本需要进到old文件夹，8月份最近的更新，另一个则是基因组版本，以及数据库来源的版本）我这次下载的是：https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.genes_vs_motifs.rankings.feather
Motif注释文件： 数据库提供了丰富的motif 和结合该motif的转录因子之间的缺失链接。该pipeline需要一个 TSV 文本文件，其中每一行都代表一个特定的注释。（这个需要重新下载，可以从这里下载相应的文件：https://resources.aertslab.org/cistarget/motif2tf/）
TF列表文件：从这里下载：https://github.com/aertslab/pySCENIC/tree/master/resources。我下载的是人的。https://github.com/aertslab/pySCENIC/blob/master/resources/hs_hgnc_tfs.txt
======命令行界面=====
下面主要参考这个帖子：https://www.jianshu.com/p/71ed43163ce1
========第一步：从Seurat对象中提取相关的信息========
此步骤的目的是存储表达矩阵和注释信息，为后面的分析生成输入文件。
get_count_from_seurat.R文件代码如下：
library(optparse)
op_list <- list(
make_option(c("-i", "--inrds"), type = "character", default = NULL, action = "store", help = "The input of Seurat RDS",metavar="rds"),
make_option(c("-d", "--ident"), type = "character", default = NULL, action = "store", help = "The sample Ident of Seurat object",metavar="idents"),
make_option(c("-s", "--size"), type = "integer", default = 1000, action = "store", help = "The sample size of Seurat object",metavar="size"),
make_option(c("-l", "--label"), type = "character", default = "out", action = "store", help = "The label of output file",metavar="label")
)
parser <- OptionParser(option_list = op_list)
opt = parse_args(parser)
library(Seurat)
obj <- readRDS(opt$inrds)
if (is.null(opt$ident)) {
Idents(obj) <- opt$ident
obj <- subset(x = obj, downsample = opt$size)
saveRDS(obj,"subset.rds")
}
if (is.null(opt$label)) {
label1 <- 'out'
}else{
label1 <- opt$label
}
write.csv(t(as.matrix(obj@assays$RNA@counts)),file = paste0(label1,'.csv'),quote=F)
write.table(obj@meta.data,'metadata_subset.xls',sep='\t',quote=F)
其中：
-i，输入Seurat对象的RDS文件
-d，随机取样分组的列名，例如groups，如果不赋值则表示不随机取样，使用全部细胞
-s，随机取样的大小，例如20，因为这里用的是pbmc_small，所以设置小一点，实际情况可能需要设置大一点
-l，输出文件的标签，默认为out。
运行代码如下：Rscript get_count_from_seurat.R -i test.rds -d groups -s 20 -l out
运行后会生成3个文件：矩阵out.csv，metadata文件metadata_subset.xls和取子集后的RDS文件subset.rds（如果不取子集，这个文件不会生成）。
=========第二步：使用python导入csv文件后生成loom文件========
此步骤是将表达矩阵转换为loom格式，因为pyscenic的输入为loom格式。
create_loom_file_from_scanpy.py 文件代码如下：
import argparse
import os, sys
import loompy as lp;
import numpy as np;
import scanpy as sc;
def main():
parser= argparse.ArgumentParser(description='make input for pySCENIC')
parser.add_argument('-i', '--input', type=str, required=True, metavar='input_csv')
args = parser.parse_args()
x=sc.read_csv(args.input)
row_attrs = {"Gene": np.array(x.var_names),}
col_attrs = {"CellID": np.array(x.obs_names)}
name = args.input.split('.')[0]
lp.create(name+'.loom',x.X.transpose(),row_attrs,col_attrs)
if __name__ == '__main__':
main()
使用方法：
-i，传入的csv矩阵文件，例如第一步得到的out.csv
python create_loom_file_from_scanpy.py -i out.csv
运行后生成out.loom文件。
==========第三步：运行SCENIC的python版本==========
SCENIC的标准流程分为3步，前面的帖子已经讲过了：第一步利用GENIE3构建转录因子与基因的调控网络，第二步利用RcisTarget验证转录因子与基因的调控网络的真实性，第三步计算AUC曲线筛选调控网络。
生成一个shell脚本，pyscenic_from_loom.sh文件代码如下：
#default value
input_loom=out.loom
n_workers=20
#help function
function usage() {
echo -e "OPTIONS:\n-i|--input_loom:\t input loom file"
echo -e "-n|--n_workers:\t working core number"
echo -e "-h|--help:\t Usage information"
exit 1
}
#get value
while getopts :i:n:h opt
do
case "$opt" in
i) input_loom="$OPTARG" ;;
n) n_workers="$OPTARG" ;;
h) usage ;;
:) echo "This option -$OPTARG requires an argument."
exit 1 ;;
?) echo "-$OPTARG is not an option"
exit 2 ;;
esac
done
#需要更改为自己的路径
tfs=hs_hgnc_tfs.txt
feather=hg19-tss-centered-10kb-7species.mc9nr.genes_vs_motifs.rankings.feather
tbl=motifs-v9-nr.hgnc-m0.001-o0.0.tbl
pyscenic=pyscenic
# grn
$pyscenic grn \
--num_workers $n_workers \
--output grn.tsv \
--method grnboost2 \
$input_loom $tfs
# cistarget
$pyscenic ctx \
grn.tsv $feather \
--annotations_fname $tbl \
--expression_mtx_fname $input_loom \
--mode "dask_multiprocessing" \
--output ctx.csv \
--num_workers $n_workers \
--mask_dropouts
# AUCell
$pyscenic aucell \
$input_loom \
ctx.csv \
--output aucell.loom \
--num_workers $n_workers
使用方法：
-i，输入的loom文件，例如第二步生成的out.loom
-n，运行线程数，设置越大越快，根据实际情况设置
./pyscenic_from_loom.sh -i out.loom -n 20
注：其中最耗时的也就是这一步，但是运行pbmc的数据也只花了5分钟左右，比如R版本已经是快太多了。
========第四步：计算RSS=====
此步骤的作用是通过计算RSS值找到细胞类型的特异TF。
calcRSS_by_scenic.R文件代码如下：
library(optparse)
op_list <- list(
make_option(c("-l", "--input_loom"), type = "character", default = NULL, action = "store", help = "The input of aucell loom file",metavar="rds"),
make_option(c("-m", "--input_meta"), type = "character", default = NULL, action = "store", help = "The metadata of Seurat object",metavar="idents"),
make_option(c("-c", "--celltype"), type = "character", default = NULL, action = "store", help = "The colname of metadata to calculate RSS",metavar="lab
el")
)
parser <- OptionParser(option_list = op_list)
opt = parse_args(parser)
library(Seurat)
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
library(dplyr)
library(KernSmooth)
library(RColorBrewer)
library(plotly)
library(BiocParallel)
loom <- open_loom(opt$input_loom)
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')
regulonAucThresholds <- get_regulon_thresholds(loom)
close_loom(loom)
meta <- read.table(opt$input_meta,sep='\t',header=T,stringsAsFactor=F)
#因为我也不想生成单个cell的文件，所以改成所有的cell_type了
cellinfo <- meta[,c(opt$celltype,"nFeature_RNA","nCount_RNA")]
colnames(cellinfo)=c('celltype', 'nGene' ,'nUMI')
cellTypes <- as.data.frame(subset(cellinfo,select = 'celltype'))
selectedResolution <- "celltype"
sub_regulonAUC <- regulonAUC
rss <- calcRSS(AUC=getAUC(sub_regulonAUC),
cellAnnotation=cellTypes[colnames(sub_regulonAUC),
selectedResolution])
rss=na.omit(rss)
try({
rssPlot <- plotRSS(rss)
save(regulonAUC,rssPlot,regulons,file='regulon_RSS.Rdata')
})
saveRDS(rss,paste0(opt$celltype,"_rss.rds"))
使用方法：
-i，第三步得到aucell.loom文件
-m，第一步得到的metadata_subset.xls文件
-c，计算RSS的分组列
Rscript calcRSS_by_scenic.R -l aucell.loom -m metadata_subset.xls -c cell_type
运行后生成regulon_RSS.Rdata和cell_type_rss.rds文件
下面的帖子，我们学习一下可视化python版本的结果。
