%prog swap agpfile

    Swap objects and components. Will add gap lines. This is often used in
    conjuction with formats.chain.fromagp() to convert between different
    coordinate systems.
546617
def swap(args):
    
    from jcvi.utils.range import range_interleave

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

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

    agpfile, = args

    agp = AGP(agpfile, nogaps=True, validate=False)
    agp.sort(key=lambda x: (x.component_id, x.component_beg))

    newagpfile = agpfile.rsplit(".", 1)[0] + ".swapped.agp"
    fw = open(newagpfile, "w")
    agp.transfer_header(fw)
    for cid, aa in groupby(agp, key=(lambda x: x.component_id)):
        aa = list(aa)
        aranges = [(x.component_id, x.component_beg, x.component_end) \
                    for x in aa]
        gaps = range_interleave(aranges)
        for a, g in zip_longest(aa, gaps):
            a.object, a.component_id = a.component_id, a.object
            a.component_beg = a.object_beg
            a.component_end = a.object_end
            print(a, file=fw)
            if not g:
                continue

            aline = [cid, 0, 0, 0]
            gseq, ga, gb = g
            cspan = gb - ga + 1
            aline += ["N", cspan, "fragment", "yes"]
            print("\t".join(str(x) for x in aline), file=fw)

    fw.close()
    # Reindex
    idxagpfile = reindex([newagpfile, "--inplace"])

    return newagpfile
552807
def fromsoap(args):
    
    p = OptionParser(fromsoap.__doc__)
    p.add_option("--type", default="nucleotide_match",
                 help="GFF feature type [default: %default]")
    p.add_option("--source", default="soap",
                help="GFF source qualifier [default: %default]")
    p.set_fixchrnames(orgn="maize")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    soapfile, = args
    pad0 = len(str(sum(1 for line in open(soapfile))))

    fw = must_open(opts.outfile, "w")
    fp = must_open(soapfile)
    for idx, line in enumerate(fp):
        if opts.fix_chr_name:
            from jcvi.utils.cbook import fixChromName
            line = fixChromName(line, orgn=opts.fix_chr_name)

        atoms = line.strip().split("\t")
        attributes = "ID=match{0};Name={1}".format(str(idx).zfill(pad0), atoms[0])
        start, end = int(atoms[8]), int(atoms[5]) + int(atoms[8]) - 1
        seqid = atoms[7]

        print("\t".join(str(x) for x in (seqid, opts.source, opts.type, \
            start, end, ".", atoms[6], ".", attributes)), file=fw)
617503
def flip(args):
    
    p = OptionParser(flip.__doc__)
    opts, args = p.parse_args(args)

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

    fastafile, = args
    outfastafile = fastafile.rsplit(".", 1)[0] + ".flipped.fasta"
    fo = open(outfastafile, "w")
    f = Fasta(fastafile, lazy=True)
    for name, rec in f.iteritems_ordered():
        tmpfasta = "a.fasta"
        fw = open(tmpfasta, "w")
        SeqIO.write([rec], fw, "fasta")
        fw.close()

        o = overlap([tmpfasta, name])
        if o.orientation == '-':
            rec.seq = rec.seq.reverse_complement()

        SeqIO.write([rec], fo, "fasta")
        os.remove(tmpfasta)
