diff --git a/jcvi/graphics/chromosome.py b/jcvi/graphics/chromosome.py index 5f590283..ef9ce000 100644 --- a/jcvi/graphics/chromosome.py +++ b/jcvi/graphics/chromosome.py @@ -6,9 +6,10 @@ from the script used in the Tang et al. PNAS 2010 paper, sigma figure. """ import sys + from itertools import groupby from math import ceil -from typing import Tuple +from typing import Optional, Tuple import numpy as np @@ -24,6 +25,7 @@ Rectangle, latex, markup, + normalize_axes, plt, savefig, set1_n, @@ -480,7 +482,7 @@ class will get assigned a unique color. `id_mappings` file is optional (if mappingfile = args[1] fig = plt.figure(1, (iopts.w, iopts.h)) - root = fig.add_axes([0, 0, 1, 1]) + root = fig.add_axes((0, 0, 1, 1)) draw_chromosomes( root, @@ -497,9 +499,7 @@ class will get assigned a unique color. `id_mappings` file is optional (if title=opts.title, ) - root.set_xlim(0, 1) - root.set_ylim(0, 1) - root.set_axis_off() + normalize_axes(root) prefix = bedfile.rsplit(".", 1)[0] figname = prefix + "." + opts.format @@ -511,14 +511,14 @@ def draw_chromosomes( bedfile, sizes, iopts, - mergedist, - winsize, - imagemap, - mappingfile=None, - gauge=False, - legend=True, - empty=False, - title=None, + mergedist: int, + winsize: int, + imagemap: bool = False, + mappingfile: Optional[str] = None, + gauge: bool = False, + legend: bool = True, + empty: bool = False, + title: Optional[str] = None, ): bed = Bed(bedfile) prefix = bedfile.rsplit(".", 1)[0] diff --git a/jcvi/projects/jcvi.py b/jcvi/projects/jcvi.py index 2d7623a7..dc38dd02 100644 --- a/jcvi/projects/jcvi.py +++ b/jcvi/projects/jcvi.py @@ -15,81 +15,10 @@ from ..assembly.kmer import draw_ks_histogram from ..compara.pedigree import Pedigree, calculate_inbreeding from ..graphics.base import normalize_axes, panel_labels, plt, savefig +from ..graphics.chromosome import draw_chromosomes from ..graphics.landscape import draw_multi_depth -def genomebuild(args): - """ - %prog genomebuild reads.histo geneticmap.matrix hic.resolution_500000.npy hic.resolution_500000.json - - Plot genome build composite figure, including: - A. Read kmer histogram - B. Genetic map concordance - C. Hi-C contact map concordance - """ - p = OptionParser(genomebuild.__doc__) - _, args, iopts = p.set_image_options(args, figsize="21x7") - - if len(args) != 4: - sys.exit(not p.print_help()) - - reads_histo, mstmap, hic_matrix, hic_json = args - - fig = plt.figure(1, (iopts.w, iopts.h)) - root = fig.add_axes((0, 0, 1, 1)) - ax1_root = fig.add_axes((0, 0, 1 / 3, 1)) - ax2_root = fig.add_axes((1 / 3, 0, 1 / 3, 1)) - ax3_root = fig.add_axes((2 / 3, 0, 1 / 3, 1)) - ax1 = fig.add_axes((1 / 3 * 0.1, 0.1, 1 / 3 * 0.8, 0.8)) - ax2 = fig.add_axes((1 / 3 * 1.1, 0.1, 1 / 3 * 0.8, 0.8)) - ax3 = fig.add_axes((1 / 3 * 2.1, 0.1, 1 / 3 * 0.8, 0.8)) - - # Panel A - logger.info("Plotting read kmer histogram") - _ = draw_ks_histogram( - ax1, - reads_histo, - method="nbinom", - coverage=0, - vmin=2, - vmax=200, - species="*S. species* ‘Variety 1’", - K=21, - maxiter=100, - peaks=False, - ) - - # Panel B - logger.info("Plotting genetic map concordance") - draw_geneticmap_heatmap(ax2_root, ax2, mstmap, 1000) - - # Panel C - logger.info("Plotting Hi-C contact map concordance") - draw_hic_heatmap( - ax3_root, - ax3, - hic_matrix, - hic_json, - contig=None, - groups_file="groups", - title="*S. species* Hi-C contact map", - vmin=1, - vmax=6, - plot_breaks=True, - ) - - labels = ( - (1 / 3 * 0.1, 0.95, "A"), - (1 / 3 * 1.1, 0.95, "B"), - (1 / 3 * 2.1, 0.95, "C"), - ) - panel_labels(root, labels) - normalize_axes([root, ax1_root, ax2_root, ax3_root]) - - image_name = "genomebuild.pdf" - savefig(image_name, dpi=iopts.dpi, iopts=iopts) - - def diversity(args): """ %prog diversity pedigree.ped VAR?_srtd.wgs.regions.bed.gz @@ -124,12 +53,8 @@ def diversity(args): A.draw(pngfile, prog="dot", args=f"-Gdpi={dpi}") logger.info("Pedigree graph written to `%s`", pngfile) - M = plt.imread(pngfile) - width, height = M.shape[1] / dpi, M.shape[0] / dpi - logger.info("Image size: %.2f x %.2f (dpi=%d)", width, height, dpi) - - # Force aspect ratio to be equal - ax1_root.imshow(M) + # Show the image as is + ax1_root.imshow(plt.imread(pngfile)) ax1_root.set_axis_off() # Panel B @@ -169,11 +94,152 @@ def diversity(args): savefig(image_name, dpi=iopts.dpi, iopts=iopts) +def landscape(args): + """ + %prog landscape features.bed athaliana.sizes Fig5.png + + Plot landscape composite figure, including: + A. Example genomic features painted on Arabidopsis genome + B. Landscape of genomic features across the genome + """ + p = OptionParser(landscape.__doc__) + _, args, iopts = p.set_image_options(args, figsize="10x7") + + if len(args) != 3: + sys.exit(not p.print_help()) + + bedfile, sizesfile, pngfile = args + + fig = plt.figure(1, (iopts.w, iopts.h)) + root = fig.add_axes((0, 0, 1, 1)) + aspect_ratio = iopts.w / iopts.h + ax1_root = fig.add_axes((0, 1 / 4, 0.4, 0.5 * aspect_ratio)) + ax2_root = fig.add_axes((0.4, 0.43, 0.54, 0.57)) + ax3_root = fig.add_axes((0.43, -0.13, 0.54, 0.57)) + + # Panel A + logger.info("Plotting example genomic features painted on Arabidopsis genome") + draw_chromosomes( + ax1_root, + bedfile, + sizesfile, + iopts=iopts, + mergedist=0, + winsize=50000, + gauge=True, + legend=True, + empty=False, + title="*Arabidopsis* genome features", + ) + + # Panel B + logger.info("Plotting landscape of genomic features across the genome") + M = plt.imread(pngfile) + width, height = M.shape[1], M.shape[0] + # Split the image into left and right parts + mid = width // 2 - 10 + left_cut = right_cut = 900 + logger.info("Image size: %dx%d", width, height) + logger.info("Splitting image at %d", mid) + + LM, RM = M[:left_cut, :mid], M[:right_cut, mid:] + logger.info("Left image size: %dx%d", LM.shape[1], LM.shape[0]) + logger.info("Right image size: %dx%d", RM.shape[1], RM.shape[0]) + + ax2_root.imshow(LM, extent=(0.4, 1, 0.5, 1), aspect="auto") + ax3_root.imshow(RM, extent=(0.4, 1, 0, 0.5), aspect="auto") + ax2_root.set_axis_off() + ax3_root.set_axis_off() + + labels = ( + (0.02, 0.95, "A"), + (0.42, 0.95, "B"), + ) + panel_labels(root, labels) + normalize_axes([root, ax1_root]) + + image_name = "landscape.pdf" + savefig(image_name, dpi=iopts.dpi, iopts=iopts) + + +def genomebuild(args): + """ + %prog genomebuild reads.histo geneticmap.matrix hic.resolution_500000.npy hic.resolution_500000.json + + Plot genome build composite figure, including: + A. Read kmer histogram + B. Genetic map concordance + C. Hi-C contact map concordance + """ + p = OptionParser(genomebuild.__doc__) + _, args, iopts = p.set_image_options(args, figsize="21x7") + + if len(args) != 4: + sys.exit(not p.print_help()) + + reads_histo, mstmap, hic_matrix, hic_json = args + + fig = plt.figure(1, (iopts.w, iopts.h)) + root = fig.add_axes((0, 0, 1, 1)) + ax1_root = fig.add_axes((0, 0, 1 / 3, 1)) + ax2_root = fig.add_axes((1 / 3, 0, 1 / 3, 1)) + ax3_root = fig.add_axes((2 / 3, 0, 1 / 3, 1)) + ax1 = fig.add_axes((1 / 3 * 0.1, 0.1, 1 / 3 * 0.8, 0.8)) + ax2 = fig.add_axes((1 / 3 * 1.1, 0.1, 1 / 3 * 0.8, 0.8)) + ax3 = fig.add_axes((1 / 3 * 2.1, 0.1, 1 / 3 * 0.8, 0.8)) + + # Panel A + logger.info("Plotting read kmer histogram") + _ = draw_ks_histogram( + ax1, + reads_histo, + method="nbinom", + coverage=0, + vmin=2, + vmax=200, + species="*S. species* ‘Variety 1’", + K=21, + maxiter=100, + peaks=False, + ) + + # Panel B + logger.info("Plotting genetic map concordance") + draw_geneticmap_heatmap(ax2_root, ax2, mstmap, 1000) + + # Panel C + logger.info("Plotting Hi-C contact map concordance") + draw_hic_heatmap( + ax3_root, + ax3, + hic_matrix, + hic_json, + contig=None, + groups_file="groups", + title="*S. species* Hi-C contact map", + vmin=1, + vmax=6, + plot_breaks=True, + ) + + labels = ( + (1 / 3 * 0.1, 0.95, "A"), + (1 / 3 * 1.1, 0.95, "B"), + (1 / 3 * 2.1, 0.95, "C"), + ) + panel_labels(root, labels) + normalize_axes([root, ax1_root, ax2_root, ax3_root]) + + image_name = "genomebuild.pdf" + savefig(image_name, dpi=iopts.dpi, iopts=iopts) + + def main(): actions = ( - ("genomebuild", "Plot genome build composite figure"), ("diversity", "Plot diversity composite figure"), + ("genomebuild", "Plot genome build composite figure"), + ("landscape", "Plot landscape composite figure"), ) p = ActionDispatcher(actions) p.dispatch(globals())