%prog covlen covfile fastafile

    Plot coverage vs length. `covfile` is two-column listing contig id and
    depth of coverage.
177840
def covlen(args):
    
    import numpy as np
    import pandas as pd
    import seaborn as sns
    from jcvi.formats.base import DictFile

    p = OptionParser(covlen.__doc__)
    p.add_option("--maxsize", default=1000000, type="int", help="Max contig size")
    p.add_option("--maxcov", default=100, type="int", help="Max contig size")
    p.add_option("--color", default='m', help="Color of the data points")
    p.add_option("--kind", default="scatter",
                 choices=("scatter", "reg", "resid", "kde", "hex"),
                 help="Kind of plot to draw")
    opts, args, iopts = p.set_image_options(args, figsize="8x8")

    if len(args) != 2:
        sys.exit(not p.print_help())

    covfile, fastafile = args
    cov = DictFile(covfile, cast=float)
    s = Sizes(fastafile)
    data = []
    maxsize, maxcov = opts.maxsize, opts.maxcov
    for ctg, size in s.iter_sizes():
        c = cov.get(ctg, 0)
        if size > maxsize:
            continue
        if c > maxcov:
            continue
        data.append((size, c))

    x, y = zip(*data)
    x = np.array(x)
    y = np.array(y)
    logging.debug("X size {0}, Y size {1}".format(x.size, y.size))

    df = pd.DataFrame()
    xlab, ylab = "Length", "Coverage of depth (X)"
    df[xlab] = x
    df[ylab] = y
    sns.jointplot(xlab, ylab, kind=opts.kind, data=df,
                  xlim=(0, maxsize), ylim=(0, maxcov),
                  stat_func=None, edgecolor="w", color=opts.color)

    figname = covfile + ".pdf"
    savefig(figname, dpi=iopts.dpi, iopts=iopts)
