Can run_pycistarget be used for non-model organisms?

[SidG13]

Hi, thanks for this great tool! I was running run_pycistarget, but realized I can't change the species from 'homo_sapiens' to my non-model, non-Ensemble species. I have made cistarget region rankings and scores databases already specific to my organism, with consensus peaks called as the regions.

I'm assuming there is technically support for non-model species, but I was just not clear how to proceed. Thanks!

[SeppeDeWinter]

Hi Sid!

You're right at the moment this function is not very non-model-organism friendly. However, the only reason why you need to provide the species name is to download the Transcription Starting Site annotations from biomart (this so promoters can be excluded from the analysis, which is optional). So you can actually run the analysis on your non-model organism, just not using this function (at this moment). I would suggest to follow the tutorials on https://pycistarget.readthedocs.io/en/latest/ for now.

Hope this helps.

We hope to make the functions a bit more non-model organism friendly in the future.

Best,

Seppe

[SidG13]

Hi, thanks for the quick response! I was following the tutorial you suggested, initially running the run_cistarget function, however I'm getting an issue where there are "No enriched motifs found for consensus_regions". I was wondering if I could get some guidance here since I'm sure there is an issue, and there are at least some enriched motifs.

Here's a bit of my regions_vs_motifs.rankings.feather file (shifted all the way to the right of the table) :

df.shape
(1956, 139151)

My motif2TF table, which is simply the motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl, subset for only JASPAR-sourced matrices.

And my regions table ("region_sets") object (just consensus peaks called from scATACseq in an earlier step):

The command I'm running is:

cistarget_dict = run_cistarget(ctx_db = '/home/administrator/sc_multiome/SCENICplus_multiome/for_cistarget/consensuspeaks.regions_vs_motifs.rankings.feather',
                                                      region_sets = region_sets,
                                                      specie = 'species',
                                                      auc_threshold = 0.005,
                                                      nes_threshold = 3.0,
                                                      rank_threshold = 0.05,
                                                      annotation = ['Direct_annot', 'Orthology_annot'],
                                                      annotation_version = 'v1nr_clust',
                                                      path_to_motif_annotations = '/home/administrator/sc_multiome/SCENICplus_multiome/for_cistarget/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0_just_jaspar.tbl',
                                                      n_cpu = 4,
                                                      _temp_dir='/home/administrator/tmp')

cistarget_results(cistarget_dict, name='consensus_regions')

The output is empty however:

Is there anything obvious I'm doing wrong? I've tried playing with the threshold parameters, but it doesn't seem to help. Thanks for your assistance!

[SeppeDeWinter]

Hi

There are two things:

1. motif enrichment analysis will not work when you use all consensus regions, you should take smaller region sets of interest containing regions which are co-accessible in your multiple celltypes of interest. This can be either topics or differentially accessible regions, or both. The reason is that you're actually testing which motifs are enriched in your regions of interest versus all the consensus regions.
2. Make sure that the motif names in the regions_vs_motifs.rankings.feather file (motifs column) match with the #motif_id column of your motif2TF table. From a first glance they don't seem to be matching. Can you check this?

Hope this helps!

Best,

Seppe

[SidG13]

Thanks again for your reply! Your suggestions do make sense, and sorry for the errors, but I fixed those two things and am still having issues, which I believe are now related to the region_sets.

I created my region_sets following the SCENIC+ tutorial, so it is a list with 3 dictionaries, each populated with pr.PyRanges objects:

>>>for key in region_sets.keys():
    print(f'{key}: {region_sets[key].keys()}')

topics_otsu: dict_keys(['Topic1', 'Topic2', 'Topic3', 'Topic4', 'Topic5', 'Topic6', 'Topic7', 'Topic8', 'Topic9', 'Topic10', 'Topic11', 'Topic12', 'Topic13', 'Topic14', 'Topic15', 'Topic16'])
topics_top_3: dict_keys(['Topic1', 'Topic2', 'Topic3', 'Topic4', 'Topic5', 'Topic6', 'Topic7', 'Topic8', 'Topic9', 'Topic10', 'Topic11', 'Topic12', 'Topic13', 'Topic14', 'Topic15', 'Topic16'])
DARs: dict_keys(['0', '1', '2', '3', '4', '5', '6', '7'])

I'm running the same command as above, now specifying new region sets and the modified motif2TF database but ran into an error when I try to run the entire region_sets list of dicts together:

2022-10-21 12:30:05,142 cisTarget    INFO     Reading cisTarget database
---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
Cell In [8], line 24
     21 from pycistarget.motif_enrichment_cistarget import *
     23 # Run, using precomputed database