79401
def insert(args):
    
    from jcvi.formats.agp import mask, bed
    from jcvi.formats.sizes import agp

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

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

    candidates, gapsbed, chrfasta, unplacedfasta = args
    refinedbed = refine([candidates, gapsbed])
    sizes = Sizes(unplacedfasta).mapping
    cbed = Bed(candidates)
    corder = cbed.order
    gbed = Bed(gapsbed)
    gorder = gbed.order

    gpbed = Bed()
    gappositions = {}  # (chr, start, end) => gapid

    fp = open(refinedbed)
    gap_to_scf = defaultdict(list)
    seen = set()
    for row in fp:
        atoms = row.split()
        if len(atoms) <= 6:
            continue
        unplaced = atoms[3]
        strand = atoms[5]
        gapid = atoms[9]
        if gapid not in seen:
            seen.add(gapid)
            gi, gb = gorder[gapid]
            gpbed.append(gb)
            gappositions[(gb.seqid, gb.start, gb.end)] = gapid
        gap_to_scf[gapid].append((unplaced, strand))

    gpbedfile = "candidate.gaps.bed"
    gpbed.print_to_file(gpbedfile, sorted=True)

    agpfile = agp([chrfasta])
    maskedagpfile = mask([agpfile, gpbedfile])
    maskedbedfile = maskedagpfile.rsplit(".", 1)[0] + ".bed"
    bed([maskedagpfile, "--outfile={0}".format(maskedbedfile)])

    mbed = Bed(maskedbedfile)
    finalbed = Bed()
    for b in mbed:
        sid = b.seqid
        key = (sid, b.start, b.end)
        if key not in gappositions:
            finalbed.add("{0}\n".format(b))
            continue

        gapid = gappositions[key]
        scfs = gap_to_scf[gapid]

        # For scaffolds placed in the same gap, sort according to positions
        scfs.sort(key=lambda x: corder[x[0]][1].start + corder[x[0]][1].end)
        for scf, strand in scfs:
            size = sizes[scf]
            finalbed.add("\t".join(str(x) for x in \
                    (scf, 0, size, sid, 1000, strand)))

    finalbedfile = "final.bed"
    finalbed.print_to_file(finalbedfile)

    # Clean-up
    toclean = [gpbedfile, agpfile, maskedagpfile, maskedbedfile]
    FileShredder(toclean)
1125315
def gaps(args):
    
    from jcvi.formats.base import DictFile
    from jcvi.apps.base import popen
    from jcvi.utils.cbook import percentage

    p = OptionParser(gaps.__doc__)
    p.add_option("--bdist", default=0, type="int",
                 help="Base pair distance [default: %default]")
    opts, args = p.parse_args(args)

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

    idsfile, frfile, gapsbed = args
    bdist = opts.bdist
    d =  DictFile(frfile, keypos=1, valuepos=2)
    bedfile = idsfile + ".bed"
    fw = open(bedfile, "w")
    fp = open(idsfile)
    total = 0
    for row in fp:
        id = row.strip()
        hit = d[id]
        tag, pos = get_tag(hit, None)
        seqid, start, end = pos
        start, end = max(start - bdist, 1), end + bdist
        print("\t".join(str(x) for x in (seqid, start - 1, end, id)), file=fw)
        total += 1
    fw.close()

    cmd = "intersectBed -a {0} -b {1} -v | wc -l".format(bedfile, gapsbed)
    not_in_gaps = popen(cmd).read()
    not_in_gaps = int(not_in_gaps)
    in_gaps = total - not_in_gaps
    print("Ids in gaps: {1}".\
            format(total, percentage(in_gaps, total)), file=sys.stderr)
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)
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)
1199766
def merge(args):
    
    p = OptionParser(merge.__doc__)
    p.add_option("--seq", default=False, action="store_true",
                 help="Print FASTA sequences at the end")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    if nargs == 1:
        listfile, = args
        fp = open(listfile)
        gffiles = [x.strip() for x in fp]
    else:
        gffiles = args

    outfile = opts.outfile

    deflines = set()
    fw = must_open(outfile, "w")
    fastarecs = {}
    for gffile in natsorted(gffiles, key=lambda x: op.basename(x)):
        logging.debug(gffile)
        fp = open(gffile)
        for row in fp:
            row = row.rstrip()
            if not row or row[0] == '#':
                if row == FastaTag:
                    break
                if row in deflines:
                    continue
                else:
                    deflines.add(row)

            print(row, file=fw)

        if not opts.seq:
            continue

        f = Fasta(gffile, lazy=True)
        for key, rec in f.iteritems_ordered():
            if key in fastarecs:
                continue
            fastarecs[key] = rec

    if opts.seq:
        print(FastaTag, file=fw)
        SeqIO.write(fastarecs.values(), fw, "fasta")

    fw.close()
553147
def write_agp(self, obj, sizes, fw=sys.stdout, gapsize=100,
                  gaptype="contig", evidence="map"):
        '''
        '''
        contigorder = [(x.contig_name, x.strand) for x in self]
        order_to_agp(obj, contigorder, sizes, fw,
                     gapsize=gapsize, gaptype=gaptype, evidence=evidence)
864831
def plotall(xargs):
    
    p = OptionParser(plotall.__doc__)
    add_allmaps_plot_options(p)
    opts, args, iopts = p.set_image_options(xargs, figsize="10x6")

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

    inputbed, = args
    pf = inputbed.rsplit(".", 1)[0]
    agpfile = pf + ".chr.agp"
    agp = AGP(agpfile)
    objects = [ob for ob, lines in agp.iter_object()]
    for seqid in natsorted(objects):
        plot(xargs + [seqid])
