Gene regulatory analysis with pySCENIC

This exercise has been largely borrowed from the pySCENIC tutorial and adapted to work with a dataset of 10X PBMCs which you have downloaded with the course data. As it takes a long time run, we suggest if you are curious in the method to run it after the course.

Related paper – SCENIC: single-cell regulatory network inference and clustering.

Download Presentation: SCENIC Overview

Once single-cell genomics data has been processed, one can dissect important relationships between observed features in their genome context. In our genome, the activation of genes is controlled in the nucleus by the RNA transcriptional machinery, which activates local (promoters) or distal cis-regulatory elements (enhancers), to control the amount of RNA produced by every gene.

Conceptually, a Gene Regulatory Network (GRN) refers to a graph representation of how certain genes that control transcription i.e. “Transcription Factors” (TF) are in charge of directly controlling the transcription rates of their target genes (cis-regulation). At the same time, such target genes once activated, can be in charge of controlling other downstream target genes (trans-regulation). Computationally, methods that infer GRNs consider the co-variation of gene and chromatin accessibility features to identify modules that could be grouped and associated simultaneously with a few TFs. A group of genes controlled by the activity of the same TF is defined as a regulon.

In addition to co-variation, several methods have recognized the gathering and injection of prior knowledge data, such as the locations where a TF bind in the genome, or previously reported TF-target gene associations, to pre-define gene-gene edges that guide the inference of GRNs that are most-supported by that type of evidence.

The purpose of this exercise is to showcase the ability to infer information about gene regulatory networks using transcriptomic data alone: by pairing expression of transcription factors (TFs) with co-expressed genes, we can predict which genes are likely being activate by a particular TF in a specific cell population.

pySCENIC

# import dependencies
import os
import numpy as np
import pandas as pd
import scanpy as sc
import loompy as lp

import seaborn as sns
import matplotlib.pyplot as plt
import glob

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:

Load dataset

adata = sc.read_10x_h5("course_data/pbmc_10k_v3_filtered_feature_bc_matrix.h5")

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/anndata/_core/anndata.py:1840: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/anndata/_core/anndata.py:1840: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")

adata.var_names_make_unique()

nCountsPerGene = np.sum(adata.X, axis=0)
nCellsPerGene = np.sum(adata.X>0, axis=0)

# Show info
print("Number of counts (in the dataset units) per gene:", nCountsPerGene.min(), " - " ,nCountsPerGene.max())
print("Number of cells in which each gene is detected:", nCellsPerGene.min(), " - " ,nCellsPerGene.max())

Number of counts (in the dataset units) per gene: 0.0  -  3567008.0
Number of cells in which each gene is detected: 0  -  11766

nCells=adata.X.shape[0]

# pySCENIC thresholds
minCountsPerGene=3*.01*nCells # 3 counts in 1% of cells
print("minCountsPerGene: ", minCountsPerGene)

minSamples=.01*nCells # 1% of cells
print("minSamples: ", minSamples)

minCountsPerGene:  353.07
minSamples:  117.69

# simply compute the number of genes per cell (computers '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.sum(
    adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
# add the total counts per cell as observations-annotation to adata
adata.obs['n_counts'] = adata.X.sum(axis=1).A1

# initial cuts
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, :]

# save a copy of the raw data
adata.raw = adata

# Total-count normalize (library-size correct) to 10,000 reads/cell
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)

# log transform the data.
sc.pp.log1p(adata)

# identify highly variable genes.
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)

# keep only highly variable genes:
adata = adata[:, adata.var['highly_variable']]

# regress out total counts per cell and the percentage of mitochondrial genes expressed
sc.pp.regress_out(adata, ['n_counts', 'percent_mito'] ) #, n_jobs=args.threads)

# scale each gene to unit variance, clip values exceeding SD 10.
sc.pp.scale(adata, max_value=10)

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/anndata/_core/anndata.py:1230: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.
  df[key] = c
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/anndata/_core/anndata.py:1230: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.
  df[key] = c

# principal component analysis
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True)

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)