953076
def stats(args):
    
    from jcvi.utils.table import tabulate

    p = OptionParser(stats.__doc__)
    p.add_option("--warn", default=False, action="store_true",
                 help="Warnings on small component spans [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(p.print_help())

    agpfile, = args

    agp = AGP(agpfile)
    gap_lengths = []
    component_lengths = []
    for a in agp:
        span = a.object_span
        if a.is_gap:
            label = a.gap_type
            gap_lengths.append((span, label))
        else:
            label = "{0}:{1}-{2}".format(a.component_id, a.component_beg, \
                   a.component_end)
            component_lengths.append((span, label))
            if opts.warn and span < 50:
                logging.error("component span too small ({0}):\n{1}".\
                    format(span, a))

    table = dict()
    for label, lengths in zip(("Gaps", "Components"),
            (gap_lengths, component_lengths)):

        if not lengths:
            table[(label, "Min")] = table[(label, "Max")] \
                                  = table[(label, "Sum")] = "n.a."
            continue

        table[(label, "Min")] = "{0} ({1})".format(*min(lengths))
        table[(label, "Max")] = "{0} ({1})".format(*max(lengths))
        table[(label, "Sum")] = sum(x[0] for x in lengths)

    print(tabulate(table), file=sys.stderr)
64284
def find_default_cogfile
      parts = Dir.pwd.split File::SEPARATOR
      i = parts.length
      while i >= 0 && !File.exists?(File.join(parts.slice(0, i) + ['Cogfile']))
        i -= 1
      end
      path = File.join(parts.slice(0, i) + ['Cogfile']) if i >= 0
      (path && File.exists?(path)) ? path : nil
    end
1146871
def cov(self, other=None, pairwise=None, bias=False, **kwargs):
        
        if other is None:
            other = self._selected_obj
            # only default unset
            pairwise = True if pairwise is None else pairwise
        other = self._shallow_copy(other)

        def _get_cov(X, Y):
            X = self._shallow_copy(X)
            Y = self._shallow_copy(Y)
            cov = libwindow.ewmcov(X._prep_values(), Y._prep_values(),
                                   self.com, int(self.adjust),
                                   int(self.ignore_na), int(self.min_periods),
                                   int(bias))
            return X._wrap_result(cov)

        return _flex_binary_moment(self._selected_obj, other._selected_obj,
                                   _get_cov, pairwise=bool(pairwise))
795029
def load_covarfile(self, file, indices=[], names=[], sample_file=False):
        """"""

        # Clean up input in case we are given some empty values
        var_indices = []
        for x in indices:
            try:
                var_indices.append(int(x) + 1)
            except:
                pass

        var_names = []
        for name in names:
            if name.strip() != "":
                var_names.append(name)

        if PhenoCovar.sex_as_covariate:
            self.covariate_labels = ["SEX"]
        file.seek(0)
        if file:
            header = file.readline().strip().split()
            line_number = 0

            # Rows ignored because family/individual missing from our pedigree data
            ignored_data = []
            # We can accept a default phenotype column if we only have 3 columns
            if len(header) == 3:
                if len(var_names) + len(var_indices) == 0:
                    var_indices.append(2)

            if header[0].upper() == "FID":
                for name in var_names:
                    if name.strip() != "":
                        try:
                            var_indices.append(header.index(name))
                        except:
                            raise InvalidSelection("The name, %s, was not found in %s" % (name, file.name))
                line_number = 1


                if len(var_indices) > 0 and max(var_indices) > (len(header)):
                    raise InvalidSelection("The index, %s, is larger than the number of entries in the file, %s:%s" % (max(var_indices), file.name, line_number))

                for index in var_indices:
                    self.covariate_labels.append(header[index])
                # Dump the second line for sample_file==True
                if sample_file:
                    file.readline()
            else:
                file.seek(0)
                if len(var_names) > 0:
                    raise MalformedInputFile("Names only work with covariate files with headers: %s" % (",".join(names), file.name))
                if len(var_indices) > 0 and max(var_indices) > (len(header)):
                    raise InvalidSelection("The index, %s, is larger than the number of entries in the file, %s:%s" % (max(var_indices), file.name, line_number))

                for i in xrange(0, len(var_indices)):
                    self.covariate_labels.append("Cov-%s" % (var_indices[i]-1))

            covar_data = numpy.empty((len(var_indices)+len(self.covariate_data), len(self.pedigree_data)))
            covar_data.fill(PhenoCovar.missing_encoding)

            # We have to be careful to keep the sex covariate data if it came from the pedigree
            if PhenoCovar.sex_as_covariate:
                covar_data[0] = self.covariate_data[0]
            self.covariate_data = covar_data

            for line in file:
                line_number += 1
                words   = line.split()
                iid = ":".join(words[0:2])

                # Indexes are 1 based...silly humans
                if len(var_indices) > 0 and max(var_indices) > (len(words)):
                    raise InvalidSelection("The index, %s, is larger than the number of entries in the file, %s:%s" % (max(var_indices), file.name, line_number))

                cidx = 0
                if PhenoCovar.sex_as_covariate:
                    cidx += 1
                for idx in var_indices:
                    cov = float(words[idx])
                    if iid not in self.pedigree_data:
                        ignored_data.append(iid)
                    else:
                        self.covariate_data[cidx, self.pedigree_data[iid]] = cov
                    cidx += 1


        for covar in self.covariate_data:
            if sum(covar !=PhenoCovar.missing_encoding ) == 0:
                raise NoMatchedPhenoCovars("No matching individuals were found in the covariate file")
751101
def qual(args):
    
    from jcvi.formats.sizes import Sizes

    p = OptionParser(qual.__doc__)
    p.add_option("--qv", default=31, type="int",
                 help="Dummy qv score for extended bases")
    p.set_outfile()
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    fastafile, = args
    sizes = Sizes(fastafile)
    qvchar = str(opts.qv)
    fw = must_open(opts.outfile, "w")
    total = 0
    for s, slen in sizes.iter_sizes():
        print(">" + s, file=fw)
        print(" ".join([qvchar] * slen), file=fw)
        total += 1
    fw.close()
    logging.debug("Written {0} records in `{1}`.".format(total, opts.outfile))
629746
def _prep_streamlined_config(config_script, work_dir):
    
    new_min_size = 100
    in_file = config_script + ".ini"
    out_file = os.path.join(work_dir, os.path.basename(in_file))
    with open(in_file) as in_handle:
        with open(out_file, "w") as out_handle:
            for line in in_handle:
                if line.startswith("minCandidateVariantSize"):
                    out_handle.write("minCandidateVariantSize = %s\n" % new_min_size)
                else:
                    out_handle.write(line)
    return out_file
1251209
def script_lines_workaround(line_info, coverage_info, count_info)
      is_repeated = lambda do |div|
        n = line_info.size / div
        break false unless line_info.size % div == 0 && n > 1
        different = false
        n.times do |i|
        
          things = (0...div).map { |j| line_info[i + j * n] }
          if things.uniq.size != 1
            different = true
            break
          end
        end

        ! different
      end

      factors = braindead_factorize(line_info.size)
      factors.each do |n|
        if is_repeated[n]
          line_info = line_info[0, line_info.size / n]
          coverage_info = coverage_info[0, coverage_info.size / n]
          count_info = count_info[0, count_info.size / n]
        end
      end if factors.size > 1   # don't even try if it's prime

      [line_info, coverage_info, count_info]
    end
728043
def parse_gcov_file(args, fobj, filename):
    
    coverage = []
    ignoring = False
    for line in fobj:
        report_fields = line.decode('utf-8', 'replace').split(':', 2)
        if len(report_fields) == 1:
            continue
        line_num = report_fields[1].strip()
        if line_num == '':
            continue

        cov_num = report_fields[0].strip()
        line_num = int(line_num)
        text = report_fields[2]
        if line_num == 0:
            continue
        if re.search(r'\bLCOV_EXCL_START\b', text):
            if ignoring:
                sys.stderr.write("Warning: %s:%d: nested LCOV_EXCL_START, "
                                 "please fix\n" % (filename, line_num))
            ignoring = True
        elif re.search(r'\bLCOV_EXCL_(STOP|END)\b', text):
            if not ignoring:
                sys.stderr.write("Warning: %s:%d: LCOV_EXCL_STOP outside of "
                                 "exclusion zone, please fix\n" % (filename,
                                                                   line_num))
            if 'LCOV_EXCL_END' in text:
                sys.stderr.write("Warning: %s:%d: LCOV_EXCL_STOP is the "
                                 "correct keyword\n" % (filename, line_num))
            ignoring = False
        if cov_num == '-':
            coverage.append(None)
        elif cov_num == '#####':
            # Avoid false positives.
            if (
                ignoring or
                any([re.search(pattern, text) for pattern in args.exclude_lines_pattern])
            ):
                coverage.append(None)
            else:
                coverage.append(0)
        elif cov_num == '=====':
            # This is indicitive of a gcov output parse
            # error.
            coverage.append(0)
        else:
            coverage.append(int(cov_num.rstrip('*')))
    return coverage
896393
def fig_lfp_corr(params, savefolders, transient=200,
                 channels=[0,3,7,11,13], Df=None,
                 mlab=True, NFFT=256, noverlap=128,
                 window=plt.mlab.window_hanning,
                 letterslist=['AB', 'CD'], data_type = 'LFP'):
    ''' 

    '''
    ana_params.set_PLOS_2column_fig_style(ratio=0.5)
    fig = plt.figure()
    fig.subplots_adjust(left=0.07, right=0.95, bottom=0.1, wspace=0.3, hspace=0.1)

    gs = gridspec.GridSpec(5, 4)


    for i, (savefolder, letters) in enumerate(zip(savefolders, letterslist)):
        # path to simulation files
        params.savefolder = os.path.join(os.path.split(params.savefolder)[0],
                                         savefolder)
        params.figures_path = os.path.join(params.savefolder, 'figures')
        params.spike_output_path = os.path.join(params.savefolder,
                                                           'processed_nest_output')
        params.networkSimParams['spike_output_path'] = params.spike_output_path


        ## Including correlations
        f = h5py.File(os.path.join(params.savefolder, ana_params.analysis_folder, data_type + ana_params.fname_psd),'r')
        freqs = f['freqs'].value
        LFP_PSD_corr = f['psd'].value
        f.close()
    
        ## Excluding correlations
        f = h5py.File(os.path.join(params.savefolder, ana_params.analysis_folder, data_type +  ana_params.fname_psd_uncorr),'r')
        freqs = f['freqs'].value
        LFP_PSD_uncorr = f['psd'].value
        f.close()
        
     
        ##################################
        ###  Single channel LFP PSDs   ###
        ##################################
    
        ax = fig.add_subplot(gs[0, (i % 2)*2]) 
        phlp.remove_axis_junk(ax)
        ax.loglog(freqs,LFP_PSD_corr[channels[0]], color='k', label='$P$')
        ax.loglog(freqs,LFP_PSD_uncorr[channels[0]],
                  color='gray' if analysis_params.bw else analysis_params.colorP,
                  lw=1,
                  label='$\tilde{P}$')
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')
        ax.text(0.80,0.82,'ch. %i' %(channels[0]+1),horizontalalignment='left',
            verticalalignment='center',
            fontsize=6, 
            transform=ax.transAxes)
        ax.yaxis.set_minor_locator(plt.NullLocator())
        ax.set_ylabel('(mV$^2$/Hz)', labelpad=0.)
        ax.set_xticks([])
        ax.set_xticklabels([])
        ax.tick_params(axis='y',which='minor',bottom='off')
        ax.set_xlim(4E0,4E2)
        ax.set_ylim(1E-8,1.5E-4)
        ax.set_yticks([1E-8,1E-6,1E-4])
        ax.set_title('power spectra')
        phlp.annotate_subplot(ax, ncols=4, nrows=5, letter=letters[0],
                           linear_offset=0.065)       
    
        ax = fig.add_subplot(gs[1, (i % 2)*2])
        phlp.remove_axis_junk(ax)
        ax.loglog(freqs,LFP_PSD_corr[channels[1]], color='k', label='corr')
        ax.loglog(freqs,LFP_PSD_uncorr[channels[1]],
                  color='gray' if analysis_params.bw else analysis_params.colorP,
                  lw=1,
                  label='uncorr')
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')
        ax.text(0.80,0.82,'ch. %i' %(channels[1]+1),horizontalalignment='left',
            verticalalignment='center',
            fontsize=6, 
            transform=ax.transAxes)
        ax.yaxis.set_minor_locator(plt.NullLocator())
        ax.set_xticks([])
        ax.set_xticklabels([])
        ax.tick_params(axis='y',which='minor',bottom='off')
        ax.set_xlim(4E0,4E2)
        ax.set_ylim(1E-8,1.5E-4)
        ax.set_yticks([1E-8,1E-6,1E-4])
        ax.set_yticklabels([])
    
    
        ax = fig.add_subplot(gs[2, (i % 2)*2])
        phlp.remove_axis_junk(ax)
        ax.loglog(freqs,LFP_PSD_corr[channels[2]], color='k', label='corr')
        ax.loglog(freqs,LFP_PSD_uncorr[channels[2]],
                  color='gray' if analysis_params.bw else analysis_params.colorP,
                  lw=1,
                  label='uncorr')   
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')
        ax.text(0.80,0.82,'ch. %i' %(channels[2]+1),horizontalalignment='left',
            verticalalignment='center',
            fontsize=6, 
            transform=ax.transAxes)
        ax.yaxis.set_minor_locator(plt.NullLocator())
        ax.set_xticks([])
        ax.set_xticklabels([])
        ax.tick_params(axis='y',which='minor',bottom='off')
        ax.set_xlim(4E0,4E2)
        ax.set_ylim(1E-8,1.5E-4)
        ax.set_yticks([1E-8,1E-6,1E-4])
        ax.set_yticklabels([])
    
    
        ax = fig.add_subplot(gs[3, (i % 2)*2]) 
        phlp.remove_axis_junk(ax)
        ax.loglog(freqs,LFP_PSD_corr[channels[3]], color='k', label='corr')
        ax.loglog(freqs,LFP_PSD_uncorr[channels[3]],
                  color='gray' if analysis_params.bw else analysis_params.colorP,
                  lw=1,
                  label='uncorr')   
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')
        ax.text(0.80,0.82,'ch. %i' %(channels[3]+1),horizontalalignment='left',
            verticalalignment='center',
            fontsize=6, 
            transform=ax.transAxes)
        ax.yaxis.set_minor_locator(plt.NullLocator())
        ax.set_xticks([])
        ax.set_xticklabels([])
        ax.tick_params(axis='y',which='minor',bottom='off')
        ax.set_xlim(4E0,4E2)
        ax.set_ylim(1E-8,1.5E-4)
        ax.set_yticks([1E-8,1E-6,1E-4])
        ax.set_yticklabels([])
    
    
        ax = fig.add_subplot(gs[4, (i % 2)*2])
        phlp.remove_axis_junk(ax)
        ax.loglog(freqs,LFP_PSD_corr[channels[4]], color='k', label='corr')
        ax.loglog(freqs,LFP_PSD_uncorr[channels[4]],
                  color='gray' if analysis_params.bw else analysis_params.colorP,
                  lw=1,
                  label='uncorr')   
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')
        ax.set_xlabel(r'$f$ (Hz)', labelpad=0.2)
        ax.text(0.80,0.82,'ch. %i' %(channels[4]+1),horizontalalignment='left',
            verticalalignment='center',
            fontsize=6, 
            transform=ax.transAxes)
        ax.yaxis.set_minor_locator(plt.NullLocator())
        ax.tick_params(axis='y',which='minor',bottom='off')
        ax.set_xlim(4E0,4E2)
        ax.set_ylim(1E-8,1.5E-4)
        ax.set_yticks([1E-8,1E-6,1E-4])
        ax.set_yticklabels([])
    
    
        ##################################
        ###  LFP PSD ratios            ###
        ##################################
    
        ax = fig.add_subplot(gs[:, (i % 2)*2 + 1]) 
        phlp.annotate_subplot(ax, ncols=4, nrows=1, letter=letters[1],
                           linear_offset=0.065)
        phlp.remove_axis_junk(ax)
        ax.set_title('power ratio')
        PSD_ratio = LFP_PSD_corr/LFP_PSD_uncorr
        
        zvec = np.r_[params.electrodeParams['z']]
        zvec = np.r_[zvec, zvec[-1] + np.diff(zvec)[-1]]
        inds = freqs >= 1  # frequencies greater than 4 Hz  
        im = ax.pcolormesh(freqs[inds], zvec+40, PSD_ratio[:, inds],
                            rasterized=False,
                            cmap=plt.get_cmap('gray_r', 12) if analysis_params.bw else plt.get_cmap('Reds', 12),
                            vmin=10**-0.25,vmax=10**2.75,norm=LogNorm())
        ax.set_xscale('log')
    
        ax.set_yticks(zvec)
        yticklabels = ['ch. %i' %i for i in np.arange(len(zvec))+1]
        ax.set_yticklabels(yticklabels)
        ax.set_xlabel(r'$f$ (Hz)',labelpad=0.2)
        plt.axis('tight')
        ax.set_xlim([4E0, 4E2])
    
        cb = phlp.colorbar(fig, ax, im,
                           width=0.05, height=0.5,
                           hoffset=-0.05, voffset=0.0)
        cb.set_label('(-)', labelpad=0.1)


 
    return fig
482615
def evidences(args):
    
    p = OptionParser(evidences.__doc__)
    p.add_option("--csv", default="hli.20170328.tred.tsv",
                 help="TRED csv output to plot")
    opts, args, iopts = p.set_image_options(args, format="pdf")

    if len(args) != 0:
        sys.exit(not p.print_help())

    format = iopts.format

    # Extract sample coverage first
    df = pd.read_csv("qc-export-MeanCoverage.csv", header=None,
                     names=["Samplekey", "MeanCoverage"], index_col=0)

    # Find coverage for HD
    xf = pd.read_csv(opts.csv, sep="\t", index_col=0)
    dp = {}
    tred = "HD"
    for sk, row in xf.iterrows():
        sk = str(sk)
        a1 = row[tred + ".1"]
        a2 = row[tred + ".2"]
        fdp = row[tred + ".FDP"]
        pdp = row[tred + ".PDP"]
        pedp = row[tred + ".PEDP"]
        dp[sk] = (a1, a2, fdp, pdp, pedp)

    # Build a consolidated dataframe
    ef = pd.DataFrame.from_dict(dp, orient="index")
    ef.columns = [tred + ".1", tred + ".2", tred + ".FDP",
                  tred + ".PDP", tred + ".PEDP"]
    ef.index.name = "SampleKey"
    mf = df.merge(ef, how="right", left_index=True, right_index=True)

    # Plot a bunch of figures
    outdir = "output"
    mkdir(outdir)
    xlim = ylim = (0, 100)
    draw_jointplot(outdir + "/A", "MeanCoverage", "HD.FDP",
                   data=mf, xlim=xlim, ylim=ylim, format=format)
    draw_jointplot(outdir + "/B", "MeanCoverage", "HD.PDP",
                   data=mf, color='g', xlim=xlim, ylim=ylim, format=format)
    draw_jointplot(outdir + "/C", "MeanCoverage", "HD.PEDP",
                   data=mf, color='m', xlim=xlim, ylim=ylim, format=format)

    xlim = (0, 50)
    draw_jointplot(outdir + "/D", "HD.2", "HD.FDP",
                   data=mf, xlim=xlim, ylim=ylim, format=format)
    draw_jointplot(outdir + "/E", "HD.2", "HD.PDP",
                   data=mf, color='g', xlim=xlim, ylim=ylim, format=format)
    draw_jointplot(outdir + "/F", "HD.2", "HD.PEDP",
                   data=mf, color='m', xlim=xlim, ylim=ylim, format=format)
511061
def fasta(args):
    
    from jcvi.formats.sizes import Sizes
    from jcvi.formats.agp import OO, build

    p = OptionParser(fasta.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    bedfile, scffasta, pmolfasta = args
    pf = bedfile.rsplit(".", 1)[0]
    bed = Bed(bedfile)
    selected = select_bed(bed)
    oo = OO()
    seen = set()
    sizes = Sizes(scffasta).mapping
    agpfile = pf + ".agp"
    agp = open(agpfile, "w")
    for b in selected:
        scf = range_parse(b.accn).seqid
        chr = b.seqid
        cs = (chr, scf)
        if cs not in seen:
            oo.add(chr, scf, sizes[scf], b.strand)
            seen.add(cs)
        else:
            logging.debug("Seen {0}, ignored.".format(cs))

    oo.write_AGP(agp, gaptype="contig")
    agp.close()
    build([agpfile, scffasta, pmolfasta])
62368
def fastq(args):
    
    from jcvi.formats.fastq import FastqLite

    p = OptionParser(fastq.__doc__)
    p.add_option("--qv", type="int",
                 help="Use generic qv value [dafault: %default]")

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    fastafile, = args
    fastqfile = fastafile.rsplit(".", 1)[0] + ".fastq"
    fastqhandle = open(fastqfile, "w")
    num_records = 0

    if opts.qv is not None:
        qv = chr(ord('!') + opts.qv)
        logging.debug("QV char '{0}' ({1})".format(qv, opts.qv))
    else:
        qv = None

    if qv:
        f = Fasta(fastafile, lazy=True)
        for name, rec in f.iteritems_ordered():
            r = FastqLite("@" + name, str(rec.seq).upper(), qv * len(rec.seq))
            print(r, file=fastqhandle)
            num_records += 1

    else:
        qualfile = get_qual(fastafile)
        for rec in iter_fasta_qual(fastafile, qualfile):
            SeqIO.write([rec], fastqhandle, "fastq")
            num_records += 1

    fastqhandle.close()
    logging.debug("A total of %d records written to `%s`" % \
            (num_records, fastqfile))
667749
def cov(args):
    
    p = OptionParser(cov.__doc__)
    p.add_option("--order",
                default="swede,kale,h165,yudal,aviso,abu,bristol,bzh",
                help="The order to plot the tracks, comma-separated")
    p.add_option("--reverse", default=False, action="store_true",
                help="Plot the order in reverse")
    p.add_option("--gauge_step", default=5000000, type="int",
                help="Step size for the base scale")
    p.add_option("--hlsuffix", default="regions.forhaibao",
                help="Suffix for the filename to be used to highlight regions")
    opts, args, iopts = p.set_image_options(args, figsize="11x8")

    if len(args) != 4:
        sys.exit(not p.print_help())

    chr1, chr2, sizesfile, datadir = args
    chr1 = chr1.split(",")
    chr2 = chr2.split(",")

    order = opts.order
    hlsuffix = opts.hlsuffix
    if order:
        order = order.split(",")
        if opts.reverse:
            order.reverse()
    sizes = Sizes(sizesfile).mapping
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    chrs = (chr1, chr2)
    chr_sizes, chr_sum_sizes, ratio = calc_ratio(chrs, sizes)
    chr_size1, chr_size2 = chr_sum_sizes
    chr_sizes1, chr_sizes2 = chr_sizes

    w1_start, w1_end = center_panel(chr1, chr_size1, ratio)
    w2_start, w2_end = center_panel(chr2, chr_size2, ratio)
    w1s = w1_start
    w2s = w2_start

    dsg = "gray"
    i = 0
    for c1, s1 in zip(chr1, chr_sizes1):
        w1 = ratio * s1
        plot_label = i == 0
        i += 1
        canvas1 = (w1s, .6, w1, .3)
        Coverage(fig, root, canvas1, c1, (0, s1), datadir,
                     order=order, gauge="top", plot_label=plot_label,
                     gauge_step=opts.gauge_step, palette=dsg,
                     cap=40, hlsuffix=hlsuffix)
        w1s += w1 + gap

    i = 0
    for c2, s2 in zip(chr2, chr_sizes2):
        w2 = ratio * s2
        plot_label = i == 0
        i += 1
        canvas2 = (w2s, .15, w2, .3)
        Coverage(fig, root, canvas2, c2, (0, s2), datadir,
                     order=order, gauge="bottom", plot_label=plot_label,
                     gauge_step=opts.gauge_step, palette=dsg,
                     cap=40, hlsuffix=hlsuffix)
        w2s += w2 + gap

    # Synteny panel
    seqidsfile = make_seqids(chrs)
    klayout = make_layout(chrs, chr_sum_sizes, ratio, template_cov)
    Karyotype(fig, root, seqidsfile, klayout, gap=gap,
              generank=False, sizes=sizes)

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    chr2 = "_".join(chr2)
    if opts.reverse:
        chr2 += ".reverse"
    image_name = chr2 + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
203673
function populateCoverage(cov) {
    cov.LOC =
    cov.SLOC =
    cov.totalFiles =
    cov.totalHits =
    cov.totalMisses =
    cov.coverage = 0;
    for (var name in cov) {
        var file = cov[name];
        if (Array.isArray(file)) {
            // Stats
            ++cov.totalFiles;
            cov.totalHits += file.totalHits = coverage(file, true);
            cov.totalMisses += file.totalMisses = coverage(file, false);
            file.totalLines = file.totalHits + file.totalMisses;
            cov.SLOC += file.SLOC = file.totalLines;
            if (!file.source) file.source = [];
            cov.LOC += file.LOC = file.source.length;
            file.coverage = (file.totalHits / file.totalLines) * 100;
            // Source
            var width = file.source.length.toString().length;
            file.source = file.source.map(function(line, i){
                ++i;
                var hits = file[i] === 0 ? 0 : (file[i] || ' ');
                if (hits === 0) {
                    hits = '\x1b[31m' + hits + '\x1b[0m';
                    line = '\x1b[41m' + line + '\x1b[0m';
                } else {
                    hits = '\x1b[32m' + hits + '\x1b[0m';
                }
                return '\n     ' + lpad(i, width) + ' | ' + hits + ' | ' + line;
            }).join('');
        }
    }
    cov.coverage = (cov.totalHits / cov.SLOC) * 100;
}
1149615
def multilineplot(args):
    
    p = OptionParser(multilineplot.__doc__)
    p.add_option("--lines",
                 help="Features to plot in lineplot [default: %default]")
    p.add_option("--colors",
                 help="List of colors matching number of input bed files")
    p.add_option("--mode", default="span", choices=("span", "count", "score"),
                 help="Accumulate feature based on [default: %default]")
    p.add_option("--binned", default=False, action="store_true",
                 help="Specify whether the input is already binned; " +
                 "if True, input files are considered to be binfiles")
    p.add_option("--ymax", type="int", help="Set Y-axis max")
    add_window_options(p)
    opts, args, iopts = p.set_image_options(args, figsize="8x5")

    if len(args) != 2:
        sys.exit(not p.print_help())

    fastafile, chr = args
    window, shift, subtract, merge = check_window_options(opts)
    linebeds = []
    colors = opts.colors
    if opts.lines:
        lines = opts.lines.split(",")
        assert len(colors) == len(lines), "Number of chosen colors must match" + \
                " number of input bed files"
        linebeds = get_beds(lines, binned=opts.binned)

    linebins = get_binfiles(linebeds, fastafile, shift, mode=opts.mode,
                            binned=opts.binned, merge=merge)

    clen = Sizes(fastafile).mapping[chr]
    nbins = get_nbins(clen, shift)

    plt.rcParams["xtick.major.size"] = 0
    plt.rcParams["ytick.major.size"] = 0
    plt.rcParams["figure.figsize"] = iopts.w, iopts.h

    fig, axarr = plt.subplots(nrows=len(lines))
    if len(linebeds) == 1:
        axarr = (axarr, )
    fig.suptitle(latex(chr), color="darkslategray")

    for i, ax in enumerate(axarr):
        lineplot(ax, [linebins[i]], nbins, chr, window, shift, \
                color="{0}{1}".format(colors[i], 'r'))

    if opts.ymax:
        ax.set_ylim(0, opts.ymax)

    plt.subplots_adjust(hspace=0.5)

    image_name = chr + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
173485
function reportCoverage(cov) {
    // Stats
    print('\n   [bold]{Test Coverage}\n');
    var sep = '   +------------------------------------------+----------+------+------+--------+',
        lastSep = '                                              +----------+------+------+--------+';
	result = sep+'\n';
    result += '   | filename                                 | coverage | LOC  | SLOC | missed |\n';
    result += sep+'\n';
    for (var name in cov) {
        var file = cov[name];
        if (Array.isArray(file)) {
            result += '   | ' + rpad(name, 40);
            result += ' | ' + lpad(file.coverage.toFixed(2), 8);
            result += ' | ' + lpad(file.LOC, 4);
            result += ' | ' + lpad(file.SLOC, 4);
            result += ' | ' + lpad(file.totalMisses, 6);
            result += ' |\n';
        }
    }
    result += sep+'\n';
    result += '     ' + rpad('', 40);
    result += ' | ' + lpad(cov.coverage.toFixed(2), 8);
    result += ' | ' + lpad(cov.LOC, 4);
    result += ' | ' + lpad(cov.SLOC, 4);
    result += ' | ' + lpad(cov.totalMisses, 6);
    result += ' |\n';
    result += lastSep;
	console.log(result);

	for (var name in cov) {
		if (name.match(file_matcher)) {
			var file = cov[name];
			var annotated = '';
			annotated += colorize('\n  [bold]{' + name + '}:');
			annotated += colorize(file.source);
			annotated += '\n';
			fs.writeFileSync('annotated/'+name, annotated);
		}
	}
}