Plotting TADs and loops

[Gonysh]

Hi,

I wanted to make sure there is no function to plot/annotate loops and contact regions over interaction matrices using fancplot, right? For example, using a .bedpe file to mark corner dots or a .bed file to plot a square marking a TAD.

Comparing hic matrices and RNA-seq works beautifully :), just missing the 2D feature annotation :(.

Thanks, Gony.

[kaukrise]

Hi, there is currently nothing in the fancplot command line tool, unfortunately. If you can use the API, have a look here for an "arc" style plotting base and also at HicPeakPlot - this currently uses FAN-C peaks files for plotting and is not properly documented, but it might point you in the right direction. An instance of HicPeakPlot can be added to a matrix plot using the add_overlay function.

When I have more time I'll consider adding something more versatile.

[irenemota]

Hi,

I am trying to use a HicPeakPlot in order to plot the loops on top of the matrices. The loop files that I have, come from HICCUPS directly, so they are in a bedpe format and when I try to do the add_overlap function in the hic object with my loop file I get the following error message:

AttributeError: 'Bedpe' object has no attribute 'compatibility'

The way that I am doing it is as follows:

sert_10k = fanc.load("../Low_coverage/Sox9_rep1_2_10k.fanc")
loops_sert_vs_sf1xy= fanc.load("../Loops/sf1xy_vs_sert/differential_loops1.bedpe")
hic= fancplot.HicPlot(sert_10k)
hic.add_overlay(loops_sert_vs_sf1xy)

Is there a way that I can import this file as a Peakinfo file from fanc or make it compatible?

Thank you very much in advance,

Irene

[kaukrise]

Hey, the overlay functionality is not well developed and therefore not documented. I have a solution for you below, but please note that it is not properly tested, cumbersome to use, and in no way optimised for speed or resources.

import fanc
import fanc.plotting as fancplot
from fanc.plotting.base_plotter import BaseOverlayPlotter
import genomic_regions as gr

from matplotlib import patches
from future.utils import string_types

class RectangleOverlay(BaseOverlayPlotter):
        def __init__(self, region_pairs, size=None, width=None, height=None,
                 **kwargs):
        BaseOverlayPlotter.__init__(self)
        self._region_pairs = region_pairs
        if size is not None:
            self._width, self._height = size, size
        else:
            self._width = width
            self._height = height

        kwargs.setdefault('linewidth', 2)
        kwargs.setdefault('edgecolor', 'r')
        kwargs.setdefault('facecolor', 'none')
        self._rect_args = kwargs
        self._patches = []

    def _plot(self, base_plot, region):
        for r1, r2 in self._region_pairs:
            r1 = gr.as_region(r1)
            r2 = gr.as_region(r2)

            rect = patches.Rectangle(xy=(r1.start, r2.start),
                                     width=self._width if self._width is not None else len(r1),
                                     height=self._height if self._height is not None else len(r2),
                                     **self._rect_args)
            base_plot.ax.add_patch(rect)
            self._patches.append(rect)

    @property
    def compatibility(self):
        return ['SquareMatrixPlot']

sert_10k = fanc.load("../Low_coverage/Sox9_rep1_2_10k.fanc")
loops_sert_vs_sf1xy = fanc.load("../Loops/sf1xy_vs_sert/differential_loops1.bedpe")
region_pairs = list(loops_sert_vs_sf1xy.region_pairs())

loop_overlay = RectangleOverlay(region_pairs)
p = fancplot.SquareMatrixPlot(sert_10k)
p.add_overlay(loop_overlay)
p.plot('chr12:21mb-22mb')  # change to something with a loop, of course)

[irenemota]

Thank you very much for the quick reply.

I have tried running it , however the loops seem to not match with the hic (maybe they are shifted?). I also loaded the same files in juicebox to make sure that the problem was not in the loop calling and there it looks fine:

Do you know where it could be the problem?

Thank you very much again,

Irene

[kaukrise]

Ah, I think I see the issue. Can you replace the following:

    def _plot(self, base_plot, region):
        for r1, r2 in self._region_pairs:
            r1 = gr.as_region(r1)
            r2 = gr.as_region(r2)

            rect = patches.Rectangle(xy=(r1.start, r2.start),
                                     width=self._width if self._width is not None else len(r1),
                                     height=self._height if self._height is not None else len(r2),
                                     **self._rect_args)
            base_plot.ax.add_patch(rect)
            self._patches.append(rect)

with this?

    def _plot(self, base_plot, region):
        for r1, r2 in self._region_pairs:
            r1 = gr.as_region(r1)
            r2 = gr.as_region(r2)

            if r1.chromosome != region[0].chromosome or r2.chromosome != region[1].chromosome:
                continue

            rect = patches.Rectangle(xy=(r1.start, r2.start),
                                     width=self._width if self._width is not None else len(r1),
                                     height=self._height if self._height is not None else len(r2),
                                     **self._rect_args)
            base_plot.ax.add_patch(rect)
            self._patches.append(rect)

[irenemota]

Thank you very much! This solved the problem

Best,

Irene