实用配景
单细胞转录组调控网络分析是单细胞转录组分析内容的高级分析之一，本文将先容SCENIC/pySCENIC的流程，详细原理和内容不睁开，主要展示代码复现流程。R的SCENIC基于AUCell，RcisTarget和GENIE3三个包举行分析，以是要先安装这些依靠包，而pySCENIC则已经封装好，直接用pip安装即可。只用SCENIC或pySCENIC也可以单独完身分析，但R语言运行起来很慢，pySCENIC可以有效提拔分析速率，还用SCENIC是由于可视化用R语言会简朴一些。
可视化部分请看这篇文章SCENIC/pySCENIC结果可视化 2022-11-08
快来看看三步完成单细胞数据调控网络流程分析和可视化!
运行情况预备
Python安装pyscenic
pip install scanpy -i https://mirrors.aliyun.com/pypi/simple/pip install loompy -i https://mirrors.aliyun.com/pypi/simple/pip install pyscenic -i https://mirrors.aliyun.com/pypi/simple/R安装SCENIC
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")BiocManager::install(c("AUCell", "RcisTarget","GENIE3"))if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")devtools::install_github("aertslab/SCENIC")代码实践
第一步，从Seurat对象中取子集并获取矩阵存为csv文件，并提取metadata信息。
此步调的目标是存储表达矩阵和解释信息，为反面的分析天生输入文件。
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$identobj <- 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 argparseimport os, sysimport 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版本，pyscenic
pyscenic的设置文件在官方提供的网站可以找到，须要三个文件：
- 1、TF列表文件：https://github.com/aertslab/pySCENIC/tree/master/resources
或
https://resources.aertslab.org/cistarget/tf_lists/
人的列表文件：https://github.com/aertslab/pySCENIC/blob/master/resources/hs_hgnc_tfs.txt
- 2、数据库feather文件，需留意版本标题，一个是 pySCENIC < 0.12.0 和ctxcore < 0.2.0软件版本须要进到old文件夹，8月份迩来的更新，另一个则是基因组版本，以及数据库泉源的版本：
https://resources.aertslab.org/cistarget/databases/
人的feather文件：
https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-10species.mc9nr.feather
- 3、motif2tf数据库tbl文件：
https://resources.aertslab.org/cistarget/motif2tf/
人的文件：
https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
以上文件选取仅供参考，实际情况大概须要选择差别文件。SCENIC的尺度流程分为3步，第一步利用GENIE3构建转录因子与基因的调控网络，第二步利用RcisTarget验证转录因子与基因的调控网络的真实性，第三步盘算AUC曲线筛选调控网络。
pyscenic_from_loom.sh文件代码如下：
#default valueinput_loom=out.loomn_workers=20#help functionfunction 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 valuewhile getopts :i:n:h optdo 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 ;; esacdone#须要更改路径tfs=hs_hgnc_tfs.txtfeather=hg19-tss-centered-10kb-7species.mc9nr.feathertbl=motifs-v9-nr.hgnc-m0.001-o0.0.tblpyscenic=/app/bin/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，运行线程数，设置越大越快，根据实际情况设置
sh pyscenic_from_loom.sh -i out.loom -n 10第四步，盘算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="label"))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)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 <- regulonAUCrss <- 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(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 groups运行后天生regulon_RSS.Rdata和groups_rss.rds文件
本文只展示了SCENIC/pySCENIC的常规分析流程，没有涉及可视化，下一篇操持先容本文运行结果的可视化。 |