sc.tl.louvain(adata,resolution=0.4)
sc.pl.umap(adata, color=['louvain'] )

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:394: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(

sc.tl.rank_genes_groups(adata, 'louvain', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

WARNING: It seems you use rank_genes_groups on the raw count data. Please logarithmize your data before calling rank_genes_groups.

Step 1: Gene regulatory network inference, and generation of co-expression modules

GRN inference using the GRNBoost2 algorithm: for this step the CLI version of SCENIC is used. This step can be deployed on an High Performance Computing system. We use the counts matrix (without log transformation or further processing) from the loom file we wrote earlier.

Output: List of adjacencies between a TF and its targets stored in ADJACENCIES_FNAME.

## this file has to be downloaded if not found
!wget -nc https://raw.githubusercontent.com/aertslab/SCENICprotocol/master/example/allTFs_hg38.txt

File ‘allTFs_hg38.txt’ already there; not retrieving.

f_tfs = "allTFs_hg38.txt"
f_loom_path_scenic = "pbmc10k_filtered_scenic.loom"

# create basic row and column attributes for the loom file:
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( f_loom_path_scenic, adata.X.transpose(), row_attrs, col_attrs)

!pyscenic grn {f_loom_path_scenic} {f_tfs} -o adj.csv --num_workers 20

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:

2024-05-28 21:54:30,489 - pyscenic.cli.pyscenic - INFO - Loading expression matrix.

2024-05-28 21:54:31,078 - pyscenic.cli.pyscenic - INFO - Inferring regulatory networks.
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:
preparing dask client
parsing input
creating dask graph
20 partitions
computing dask graph
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/distributed/client.py:3108: UserWarning: Sending large graph of size 21.19 MiB.
This may cause some slowdown.
Consider scattering data ahead of time and using futures.
  warnings.warn(
not shutting down client, client was created externally
finished

2024-05-28 21:58:19,684 - pyscenic.cli.pyscenic - INFO - Writing results to file.

adjacencies = pd.read_csv("adj.csv", index_col=False)

Visualize the distribution of weights for general inspection of the quantiles and thresholds obtained from pyscenic. As provided by the pyscenic grn step, the importance scores follow a unimodal distribution, with negative/positive values indicating TF-gene associations with less/more importance, respectively. From the right-tail of this distribution, we can recover the most relevant interactions between TFs and potential target genes, supported by gene expression values and the analysis done by pyscenic.

adjacencies.head()

	TF	target	importance
0	CEBPD	VCAN	33.587173
1	ZEB2	LTB	33.086506
2	KLF4	VCAN	29.961844
3	CEBPD	SRGN	29.306472
4	MEF2C	HLA-DRA	28.519784

plt.hist(np.log10(adjacencies["importance"]), bins=50)
plt.xlim([-10, 10])

As targets genes have DNA motifs at promoters (sequence specific DNA motifs), those can be used to link TFs to target genes. Next, we use an annotation of TF associations to Transcription Start Sites (TSSs) to refine this annotation.

Step 2-3: Regulon prediction aka cisTarget from CLI

For this step the CLI version of SCENIC is used. This step can be deployed on an High Performance Computing system.

Output: List of adjacencies between a TF and its targets stored in MOTIFS_FNAME.

Download TSS annotations precalculated by Aerts’s lab:

!wget -nc https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather

File ‘hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather’ already there; not retrieving.

# ranking databases
db_glob = "*feather"
db_names = " ".join(glob.glob(db_glob))

db_names

'hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather'

Download a catalog of motif-to-TF associations

!wget -nc https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl

File ‘motifs-v9-nr.hgnc-m0.001-o0.0.tbl’ already there; not retrieving.

# motif databases
motif_path = "motifs-v9-nr.hgnc-m0.001-o0.0.tbl"

!pyscenic ctx adj.csv \
    'hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather' \
    --annotations_fname "motifs-v9-nr.hgnc-m0.001-o0.0.tbl" \
    --expression_mtx_fname "pbmc10k_filtered_scenic.loom" \
    --output reg.csv \
    --mask_dropouts \
    --num_workers 3 > pyscenic_ctx_stdout.txt

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:

2024-05-28 21:58:28,417 - pyscenic.cli.pyscenic - INFO - Creating modules.

2024-05-28 21:58:28,483 - pyscenic.cli.pyscenic - INFO - Loading expression matrix.

2024-05-28 21:58:29,374 - pyscenic.utils - INFO - Calculating Pearson correlations.

2024-05-28 21:58:29,393 - pyscenic.utils - WARNING - Note on correlation calculation: the default behaviour for calculating the correlations has changed after pySCENIC verion 0.9.16. Previously, the default was to calculate the correlation between a TF and target gene using only cells with non-zero expression values (mask_dropouts=True). The current default is now to use all cells to match the behavior of the R verision of SCENIC. The original settings can be retained by setting 'rho_mask_dropouts=True' in the modules_from_adjacencies function, or '--mask_dropouts' from the CLI.
    Dropout masking is currently set to [True].

2024-05-28 21:58:36,405 - pyscenic.utils - INFO - Creating modules.

2024-05-28 21:58:43,251 - pyscenic.cli.pyscenic - INFO - Loading databases.

2024-05-28 21:58:43,507 - pyscenic.cli.pyscenic - INFO - Calculating regulons.

2024-05-28 21:58:43,507 - pyscenic.prune - INFO - Using 3 workers.

2024-05-28 21:58:43,507 - pyscenic.prune - INFO - Using 3 workers.

2024-05-28 21:58:46,654 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): database loaded in memory.

2024-05-28 21:58:46,654 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): database loaded in memory.

2024-05-28 21:58:46,657 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): database loaded in memory.

