%prog plot workdir sample chr1,chr2

    Plot some chromosomes for visual proof. Separate multiple chromosomes with
    comma. Must contain folder workdir/sample-cn/.
345955
def plot(args):
    
    from jcvi.graphics.base import savefig

    p = OptionParser(plot.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="8x7", format="png")

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

    workdir, sample_key, chrs = args
    chrs = chrs.split(",")
    hmm = CopyNumberHMM(workdir=workdir)
    hmm.plot(sample_key, chrs=chrs)

    image_name = sample_key + "_cn." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
673621
def _goleft_indexcov(bam_file, data, out_dir):
    
    if not dd.get_coverage_interval(data) == "genome":
        return []
    out_dir = utils.safe_makedir(os.path.join(out_dir, "indexcov"))
    out_files = [os.path.join(out_dir, "%s-indexcov.%s" % (dd.get_sample_name(data), ext))
                 for ext in ["roc", "ped", "bed.gz"]]
    if not utils.file_uptodate(out_files[-1], bam_file):
        with transaction.tx_tmpdir(data) as tmp_dir:
            tmp_dir = utils.safe_makedir(os.path.join(tmp_dir, dd.get_sample_name(data)))
            gender_chroms = [x.name for x in ref.file_contigs(dd.get_ref_file(data)) if chromhacks.is_sex(x.name)]
            gender_args = "--sex %s" % (",".join(gender_chroms)) if gender_chroms else ""
            cmd = "goleft indexcov --directory {tmp_dir} {gender_args} -- {bam_file}"
            try:
                do.run(cmd.format(**locals()), "QC: goleft indexcov")
            except subprocess.CalledProcessError as msg:
                if not ("indexcov: no usable" in str(msg) or
                        ("indexcov: expected" in str(msg) and "sex chromosomes, found:" in str(msg))):
                    raise
            for out_file in out_files:
                orig_file = os.path.join(tmp_dir, os.path.basename(out_file))
                if utils.file_exists(orig_file):
                    utils.copy_plus(orig_file, out_file)
    # MultiQC needs non-gzipped/BED inputs so unpack the file
    out_bed = out_files[-1].replace(".bed.gz", ".tsv")
    if utils.file_exists(out_files[-1]) and not utils.file_exists(out_bed):
        with transaction.file_transaction(data, out_bed) as tx_out_bed:
            cmd = "gunzip -c %s > %s" % (out_files[-1], tx_out_bed)
            do.run(cmd, "Unpack indexcov BED file")
    out_files[-1] = out_bed
    return [x for x in out_files if utils.file_exists(x)]