161490
def bed(args):
    
    from jcvi.formats.obo import validate_term

    p = OptionParser(bed.__doc__)
    p.add_option("--gaps", default=False, action="store_true",
            help="Only print bed lines for gaps [default: %default]")
    p.add_option("--nogaps", default=False, action="store_true",
            help="Do not print bed lines for gaps [default: %default]")
    p.add_option("--bed12", default=False, action="store_true",
            help="Produce bed12 formatted output [default: %default]")
    p.add_option("--component", default=False, action="store_true",
            help="Generate bed file for components [default: %default]")
    p.set_outfile()
    g1 = OptionGroup(p, "GFF specific parameters",
            "Note: If not specified, output will be in `bed` format")
    g1.add_option("--gff", default=False, action="store_true",
            help="Produce gff3 formatted output. By default, ignores " +\
                 "AGP gap lines. [default: %default]")
    g1.add_option("--source", default="MGSC",
            help="Specify a gff3 source [default: `%default`]")
    g1.add_option("--feature", default="golden_path_fragment",
            help="Specify a gff3 feature type [default: `%default`]")
    p.add_option_group(g1)
    p.set_SO_opts()

    opts, args = p.parse_args(args)

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

    if opts.component:
        opts.nogaps = True

    # If output format is gff3 and 'verifySO' option is invoked, validate the SO term
    if opts.gff and opts.verifySO:
        validate_term(opts.feature, method=opts.verifySO)

    agpfile, = args
    agp = AGP(agpfile)
    fw = must_open(opts.outfile, "w")
    if opts.gff:
        print("##gff-version 3", file=fw)

    for a in agp:
        if opts.nogaps and a.is_gap:
            continue
        if opts.gaps and not a.is_gap:
            continue
        if opts.bed12:
            print(a.bed12line, file=fw)
        elif opts.gff:
            print(a.gffline(gff_source=opts.source, gff_feat_type=opts.feature), file=fw)
        elif opts.component:
            name = "{0}:{1}-{2}".\
                        format(a.component_id, a.component_beg, a.component_end)
            print("\t".join(str(x) for x in (a.component_id, a.component_beg - 1,
                                a.component_end, name,
                                a.component_type, a.orientation)), file=fw)
        else:
            print(a.bedline, file=fw)
    fw.close()

    return fw.name
160664
def refine(args):
    
    p = OptionParser(refine.__doc__)
    p.add_option("--closest", default=False, action="store_true",
                 help="In case of no gaps, use closest [default: %default]")
    opts, args = p.parse_args(args)

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

    breakpointsbed, gapsbed = args
    ncols = len(open(breakpointsbed).next().split())
    logging.debug("File {0} contains {1} columns.".format(breakpointsbed, ncols))
    cmd = "intersectBed -wao -a {0} -b {1}".format(breakpointsbed, gapsbed)

    pf = "{0}.{1}".format(breakpointsbed.split(".")[0], gapsbed.split(".")[0])
    ingapsbed = pf + ".bed"
    sh(cmd, outfile=ingapsbed)

    fp = open(ingapsbed)
    data = [x.split() for x in fp]

    nogapsbed = pf + ".nogaps.bed"
    largestgapsbed = pf + ".largestgaps.bed"
    nogapsfw = open(nogapsbed, "w")
    largestgapsfw = open(largestgapsbed, "w")
    for b, gaps in groupby(data, key=lambda x: x[:ncols]):
        gaps = list(gaps)
        gap = gaps[0]
        if len(gaps) == 1 and gap[-1] == "0":
            assert gap[-3] == "."
            print("\t".join(b), file=nogapsfw)
            continue

        gaps = [(int(x[-1]), x) for x in gaps]
        maxgap = max(gaps)[1]
        print("\t".join(maxgap), file=largestgapsfw)

    nogapsfw.close()
    largestgapsfw.close()
    beds = [largestgapsbed]
    toclean = [nogapsbed, largestgapsbed]

    if opts.closest:
        closestgapsbed = pf + ".closestgaps.bed"
        cmd = "closestBed -a {0} -b {1} -d".format(nogapsbed, gapsbed)
        sh(cmd, outfile=closestgapsbed)
        beds += [closestgapsbed]
        toclean += [closestgapsbed]
    else:
        pointbed = pf + ".point.bed"
        pbed = Bed()
        bed = Bed(nogapsbed)
        for b in bed:
            pos = (b.start + b.end) / 2
            b.start, b.end = pos, pos
            pbed.append(b)
        pbed.print_to_file(pointbed)
        beds += [pointbed]
        toclean += [pointbed]

    refinedbed = pf + ".refined.bed"
    FileMerger(beds, outfile=refinedbed).merge()

    # Clean-up
    FileShredder(toclean)

    return refinedbed
