Edinburgh-Genome-Foundry / DnaFeaturesViewer

:eye: Python library to plot DNA sequence features (e.g. from Genbank files)
https://edinburgh-genome-foundry.github.io/DnaFeaturesViewer/
MIT License
584 stars 90 forks source link

Zoom in capability? #30

Closed mortonjt closed 4 years ago

mortonjt commented 4 years ago

This is an amazing library - I love how tightly it is integrated with matplotlib.

I've got a genome that I'm trying to perform more in depth analysis (example code is below using the ncov2019 genome)

from dna_features_viewer import BiopythonTranslator
from Bio import SeqIO
import matplotlib.pyplot as plt

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(
    5, 1, figsize=(15, 10), sharex=True, gridspec_kw={"height_ratios": [5, 1, 1, 1, 1]}
)
record = SeqIO.read("sequence.gb", "genbank")
graphic_record = BiopythonTranslator().translate_record(record)
graphic_record.plot(ax=ax1, with_ruler=False, strand_in_label_threshold=4)

gc = lambda s: 100.0 * len([c for c in s if c in "GC"]) / 50
xx = np.arange(len(record.seq) - 50)
yy = [gc(record.seq[x : x + 50]) for x in xx]
ax2.fill_between(xx + 25, yy, alpha=0.3)
ax2.set_ylim(bottom=0)
ax2.set_ylabel("GC(%)")

ax3.plot(filtY[:, 0])
ax3.set_ylabel('(R/Y)')
ax4.plot(filtY[:, 1])
ax4.set_ylabel('(A/G)')
ax5.plot(filtY[:, 2])
ax5.set_ylabel('(C/T)')

So far, I'm able to generate a beautiful genome plot

Screenshot 2020-02-14 12 04 18

But when I try to zoom in with

for ax in [ax1, ax2, ax3, ax4, ax5]:
    ax.set_xlim([24000, 25000])

The plot seems to implode

Screenshot 2020-02-14 12 04 26

Not sure what is going on here (may require a closer look at the internals) - but these sorts of features would be incredibly useful to have ready out-of-the-box for genome wide analysis.

EDIT: the problem seems to be arising due to Jupyter. When I save these plots to a file, it seems to be ok

image

Although, this is still problematic - if I zoom in to a 50bp region, I get

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/miniconda3/envs/gert/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

~/miniconda3/envs/gert/lib/python3.7/site-packages/IPython/core/pylabtools.py in <lambda>(fig)
    242 
    243     if 'png' in formats:
--> 244         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    245     if 'retina' in formats or 'png2x' in formats:
    246         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

~/miniconda3/envs/gert/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, **kwargs)
    126 
    127     bytes_io = BytesIO()
--> 128     fig.canvas.print_figure(bytes_io, **kw)
    129     data = bytes_io.getvalue()
    130     if fmt == 'svg':

~/miniconda3/envs/gert/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, **kwargs)
   2087                     orientation=orientation,
   2088                     bbox_inches_restore=_bbox_inches_restore,
-> 2089                     **kwargs)
   2090             finally:
   2091                 if bbox_inches and restore_bbox:

~/miniconda3/envs/gert/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in print_png(self, filename_or_obj, metadata, pil_kwargs, *args, **kwargs)
    525 
    526         else:
--> 527             FigureCanvasAgg.draw(self)
    528             renderer = self.get_renderer()
    529             with cbook._setattr_cm(renderer, dpi=self.figure.dpi), \

~/miniconda3/envs/gert/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in draw(self)
    384         Draw the figure using the renderer.
    385         """
--> 386         self.renderer = self.get_renderer(cleared=True)
    387         with RendererAgg.lock:
    388             self.figure.draw(self.renderer)

~/miniconda3/envs/gert/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_renderer(self, cleared)
    397                           and getattr(self, "_lastKey", None) == key)
    398         if not reuse_renderer:
--> 399             self.renderer = RendererAgg(w, h, self.figure.dpi)
    400             self._lastKey = key
    401         elif cleared:

~/miniconda3/envs/gert/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in __init__(self, width, height, dpi)
     84         self.width = width
     85         self.height = height
---> 86         self._renderer = _RendererAgg(int(width), int(height), dpi)
     87         self._filter_renderers = []
     88 

ValueError: Image size of 248267x575 pixels is too large. It must be less than 2^16 in each direction.

<Figure size 1080x720 with 4 Axes>
Zulko commented 4 years ago

To zoom in, use .crop(). See the examples in the Readme

On Fri, 14 Feb 2020, 17:08 Jamie Morton, notifications@github.com wrote:

This is an amazing library - I love how tightly it is integrated with matplotlib.

I've got a genome that I'm trying to perform more in depth analysis (example code is below using the ncov2019 genome https://www.ncbi.nlm.nih.gov/nuccore/MN908947)

from dna_features_viewer import BiopythonTranslatorfrom Bio import SeqIOimport matplotlib.pyplot as plt

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots( 5, 1, figsize=(15, 10), sharex=True, gridspec_kw={"height_ratios": [5, 1, 1, 1, 1]} ) record = SeqIO.read("sequence.gb", "genbank") graphic_record = BiopythonTranslator().translate_record(record) graphic_record.plot(ax=ax1, with_ruler=False, strand_in_label_threshold=4)

gc = lambda s: 100.0 * len([c for c in s if c in "GC"]) / 50 xx = np.arange(len(record.seq) - 50) yy = [gc(record.seq[x : x + 50]) for x in xx] ax2.fill_between(xx + 25, yy, alpha=0.3) ax2.set_ylim(bottom=0) ax2.set_ylabel("GC(%)")

ax3.plot(filtY[:, 0]) ax3.set_ylabel('(R/Y)') ax4.plot(filtY[:, 1]) ax4.set_ylabel('(A/G)') ax5.plot(filtY[:, 2]) ax5.set_ylabel('(C/T)')

So far, I'm able to generate a beautiful genome plot [image: Screenshot 2020-02-14 12 04 18] https://user-images.githubusercontent.com/4184797/74551947-33841d00-4f22-11ea-866e-ab00e810d73f.png

But when I try to zoom in with

for ax in [ax1, ax2, ax3, ax4, ax5]: ax.set_xlim([24000, 25000])

The plot seems to implode

[image: Screenshot 2020-02-14 12 04 26] https://user-images.githubusercontent.com/4184797/74551941-2f57ff80-4f22-11ea-93bd-56f9fe3802ab.png

Not sure what is going on here (may require a closer look at the internals) - but these sorts of features would be incredibly useful to have ready out-of-the-box for genome wide analysis.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/Edinburgh-Genome-Foundry/DnaFeaturesViewer/issues/30?email_source=notifications&email_token=AAMR4KRI2X2IYNLZRTAVJGDRC3FYXA5CNFSM4KVMFGM2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4INUDNBQ, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAMR4KUYQMIZOE2DX7PFP23RC3FYXANCNFSM4KVMFGMQ .

mortonjt commented 4 years ago

Ok, that fixed it. Thanks!