1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
|
import sys
import click
import numpy as np
from .. import util
from ..api import Cooler
from . import cli
MAX_MATRIX_SIZE_FILE = int(1e8)
MAX_MATRIX_SIZE_INTERACTIVE = int(1e7)
def get_matrix_size(c, row_region, col_region):
nrows = c.extent(row_region)[1] - c.extent(row_region)[0]
ncols = c.extent(col_region)[1] - c.extent(col_region)[0]
return ncols * nrows
def load_matrix(c, row_region, col_region, field, balanced, scale):
mat = c.matrix(balance=balanced, field=field).fetch(row_region, col_region)
if scale == "log2":
mat = np.log2(mat)
elif scale == "log10":
mat = np.log10(mat)
return mat
def interactive(ax, c, row_chrom, col_chrom, field, balanced, scale): # pragma: no cover
import matplotlib.pyplot as plt
# The code is heavily insired by
# https://gist.github.com/mdboom/048aa35df685fe694330764894f0e40a
def get_extent(ax):
xstart, ystart, xdelta, ydelta = ax.viewLim.bounds
xend = xstart + xdelta
yend = ystart + ydelta
return xstart, xend, ystart, yend
def round_trim_extent(extent, binsize, row_chrom_len, col_chrom_len):
xstart = int(np.floor(extent[0] / binsize) * binsize)
xend = int(np.ceil(extent[1] / binsize) * binsize)
ystart = int(np.floor(extent[3] / binsize) * binsize)
yend = int(np.ceil(extent[2] / binsize) * binsize)
xstart = max(0, xstart)
ystart = max(0, ystart)
# For now, don't let users to request the last bin, b/c its end
# lies outside of the genome
xend = min(xend, int(np.floor(col_chrom_len / binsize) * binsize))
yend = min(yend, int(np.floor(row_chrom_len / binsize) * binsize))
return xstart, xend, yend, ystart
def update_heatmap(event):
ax.set_autoscale_on(False) # Otherwise, infinite loop
extent = get_extent(ax)
extent = round_trim_extent(extent, binsize, row_chrom_len, col_chrom_len)
if extent == plotstate["prev_extent"]:
return
plotstate["prev_extent"] = extent
new_col_region = col_chrom, int(extent[0]), int(extent[1])
new_row_region = row_chrom, int(extent[3]), int(extent[2])
im = ax.images[-1]
nelem = get_matrix_size(c, new_row_region, new_col_region)
if nelem >= MAX_MATRIX_SIZE_INTERACTIVE:
# requested area too large
im.set_data(np.ones(1)[:, None] * np.nan)
if not plotstate["placeholders"]:
box, = plt.plot(
[0, col_chrom_len, col_chrom_len, 0, 0, col_chrom_len],
[0, row_chrom_len, 0, 0, row_chrom_len, row_chrom_len],
c="k",
lw=0.5,
)
txt = plt.text(
0.5,
0.5,
"The requested region is too large\n"
"to display at this resolution.",
horizontalalignment="center",
verticalalignment="center",
transform=ax.transAxes,
)
plotstate["placeholders"] = [box, txt]
else:
# remove placeholders if any and update
while plotstate["placeholders"]:
plotstate["placeholders"].pop().remove()
im.set_data(
load_matrix(c, new_row_region, new_col_region, field, balanced, scale)
)
im.set_extent(extent)
ax.figure.canvas.draw_idle()
binsize = c.info["bin-size"]
row_chrom_len = c.chromsizes[row_chrom]
col_chrom_len = c.chromsizes[col_chrom]
plotstate = {"placeholders": [], "prev_extent": get_extent(plt.gca())}
plt.gcf().canvas.mpl_connect("button_release_event", update_heatmap)
plt.show()
@cli.command()
@click.argument(
"cool_uri",
metavar="COOL_PATH"
)
@click.argument(
"range",
type=str
)
@click.option(
"--range2", "-r2",
type=str,
help="The coordinates of a genomic region shown along the column dimension. "
"If omitted, the column range is the same as the row range. "
"Use to display asymmetric matrices or trans interactions.",
)
@click.option(
"--balanced", "-b",
is_flag=True,
default=False,
help="Show the balanced contact matrix. "
"If not provided, display the unbalanced counts.",
)
@click.option(
"--out", "-o",
help="Save the image of the contact matrix to a file. "
"If not specified, the matrix is displayed in an interactive window. "
"The figure format is deduced from the extension of the file, "
"the supported formats are png, jpg, svg, pdf, ps and eps.",
)
@click.option(
"--dpi",
type=int,
help="The DPI of the figure, if saving to a file"
)
@click.option(
"--scale", "-s",
type=click.Choice(["linear", "log2", "log10"]),
help="Scale transformation of the colormap: linear, log2 or log10. "
"Default is log10.",
default="log10",
)
@click.option(
"--force", "-f",
is_flag=True,
default=False,
help="Force display very large matrices (>=10^8 pixels). "
"Use at your own risk as it may cause performance issues.",
)
@click.option(
"--zmin",
type=float,
help="The minimal value of the color scale. Units must match those of the colormap scale. "
"To provide a negative value use a equal sign and quotes, e.g. -zmin='-0.5'",
)
@click.option(
"--zmax",
type=float,
help="The maximal value of the color scale. Units must match those of the colormap scale. "
"To provide a negative value use a equal sign and quotes, e.g. -zmax='-0.5'",
)
@click.option(
"--cmap",
default="YlOrRd",
help="The colormap used to display the contact matrix. "
"See the full list at http://matplotlib.org/examples/color/colormaps_reference.html",
)
@click.option(
"--field",
default="count",
show_default=True,
help="Pixel values to display."
)
def show(
cool_uri, range, range2, balanced, out, dpi, scale, force, zmin, zmax, cmap, field
):
"""
Display and browse a cooler in matplotlib.
COOL_PATH : Path to a COOL file or Cooler URI.
RANGE : The coordinates of the genomic region to display, in UCSC notation.
Example: chr1:10,000,000-11,000,000
"""
try:
import matplotlib as mpl
if out is not None:
mpl.use("Agg")
import matplotlib.pyplot as plt
except ImportError: # pragma: no cover
print("Install matplotlib to use cooler show", file=sys.stderr)
sys.exit(1)
c = Cooler(cool_uri)
chromsizes = c.chromsizes
row_region = range
col_region = row_region if range2 is None else range2
row_chrom, row_lo, row_hi = util.parse_region(row_region, chromsizes)
col_chrom, col_lo, col_hi = util.parse_region(col_region, chromsizes)
if (
get_matrix_size(c, row_region, col_region) >= MAX_MATRIX_SIZE_FILE
) and not force:
print(
"The matrix of the selected region is too large. "
"Try using lower resolution, selecting a smaller region, or use "
"the '--force' flag to override this safety limit.",
file=sys.stderr,
)
sys.exit(1)
plt.figure(figsize=(11, 10))
plt.get_current_fig_manager().set_window_title("Contact matrix")
plt.title("")
plt.imshow(
load_matrix(c, row_region, col_region, field, balanced, scale),
interpolation="none",
extent=[col_lo, col_hi, row_hi, row_lo],
vmin=zmin,
vmax=zmax,
cmap=cmap,
)
# If plotting into a file, plot and quit
plt.ylabel(f"{row_chrom} coordinate")
plt.xlabel(f"{col_chrom} coordinate")
cb = plt.colorbar()
cb.set_label(
{
"linear": "relative contact frequency",
"log2": "log 2 ( relative contact frequency )",
"log10": "log 10 ( relative contact frequency )",
}[scale]
)
if out:
plt.savefig(out, dpi=dpi)
else:
interactive(plt.gca(), c, row_chrom, col_chrom, field, balanced, scale)
|