適用背景
單細(xì)胞轉(zhuǎn)錄組調(diào)控網(wǎng)絡(luò)分析是單細(xì)胞轉(zhuǎn)錄組分析內(nèi)容的高級(jí)分析之一，本文將介紹SCENIC/pySCENIC的流程，具體原理和內(nèi)容不展開憨愉，主要展示代碼復(fù)現(xiàn)流程。R的SCENIC基于AUCell卿捎，RcisTarget和GENIE3三個(gè)包進(jìn)行分析配紫，所以要先安裝這些依賴包，而pySCENIC則已經(jīng)封裝好午阵，直接用pip安裝即可躺孝。只用SCENIC或pySCENIC也可以單獨(dú)完成分析，但R語(yǔ)言運(yùn)行起來(lái)很慢底桂，pySCENIC可以有效提升分析速度植袍，還用SCENIC是因?yàn)榭梢暬肦語(yǔ)言會(huì)簡(jiǎn)單一些。
可視化部分請(qǐng)看這篇文章SCENIC/pySCENIC結(jié)果可視化 2022-11-08
快來(lái)看看三步完成單細(xì)胞數(shù)據(jù)調(diào)控網(wǎng)絡(luò)流程分析和可視化!
運(yùn)行環(huán)境準(zhǔn)備
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")
代碼實(shí)踐
第一步籽懦，從Seurat對(duì)象中取子集并獲取矩陣存為csv文件于个，并提取metadata信息。
此步驟的目的是存儲(chǔ)表達(dá)矩陣和注釋信息暮顺，為后面的分析生成輸入文件厅篓。
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對(duì)象的RDS文件
-d，隨機(jī)取樣分組的列名羽氮，例如groups或链，如果不賦值則表示不隨機(jī)取樣，使用全部細(xì)胞
-s档押，隨機(jī)取樣的大小澳盐，例如20，因?yàn)檫@里用的是pbmc_small汇荐，所以設(shè)置小一點(diǎn)洞就，實(shí)際情況可能需要設(shè)置大一點(diǎn)
-l，輸出文件的標(biāo)簽掀淘，默認(rèn)為out旬蟋。
Rscript get_count_from_seurat.R -i test.rds -d groups -s 20 -l out
運(yùn)行后會(huì)生成3個(gè)文件：矩陣out.csv，metadata文件metadata_subset.xls和取子集后的RDS文件subset.rds（如果不取子集革娄，這個(gè)文件不會(huì)生成）倾贰。
第二步，使用python導(dǎo)入csv文件后生成loom文件
此步驟是將表達(dá)矩陣轉(zhuǎn)換為loom格式拦惋，因?yàn)閜yscenic的輸入為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
運(yùn)行后生成out.loom文件首尼。
第三步，運(yùn)行SCENIC的python版本言秸，pyscenic
pyscenic的配置文件在官方提供的網(wǎng)站可以找到软能，需要三個(gè)文件：
- 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举畸、數(shù)據(jù)庫(kù)feather文件查排，需注意版本問(wèn)題，一個(gè)是 pySCENIC < 0.12.0 和ctxcore < 0.2.0軟件版本需要進(jìn)到old文件夾抄沮，8月份最近的更新跋核，另一個(gè)則是基因組版本，以及數(shù)據(jù)庫(kù)來(lái)源的版本：
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
cis Target databases - 3叛买、motif2tf數(shù)據(jù)庫(kù)tbl文件：
https://resources.aertslab.org/cistarget/motif2tf/
人的文件：
https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
以上文件選取僅供參考砂代，實(shí)際情況可能需要選擇不同文件。SCENIC的標(biāo)準(zhǔn)流程分為3步率挣，第一步利用GENIE3構(gòu)建轉(zhuǎn)錄因子與基因的調(diào)控網(wǎng)絡(luò)泊藕，第二步利用RcisTarget驗(yàn)證轉(zhuǎn)錄因子與基因的調(diào)控網(wǎng)絡(luò)的真實(shí)性，第三步計(jì)算AUC曲線篩選調(diào)控網(wǎng)絡(luò)。
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.feather
tbl=motifs-v9-nr.hgnc-m0.001-o0.0.tbl
pyscenic=/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蛾茉，運(yùn)行線程數(shù)讼呢，設(shè)置越大越快，根據(jù)實(shí)際情況設(shè)置
sh pyscenic_from_loom.sh -i out.loom -n 10
第四步谦炬，計(jì)算RSS
此步驟的作用是通過(guò)計(jì)算RSS值找到細(xì)胞類型的特異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)
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(celltype,"_rss.rds"))
使用示例：
-i，第三步得到aucell.loom文件
-m键思，第一步得到的metadata_subset.xls文件
-c础爬，計(jì)算RSS的分組列
Rscript calcRSS_by_scenic.R -l aucell.loom -m metadata_subset.xls -c groups
運(yùn)行后生成regulon_RSS.Rdata和groups_rss.rds文件
