def multihistogram(args):
    "\n    %prog multihistogram *.histogram species\n\n    Plot the histogram based on a set of K-mer hisotograms. The method is based\n    on Star et al.'s method (Atlantic Cod genome paper).\n    "
    p = OptionParser(multihistogram.__doc__)
    p.add_option('--kmin', default=15, type='int', help='Minimum K-mer size, inclusive')
    p.add_option('--kmax', default=30, type='int', help='Maximum K-mer size, inclusive')
    p.add_option('--vmin', default=2, type='int', help='Minimum value, inclusive')
    p.add_option('--vmax', default=100, type='int', help='Maximum value, inclusive')
    (opts, args, iopts) = p.set_image_options(args, figsize='10x5', dpi=300)
    if (len(args) < 1):
        sys.exit((not p.print_help()))
    histfiles = args[:(- 1)]
    species = args[(- 1)]
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])
    A = fig.add_axes([0.08, 0.12, 0.38, 0.76])
    B = fig.add_axes([0.58, 0.12, 0.38, 0.76])
    lines = []
    legends = []
    genomesizes = []
    for histfile in histfiles:
        ks = KmerSpectrum(histfile)
        (x, y) = ks.get_xy(opts.vmin, opts.vmax)
        K = get_number(op.basename(histfile).split('.')[0].split('-')[(- 1)])
        if (not (opts.kmin <= K <= opts.kmax)):
            continue
        (line,) = A.plot(x, y, '-', lw=1)
        lines.append(line)
        legends.append('K = {0}'.format(K))
        ks.analyze(K=K)
        genomesizes.append((K, (ks.genomesize / 1000000.0)))
    leg = A.legend(lines, legends, shadow=True, fancybox=True)
    leg.get_frame().set_alpha(0.5)
    title = '{0} genome K-mer histogram'.format(species)
    A.set_title(markup(title))
    (xlabel, ylabel) = ('Coverage (X)', 'Counts')
    A.set_xlabel(xlabel)
    A.set_ylabel(ylabel)
    set_human_axis(A)
    title = '{0} genome size estimate'.format(species)
    B.set_title(markup(title))
    (x, y) = zip(*genomesizes)
    B.plot(x, y, 'ko', mfc='w')
    t = np.linspace((opts.kmin - 0.5), (opts.kmax + 0.5), 100)
    p = np.poly1d(np.polyfit(x, y, 2))
    B.plot(t, p(t), 'r:')
    (xlabel, ylabel) = ('K-mer size', 'Estimated genome size (Mb)')
    B.set_xlabel(xlabel)
    B.set_ylabel(ylabel)
    set_ticklabels_arial(B)
    labels = ((0.04, 0.96, 'A'), (0.54, 0.96, 'B'))
    panel_labels(root, labels)
    normalize_axes(root)
    imagename = (species + '.multiK.pdf')
    savefig(imagename, dpi=iopts.dpi, iopts=iopts)
%prog multihistogram *.histogram species

Plot the histogram based on a set of K-mer hisotograms. The method is based
on Star et al.'s method (Atlantic Cod genome paper).
def multihistogram(args):
    "\n    %prog multihistogram *.histogram species\n\n    Plot the histogram based on a set of K-mer hisotograms. The method is based\n    on Star et al.'s method (Atlantic Cod genome paper).\n    "
    p = OptionParser(multihistogram.__doc__)
    p.add_option('--kmin', default=15, type='int', help='Minimum K-mer size, inclusive')
    p.add_option('--kmax', default=30, type='int', help='Maximum K-mer size, inclusive')
    p.add_option('--vmin', default=2, type='int', help='Minimum value, inclusive')
    p.add_option('--vmax', default=100, type='int', help='Maximum value, inclusive')
    (opts, args, iopts) = p.set_image_options(args, figsize='10x5', dpi=300)
    if (len(args) < 1):
        sys.exit((not p.print_help()))
    histfiles = args[:(- 1)]
    species = args[(- 1)]
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])
    A = fig.add_axes([0.08, 0.12, 0.38, 0.76])
    B = fig.add_axes([0.58, 0.12, 0.38, 0.76])
    lines = []
    legends = []
    genomesizes = []
    for histfile in histfiles:
        ks = KmerSpectrum(histfile)
        (x, y) = ks.get_xy(opts.vmin, opts.vmax)
        K = get_number(op.basename(histfile).split('.')[0].split('-')[(- 1)])
        if (not (opts.kmin <= K <= opts.kmax)):
            continue
        (line,) = A.plot(x, y, '-', lw=1)
        lines.append(line)
        legends.append('K = {0}'.format(K))
        ks.analyze(K=K)
        genomesizes.append((K, (ks.genomesize / 1000000.0)))
    leg = A.legend(lines, legends, shadow=True, fancybox=True)
    leg.get_frame().set_alpha(0.5)
    title = '{0} genome K-mer histogram'.format(species)
    A.set_title(markup(title))
    (xlabel, ylabel) = ('Coverage (X)', 'Counts')
    A.set_xlabel(xlabel)
    A.set_ylabel(ylabel)
    set_human_axis(A)
    title = '{0} genome size estimate'.format(species)
    B.set_title(markup(title))
    (x, y) = zip(*genomesizes)
    B.plot(x, y, 'ko', mfc='w')
    t = np.linspace((opts.kmin - 0.5), (opts.kmax + 0.5), 100)
    p = np.poly1d(np.polyfit(x, y, 2))
    B.plot(t, p(t), 'r:')
    (xlabel, ylabel) = ('K-mer size', 'Estimated genome size (Mb)')
    B.set_xlabel(xlabel)
    B.set_ylabel(ylabel)
    set_ticklabels_arial(B)
    labels = ((0.04, 0.96, 'A'), (0.54, 0.96, 'B'))
    panel_labels(root, labels)
    normalize_axes(root)
    imagename = (species + '.multiK.pdf')
    savefig(imagename, dpi=iopts.dpi, iopts=iopts)