41894
def depth(args):
    
    from jcvi.utils.range import range_depth

    p = OptionParser(depth.__doc__)
    p.add_option("--depthfile",
                 help="Generate file with gene and depth [default: %default]")
    p.add_option("--histogram", default=False, action="store_true",
                 help="Plot histograms in PDF")
    p.add_option("--xmax", type="int", help="x-axis maximum to display in plot")
    p.add_option("--title", default=None, help="Title to display in plot")
    p.add_option("--quota", help="Force to use this quota, e.g. 1:1, 1:2 ...")
    p.set_beds()

    opts, args = p.parse_args(args)

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

    anchorfile, = args
    qbed, sbed, qorder, sorder, is_self = check_beds(anchorfile, p, opts)
    depthfile = opts.depthfile
    ac = AnchorFile(anchorfile)
    qranges = []
    sranges = []
    blocks = ac.blocks
    for ib in blocks:
        q, s, t = zip(*ib)
        q = [qorder[x] for x in q]
        s = [sorder[x] for x in s]
        qrange = (min(q)[0], max(q)[0])
        srange = (min(s)[0], max(s)[0])
        qranges.append(qrange)
        sranges.append(srange)
        if is_self:
            qranges.append(srange)

    qgenome = op.basename(qbed.filename).split(".")[0]
    sgenome = op.basename(sbed.filename).split(".")[0]
    qtag = "Genome {0} depths".format(qgenome)
    print("{}:".format(qtag), file=sys.stderr)
    dsq, details = range_depth(qranges, len(qbed))
    if depthfile:
        fw = open(depthfile, "w")
        write_details(fw, details, qbed)

    if is_self:
        return

    stag = "Genome {0} depths".format(sgenome)
    print("{}:".format(stag), file=sys.stderr)
    dss, details = range_depth(sranges, len(sbed))
    if depthfile:
        write_details(fw, details, sbed)
        fw.close()
        logging.debug("Depth written to `{0}`.".format(depthfile))

    if not opts.histogram:
        return

    from jcvi.graphics.base import plt, quickplot_ax, savefig, normalize_axes

    # Plot two histograms one for query genome, one for subject genome
    plt.figure(1, (6, 3))
    f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)

    xmax = opts.xmax or max(4, max(dsq.keys() + dss.keys()))
    if opts.quota:
        speak, qpeak = opts.quota.split(":")
        qpeak, speak = int(qpeak), int(speak)
    else:
        qpeak = find_peak(dsq)
        speak = find_peak(dss)

    qtag = "# of {} blocks per {} gene".format(sgenome, qgenome)
    stag = "# of {} blocks per {} gene".format(qgenome, sgenome)
    quickplot_ax(ax1, dss, 0, xmax, stag, ylabel="Percentage of genome",
                 highlight=range(1, speak + 1))
    quickplot_ax(ax2, dsq, 0, xmax, qtag, ylabel=None,
                 highlight=range(1, qpeak + 1))

    title = opts.title or "{} vs {} syntenic depths\n{}:{} pattern"\
                    .format(qgenome, sgenome, speak, qpeak)
    root = f.add_axes([0, 0, 1, 1])
    vs, pattern = title.split('\n')
    root.text(.5, .97, vs, ha="center", va="center", color="darkslategray")
    root.text(.5, .925, pattern, ha="center", va="center",
                                 color="tomato", size=16)
    print(title, file=sys.stderr)

    normalize_axes(root)

    pf = anchorfile.rsplit(".", 1)[0] + ".depth"
    image_name = pf + ".pdf"
    savefig(image_name)