2024-05-28 21:58:46,657 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): database loaded in memory.

2024-05-28 21:58:46,718 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): database loaded in memory.

2024-05-28 21:58:46,718 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): database loaded in memory.

2024-05-28 21:58:47,490 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): motif annotations loaded in memory.

2024-05-28 21:58:47,490 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): motif annotations loaded in memory.

2024-05-28 21:58:47,491 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): motif annotations loaded in memory.

2024-05-28 21:58:47,491 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): motif annotations loaded in memory.

2024-05-28 21:58:47,491 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): motif annotations loaded in memory.

2024-05-28 21:58:47,491 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): motif annotations loaded in memory.

2024-05-28 21:59:02,767 - pyscenic.transform - WARNING - Less than 80% of the genes in KLF12 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 21:59:04,752 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for KLF12 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 21:59:28,614 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for EGR1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 21:59:34,853 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for MYC could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 21:59:36,692 - pyscenic.transform - WARNING - Less than 80% of the genes in PBX4 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 21:59:45,671 - pyscenic.transform - WARNING - Less than 80% of the genes in RPS6KA5 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:00:25,006 - pyscenic.transform - WARNING - Less than 80% of the genes in ZNF595 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:00:37,943 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for MYC could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:00:47,258 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for PBX4 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:00:50,970 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for ZNF367 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:00:51,277 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for ZNF527 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:01:07,842 - pyscenic.transform - WARNING - Less than 80% of the genes in BNC2 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:01:18,465 - pyscenic.transform - WARNING - Less than 80% of the genes in CYB5R1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:01:22,986 - pyscenic.transform - WARNING - Less than 80% of the genes in EBF1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:01:30,579 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for ZNF595 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:01:30,774 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for ZNF83 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:02:02,012 - pyscenic.transform - WARNING - Less than 80% of the genes in KLF12 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:02:18,349 - pyscenic.transform - WARNING - Less than 80% of the genes in MYC could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:02:23,860 - pyscenic.transform - WARNING - Less than 80% of the genes in NPDC1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:03:15,944 - pyscenic.transform - WARNING - Less than 80% of the genes in ZEB1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:03:16,218 - pyscenic.transform - WARNING - Less than 80% of the genes in ZNF527 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:03:26,272 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): All regulons derived.

2024-05-28 22:03:26,272 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): All regulons derived.