836136
void swap(int element1, int element2) {
		int ar1[] = m_buffer[element1 >> m_blockPower];
		int ar2[] = m_buffer[element2 >> m_blockPower];
		int ind1 = element1 & m_blockMask;
		int ind2 = element2 & m_blockMask;
		for (int i = 0; i < m_stride; i++) {
			int tmp = ar1[ind1 + i];
			ar1[ind1 + i] = ar2[ind2 + i];
			ar2[ind2 + i] = tmp;
		}
	}
204639
def frombed(args):
    
    p = OptionParser(frombed.__doc__)
    p.add_option("--type", default="match",
                 help="GFF feature type [default: %default]")
    p.add_option("--source", default="default",
                help="GFF source qualifier [default: %default]")
    opts, args = p.parse_args(args)

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

    bedfile, = args
    bed = Bed(bedfile)

    for b in bed:
        print(b.gffline(type=opts.type, source=opts.source))
70351
def fromagp(args):
    
    from jcvi.formats.agp import AGP
    from jcvi.formats.sizes import Sizes

    p = OptionParser(fromagp.__doc__)
    p.add_option("--novalidate", default=False, action="store_true",
                 help="Do not validate AGP")
    opts, args = p.parse_args(args)

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

    agpfile, componentfasta, objectfasta = args
    chainfile = agpfile.rsplit(".", 1)[0] + ".chain"
    fw = open(chainfile, "w")
    agp = AGP(agpfile, validate=(not opts.novalidate))
    componentsizes = Sizes(componentfasta).mapping
    objectsizes = Sizes(objectfasta).mapping
    chain = "chain"
    score = 1000
    tStrand = "+"
    id = 0
    for a in agp:
        if a.is_gap:
            continue

        tName = a.component_id
        tSize = componentsizes[tName]
        tStart = a.component_beg
        tEnd = a.component_end
        tStart -= 1

        qName = a.object
        qSize = objectsizes[qName]
        qStrand = "-" if a.orientation == "-" else "+"
        qStart = a.object_beg
        qEnd = a.object_end
        if qStrand == '-':
            _qStart = qSize - qEnd + 1
            _qEnd = qSize - qStart + 1
            qStart, qEnd = _qStart, _qEnd
        qStart -= 1

        id += 1
        size = a.object_span
        headerline = "\t".join(str(x) for x in (
             chain, score, tName, tSize, tStrand, tStart,
             tEnd, qName, qSize, qStrand, qStart, qEnd, id
        ))
        alignmentline = size
        print(headerline, file=fw)
        print(alignmentline, file=fw)
        print(file=fw)

    fw.close()
    logging.debug("File written to `{0}`.".format(chainfile))
1134200
def conv(arg, source, target, filename, section):
    
    logging.info(_('This is Binary Conversion mode.'))
    section = section.encode('utf-8')
    if source == 'sec':
        arg = open(arg, 'rb')
    if source == 'sec':
        kwargs = dict(section_name=section)
    else:
        kwargs = dict()
    result = Converter.uni_from(source, arg, **kwargs).uni_to(target)
    if result:
        if filename:
            logging.info(
                _('Writing shellcode to the file: %s'), filename)
            mode = 'wb' if target == 'bin' else 'w'
            with open(filename, mode) as output:
                output.write(result)
        else:
            print(result)
    else:
        logging.error(_('Failed.'))
    if source == 'sec':
        arg.close()
    return 0
1044105
def fromaligns(args):
    
    p = OptionParser(fromaligns.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    alignsfile, = args
    fp = must_open(alignsfile)
    fw = must_open(opts.outfile, "w")
    for row in fp:
        if row.startswith("## Alignment"):
            print("###", file=fw)
            continue
        if row[0] == '#' or not row.strip():
            continue
        atoms = row.split(':')[-1].split()
        print("\t".join(atoms[:2]), file=fw)
    fw.close()