---> 24 cistarget_dict = run_cistarget(ctx_db = '/home/administrator/Desktop/sc_multiome/SCENICplus_multiome/data_manipulation/for_cistarget/consensuspeaks.regions_vs_motifs.rankings.feather',
     25                                                       region_sets = region_sets,
     26                                                       specie = 'species',
     27                                                       auc_threshold = 0.005,
     28                                                       nes_threshold = 3.0,
     29                                                       rank_threshold = 0.05,
     30                                                       annotation = ['Direct_annot', 'Orthology_annot'],
     31                                                       annotation_version = 'v1nr_clust',
     32                                                       path_to_motif_annotations = '/home/administrator/Desktop/sc_multiome/SCENICplus_multiome/data_manipulation/for_cistarget/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0_just_jaspar2.tbl',
     33                                                       n_cpu = 4,
     34                                                       _temp_dir='/home/administrator/Desktop/tmp')

File ~/anaconda2/envs/py38/lib/python3.8/site-packages/pycistarget/motif_enrichment_cistarget.py:518, in run_cistarget(ctx_db, region_sets, specie, name, fraction_overlap, auc_threshold, nes_threshold, rank_threshold, path_to_motif_annotations, annotation_version, annotation, motif_similarity_fdr, orthologous_identity_threshold, n_cpu, motifs_to_use, **kwargs)
    516 # Load database
    517 if isinstance(ctx_db, str):
--> 518     ctx_db = cisTargetDatabase(ctx_db,
    519                          region_sets,
    520                          name = name,
    521                          fraction_overlap = fraction_overlap)
    523 # Run cistarget analysis in parallel
    524 if n_cpu > 1:

File ~/anaconda2/envs/py38/lib/python3.8/site-packages/pycistarget/motif_enrichment_cistarget.py:67, in cisTargetDatabase.__init__(self, fname, region_sets, name, fraction_overlap)
     48 def __init__(self, 
     49             fname: str,
     50             region_sets: Union[Dict[str, pr.PyRanges], pr.PyRanges] = None,
     51             name: str = None,
     52             fraction_overlap: float = 0.4):
     53     """
     54     Initialize cisTargetDatabase
     55     
   (...)
     65         Minimal overlap between query and regions in the database for the mapping.     
     66     """
---> 67     self.regions_to_db, self.db_rankings, self.total_regions = self.load_db(fname,
     68                                                       region_sets,
     69                                                       name,
     70                                                       fraction_overlap)

File ~/anaconda2/envs/py38/lib/python3.8/site-packages/pycistarget/motif_enrichment_cistarget.py:119, in cisTargetDatabase.load_db(self, fname, region_sets, name, fraction_overlap)
    117 if region_sets is not None:
    118     if type(region_sets) == dict:
--> 119         target_to_db_dict = {x: target_to_query(region_sets[x], list(db_regions), fraction_overlap = fraction_overlap) for x in region_sets.keys()}
    120         target_regions_in_db = list(set(sum([target_to_db_dict[x]['Query'].tolist() for x in target_to_db_dict.keys()],[])))
    121     elif type(region_sets) == pr.PyRanges:

File ~/anaconda2/envs/py38/lib/python3.8/site-packages/pycistarget/motif_enrichment_cistarget.py:119, in <dictcomp>(.0)
    117 if region_sets is not None:
    118     if type(region_sets) == dict:
--> 119         target_to_db_dict = {x: target_to_query(region_sets[x], list(db_regions), fraction_overlap = fraction_overlap) for x in region_sets.keys()}
    120         target_regions_in_db = list(set(sum([target_to_db_dict[x]['Query'].tolist() for x in target_to_db_dict.keys()],[])))
    121     elif type(region_sets) == pr.PyRanges:

File ~/anaconda2/envs/py38/lib/python3.8/site-packages/pycistarget/utils.py:267, in target_to_query(target, query, fraction_overlap)
    264 if isinstance(query, pr.PyRanges):
    265     query_pr=query
--> 267 join_pr = target_pr.join(query_pr, report_overlap = True)
    268 join_pr.Overlap_query =  join_pr.Overlap/(join_pr.End_b - join_pr.Start_b)
    269 join_pr.Overlap_target =  join_pr.Overlap/(join_pr.End - join_pr.Start)

UnboundLocalError: local variable 'target_pr' referenced before assignment

The function works when I run one dictionary element at a time, e.g., region_sets=region_sets['DARs']. I also made sure that all values in the dictionary are populated (there is a Pyranges set for every topic and DAR).

My understanding of the utils.py script is that my region_sets are not being read in properly. target_pr should be given the value of target as it's a set of pr.Pyranges.

Thanks very much for your help with this, it's a fantastic resource!

[SeppeDeWinter]

Hi @SidG13

This function from pycistarget is indeed not designed to work with a dictionary of dictionaries. As you suggested yourself you should run the function for one dictionary element at a time (this is also what is done internally in the wrapper function in scenicplus).

However as of commit 22735ca with code from PR #56 you should be able to use the wrapper function now with your non-model organism.

Best,

Seppe