注：pyscenic发表在Nature protocol上，题为：A scalable SCENIC workflow for single-cell gene regulatory network analysis
在原文中有完整的环境配置流程和示例代码，完全可以自己无障碍完成示例代码的运行
但是!我不会python。。。。。因此在我完全学会python之前，需要通过R的代码进行必要的前期的文件处理！
我们先展示一下protocol上的示例代码（使用scanpy处理）
#1、Perform cleaning and quality control on the downloaded expression data set. We use the Scanpy
#(https://scanpy.readthedocs.io/en/latest/) package in this example and run the following statements
#in a Python interpreter session
import scanpy as sc
adata = sc.read_10x_mtx('filtered_feature_bc_matrix/',var_names='gene_symbols')
adata.var_names_make_unique()
# compute the number of genes per cell (computes ‘n_genes' column)
sc.pp.filter_cells(adata, min_genes=0)
# mito and genes/counts cuts
mito_genes = adata.var_names.str.startswith('MT-')
# for each cell compute fraction of counts in mito genes vs. all genes
adata.obs['percent_mito'] = np.ravel(np.sum(
adata[:, mito_genes].X, axis=1)) / np.ravel(np.sum(adata.X,
axis=1))
# add the total counts per cell as observations-annotation to adata
adata.obs['n_counts'] = np.ravel(adata.X.sum(axis=1))
#2、Carry out the filtering steps with basic thresholds for genes and cells
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs['n_genes'] < 4000, :]
adata = adata[adata.obs['percent_mito'] < 0.15, :]
#3、Create a loom file with the filtered data, to be used in downstream analyses
import loompy as lp
row_attrs = {
"Gene": np.array(adata.var_names),
col_attrs = {
"CellID": np.array(adata.obs_names),
"nGene": np.array(np.sum(adata.X.transpose()>0, axis=0)).flatten(),
"nUMI": np.array(np.sum(adata.X.transpose(),axis=0)).flatten(),
lp.create("PBMC10k_filtered.loom",adata.X.transpose(),row_attrs,
col_attrs)
总结一下，基本流程是将数据读入后，先计算了mito基因的比例，然后筛选了细胞和基因以及mito基因的比例，最后再保存为loom文件。
因此如果我们在seurat中将数据按照我们自己的标准处理，然后保存 成loom应该也可以。
library(Seurat)
x<-Read10X("../pSCENIC/cisdatabase/filtered_feature_bc_matrix/")
dim(x)
x<-CreateSeuratObject(x,min.cells = 3,min.features = 200)
x[["mito"]]<-PercentageFeatureSet(x, pattern = "^MT-")
x<-subset(x,subset = nFeature_RNA < 4000 & mito<15)
#注意，这里不管转不转置，都会报下面的错
create(filename = "data.loom",data = t(as.matrix(x@assays$RNA@counts)))
报错信息：
Traceback (most recent call last):
File "/home/whq/miniconda3/envs/scenic_protocol/bin/pyscenic", line 8, in <module>
sys.exit(main())
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/pyscenic/cli/pyscenic.py", line 421, in main
args.func(args)
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/pyscenic/cli/pyscenic.py", line 49, in find_adjacencies_command
args.gene_attribute)
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/pyscenic/cli/utils.py", line 116, in load_exp_matrix
return load_exp_matrix_as_loom(fname, return_sparse, attribute_name_cell_id, attribute_name_gene)
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/pyscenic/cli/utils.py", line 69, in load_exp_matrix_as_loom
with lp.connect(fname,mode='r',validate=False) as ds:
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/loompy/loompy.py", line 1389, in connect
return LoomConnection(filename, mode, validate=validate)
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/loompy/loompy.py", line 84, in init
self._file = h5py.File(filename, mode)
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/h5py/_hl/files.py", line 427, in init
swmr=swmr)
File "/home/whq/miniconda3/envs/scenic_protocol/lib/python3.6/site-packages/h5py/_hl/files.py", line 190, in make_fid
fid = h5f.open(name, flags, fapl=fapl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 96, in h5py.h5f.open
OSError: Unable to open file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
所以我们使用R得到的loom应该与python得到的loom是不一样的！
因此我们只能用planB，先用seurat将矩阵筛选，然后输出成csv，用python读入，然后打包成loom
此处为了快读快写快运行，我就用了seurat中自带的示例对象进行演示
#注意矩阵一定要转置，不然会报错
write.csv(t(as.matrix(pbmc_small@assays$RNA@counts)),file = "sample.csv")
#以下为python代码
import loompy as lp
import numpy as np
import scanpy as sc
x=sc.read_csv("../../code/sample.csv")
row_attrs = {
"Gene": np.array(x.var_names),
col_attrs = {
"CellID": np.array(x.obs_names)
lp.create("sample.loom",x.X.transpose(),row_attrs,
col_attrs)
这样，我们在后面接pyscenic grn就可以了
pyscenic grn \
--num_workers 20 \
--output adj.sample.tsv \
--method grnboost2 \
sample1.loom \
hs_hgnc_tfs.txt
最后可以正常输出了,就是临界矩阵转化的module
TF target importance
ODC1 NGFRAP1 274.24821926691243
ODC1 SPARC 208.2001864783341
ODC1 PF4 206.9865692985902
ODC1 SDPR 191.50108534167583
ODC1 PPBP 157.22866071698584
ODC1 CLU 140.29117587982287
ODC1 NCOA4 138.34975956426598
ODC1 GP9 115.40677095419836
ODC1 RUFY1 113.59279900250596
然后再跑cistarget
pyscenic ctx \
adj.sample.tsv \
hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather \
--annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl \
--expression_mtx_fname sample1.loom \
--mode "dask_multiprocessing" \
--output reg.csv \
--num_workers 3 \
--mask_dropouts