coverageplot does not show gene body, exon, etc

[mervedede]

Hello,

I have been trying to use the coverage_plot function to generate integrated multiome plots for my scplus object that I generated for my dataset. For some of my regulon hits, I was able to create the coverage plots properly but for some regulons, I have noticed that the gene body/exon section at the bottom of the coverageplot appears empty.

I have checked the gtf file and pr_interact and could not figure out why it appears empty for some regulons that have passed the threshold. I am pasting the images for a case that works and a case that fails.

works:

fails:

[SeppeDeWinter]

Hi @mervedede

Thanks for using SCENIC+.
I transferred your question to the discussion section.

For the ones where it appears empty, can you check whether there are any genes located at those specific chromosomal positions in the GTF file you are using?

Best,

Seppe

[mervedede]

Hi Seppe,

Thanks for developing this tool and your answer. Yes, I have checked that the gtf file contains genes in this region. For example, for my specific example above, the region I was trying to plot has the following genes (and my gene of interest NFKB1):

[SeppeDeWinter]

Hi @mervedede

hmm... not sure what is going on. The code for this plotting is a bit buggy though ... Would you mind sharing the GTF file + region of interest, I'll do some troubleshooting.

Best,

Seppe

[mervedede]

Hi Seppe,

Thank you for your help. The gtf file I am using can be downloaded from here: "https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.annotation.gtf.gz". The example region that does not work for me above is "chr4:102430447-102736464".

Best,
Merve

[SeppeDeWinter]

This commit should fix your issue: e4bdd9f

Best,

Seppe

[jklvrt]

Hi @SeppeDeWinter,

Having this same problem now. For some ranges (which I'm sure contain genes) no genes get shown at all... for other ranges only 1 gene shown when there are several more (protein-coding) genes in this window... tried the fix (deleting/#-ing those 3 lines in coverageplot.py) but nothing changes. As always, thanks in advance :)

Also can't get any of the gene violin plots to show, error looks something like:::

---> 17 scenicplus.plotting.coverageplot.coverage_plot(scplus_obj, bw_dict = bw_dict, region = 'chr4:89600344-89848315',  meta_data_key = 'cell_type2', pr_consensus_bed = consensus_peaks, pr_gtf = gene_info, genes_violin_plot = 'SNCA')

File [~/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/scenicplus/plotting/coverageplot.py:391](http://localhost:8888/lab/tree/DATA2/home/jonathan/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/scenicplus/plotting/coverageplot.py#line=390), in coverage_plot(SCENICPLUS_obj, bw_dict, region, genes_violin_plot, genes_arcs, gene_height, exon_height, meta_data_key, pr_consensus_bed, region_bed_height, pr_gtf, pr_interact, bw_ymax, color_dict, cmap, plot_order, figsize, fontsize_dict, gene_label_offset, arc_rad, arc_lw, cmap_violinplots, violinplots_means_color, violinplots_edge_color, violoinplots_alpha, width_ratios_dict, height_ratios_dict, sort_vln_plots, add_custom_ax)
    388 # get expression values for gene[/genes](http://localhost:8888/genes)
    389 idx_gene = list(SCENICPLUS_obj.gene_names).index(genes_violin_plot) if type(
    390     genes_violin_plot) == str else [list(SCENICPLUS_obj.gene_names).index(g) for g in genes_violin_plot]
--> 391 expr_vals = mtx[:, idx_gene]
    393 norm = matplotlib.colors.Normalize(vmin=0, vmax=len(plot_order))
    394 mapper = cm.ScalarMappable(norm=norm, cmap=cmap_violinplots)

File [~/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/pandas/core/frame.py:3805](http://localhost:8888/lab/tree/DATA2/home/jonathan/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/pandas/core/frame.py#line=3804), in DataFrame.__getitem__(self, key)
   3803 if self.columns.nlevels > 1:
   3804     return self._getitem_multilevel(key)
-> 3805 indexer = self.columns.get_loc(key)
   3806 if is_integer(indexer):
   3807     indexer = [indexer]

File [~/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/pandas/core/indexes/base.py:3807](http://localhost:8888/lab/tree/DATA2/home/jonathan/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/pandas/core/indexes/base.py#line=3806), in Index.get_loc(self, key, method, tolerance)
   3802         raise KeyError(key) from err
   3803     except TypeError:
   3804         # If we have a listlike key, _check_indexing_error will raise
   3805         #  InvalidIndexError. Otherwise we fall through and re-raise
   3806         #  the TypeError.
-> 3807         self._check_indexing_error(key)
   3808         raise
   3810 # GH#42269

File [~/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/pandas/core/indexes/base.py:5963](http://localhost:8888/lab/tree/DATA2/home/jonathan/DATA2/home/jonathan/miniconda3/envs/scenicplus/lib/python3.11/site-packages/pandas/core/indexes/base.py#line=5962), in Index._check_indexing_error(self, key)
   5959 def _check_indexing_error(self, key):
   5960     if not is_scalar(key):
   5961         # if key is not a scalar, directly raise an error (the code below
   5962         # would convert to numpy arrays and raise later any way) - GH29926
-> 5963         raise InvalidIndexError(key)

InvalidIndexError: (slice(None, None, None), 7500)

[yojetsharma]

Am I doing this right?Given such a sparse documentation I am finding it hard to understand each parameter that is being used to generate a coverage plot.

from scenicplus.plotting.coverageplot import coverage_plot
# Example data
bw_dict = {'Glioblast': '~/consensus_peak_calling/pseudobulk_bw_files/Glioblast.bw'}
region = 'chr19:8423850-8424350'  # Example genomic region

# Optional gene list for violin plots and arcs
#genes_violin_plot = ['DLK1', 'YY1']
#genes_arcs = ['DLK1']

# Optional metadata, pyranges, etc. - adjust as needed for your data
meta_data_key = 'broader_annotation'  # If available in SCENICPLUS_obj.metadata_cell
pr_consensus_bed = "~/consensus_peak_calling/consensus_regions.bed"  # Provide PyRanges object if available
pr_gtf = pr_annotation  # Provide PyRanges object from pycistopic tutorial?
pr_interact = None  # not sure where to provide this from 

# Plot
coverage_plot(
    SCENICPLUS_obj=scplus_obj,     # Your SCENICPLUS object
    bw_dict=bw_dict,               # BigWig files for glioblast
    region=region,                 # Genomic region of interest
    genes_violin_plot=None,   # Genes to include in violin plots
    genes_arcs=None,         # Genes to plot region-to-gene arcs
    meta_data_key='broader_annotation',   # Metadata key for annotations
    pr_consensus_bed=None,     # PyRanges object for peaks, if available
    pr_gtf=None,                 # Gene structure PyRanges object, if available
    pr_interact=None,       # Region-to-gene interaction PyRanges object
    figsize=(6, 8),                # Figure size
    cmap='tab20'                   # Colormap for coverage plots
)