805314
def lms(args):
    
    from random import randint
    from jcvi.graphics.chromosome import HorizontalChromosome

    p = OptionParser(lms.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="6x6", dpi=300)

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    # Panel A
    w, h = .7, .35
    ax = fig.add_axes([.15, .6, w, h])

    xdata = [x + randint(-3, 3) for x in range(10, 110, 10)]
    ydata = [x + randint(-3, 3) for x in range(10, 110, 10)]
    ydata[3:7] = ydata[3:7][::-1]
    xydata = zip(xdata, ydata)
    lis = xydata[:3] + [xydata[4]] + xydata[7:]
    lds = xydata[3:7]
    xlis, ylis = zip(*lis)
    xlds, ylds = zip(*lds)
    ax.plot(xlis, ylis, "r-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xlds, ylds, "g-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xdata, ydata, "k.", mec="k", mfc="w", mew=3, ms=12)
    HorizontalChromosome(root, .15, .15 + w, .57, height=.02, lw=2)
    root.text(.15 + w / 2, .55, "Chromosome location (bp)", ha="center", va="top")

    ax.text(80, 30, "LIS = 7", color="r", ha="center", va="center")
    ax.text(80, 20, "LDS = 4", color="g", ha="center", va="center")
    ax.text(80, 10, "LMS = $max$(LIS, LDS) = 7", ha="center", va="center")
    normalize_lms_axis(ax, xlim=110, ylim=110)

    # Panel B
    w = .37
    p = (0, 45, 75, 110)
    ax = fig.add_axes([.1, .12, w, h])
    xdata = [x for x in range(10, 110, 10)]
    ydata = ydata_orig = [x for x in range(10, 110, 10)]
    ydata = ydata[:4] + ydata[7:] + ydata[4:7][::-1]
    xydata = zip(xdata, ydata)
    lis = xydata[:7]
    xlis, ylis = zip(*lis)
    ax.plot(xlis, ylis, "r-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xdata, ydata, "k.", mec="k", mfc="w", mew=3, ms=12)
    ax.vlines(p, 0, 110, colors="beige", lw=3)
    normalize_lms_axis(ax, xlim=110, ylim=110)
    patch = [.1 + w * x / 110. for x in p]
    HorizontalChromosome(root, .1, .1 + w, .09, patch=patch,
                         height=.02, lw=2)
    scaffolds = ("a", "b", "c")
    for i, s in enumerate(scaffolds):
        xx = (patch[i] + patch[i + 1]) / 2
        root.text(xx, .09, s, va="center", ha="center")
    root.text(.1 + w / 2, .04, "LMS($a||b||c$) = 7", ha="center")

    # Panel C
    ax = fig.add_axes([.6, .12, w, h])
    patch = [.6 + w * x / 110. for x in p]
    ydata = ydata_orig
    ax.plot(xdata, ydata, "r-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xdata, ydata, "k.", mec="k", mfc="w", mew=3, ms=12)
    ax.vlines(p, [0], [110], colors="beige", lw=3)
    normalize_lms_axis(ax, xlim=110, ylim=110)
    HorizontalChromosome(root, .6, .6 + w, .09, patch=patch,
                         height=.02, lw=2)
    scaffolds = ("a", "-c", "b")
    for i, s in enumerate(scaffolds):
        xx = (patch[i] + patch[i + 1]) / 2
        root.text(xx, .09, s, va="center", ha="center")
    root.text(.6 + w / 2, .04, "LMS($a||-c||b$) = 10", ha="center")

    labels = ((.05, .95, 'A'), (.05, .48, 'B'), (.55, .48, 'C'))
    panel_labels(root, labels)

    normalize_axes(root)

    pf = "lms"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
689777
def main(argv=None):
    

    if argv is None:
        argv = sys.argv

    # setup command line parser
    parser = U.OptionParser(version="%prog version: $Id$",
                            usage=usage,
                            description=globals()["__doc__"])

    group = U.OptionGroup(parser, "count-specific options")

    parser.add_option("--wide-format-cell-counts", dest="wide_format_cell_counts",
                      action="store_true",
                      default=False,
                      help=("output the cell counts in a wide format "
                            "(rows=genes, columns=cells)"))

    parser.add_option_group(group)

    # add common options (-h/--help, ...) and parse command line
    (options, args) = U.Start(parser, argv=argv, add_group_dedup_options=False)

    options.per_gene = True  # hardcodes counting to per-gene only

    U.validateSamOptions(options, group=False)

    if options.random_seed:
        np.random.seed(options.random_seed)

    if options.stdin != sys.stdin:
        in_name = options.stdin.name
        options.stdin.close()
    else:
        raise ValueError("Input on standard in not currently supported")

    if options.in_sam:
        in_mode = "r"
    else:
        in_mode = "rb"

    infile = pysam.Samfile(in_name, in_mode)

    # write out to tempfile and then sort to stdout
    tmpfilename = U.getTempFilename(dir=options.tmpdir)
    tmpfile = U.openFile(tmpfilename, mode="w")

    nInput, nOutput, input_reads = 0, 0, 0

    gene_tag = options.gene_tag
    metacontig2contig = None

    if options.chrom:
        inreads = infile.fetch(reference=options.chrom)
    else:
        if options.gene_transcript_map:
            metacontig2contig = sam_methods.getMetaContig2contig(
                infile, options.gene_transcript_map)
            metatag = "MC"
            inreads = sam_methods.metafetcher(infile, metacontig2contig, metatag)
            gene_tag = metatag
        else:
            inreads = infile.fetch()

    bundle_iterator = sam_methods.get_bundles(
        options,
        only_count_reads=True,
        metacontig_contig=metacontig2contig)

    # set up UMIClusterer functor with methods specific to
    # specified options.method

    processor = network.UMIClusterer(options.method)

    for bundle, key, status in bundle_iterator(inreads):
        if status == "single_read":
            continue

        gene, cell = key

        umis = bundle.keys()
        counts = {umi: bundle[umi]["count"] for umi in umis}

        nInput += sum(counts.values())

        while nInput >= input_reads + 1000000:
            input_reads += 1000000
            U.info("Parsed %i input reads" % input_reads)

        # group the umis
        groups = processor(
            counts,
            threshold=options.threshold)

        gene_count = len(groups)

        if options.per_cell:
            tmpfile.write("%s\n" % "\t".join((gene, cell.decode(), str(gene_count))))
        else:
            tmpfile.write("%s\n" % "\t".join((gene, str(gene_count))))

        nOutput += gene_count

    tmpfile.close()

    if options.per_cell:

        gene_counts_dict = {}

        with U.openFile(tmpfilename, mode="r") as inf:
            genes = set()
            cells = set()
            for line in inf:
                gene, cell, gene_count = line.strip().split("\t")
                genes.add(gene)
                cells.add(cell)

                if gene not in gene_counts_dict:
                    gene_counts_dict[gene] = {}

                gene_counts_dict[gene][cell] = gene_count

        if options.wide_format_cell_counts:  # write out in wide format

            options.stdout.write(
                "%s\t%s\n" % ("gene", "\t".join(sorted(cells))))

            for gene in sorted(genes):
                counts = []
                for cell in sorted(cells):
                    if cell in gene_counts_dict[gene]:
                        counts.append(gene_counts_dict[gene][cell])
                    else:
                        counts.append(0)
                options.stdout.write(
                    "%s\t%s\n" % (gene, "\t".join(map(str, counts))))

        else:  # write out in long format
            options.stdout.write("%s\t%s\t%s\n" % ("gene", "cell", "count"))
            for gene in sorted(genes):
                for cell in sorted(list(gene_counts_dict[gene].keys())):
                    options.stdout.write("%s\t%s\t%s\n" % (
                        gene, cell, gene_counts_dict[gene][cell]))
    else:
        options.stdout.write("%s\t%s\n" % ("gene", "count"))

        with U.openFile(tmpfilename, mode="r") as inf:
            for line in inf:
                options.stdout.write(line)

    os.unlink(tmpfilename)

    # output reads events and benchmark information.
    for event in bundle_iterator.read_events.most_common():
        U.info("%s: %s" % (event[0], event[1]))

    U.info("Number of (post deduplication) reads counted: %i" % nOutput)

    U.Stop()
151202
def main(argv=None):
    

    if argv is None:
        argv = sys.argv

    # setup command line parser
    parser = U.OptionParser(version="%prog version: $Id$",
                            usage=usage,
                            description=globals()["__doc__"])

    group = U.OptionGroup(parser, "count_tab-specific options")

    group.add_option("--barcode-separator", dest="bc_sep",
                     type="string", help="separator between read id and UMI "
                     " and (optionally) the cell barcode", default="_")

    group.add_option("--per-cell", dest="per_cell",
                     action="store_true",
                     help="Readname includes cell barcode as well as UMI in "
                     "format: read[sep]UMI[sep]CB")

    parser.add_option_group(group)

    # add common options (-h/--help, ...) and parse command line
    (options, args) = U.Start(parser, argv=argv, add_group_dedup_options=False,
                              add_sam_options=False)

    nInput, nOutput = 0, 0

    # set the method with which to extract umis from reads
    if options.per_cell:
        bc_getter = partial(
            sam_methods.get_cell_umi_read_string, sep=options.bc_sep)
    else:
        bc_getter = partial(
            sam_methods.get_umi_read_string, sep=options.bc_sep)

    if options.per_cell:
        options.stdout.write("%s\t%s\t%s\n" % ("cell", "gene", "count"))
    else:
        options.stdout.write("%s\t%s\n" % ("gene", "count"))

    # set up UMIClusterer functor with methods specific to
    # specified options.method
    processor = network.UMIClusterer(options.method)

    for gene, counts in sam_methods.get_gene_count_tab(
            options.stdin,
            bc_getter=bc_getter):

        for cell in counts.keys():
            umis = counts[cell].keys()

            nInput += sum(counts[cell].values())

            # group the umis
            groups = processor(
                counts[cell],
                threshold=options.threshold)

            gene_count = len(groups)
            if options.per_cell:
                options.stdout.write("%s\t%s\t%i\n" % (cell, gene, gene_count))
            else:
                options.stdout.write("%s\t%i\n" % (gene, gene_count))
                nOutput += gene_count

    U.info("Number of reads counted: %i" % nOutput)

    U.Stop()
999395
def prepare(args):
    
    from jcvi.utils.table import write_csv
    from jcvi.formats.base import write_file
    from jcvi.formats.fastq import guessoffset, readlen

    p = OptionParser(prepare.__doc__ + FastqNamings)
    p.add_option("--corr", default=False, action="store_true",
                 help="Extra parameters for corrected data [default: %default]")
    p.add_option("--norun", default=False, action="store_true",
                 help="Don't write `run.sh` script [default: %default]")
    p.add_option("--ploidy", default="2", choices=("1", "2"),
                 help="Ploidy [default: %default]")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    organism_name = args[0]
    project_name = "".join(x[0] for x in organism_name.split()).upper()
    fnames = sorted(glob("*.fastq*") if len(args) == 1 else args[1:])
    for x in fnames:
        assert op.exists(x), "File `{0}` not found.".format(x)

    groupheader = "group_name library_name file_name".split()
    libheader = "library_name project_name organism_name type paired "\
        "frag_size frag_stddev insert_size insert_stddev read_orientation "\
        "genomic_start genomic_end".split()
    groups_33 = []
    groups_64 = []
    libs = []
    for file_name in fnames:
        offset = guessoffset([file_name])
        group_name = op.basename(file_name).split(".")[0]
        library_name = "-".join(group_name.split("-")[:2])

        # Handle paired files and convert to wildcard
        if ".1." in file_name:
            file_name = file_name.replace(".1.", ".?.")
        elif ".2." in file_name:
            continue

        groupscontents = groups_64 if offset == 64 else groups_33
        groupscontents.append((group_name, library_name, file_name))
        if library_name not in libs:
            libs.append(library_name)

    libcontents = []
    for library_name in libs:
        L = Library(library_name)
        size = L.size
        stddev = L.stddev
        type = L.type
        paired = L.paired
        read_orientation = L.read_orientation

        size = size or ""
        stddev = stddev or ""
        frag_size = size if type == "fragment" else ""
        frag_stddev = stddev if type == "fragment" else ""
        insert_size = size if type != "fragment" else ""
        insert_stddev = stddev if type != "fragment" else ""
        genomic_start, genomic_end = "", ""
        libcontents.append((library_name, project_name, organism_name, type, \
            paired, frag_size, frag_stddev, insert_size, insert_stddev, \
            read_orientation, genomic_start, genomic_end))

    for groups, csvfile in ((groups_33, "in_groups_33.csv"), \
                            (groups_64, "in_groups_64.csv"), \
                            (groups_33 + groups_64, "in_groups.csv")):
        if not groups:
            continue
        write_csv(groupheader, groups, filename=csvfile, tee=True)
        logging.debug("`{0}` created (# of groups = {1}).".\
            format(csvfile, len(groups)))

    write_csv(libheader, libcontents, filename="in_libs.csv", tee=True)
    logging.debug("`in_libs.csv` created (# of libs = {0}).".\
        format(len(libcontents)))

    runfile = "run.sh"

    # ALLPATHS stalls on reads over 250bp <https://www.biostars.org/p/122091/>
    max_rd_len = max(readlen([f]) for f in fnames)
    extra = "CLOSE_UNIPATH_GAPS=False " if max_rd_len > 200 else ""
    if opts.corr:
        extra += "FE_NUM_CYCLES=1 EC_K=28 FE_QUAL_CEIL_RADIUS=0"
        extra += " REMOVE_DODGY_READS_FRAG=False FE_MAX_KMER_FREQ_TO_MARK=1"

    if not opts.norun:
        contents = ALLPATHSRUN.format(opts.ploidy, opts.cpus, extra)
        write_file(runfile, contents)
822718
def Main():
  
  argument_parser = argparse.ArgumentParser(description=(
      'Plots memory usage from profiling data.'))

  argument_parser.add_argument(
      '--output', dest='output_file', type=str, help=(
          'path of the output file to write the graph to instead of using '
          'interactive mode. The output format deduced from the extension '
          'of the filename.'))

  argument_parser.add_argument(
      'profile_path', type=str, help=(
          'path to the directory containing the profiling data.'))

  options = argument_parser.parse_args()

  if not os.path.isdir(options.profile_path):
    print('No such directory: {0:s}'.format(options.profile_path))
    return False

  names = ['time', 'queued', 'processing', 'to_merge', 'abandoned', 'total']

  glob_expression = os.path.join(options.profile_path, 'task_queue-*.csv.gz')
  for csv_file_name in glob.glob(glob_expression):
    data = numpy.genfromtxt(
        csv_file_name, delimiter='\t', dtype=None, encoding='utf-8',
        names=names, skip_header=1)

    pyplot.plot(data['time'], data['queued'], label='queued')
    pyplot.plot(data['time'], data['processing'], label='processing')
    pyplot.plot(data['time'], data['to_merge'], label='to merge')
    pyplot.plot(data['time'], data['abandoned'], label='abandoned')

  pyplot.title('Number of tasks over time')

  pyplot.xlabel('Time')
  pyplot.xscale('linear')

  pyplot.ylabel('Number of tasks')
  pyplot.yscale('linear')

  pyplot.legend()

  if options.output_file:
    pyplot.savefig(options.output_file)
  else:
    pyplot.show()

  return True
738049
def p_statement_plot(p):
    
    p[0] = make_sentence('PLOT',
                         make_typecast(TYPE.ubyte, p[2], p.lineno(3)),
                         make_typecast(TYPE.ubyte, p[4], p.lineno(3)))
29086
def get_chrs_from_bam(self, file_name):
        
        x = subprocess.check_output(self.tools.samtools + " view -H " + file_name + " | grep '^@SQ' | cut -f2| sed s'/SN://'", shell=True)
        # Chromosomes will be separated by newlines; split into list to return
        return x.split()
1356054
def main():
    parser = argparse.ArgumentParser(description='Shuffle columns in a CSV file')
    parser.add_argument(metavar="FILE", dest='input_file', type=argparse.FileType('r'), nargs='?',
                        default=sys.stdin, help='Input CSV file. If omitted, read standard input.')
    parser.add_argument('-s', '--sattolo',
                        action='store_const', const=shuffle_sattolo,
                        dest='shuffle', default=random.shuffle,
                        help="Use Sattolo's shuffle algorithm.")
    col_group = parser.add_mutually_exclusive_group()
    col_group.add_argument('-c', '--columns', type=column_list,
                           help='Comma-separated list of columns to include.')
    col_group.add_argument('-C', '--no-columns', type=column_list,
                           help='Comma-separated list of columns to exclude.')
    delim_group = parser.add_mutually_exclusive_group()
    delim_group.add_argument('-d', '--delimiter', type=str, default=',',
                             help='Input column delimiter.')
    delim_group.add_argument('-t', '--tabbed', dest='delimiter',
                             action='store_const', const='\t',
                             help='Delimit input with tabs.')
    parser.add_argument('-q', '--quotechar', type=str, default='"',
                        help='Quote character.')
    parser.add_argument('-o', '--output-delimiter', type=str, default=',',
                        help='Output column delimiter.')
    parser.add_argument('-v', '--version', action='version',
                        version='%(prog)s {version}'.format(version=__version__))
    args = parser.parse_args()

    reader = csv.reader(args.input_file, delimiter=args.delimiter, quotechar=args.quotechar)

    """"""
    headers = next(reader)

    """Create a matrix of lists of columns"""
    table = []
    for c in range(len(headers)):
        table.append([])
    for row in reader:
        for c in range(len(headers)):
            table[c].append(row[c])

    cols = args.columns
    if args.no_columns:
        """If columns to exclude are provided, get a list of all other columns"""
        cols = list(set(range(len(headers))) - set(args.no_columns))
    elif not cols:
        """If no columns are provided all columns will be shuffled"""
        cols = range(len(headers))

    for c in cols:
        if c > len(headers):
            sys.stderr.write('Invalid column {0}. Last column is {1}.\n'.format(c, len(headers)))
            exit(1)
        args.shuffle(table[c - 1])

    """Transpose the matrix"""
    table = zip(*table)

    writer = csv.writer(sys.stdout, delimiter=args.output_delimiter)
    writer.writerow(headers)
    for row in table:
        writer.writerow(row)
475105
def _prep_cnv_file(cns_file, svcaller, work_dir, data):
    
    in_file = cns_file
    out_file = os.path.join(work_dir, "%s-%s-prep.csv" % (utils.splitext_plus(os.path.basename(in_file))[0],
                                                          svcaller))
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(data, out_file) as tx_out_file:
            with open(in_file) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    reader = csv.reader(in_handle, dialect="excel-tab")
                    writer = csv.writer(out_handle)
                    writer.writerow(["chrom", "start", "end", "num.mark", "seg.mean"])
                    header = next(reader)
                    for line in reader:
                        cur = dict(zip(header, line))
                        if chromhacks.is_autosomal(cur["chromosome"]):
                            writer.writerow([_to_ucsc_style(cur["chromosome"]), cur["start"],
                                             cur["end"], cur["probes"], cur["log2"]])
    return out_file
49413
def summary(args):
    
    from jcvi.formats.blast import AlignStats

    p = OptionParser(summary.__doc__)
    p.add_option("-s", dest="single", default=False, action="store_true",
            help="provide stats per reference seq")

    opts, args = p.parse_args(args)

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

    coordsfile, = args
    alignstats = get_stats(coordsfile)
    alignstats.print_stats()
172114
def plots_from_files(imspaths, figsize=(10,5), rows=1, titles=None, maintitle=None):
    
    f = plt.figure(figsize=figsize)
    if maintitle is not None: plt.suptitle(maintitle, fontsize=16)
    for i in range(len(imspaths)):
        sp = f.add_subplot(rows, ceildiv(len(imspaths), rows), i+1)
        sp.axis('Off')
        if titles is not None: sp.set_title(titles[i], fontsize=16)
        img = plt.imread(imspaths[i])
        plt.imshow(img)
1281403
def _get_purecn_files(paired, work_dir, require_exist=False):
    
    out_base = os.path.join(work_dir, "%s-purecn" % (dd.get_sample_name(paired.tumor_data)))
    out = {"plot": {}}
    all_files = []
    for plot in ["chromosomes", "local_optima", "segmentation", "summary"]:
        if plot == "summary":
            cur_file = "%s.pdf" % out_base
        else:
            cur_file = "%s_%s.pdf" % (out_base, plot)
        if not require_exist or os.path.exists(cur_file):
            out["plot"][plot] = cur_file
            all_files.append(os.path.basename(cur_file))
    for key, ext in [["hetsummary", ".csv"], ["dnacopy", "_dnacopy.seg"], ["genes", "_genes.csv"],
                     ["log", ".log"], ["loh", "_loh.csv"], ["rds", ".rds"],
                     ["variants", "_variants.csv"]]:
        cur_file = "%s%s" % (out_base, ext)
        if not require_exist or os.path.exists(cur_file):
            out[key] = cur_file
            all_files.append(os.path.basename(cur_file))
    return out_base, out, all_files
965960
def main(args=sys.argv[1:]):
    
    parser = argparse.ArgumentParser(description="""Aggregate results of
            nestly runs""")
    subparsers = parser.add_subparsers()
    delim_parser = subparsers.add_parser('delim', help="""Combine control files
            with delimited files.""")
    delim_parser.set_defaults(func=delim)
    key_group = delim_parser.add_mutually_exclusive_group()
    key_group.add_argument('-k', '--keys', help="""Comma separated list of
            keys from the JSON file to include [default: all keys]""",
            type=comma_separated_values)
    key_group.add_argument('-x', '--exclude-keys', help="""Comma separated
            list of keys from the JSON file not to include [default:
            %(default)s]""", type=comma_separated_values)
    delim_parser.add_argument('-m', '--missing-action', choices=('fail',
        'warn'), help="""Action to take when a file is missing [default:
        %(default)s]""", default='fail')
    delim_parser.add_argument('file_template', help="""Template for the
            delimited file to read in each directory [e.g. '{run_id}.csv']""")
    delim_parser.add_argument('control_files', metavar="control.json",
            help="""Control files""", nargs="*")
    delim_parser.add_argument('-d', '--directory', help="""Run on all control
            files under %(metavar)s. May be used in place of specifying control
            files.""", metavar='DIR')
    delim_parser.add_argument('-s', '--separator', default=DEFAULT_SEP,
            help="""Separator [default: %(default)s]""")
    delim_parser.add_argument('-t', '--tab', action='store_const',
            dest='separator', const='\t', help="""Files are tab-separated""")
    delim_parser.add_argument('-o', '--output', default=sys.stdout,
        type=argparse.FileType('w'), help="""Output file [default: stdout]""")

    arguments = parser.parse_args()

    arguments.func(arguments)
939841
def _remove_haplotype_chroms(in_file, data):
    
    larger_chroms = set(_get_larger_chroms(dd.get_ref_file(data)))
    out_file = os.path.join(_sv_workdir(data), "%s-chromfilter%s" % utils.splitext_plus(os.path.basename(in_file)))
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            with open(in_file) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    for line in in_handle:
                        if ((line.find("chromosome") >= 0 and line.find("start") >= 0)
                              or line.split()[0] in larger_chroms):
                            out_handle.write(line)
    return out_file