auc_matrix is a zero matrix

[linanzhang]

Hello,

I followed strictly the pipeline. Every step returns a reasonable output, except that

auc_mtx = aucell(ex_matrix, regulons, num_workers=4)

always returns a completely zero matrix.

I've tried both v9 and v10 resources/databases, both mm and hg datasets, and both linux and ios. Where could be wrong?

Thank you.

[dmalzl]

I have the exact same issue but only if I set noweights = True. So would also appreciate an answer from the devs here

[dmalzl]

did you follow the tutorial with the command line interface. I needed a more interactive approach and used this notebook from their git repo as blueprint the code I ended up using looks as follows (note this is for human):

1. Download needed resources

wget -O ../scenic_resource/hs_hgnc_curated_tfs.txt https://raw.githubusercontent.com/aertslab/pySCENIC/master/resources/hs_hgnc_tfs.txt
wget -O ../scenic_resource/motifs-v9-nr.hgnc-m0.001-o0.0.tbl https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
wget -O ../scenic_resource/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather 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

2. run SCENIC

# import all the necessary functionality
import anndata as ad

from arboreto.utils import load_tf_names
from arboreto.algo import grnboost2

from ctxcore.rnkdb import FeatherRankingDatabase as RankingDatabase
from pyscenic.utils import modules_from_adjacencies
from pyscenic.prune import prune2df, df2regulons
from pyscenic.aucell import aucell

# read data from AnnData object
# and make sure the layer you want is active
adata = ad.read_h5ad('/path/to/anndata/file.h5ad')
adata.X = adata.layers['counts'].copy()

# read TF names from previously downloaded file
with open('../scenic_resource/hs_hgnc_tfs.txt', 'r') as tf_file:
    tf_names = [line.rstrip() for line in tf_file]


# run grnboost on your local machine specifying the desired number of cores
from distributed import Client, LocalCluster
client = Client(
    LocalCluster(
        name = 'grn_call',
        n_workers = 16,
        threads_per_worker = 1
    )
)

expr_data = adata.to_df()
adjacencies = grnboost2(
    expression_data = expr_data,
    tf_names = tf_names,
    client_or_address = client,
    verbose = True
)
adjacencies.to_csv(
    '../scenic_results/adj.tsv',
    sep = '\t'
)
    
client.close()
del client

# compute regulons from adjacencies
inferred_modules = modules_from_adjacencies(
    adjacencies,
    adata.to_df(),
    # according to pyscenic warning
    # this masks all 0 expression entries for a particular gene for the correlation computation
    rho_mask_dropouts = True
)

# read cisTarget database for subsequent pruning of the regulons
cistarget_db = RankingDatabase(
    '../scenic_resource/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather',
    'hg38__refseq-r80__10kb_up_and_down_tss.mc9nr'
)

# prune modules based on regulatory sequence occurrence up/downstream  of gene
# note here I use 8 less cores because each worker maintains its own copy of the database
# which bloats the memory (8 workers requried something around 60 GB of RAM for me)
pruned_modules = prune2df(
    [cistarget_db], 
    inferred_modules, 
    '../scenic_resource/motifs-v9-nr.hgnc-m0.001-o0.0.tbl',
    client_or_address = 'custom_multiprocessing',
    num_workers = 8
)

# here is where it went wrong for me
# pruned_modules contains columns that contain lists of genes
# writing this to a file and then reading it back in does not reconstitute this list but keeps it as a string
# which is silently ignored by the AUCell algorithm and will just give you a warning that less than 80% of genes
# of the regulon are found in the expression data, resulting in an all-zero matrix
# to avoid this you can do something like this when reading back in 
# (note `eval` is generally regarded as unsafe but I used it here for simplicity and because I am lazy)
def read_pruned_modules(filename):
#     df = pd.read_csv(
#         filename, 
#         sep = '\t', 
#         index_col = [0, 1],
#         header = [0, 1]
#     ) 
#     # needs to be reconstituted to actual list and set
#     # otherwise aucell won't work
#     for column in ['Context', 'TargetGenes']:
#         col_data = df.loc[:, ('Enrichment', column)]
#         df.loc[:, ('Enrichment', column)] = col_data.apply(eval)
        
#     return df

pruned_modules.to_csv(
    '../scenic_results/pruned_modules.tsv',
    sep = '\t'
)

# prune2df returns a DataFrame object but the AUCell algorithm requires a list of Regulon objects
# so we have to convert it to this format using this function
pruned_regulons = df2regulons(
    pruned_modules
)

# finally compute AUCell matrix setting `noweights = True`
# this is necessary because it will otherwise factor in gene expression as weights
# which, at least in my case` results in unusable AUCell scores (some are negative 
# some are very positive but none are comparable) so setting this parameter ensures that 
# you get values in the expected 0 - 1 range (since AUCell is based on a custom definition of the ROC)
auc_mtx = aucell(
    adata.to_df(), 
    pruned_regulons, 
    num_workers = 16,
    noweights = True
)

[dmalzl]

btw. this is pySCENIC 0.12.1 but you have to install it from the git source because the PyPI version has a bug that is fixed in the git but not propagated to PyPI because they did not bump the version and in turn probably could not upload the package to PyPI because you can't overwrite existing packages with the same version (at least I don't know how)

[linanzhang]

Hi dmalzl,

Thank you for the reply. I've figured out how the issue is from. After running

regulons = df2regulons(df)

I should have

Regulon(name='ARID3A(+)', gene2weight=frozendict.frozendict({'ARID3A': 1.0, 'TMEM176B': 1.0, ...

However, instead I had

Regulon(name='ARID3A(+)', gene2weight=frozendict.frozendict({'[': 1.0, '(': 1.0, "'": 1.0, 'A': 1.0, 'R': 1.0, 'I': 1.0, 'D': 1.0, '3': 1.0, ...

There must be something wrong with the string split in df2regulons.

[linanzhang]

@aertslab/pySCENIC#509

[dmalzl]

Yeah. Same thing happened to me when I read in the table after saving to file. I think it is because df2regulons expects a list of tuples but pandas does not automatically evaluate the saved string to a list. So i had to do this manually. Glad my answer though :)