2024-05-28 22:03:26,413 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): Done.

2024-05-28 22:03:26,413 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(3): Done.

2024-05-28 22:03:51,601 - pyscenic.transform - WARNING - Less than 80% of the genes in EGR1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:03:58,267 - pyscenic.transform - WARNING - Less than 80% of the genes in ESR2 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:04:12,091 - pyscenic.transform - WARNING - Less than 80% of the genes in Regulon for EBF1 could be mapped to hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings. Skipping this module.

2024-05-28 22:04:16,918 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): All regulons derived.

2024-05-28 22:04:16,918 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): All regulons derived.

2024-05-28 22:04:16,977 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): Done.

2024-05-28 22:04:16,977 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(2): Done.

2024-05-28 22:04:46,928 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): All regulons derived.

2024-05-28 22:04:46,928 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): All regulons derived.

2024-05-28 22:04:47,025 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): Done.

2024-05-28 22:04:47,025 - pyscenic.prune - INFO - Worker hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings(1): Done.

2024-05-28 22:04:47,095 - pyscenic.cli.pyscenic - INFO - Writing results to file.

To explore the candidates reported, it is recommended as a rule of thumb to explore the output by the ranking of relative contribution, or by top-quantile thresholds defined visually to obtain a high signal-to-noise ratio.

Define custom quantiles for further exploration

import numpy as np

n_genes_detected_per_cell = np.sum(adata.X > 0, axis=1)
percentiles = pd.Series(n_genes_detected_per_cell.flatten()).quantile(
    [0.01, 0.05, 0.10, 0.50, 1]
)
print(percentiles)

0.01     136.0
0.05     155.0
0.10     169.0
0.50     237.0
1.00    1017.0
dtype: float64

The histogram below indicates the distribution of genes detected per cell. This visualization is convenient to define the parameter –auc_threshold in the next step. Specifically, the default parameter of –auc_threshold is 0.05, which in this plot would result in the selection of 144 genes, to be used as a reference per cell for AUCell calculations. The modification of this parameter affects the estimation of AUC values calculated by AUCell.

fig, ax = plt.subplots(1, 1, figsize=(8, 5), dpi=100)
sns.distplot(n_genes_detected_per_cell, norm_hist=False, kde=False, bins="fd")
for i, x in enumerate(percentiles):
    fig.gca().axvline(x=x, ymin=0, ymax=1, color="red")
    ax.text(
        x=x,
        y=ax.get_ylim()[1],
        s=f"{int(x)} ({percentiles.index.values[i]*100}%)",
        color="red",
        rotation=30,
        size="x-small",
        rotation_mode="anchor",
    )
ax.set_xlabel("# of genes")
ax.set_ylabel("# of cells")
fig.tight_layout()

/tmp/ipykernel_279481/1384461256.py:2: UserWarning: 

`distplot` is a deprecated function and will be removed in seaborn v0.14.0.

Please adapt your code to use either `displot` (a figure-level function with
similar flexibility) or `histplot` (an axes-level function for histograms).

For a guide to updating your code to use the new functions, please see
https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751

  sns.distplot(n_genes_detected_per_cell, norm_hist=False, kde=False, bins="fd")

This step will use TFs to calculate Area Under the Curve scores, that summarize how well the gene expression observed in each cell can be associated by the regulation of target genes regulatred by the mentioned TFs.

Using the above-generated matrix of cell x TFs and those scores, we can calculate a new embedding using only those.

Step 4: Cellular enrichment (aka AUCell) from CLI

It is important to check that most cells have a substantial fraction of expressed/detected genes in the calculation of the AUC. The following histogram gives an idea of the distribution and allows selection of an appropriate threshold. In this plot, a few thresholds are highlighted, with the number of genes selected shown in red text and the corresponding percentile in parentheses). See the relevant section in the R tutorial for more information.

By using the default setting for –auc_threshold of 0.05, we see that 1192 genes are selected for the rankings based on the plot below.

