Motif enrichment analysis using pycistarget

[mtrebelo]

Hi!

I am using the 'Tutorial: Mix of melanoma cell lines' with a zebrafish dataset.

Until now, I have been able to run everything well until 'scATAC-seq preprocessing using pycisTopic', where I got the region_bin_topics_otsu, region_bin_topics_top3k and markers_dict.

Then, I ran this:

import pyranges as pr
from pycistarget.utils import region_names_to_coordinates
region_sets = {}
region_sets['topics_otsu'] = {}
region_sets['topics_top_3'] = {}
region_sets['DARs'] = {}

for topic in region_bin_topics_otsu.keys():
    regions = region_bin_topics_otsu[topic].index[region_bin_topics_otsu[topic].index.str.startswith('chr')] #only keep regions on known chromosomes
    if len(regions) > 0:
        region_sets['topics_otsu'][topic] = pr.PyRanges(region_names_to_coordinates(regions))

for topic in region_bin_topics_top3k.keys():
    regions = region_bin_topics_top3k[topic].index[region_bin_topics_top3k[topic].index.str.startswith('chr')] #only keep regions on known chromosomes
    if len(regions) > 0:
        region_sets['topics_top_3'][topic] = pr.PyRanges(region_names_to_coordinates(regions))
        
        
for DAR in markers_dict.keys():
    regions = markers_dict[DAR].index[markers_dict[DAR].index.str.startswith('chr')] #only keep regions on known chromosomes
    if len(regions) > 0:
        region_sets['DARs'][DAR] = pr.PyRanges(region_names_to_coordinates(regions))

I created the .feather files with the following command /create_cistarget_motif_databases.py -f ../peaks_48hpf1.fasta -M ../motifs_cb_format/ -m ../motifs.lst -o peaks_48hpf1 -t 36 and took the motifs from Jaspar and transformed them into the cluster-buster format.

For the .tbl file, I took the motifs-v10-nr.hgnc-m0.00001-o0.0.tbl file from https://resources.aertslab.org/cistarget/motif_collections/v10nr_clust_public/snapshots/motifs-v10-nr.hgnc-m0.00001-o0.0.tbl and with a database for ortholog genes, changed the gene_name from human to zebrafish. Here is a screeshot of the .tbl file:

I created the costume_annot database with my .gtf file and here is a screenshot:

After running this code:

from scenicplus.wrappers.run_pycistarget import run_pycistarget
run_pycistarget(
    region_sets = region_sets,
    species = 'custom',
    custom_annot=custom_annotation,
    save_path = os.path.join(work_dir, 'motifs'),
    ctx_db_path = rankings_db,
    dem_db_path = scores_db,
    path_to_motif_annotations = motif_annotation,
    run_without_promoters = True,
    n_cpu = 8,
    _temp_dir = os.path.join(tmp_dir, 'ray_spill'),
    annotation_version = 'v10',
    )

this is the error:

2023-06-02 13:03:33,907 pycisTarget_wrapper INFO     /work/project/ladcol_004/scSLAM/Teresa_workspace/scenicplus/2_motif_enrichment/motifs folder already exists.
2023-06-02 13:03:34,121 pycisTarget_wrapper INFO     Loading cisTarget database for topics_otsu
2023-06-02 13:03:34,123 cisTarget    INFO     Reading cisTarget database
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[18], line 2
      1 from scenicplus.wrappers.run_pycistarget import run_pycistarget
----> 2 run_pycistarget(
      3     region_sets = region_sets,
      4     species = 'custom',
      5     custom_annot=custom_annotation,
      6     save_path = os.path.join(work_dir, 'scSLAM/Teresa_workspace/scenicplus/2_motif_enrichment/motifs'),
      7     ctx_db_path = rankings_db,
      8     dem_db_path = scores_db,
      9     path_to_motif_annotations = motif_annotation,
     10     run_without_promoters = True,
     11     n_cpu = 8,
     12     _temp_dir = os.path.join(tmp_dir, 'ray_spill'),
     13     annotation_version = 'v10',
     14     )

File /work/project/ladcol_015/scenicplus/src/scenicplus/wrappers/run_pycistarget.py:182, in run_pycistarget(region_sets, species, save_path, custom_annot, save_partial, ctx_db_path, dem_db_path, run_without_promoters, biomart_host, promoter_space, ctx_auc_threshold, ctx_nes_threshold, ctx_rank_threshold, dem_log2fc_thr, dem_motif_hit_thr, dem_max_bg_regions, annotation, motif_similarity_fdr, path_to_motif_annotations, annotation_version, n_cpu, _temp_dir, exclude_motifs, exclude_collection, **kwargs)
    180 ## CISTARGET
    181 regions = region_sets[key]
--> 182 ctx_db = cisTargetDatabase(ctx_db_path, regions)  
    183 if exclude_motifs is not None:
    184     out = pd.read_csv(exclude_motifs, header=None).iloc[:,0].tolist()

File ~/miniconda3/envs/scenic+/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 ~/miniconda3/envs/scenic+/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 ~/miniconda3/envs/scenic+/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 ~/miniconda3/envs/scenic+/lib/python3.8/site-packages/pycistarget/utils.py:283, in target_to_query(target, query, fraction_overlap)
    280     query_pr=query
    282 join_pr = target_pr.join(query_pr, report_overlap = True)
--> 283 join_pr.Overlap_query =  join_pr.Overlap/(join_pr.End_b - join_pr.Start_b)
    284 join_pr.Overlap_target =  join_pr.Overlap/(join_pr.End - join_pr.Start)
    285 join_pr = join_pr[(join_pr.Overlap_query > fraction_overlap) | (join_pr.Overlap_target > fraction_overlap)]

File ~/miniconda3/envs/scenic+/lib/python3.8/site-packages/pyranges/pyranges.py:269, in PyRanges.__getattr__(self, name)
    244 """Return column.
    245 
    246 Parameters
   (...)
    264 Name: Start, dtype: int32
    265 """
    267 from pyranges.methods.attr import _getattr
--> 269 return _getattr(self, name)

File ~/miniconda3/envs/scenic+/lib/python3.8/site-packages/pyranges/methods/attr.py:67, in _getattr(self, name)
     65     return pd.concat([df[name] for df in self.values()])
     66 else:
---> 67     raise AttributeError("PyRanges object has no attribute", name)

AttributeError: ('PyRanges object has no attribute', 'Overlap')

At this point I am not sure how to bypass it and I would be very greatful if you could help me out.

Best, Maria

[YaoLi3]

for one thing maybe the strands should be formatted as 1 and -1 instead of "+", "-" in your costume_annot database