前面我們用自帶例子和人類的pbmc數(shù)據(jù)測試了SCENIC的R版本的安裝和使用博投。但是，通常情況下runGenie3這一步盯蝴，稍微大一點(diǎn)點(diǎn)的數(shù)據(jù)集毅哗，運(yùn)行時(shí)間都很長。比如pbmc這個數(shù)據(jù)捧挺，還不算太大虑绵，這一步都運(yùn)行了近15個小時(shí)? 。
所以闽烙，我們今天測試一下pySCENIC翅睛。pySCENIC是通過python實(shí)現(xiàn)的一個快速的SCENIC pipeline，pySCENIC支持多線程運(yùn)行黑竞，如果你的服務(wù)器條件滿足的話捕发，可以大大的節(jié)省我們的運(yùn)算時(shí)間。
========利用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
========下載數(shù)據(jù)庫所需要的文件======
數(shù)據(jù)庫根據(jù)調(diào)節(jié)特征（即轉(zhuǎn)錄因子）對用戶感興趣的物種的整個基因組 進(jìn)行 排名扎酷。排名數(shù)據(jù)庫通常以feather格式存儲，可以從cisTargetDBs下載遏匆。和前面下載的一樣法挨。（需要注意的是：需注意版本問題谁榜，一個是 pySCENIC < 0.12.0 和ctxcore < 0.2.0軟件版本需要進(jìn)到old文件夾，8月份最近的更新凡纳，另一個則是基因組版本窃植，以及數(shù)據(jù)庫來源的版本）我這次下載的是：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注釋文件： 數(shù)據(jù)庫提供了豐富的motif 和結(jié)合該motif的轉(zhuǎn)錄因子之間的缺失鏈接。該pipeline需要一個 TSV 文本文件荐糜，其中每一行都代表一個特定的注釋巷怜。（這個需要重新下載，可以從這里下載相應(yīng)的文件：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
======命令行界面=====
下面主要參考這個帖子：http://www.reibang.com/p/71ed43163ce1
========第一步：從Seurat對象中提取相關(guān)的信息========
此步驟的目的是存儲表達(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對象的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。
運(yùn)行代碼如下：Rscript get_count_from_seurat.R -i test.rds -d groups -s 20 -l out
運(yùn)行后會生成3個文件：矩陣out.csv裙戏，metadata文件metadata_subset.xls和取子集后的RDS文件subset.rds（如果不取子集乘凸，這個文件不會生成）。
=========第二步：使用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版本==========
SCENIC的標(biāo)準(zhǔn)流程分為3步，前面的帖子已經(jīng)講過了：第一步利用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ò)。
生成一個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，運(yùn)行線程數(shù)量瓜，設(shè)置越大越快司恳，根據(jù)實(shí)際情況設(shè)置
./pyscenic_from_loom.sh -i out.loom -n 20
注：其中最耗時(shí)的也就是這一步，但是運(yùn)行pbmc的數(shù)據(jù)也只花了5分鐘左右绍傲，比如R版本已經(jīng)是快太多了扔傅。
========第四步：計(jì)算RSS=====
此步驟的作用是通過計(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)
#因?yàn)槲乙膊幌肷蓡蝹€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杠纵，計(jì)算RSS的分組列
Rscript calcRSS_by_scenic.R -l aucell.loom -m metadata_subset.xls -c cell_type
運(yùn)行后生成regulon_RSS.Rdata和cell_type_rss.rds文件
下面的帖子荠耽，我們學(xué)習(xí)一下可視化python版本的結(jié)果。