!pyscenic aucell "pbmc10k_filtered_scenic.loom" \
    "reg.csv" \
    --output "auc.csv" \
    --num_workers 3 > "pyscenic_aucell_stdout.txt"

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_to_dna(twobit: int, size: int) -> str:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def dna_to_twobit(dna: str) -> int:
/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def twobit_1hamming(twobit: int, size: int) -> List[int]:

2024-05-28 22:04:55,503 - pyscenic.cli.pyscenic - INFO - Loading expression matrix.

2024-05-28 22:04:56,091 - pyscenic.cli.pyscenic - INFO - Loading gene signatures.

2024-05-28 22:04:57,373 - pyscenic.cli.pyscenic - INFO - Calculating cellular enrichment.

2024-05-28 22:05:05,478 - pyscenic.cli.pyscenic - INFO - Writing results to file.

auc_df = pd.read_csv("auc.csv", index_col=0)

auc_df.head()

	ARID5B(+)	BACH2(+)	BATF(+)	BCL11A(+)	BCL6(+)	CEBPB(+)	CEBPD(+)	CREB5(+)	E2F1(+)	E2F2(+)	...	STAT1(+)	TAL1(+)	TBX21(+)	TCF4(+)	TCF7(+)	TCF7L2(+)	USF2(+)	ZNF595(+)	ZNF831(+)	ZNF85(+)
Cell
AAACCCAAGCGCCCAT-1	0.035164	0.000000	0.047850	0.013579	0.024272	0.018955	0.016766	0.026724	0.013183	0.036084	...	0.034538	0.017201	0.041977	0.000000	0.053187	0.015954	0.024688	0.035599	0.035599	0.000000
AAACCCACAGAGTTGG-1	0.004977	0.010264	0.023925	0.013579	0.156211	0.087353	0.074413	0.095586	0.005611	0.101618	...	0.049658	0.030772	0.017486	0.000000	0.000000	0.066439	0.006935	0.000000	0.000000	0.000000
AAACCCACAGGTATGG-1	0.015828	0.004716	0.176838	0.003637	0.000000	0.008990	0.030588	0.008107	0.023693	0.000000	...	0.011459	0.015724	0.082923	0.000000	0.000000	0.011178	0.000000	0.000000	0.049083	0.000000
AAACCCACATAGTCAC-1	0.087623	0.046325	0.000000	0.083342	0.000000	0.002954	0.003111	0.000000	0.013093	0.006796	...	0.000000	0.008041	0.011048	0.066163	0.000000	0.002572	0.034951	0.000000	0.026969	0.101942
AAACCCACATCCAATG-1	0.016399	0.046325	0.079404	0.006239	0.000871	0.009863	0.017225	0.005705	0.023158	0.000000	...	0.011141	0.005931	0.072210	0.012945	0.024905	0.021464	0.001664	0.000000	0.056634	0.037621

5 rows × 59 columns

ad_auc_mtx = sc.AnnData(auc_df)
sc.pp.neighbors(ad_auc_mtx, n_neighbors=10, metric="correlation")
sc.tl.umap(ad_auc_mtx)

WARNING: You’re trying to run this on 59 dimensions of `.X`, if you really want this, set `use_rep='X'`.
         Falling back to preprocessing with `sc.pp.pca` and default params.

adata.obsm["X_umap_aucell"] = ad_auc_mtx.obsm["X_umap"]

sc.pl.embedding(adata, basis="X_umap_aucell", color="louvain")

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:394: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(

sc.tl.louvain(ad_auc_mtx, resolution=0.2)
sc.pl.embedding(ad_auc_mtx, basis="X_umap", color='louvain')

/home/alex/anaconda3/envs/sctp/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:394: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(

sc.tl.rank_genes_groups(ad_auc_mtx, 'louvain', method='t-test')

marker_genes = pd.DataFrame(ad_auc_mtx.uns["rank_genes_groups"]["names"])
marker_genes.columns = ["Cluster" + str(x) for x in range(0, len(ad_auc_mtx.obs["louvain"].unique()))]
marker_genes.head()

	Cluster0	Cluster1	Cluster2	Cluster3	Cluster4	Cluster5	Cluster6
0	PRDM1(+)	SPI1(+)	SPIB(+)	MEOX1(+)	TCF7L2(+)	TAL1(+)	CEBPB(+)
1	EOMES(+)	NFE2(+)	BCL11A(+)	LEF1(+)	SPI1(+)	GATA1(+)	SPI1(+)
2	TBX21(+)	CEBPD(+)	PAX5(+)	TCF7(+)	CEBPB(+)	MAFG(+)	NFE2(+)
3	RUNX3(+)	CEBPB(+)	IRF8(+)	KLF12(+)	MAFB(+)	NFE2(+)	GATA3(+)
4	KLF12(+)	MAFB(+)	IRF4(+)	BACH2(+)	E2F1(+)	E2F1(+)	CEBPD(+)

sc.pl.embedding(ad_auc_mtx, basis="X_umap", color=['SPI1(+)',
                                                  'SPIB(+)', 'MEOX1(+)', 'TCF7L2(+)'], ncols=2)

regulon_df = pd.read_csv("reg.csv", index_col=0)

regulon_df.head()

	Unnamed: 1	Enrichment	Enrichment.1	Enrichment.2	Enrichment.3	Enrichment.4	Enrichment.5	Enrichment.6	Enrichment.7
NaN	NaN	AUC	NES	MotifSimilarityQvalue	OrthologousIdentity	Annotation	Context	TargetGenes	RankAtMax
TF	MotifID	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
ARID5B	cisbp__M0110	0.09194977843426884	3.4110817264013087	0.0	0.880471	gene is orthologous to ENSMUSG00000019947 in M...	frozenset({'weight>75.0%', 'activating', 'hg38...	[('ID3', 0.8384963869548077), ('KAT6B', 0.3927...	1113
ARID5B	cisbp__M0107	0.0899019739492413	3.2853287808894263	0.0	0.880471	gene is orthologous to ENSMUSG00000019947 in M...	frozenset({'weight>75.0%', 'activating', 'hg38...	[('NXPH4', 0.883027200867596), ('TRIB2', 0.556...	818
BATF	dbcorrdb__RCOR1__ENCSR000ECM_1__m1	0.09647906429927953	3.295409378378749	0.000478	1.0	gene is annotated for similar motif factorbook...	frozenset({'weight>75.0%', 'activating', 'hg38...	[('NOL4L', 0.4502441398393177), ('LLGL2', 0.40...	1072

regulon_df.shape

(4115, 9)

regulon_df.loc["SPIB"]

	Unnamed: 1	Enrichment	Enrichment.1	Enrichment.2	Enrichment.3	Enrichment.4	Enrichment.5	Enrichment.6	Enrichment.7
SPIB	dbcorrdb__IRF1__ENCSR000EGK_1__m1	0.07806698870214379	3.90354587822946	0.000965	0.82397	motif similar to hocomoco__SPIB_MOUSE.H11MO.0....	frozenset({'weight>75.0%', 'activating', 'hg38...	[('FCGR2B', 1.869718459857277), ('STX7', 6.892...	4932
SPIB	cisbp__M6313	0.07676953171783304	3.77165735433781	0.000466	0.82397	motif similar to hocomoco__SPIB_MOUSE.H11MO.0....	frozenset({'weight>75.0%', 'activating', 'hg38...	[('PNOC', 7.0142442259089695), ('PEG10', 1.095...	856
SPIB	dbcorrdb__TBL1XR1__ENCSR000DYZ_1__m1	0.07155724779432313	3.241820572780774	3.1e-05	1.0	motif similar to hocomoco__SPIB_HUMAN.H11MO.0....	frozenset({'weight>75.0%', 'activating', 'hg38...	[('BTK', 3.20425594834044), ('CYB561A3', 9.741...	2320
SPIB	cisbp__M4489	0.08201924228512117	4.305298612545568	1e-06	1.0	gene is annotated for similar motif homer__AAA...	frozenset({'weight>75.0%', 'activating', 'hg38...	[('SSBP2', 0.869120775272342), ('ARID5B', 2.30...	1221
SPIB	hocomoco__SPIB_MOUSE.H11MO.0.A	0.0760484450477065	3.698357770867258	0.0	0.839695	gene is orthologous to ENSMUSG00000008193 in M...	frozenset({'weight>75.0%', 'activating', 'hg38...	[('FCGR2B', 1.869718459857277), ('TMEM156', 4....	4732
...	...	...	...	...	...	...	...	...	...
SPIB	taipale__Spic_DBD_AAAAAGMGGAAGTA	0.07192604342454488	4.096439803073309	0.0	1.0	gene is annotated for similar motif taipale__S...	frozenset({'activating', 'hg38__refseq-r80__10...	[('ALOX5AP', 0.2436475904030507), ('DOPEY2', 4...	3194
SPIB	hocomoco__IRF8_HUMAN.H11MO.0.B	0.08214495894116133	5.281967790108696	0.0	1.0	motif similar to hocomoco__SPIB_HUMAN.H11MO.0....	frozenset({'activating', 'hg38__refseq-r80__10...	[('ALOX5AP', 0.2436475904030507), ('P2RY14', 1...	3184
SPIB	dbcorrdb__SPI1__ENCSR000BIJ_1__m1	0.0709653406500041	3.9849857036056426	0.0	1.0	motif similar to hocomoco__SPIB_HUMAN.H11MO.0....	frozenset({'activating', 'hg38__refseq-r80__10...	[('KCNG1', 0.97964022364969), ('CXXC5', 8.0650...	581
SPIB	dbcorrdb__EP300__ENCSR000DZG_1__m1	0.07342458385785497	4.270290112384945	0.0	1.0	motif similar to hocomoco__SPIB_HUMAN.H11MO.0....	frozenset({'activating', 'hg38__refseq-r80__10...	[('TCF4', 3.091736175722826), ('MEF2C', 7.5755...	1982
SPIB	hocomoco__IRF4_HUMAN.H11MO.0.A	0.08185655324008216	5.24850895294599	0.0	1.0	motif similar to hocomoco__SPIB_HUMAN.H11MO.0....	frozenset({'activating', 'hg38__refseq-r80__10...	[('MYO1C', 1.9761547224955385), ('TPM2', 0.319...	993

143 rows × 9 columns

regulon_df.loc["SPIB"].iloc[0]["Enrichment.6"][:1000]

"[('FCGR2B', 1.869718459857277), ('STX7', 6.892725862395004), ('P2RY14', 1.8778070177222), ('TMEM156', 4.272053543910822), ('BLNK', 5.738561565806354), ('CD38', 0.5724236819241179), ('CXCR4', 1.0237407733006956), ('RIMKLB', 3.814839002835818), ('FCRL3', 1.6683434765092275), ('LINC01215', 2.4484182394885634), ('HLA-DMB', 5.275150829219921), ('BCL11A', 5.323875629596565), ('HLA-DQA1', 12.631876848026147), ('TCF4', 3.091736175722826), ('CD40', 5.98234166462111), ('DOPEY2', 4.554553673574461), ('AFF3', 3.926454150323302), ('CXXC5', 8.0650414459558), ('LPAR5', 1.3386477473776452), ('KLF8', 1.6187543679806742), ('SLC44A2', 0.3937692393269603), ('SHISA8', 0.9218564482844408), ('FCRL5', 2.534706255811908), ('IKZF3', 1.4164772622710846), ('MPP6', 0.7198674632752087), ('PLEKHG1', 3.1788634402443443), ('FMNL3', 0.4906561541272032), ('PKIG', 3.768067885404645), ('CCR6', 2.1375793714457263), ('PLAC8', 3.411256052315926), ('BLK', 7.363519956872278), ('SETBP1', 5.3370810233327015), ('HLA-DQA2', 4.0790"

sc.pl.embedding(adata, basis="X_umap_aucell", use_raw=False,
                color=['FCGR2B', 'STX7', 'P2RY14', 'TMEM156'], ncols=2)