%prog bed anchorsfile

    Convert ANCHORS file to BED format.
556347
def bed(args):
    
    from collections import defaultdict
    from jcvi.compara.synteny import AnchorFile, check_beds
    from jcvi.formats.bed import Bed
    from jcvi.formats.base import get_number

    p = OptionParser(bed.__doc__)
    p.add_option("--switch", default=False, action="store_true",
                 help="Switch reference and aligned map elements")
    p.add_option("--scale", type="float",
                 help="Scale the aligned map distance by factor")
    p.set_beds()
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    anchorsfile, = args
    switch = opts.switch
    scale = opts.scale
    ac = AnchorFile(anchorsfile)
    pairs = defaultdict(list)
    for a, b, block_id in ac.iter_pairs():
        pairs[a].append(b)

    qbed, sbed, qorder, sorder, is_self = check_beds(anchorsfile, p, opts)
    bd = Bed()
    for q in qbed:
        qseqid, qstart, qend, qaccn = q.seqid, q.start, q.end, q.accn
        if qaccn not in pairs:
            continue
        for s in pairs[qaccn]:
            si, s = sorder[s]
            sseqid, sstart, send, saccn = s.seqid, s.start, s.end, s.accn
        if switch:
            qseqid, sseqid = sseqid, qseqid
            qstart, sstart = sstart, qstart
            qend, send = send, qend
            qaccn, saccn = saccn, qaccn
        if scale:
            sstart /= scale
        try:
            newsseqid = get_number(sseqid)
        except ValueError:
            raise ValueError("`{0}` is on `{1}` with no number to extract".\
                                format(saccn, sseqid))
        bedline = "\t".join(str(x) for x in (qseqid, qstart - 1, qend,
                            "{0}:{1}".format(newsseqid, sstart)))
        bd.add(bedline)

    bd.print_to_file(filename=opts.outfile, sorted=True)
831923
def bed(args):
    
    from jcvi.formats.bed import sort

    p = OptionParser(bed.__doc__)
    p.add_option("--blockonly", default=False, action="store_true",
                 help="Only print out large blocks, not fragments [default: %default]")
    p.add_option("--point", default=False, action="store_true",
                 help="Print accesssion as single point instead of interval")
    p.add_option("--scale", type="float",
                 help="Scale the OM distance by factor")
    p.add_option("--switch", default=False, action="store_true",
                 help="Switch reference and aligned map elements [default: %default]")
    p.add_option("--nosort", default=False, action="store_true",
                 help="Do not sort bed [default: %default]")
    opts, args = p.parse_args(args)

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

    xmlfile, = args
    bedfile = xmlfile.rsplit(".", 1)[0] + ".bed"

    om = OpticalMap(xmlfile)
    om.write_bed(bedfile, point=opts.point, scale=opts.scale,
                          blockonly=opts.blockonly, switch=opts.switch)

    if not opts.nosort:
        sort([bedfile, "--inplace"])
6563
def alignextend(args):
    
    p = OptionParser(alignextend.__doc__)
    p.add_option("--len", default=100, type="int",
                 help="Extend to this length")
    p.add_option("--qv", default=31, type="int",
                 help="Dummy qv score for extended bases")
    p.add_option("--bedonly", default=False, action="store_true",
                 help="Only generate bed files, no FASTA")
    p.set_bedpe()
    opts, args = p.parse_args(args)

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

    bedpe, ref = args
    qvchar = chr(opts.qv + 33)
    pf = bedpe.split(".")[0]

    filtered = bedpe + ".filtered"
    if need_update(bedpe, filtered):
        filter_bedpe(bedpe, filtered, ref, rc=opts.rc,
                     minlen=opts.minlen, maxlen=opts.maxlen, rlen=opts.rlen)

    rmdup = filtered + ".filtered.sorted.rmdup"
    if need_update(filtered, rmdup):
        rmdup_bedpe(filtered, rmdup, dupwiggle=opts.dup)

    if opts.bedonly:
        return

    bed1, bed2 = pf + ".1e.bed", pf + ".2e.bed"
    if need_update(rmdup, (bed1, bed2)):
        sh("cut -f1-3,7-9 {0}".format(rmdup), outfile=bed1)
        sh("cut -f4-6,7-8,10 {0}".format(rmdup), outfile=bed2)

    sfa1, sfa2 = pf + ".1e.sfa", pf + ".2e.sfa"
    if need_update((bed1, bed2, ref), (sfa1, sfa2)):
        for bed in (bed1, bed2):
            fastaFromBed(bed, ref, name=True, tab=True, stranded=True)

    fq1, fq2 = pf + ".1e.fq", pf + ".2e.fq"
    if need_update((sfa1, sfa2), (fq1, fq2)):
        for sfa in (sfa1, sfa2):
            sfa_to_fq(sfa, qvchar)
810306
def bed_to_bedpe(bedfile, bedpefile, pairsbedfile=None, matesfile=None, ca=False, strand=False):
    
    fp = must_open(bedfile)
    fw = must_open(bedpefile, "w")
    if pairsbedfile:
        fwpairs = must_open(pairsbedfile, "w")

    clones = defaultdict(list)
    for row in fp:
        b = BedLine(row)
        name = b.accn
        clonename = clone_name(name, ca=ca)
        clones[clonename].append(b)

    if matesfile:
        fp = open(matesfile)
        libraryline = next(fp)
        # 'library bes     37896   126916'
        lib, name, smin, smax = libraryline.split()
        assert lib == "library"
        smin, smax = int(smin), int(smax)
        logging.debug("Happy mates for lib {0} fall between {1} - {2}".\
                      format(name, smin, smax))

    nbedpe = 0
    nspan = 0
    for clonename, blines in clones.items():
        nlines = len(blines)
        if nlines == 2:
            a, b = blines
            aseqid, astart, aend = a.seqid, a.start, a.end
            bseqid, bstart, bend = b.seqid, b.start, b.end
            outcols = [aseqid, astart - 1, aend, bseqid, bstart - 1, bend, clonename]
            if strand:
                outcols.extend([0, a.strand, b.strand])
            print("\t".join(str(x) for x in outcols), file=fw)
            nbedpe += 1
        elif nlines == 1:
            a, = blines
            aseqid, astart, aend = a.seqid, a.start, a.end
            bseqid, bstart, bend = 0, 0, 0
        else:  # More than two lines per pair
            pass

        if pairsbedfile:
            start = min(astart, bstart) if bstart > 0 else astart
            end = max(aend, bend) if bend > 0 else aend
            if aseqid != bseqid:
                continue

            span = end - start + 1
            if (not matesfile) or (smin <= span <= smax):
                print("\t".join(str(x) for x in \
                        (aseqid, start - 1, end, clonename)), file=fwpairs)
                nspan += 1

    fw.close()
    logging.debug("A total of {0} bedpe written to `{1}`.".\
                  format(nbedpe, bedpefile))
    if pairsbedfile:
        fwpairs.close()
        logging.debug("A total of {0} spans written to `{1}`.".\
                      format(nspan, pairsbedfile))
967099
def bed(args):
    
    from jcvi.formats.bed import sort as bed_sort

    p = OptionParser(bed.__doc__)
    p.add_option("--swap", default=False, action="store_true",
                 help="Write query positions [default: %default]")
    p.add_option("--both", default=False, action="store_true",
                 help="Generate one line for each of query and subject")

    opts, args = p.parse_args(args)

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

    blastfile, = args
    positive = (not opts.swap) or opts.both
    negative = opts.swap or opts.both

    fp = must_open(blastfile)
    bedfile =  "{0}.bed".format(blastfile.rsplit(".", 1)[0]) \
            if blastfile.endswith(".blast") \
            else "{0}.bed".format(blastfile)
    fw = open(bedfile, "w")
    for row in fp:
        b = BlastLine(row)
        if positive:
            print(b.bedline, file=fw)
        if negative:
            print(b.swapped.bedline, file=fw)

    logging.debug("File written to `{0}`.".format(bedfile))
    fw.close()
    bed_sort([bedfile, "-i"])

    return bedfile
723000
def bed(args):
    
    p = OptionParser(bed.__doc__)
    opts, args = p.parse_args(args)

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

    idsfile, = args
    ids = set(x.strip() for x in open(idsfile))
    data = get_bed_from_phytozome(list(ids))

    pf = idsfile.rsplit(".", 1)[0]
    bedfile = pf + ".bed"
    fw = open(bedfile, "w")
    for i, row in enumerate(data):
        row = row.strip()
        if row == "":
            continue

        print(row, file=fw)

    logging.debug("A total of {0} records written to `{1}`.".format(i + 1, bedfile))
1094557
def offdiag(args):
    
    p = OptionParser(offdiag.__doc__)
    p.set_beds()
    opts, args = p.parse_args(args)

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

    anchorsfile, = args
    qbed, sbed, qorder, sorder, is_self = check_beds(anchorsfile, p, opts)

    fp = open(anchorsfile)
    pf = "-".join(anchorsfile.split(".")[:2])
    header = "Block-id|Napus|Diploid|Napus-chr|Diploid-chr|RBH?".split("|")
    print("\t".join(header))
    i = -1
    for row in fp:
        if row[0] == '#':
            i += 1
            continue
        q, s, score = row.split()
        rbh = 'no' if score[-1] == 'L' else 'yes'
        qi, qq = qorder[q]
        si, ss = sorder[s]
        oqseqid = qseqid = qq.seqid
        osseqid = sseqid = ss.seqid
        sseqid = sseqid.split("_")[0][-3:]
        if qseqid[0] == 'A':
            qseqid = qseqid[-3:]       # A09 => A09
        elif qseqid[0] == 'C':
            qseqid = 'C0' + qseqid[-1]  # C9 => C09
        else:
            continue
        if qseqid == sseqid or sseqid[-2:] == 'nn':
            continue
        block_id = pf + "-block-{0}".format(i)
        print("\t".join((block_id, q, s, oqseqid, osseqid, rbh)))
945580
def format(args):
    
    from jcvi.formats.obo import load_GODag, validate_term

    valid_multiparent_ops = ["split", "merge"]

    p = OptionParser(format.__doc__)

    g1 = OptionGroup(p, "Parameter(s) used to modify GFF attributes (9th column)")
    g1.add_option("--name", help="Add Name attribute from two-column file [default: %default]")
    g1.add_option("--note", help="Add Note attribute from two-column file [default: %default]")
    g1.add_option("--add_attribute", dest="attrib_files", help="Add new attribute(s) " + \
                "from two-column file(s); attribute name comes from filename; " + \
                "accepts comma-separated list of files [default: %default]")
    g1.add_option("--add_dbxref", dest="dbxref_files", help="Add new Dbxref value(s) (DBTAG:ID) " + \
                "from two-column file(s). DBTAG comes from filename, ID comes from 2nd column; " + \
                "accepts comma-separated list of files [default: %default]")
    g1.add_option("--nostrict", default=False, action="store_true",
                 help="Disable strict parsing of GFF file and/or mapping file [default: %default]")
    g1.add_option("--remove_attr", dest="remove_attrs", help="Specify attributes to remove; " + \
                "accepts comma-separated list of attribute names [default: %default]")
    g1.add_option("--copy_id_attr_to_name", default=False, action="store_true",
                 help="Copy `ID` attribute value to `Name`, when `Name` is not defined")
    g1.add_option("--invent_name_attr", default=False, action="store_true",
                 help="Invent `Name` attribute for 2nd level child features; " + \
                "Formatted like PARENT:FEAT_TYPE:FEAT_INDEX")
    g1.add_option("--no_keep_attr_order", default=False, action="store_true",
                 help="Do not maintain attribute order [default: %default]")
    p.add_option_group(g1)

    g2 = OptionGroup(p, "Parameter(s) used to modify content within columns 1-8")
    g2.add_option("--seqid", help="Switch seqid from two-column file. If not" + \
                " a file, value will globally replace GFF seqid [default: %default]")
    g2.add_option("--source", help="Switch GFF source from two-column file. If not" + \
                " a file, value will globally replace GFF source [default: %default]")
    g2.add_option("--type", help="Switch GFF feature type from two-column file. If not" + \
                " a file, value will globally replace GFF type [default: %default]")
    g2.add_option("--fixphase", default=False, action="store_true",
                 help="Change phase 1<->2, 2<->1 [default: %default]")
    p.add_option_group(g2)

    g3 = OptionGroup(p, "Other parameter(s) to perform manipulations to the GFF " + \
                 "file content")
    g3.add_option("--unique", default=False, action="store_true",
                 help="Make IDs unique [default: %default]")
    g3.add_option("--chaindup", default=None, dest="duptype",
                 help="Chain duplicate features of a particular GFF3 `type`," + \
                      " sharing the same ID attribute [default: %default]")
    g3.add_option("--multiparents", default=None, choices=valid_multiparent_ops,
                 help="Split/merge identical features (same `seqid`, `source`, `type` " + \
                 "`coord-range`, `strand`, `phase`) mapping to multiple parents " + \
                 "[default: %default]")
    g3.add_option("--remove_feats", help="Comma separated list of features to remove by type" + \
                " [default: %default]")
    g3.add_option("--remove_feats_by_ID", help="List of features to remove by ID;" + \
                " accepts comma-separated list or list file [default: %default]")
    g3.add_option("--gsac", default=False, action="store_true",
                 help="Fix GSAC GFF3 file attributes [default: %default]")
    g3.add_option("--invent_protein_feat", default=False, action="store_true",
                 help="Invent a protein feature span (chain CDS feats)")
    g3.add_option("--process_ftype", default=None, type="str",
                 help="Specify feature types to process; "
                 "accepts comma-separated list of feature types [default: %default]")
    g3.add_option("--gff3", default=False, action="store_true",
                 help="Print output in GFF3 format [default: %default]")
    g3.add_option("--make_gff_store", default=False, action="store_true",
                 help="Store entire GFF file in memory during first iteration [default: %default]")
    p.add_option_group(g3)

    p.set_outfile()
    p.set_SO_opts()

    opts, args = p.parse_args(args)

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

    gffile, = args
    mapfile = opts.seqid
    names = opts.name
    note = opts.note
    source = opts.source
    ftype = opts.type
    attrib_files = opts.attrib_files.split(",") if opts.attrib_files else None
    dbxref_files = opts.dbxref_files.split(",") if opts.dbxref_files else None
    remove_attrs = opts.remove_attrs.split(",") if opts.remove_attrs else None
    process_ftype = opts.process_ftype.split(",") if opts.process_ftype else None
    gsac = opts.gsac
    assert not (opts.unique and opts.duptype), \
        "Cannot use `--unique` and `--chaindup` together"
    assert not(opts.type and opts.duptype), \
        "Cannot use `--type` and `--chaindup` together"
    unique = opts.unique
    duptype = opts.duptype
    fixphase = opts.fixphase
    phaseT = {"1":"2", "2":"1"}
    remove_feats = opts.remove_feats.split(",") if opts.remove_feats else None
    remove_feats_by_ID = None
    if opts.remove_feats_by_ID:
        remove_feats_by_ID = LineFile(opts.remove_feats_by_ID, load=True).lines \
                if op.isfile(opts.remove_feats_by_ID) else \
                opts.remove_feats_by_ID.split(",")
    strict = False if opts.nostrict else True
    make_gff_store = True if gffile in ("-", "stdin") else opts.make_gff_store
    assert not (opts.copy_id_attr_to_name and opts.invent_name_attr), \
        "Cannot use `--copy_id_attr_to_name` and `--invent_name_attr` together"
    copy_id_attr_to_name = opts.copy_id_attr_to_name
    invent_name_attr = opts.invent_name_attr
    invent_protein_feat = opts.invent_protein_feat
    compute_signature = False

    outfile = opts.outfile

    mapping = None
    mod_attrs = set()
    if mapfile and op.isfile(mapfile):
        mapping = DictFile(mapfile, delimiter="\t", strict=strict)
        mod_attrs.add("ID")
    if note:
        note = DictFile(note, delimiter="\t", strict=strict)
        mod_attrs.add("Note")
    if source and op.isfile(source):
        source = DictFile(source, delimiter="\t", strict=strict)
    if ftype and op.isfile(ftype):
        ftype = DictFile(ftype, delimiter="\t", strict=strict)
    if names:
        names = DictFile(names, delimiter="\t", strict=strict)
        mod_attrs.add("Name")

    if attrib_files:
        attr_values = {}
        for fn in attrib_files:
            attr_name = op.basename(fn).rsplit(".", 1)[0]
            if attr_name not in reserved_gff_attributes:
                attr_name = attr_name.lower()
            attr_values[attr_name] = DictFile(fn, delimiter="\t", strict=strict)
            mod_attrs.add(attr_name)
    if dbxref_files:
        dbxref_values = {}
        for fn in dbxref_files:
            dbtag = op.basename(fn).rsplit(".", 1)[0]
            dbxref_values[dbtag] = DictFile(fn, delimiter="\t", strict=strict)
        mod_attrs.add("Dbxref")

    if remove_attrs:
        mod_remove_attrs = []
        for remove_attr in remove_attrs:
            if remove_attr in mod_attrs:
                mod_remove_attrs.append(remove_attr)

        if mod_remove_attrs:
            logging.error("Attributes `{0}` cannot be removed and modified".format( \
                    ",".join(mod_remove_attrs)))
            sys.exit()

    if gsac:  # setting gsac will force IDs to be unique
        unique = True
        notes = {}

    remove = set()
    if unique or duptype or remove_feats or remove_feats_by_ID \
            or opts.multiparents == "merge" or invent_name_attr or make_gff_store \
            or invent_protein_feat:
        if unique:
            dupcounts = defaultdict(int)
            seen = defaultdict(int)
            newparentid = {}
        elif duptype:
            dupranges = AutoVivification()
            skip = defaultdict(int)
        if opts.multiparents == "merge":
            merge_feats = AutoVivification()
        if invent_name_attr:
            ft = GffFeatureTracker()
        elif copy_id_attr_to_name:
            pass
        if invent_protein_feat:
            cds_track = {}
        if opts.multiparents == "merge" or invent_name_attr:
            make_gff_store = compute_signature = True
        gff = Gff(gffile, keep_attr_order=(not opts.no_keep_attr_order), \
                make_gff_store=make_gff_store, compute_signature=compute_signature, \
                strict=strict)
        for g in gff:
            if process_ftype and g.type not in process_ftype:
                continue
            id = g.accn
            if remove_feats and g.type in remove_feats:
                remove.add(id)
            if remove_feats_by_ID and id in remove_feats_by_ID:
                remove.add(id)
            if unique:
                dupcounts[id] += 1
            elif duptype and g.type == duptype:
                dupranges[g.seqid][id][g.idx] = (g.start, g.end)
            if opts.multiparents == "merge" and g.type != "CDS": #don't merge CDS
                pp = g.get_attr("Parent", first=False)
                if pp and len(pp) > 0:
                    for parent in pp:
                        if parent not in remove:
                            sig = g.sign
                            if sig not in merge_feats:
                                merge_feats[sig]['parents'] = []
                            if parent not in merge_feats[sig]['parents']:
                                merge_feats[sig]['parents'].append(parent)
            if invent_name_attr:
                parent, iso = atg_name(g.get_attr("Parent"), retval="locus,iso")
                if not parent:
                    parent = g.get_attr("Parent")
                ft.track(parent, g)
            if invent_protein_feat:
                if g.type == 'CDS':
                    cds_parent = g.get_attr("Parent")
                    if cds_parent not in cds_track:
                        cds_track[cds_parent] = []
                    cds_track[cds_parent].append((g.start, g.end))

    if opts.verifySO:
        so = load_GODag()
        valid_soterm = {}

    fw = must_open(outfile, "w")
    if not make_gff_store:
        gff = Gff(gffile, keep_attr_order=(not opts.no_keep_attr_order), \
                strict=strict)
    for g in gff:
        if process_ftype and g.type not in process_ftype:
            print(g, file=fw)
            continue

        id = g.accn

        if opts.multiparents == "merge" and g.type != "CDS": #don't merge CDS
            sig = g.sign
            if len(merge_feats[sig]['parents']) > 1:
                if 'candidate' not in merge_feats[sig]:
                    merge_feats[sig]['candidate'] = id
                    g.set_attr("Parent", merge_feats[sig]['parents'])
                else:
                    continue

        if remove_feats or remove_feats_by_ID:
            if id in remove:
                continue
            else:
                if "Parent" in g.attributes:
                    keep, parent = [], g.get_attr("Parent", first=False)
                    for i, pid in enumerate(parent):
                        if pid not in remove:
                            keep.append(parent[i])
                        else:
                            remove.add(id)
                    if len(keep) == 0:
                        continue
                    parent = g.set_attr("Parent", keep)

        if remove_attrs:
            for remove_attr in remove_attrs:
                if remove_attr in g.attributes:
                    g.set_attr(remove_attr, None)

        if opts.verifySO:
            if g.type not in valid_soterm:
                valid_soterm[g.type] = validate_term(g.type, so=so, method=opts.verifySO)
            ntype = valid_soterm[g.type]
            if ntype and g.type != ntype:
                g.type = ntype

        origid = g.seqid
        if fixphase:
            phase = g.phase
            g.phase = phaseT.get(phase, phase)

        if mapfile:
            if isinstance(mapping, dict):
                if origid in mapping:
                    g.seqid = mapping[origid]
                else:
                    logging.error("{0} not found in `{1}`. ID unchanged.".\
                            format(origid, mapfile))
            else:
                g.seqid = mapfile

        if source:
            if isinstance(source, dict) and g.source in source:
                g.source = source[g.source]
            else:
                g.source = source

        if names:
            if id in names:
                g.set_attr("Name", names[id])

        if note:
            name = g.get_attr("Name")
            tag = None
            if id in note:
                tag = note[id]
            elif name and name in note:
                tag = note[name]

            if tag:
                g.set_attr("Note", tag, update=False)

        if attrib_files:
            for attr_name in attr_values:
                name = g.get_attr("Name")
                if id in attr_values[attr_name]:
                    g.set_attr(attr_name, attr_values[attr_name][id])
                elif name and name in attr_values[attr_name]:
                    g.set_attr(attr_name, attr_values[attr_name][name])

        if dbxref_files:
            for dbtag in dbxref_values:
                if id in dbxref_values[dbtag]:
                    g.set_attr("Dbxref", dbxref_values[dbtag][id], dbtag=dbtag, append=True)

        if unique:
            if dupcounts[id] > 1:
                seen[id] += 1
                old_id = id
                id = "{0}-{1}".format(old_id, seen[old_id])
                newparentid[old_id] = id
                g.set_attr("ID", id)

            if "Parent" in g.attributes:
                parent = g.attributes["Parent"][0]
                if dupcounts[parent] > 1:
                    g.set_attr("Parent", newparentid[parent])

        if duptype:
            if duptype == g.type and len(dupranges[g.seqid][id]) > 1:
                p = sorted(dupranges[g.seqid][id])
                s, e = dupranges[g.seqid][id][p[0]][0:2]  # get coords of first encountered feature
                if g.start == s and g.end == e and p[0] == g.idx:
                    r = [dupranges[g.seqid][id][x] for x in dupranges[g.seqid][id]]
                    g.start, g.end = range_minmax(r)
                else:
                    skip[(g.seqid, g.idx, id, g.start, g.end)] = 1

        if gsac and g.type == "gene":
            notes[id] = g.attributes["Name"]

        if ftype:
            if isinstance(ftype, dict) and g.type in ftype:
                g.type = ftype[g.type]
            else:
                g.type = ftype

        if invent_name_attr:
            ft.store_symbol(g)
            if re.search(ft.ftype, g.type):
                parent, iso = atg_name(g.get_attr("Parent"), retval="locus,iso")
                if not parent:
                    parent = g.get_attr("Parent")
                if parent in ft.tracker:
                    fidx = ft.feat_index(parent, g.type, g.strand, (g.start, g.end, g.sign))
                    symbol = ft.get_symbol(parent)
                    attr = "ID" if symbol == parent else "Name"
                    g.set_attr(attr, "{0}:{1}:{2}".format(symbol, g.type, fidx + 1))
                    if opts.multiparents == "merge" and attr == "Name":
                        g.set_attr("ID", "{0}:{1}:{2}".format(parent, g.type, fidx + 1))
        elif copy_id_attr_to_name:
            if "Name" not in g.attributes.keys():
                g.set_attr("Name", g.get_attr("ID"))

        protein_feat = None
        if invent_protein_feat:
            if g.type == 'mRNA':
                if id in cds_track:
                    pstart, pstop = range_minmax(cds_track[id])
                    protein_feat = GffLine("\t".join(str(x) for x in [g.seqid, g.source, "protein", pstart, pstop, \
                            ".", g.strand, ".", "ID={0}-Protein;Name={0};Derives_from={0}".format(id)]))
            elif g.type == 'CDS':
                parent = g.get_attr("Parent")
                if parent in cds_track:
                    _parent = [parent, "{0}-Protein".format(parent)]
                    g.set_attr("Parent", _parent)

        pp = g.get_attr("Parent", first=False)
        if opts.multiparents == "split" and (pp and len(pp) > 1) and g.type != "CDS":  # separate features with multiple parents
            id = g.get_attr("ID")
            for i, parent in enumerate(pp):
                if id: g.set_attr("ID", "{0}-{1}".format(id, i + 1))
                g.set_attr("Parent", parent, update=True, urlquote=True)
                if gsac:
                    fix_gsac(g, notes)
                print(g, file=fw)
        else:
            if g.gff3 and not opts.gff3:
                opts.gff3 = True
            g.update_attributes(gff3=opts.gff3)
            if gsac:
                fix_gsac(g, notes)
            if duptype == g.type and skip[(g.seqid, g.idx, id, g.start, g.end)] == 1:
                continue
            print(g, file=fw)
            if g.type == 'mRNA' and invent_protein_feat:
                print(protein_feat, file=fw)

    fw.close()
727207
def get_minibam_bed(bamfile, bedfile, minibam=None):
    
    pf = op.basename(bedfile).split(".")[0]
    minibamfile = minibam or op.basename(bamfile).replace(".bam", ".{}.bam".format(pf))
    minisamfile = minibam.replace(".bam", ".sam")
    baifile = minibamfile + ".bai"
    if op.exists(baifile):
        sh("rm {}".format(baifile))

    cmd = "samtools view -H {} > {}".format(bamfile, minisamfile)
    sh(cmd)

    cmd = "cat {}".format(bedfile)
    cmd += " | perl -lane 'print \"$F[0]:$F[1]-$F[2]\"'"
    cmd += " | xargs -n1 -t -I \{\}"
    cmd += " samtools view {}".format(bamfile)
    cmd += " \{\} >> " + minisamfile
    sh(cmd)

    cmd = "samtools view {} -b".format(minisamfile)
    cmd += " | samtools sort -"
    cmd += " -o {0}".format(minibamfile)

    sh(cmd)
    sh("samtools index {0}".format(minibamfile))
    return minibamfile
771881
def agp(args):
    
    p = OptionParser(agp.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    odir, contigsfasta = args
    fwagp = must_open(opts.outfile, 'w')
    orderingfiles = natsorted(iglob(odir, "*.ordering"))
    sizes = Sizes(contigsfasta).mapping
    contigs = set(sizes.keys())
    anchored = set()

    for ofile in orderingfiles:
        co = ContigOrdering(ofile)
        anchored |= set([x.contig_name for x in co])
        obj = op.basename(ofile).split('.')[0]
        co.write_agp(obj, sizes, fwagp)

    singletons = contigs - anchored
    logging.debug('Anchored: {}, Singletons: {}'.
                  format(len(anchored), len(singletons)))

    for s in natsorted(singletons):
        order_to_agp(s, [(s, "?")], sizes, fwagp)
690139
def collinear(args):
    
    p = OptionParser(collinear.__doc__)
    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)

    af = AnchorFile(anchorfile)
    newanchorfile = anchorfile.rsplit(".", 1)[0] + ".collinear.anchors"
    fw = open(newanchorfile, "w")

    blocks = af.blocks
    for block in blocks:
        print("#" * 3, file=fw)
        iblock = []
        for q, s, score in block:
            qi, q = qorder[q]
            si, s = sorder[s]
            score = int(long(score))
            iblock.append([qi, si, score])

        block = get_collinear(iblock)

        for q, s, score in block:
            q = qbed[q].accn
            s = sbed[s].accn
            print("\t".join((q, s, str(score))), file=fw)

    fw.close()
875469
def clr(args):
    
    p = OptionParser(clr.__doc__)
    opts, args = p.parse_args(args)

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

    blastfile = args[0]
    fastafiles = args[1:]

    sizes = {}
    for fa in fastafiles:
        f = Fasta(fa)
        sizes.update(f.itersizes())

    b = Blast(blastfile)
    for query, hits in b.iter_hits():

        qsize = sizes[query]
        vectors = list((x.qstart, x.qstop) for x in hits)
        vmin, vmax = range_minmax(vectors)

        left_size = vmin - 1
        right_size = qsize - vmax

        if left_size > right_size:
            clr_start, clr_end = 0, vmin
        else:
            clr_start, clr_end = vmax, qsize

        print("\t".join(str(x) for x in (query, clr_start, clr_end)))
        del sizes[query]

    for q, size in sorted(sizes.items()):
        print("\t".join(str(x) for x in (q, 0, size)))
474353
def mergebed(args):
    
    p = OptionParser(mergebed.__doc__)
    p.add_option("-w", "--weightsfile", default="weights.txt",
                 help="Write weights to file")
    p.set_outfile("out.bed")
    opts, args = p.parse_args(args)

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

    maps = args
    outfile = opts.outfile
    fp = must_open(maps)
    b = Bed()
    mapnames = set()
    for row in fp:
        mapname = filename_to_mapname(fp.filename())
        mapnames.add(mapname)
        try:
            m = BedLine(row)
            m.accn = "{0}-{1}".format(mapname, m.accn)
            m.extra = ["{0}:{1}".format(m.seqid, m.start)]
            b.append(m)
        except (IndexError, ValueError):  # header or mal-formed line
            continue

    b.print_to_file(filename=outfile, sorted=True)
    logging.debug("A total of {0} markers written to `{1}`.".\
                        format(len(b), outfile))

    assert len(maps) == len(mapnames), "You have a collision in map names"
    write_weightsfile(mapnames, weightsfile=opts.weightsfile)
280236
def closest(args):
    
    p = OptionParser(closest.__doc__)
    p.add_option("--om", default=False, action="store_true",
                 help="The bedfile is OM blocks [default: %default]")
    opts, args = p.parse_args(args)

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

    candidates, gapsbed, fastafile = args
    sizes = Sizes(fastafile).mapping
    bed = Bed(candidates)
    ranges = []
    for b in bed:
        r = range_parse(b.accn) if opts.om else b
        ranges.append([r.seqid, r.start, r.end])

    gapsbed = Bed(gapsbed)
    granges = [(x.seqid, x.start, x.end) for x in gapsbed]

    ranges = range_merge(ranges)
    for r in ranges:
        a = range_closest(granges, r)
        b = range_closest(granges, r, left=False)
        seqid = r[0]

        if a is not None and a[0] != seqid:
            a = None
        if b is not None and b[0] != seqid:
            b = None

        mmin = 1 if a is None else a[1]
        mmax = sizes[seqid] if b is None else b[2]

        print("\t".join(str(x) for x in (seqid, mmin - 1, mmax)))
587243
def _input_to_bed(theta_input, work_dir, get_coords, headers):
    
    theta_bed = os.path.join(work_dir, "%s.bed" % os.path.splitext(os.path.basename(theta_input))[0])
    with open(theta_input) as in_handle:
        with open(theta_bed, "w") as out_handle:
            for line in in_handle:
                if not line.startswith(headers):
                    chrom, start, end = get_coords(line)
                    out_handle.write("\t".join([chrom, start, end]) + "\n")
    return theta_bed
1204592
def batch(args):
    
    from jcvi.apps.grid import MakeManager

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

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

    cdsfile = args[0]
    anchors = args[1:]
    workdirs = [".".join(op.basename(x).split(".")[:2]) for x in anchors]
    for wd in workdirs:
        mkdir(wd)

    mm = MakeManager()
    for wd, ac in zip(workdirs, anchors):
        pairscdsfile = wd + ".cds.fasta"
        cmd = "python -m jcvi.apps.ks prepare {} {} -o {}".\
                format(ac, cdsfile, pairscdsfile)
        mm.add((ac, cdsfile), pairscdsfile, cmd)
        ksfile = wd + ".ks"
        cmd = "python -m jcvi.apps.ks calc {} -o {} --workdir {}".\
                format(pairscdsfile, ksfile, wd)
        mm.add(pairscdsfile, ksfile, cmd)
    mm.write()
713481
def chain(args):
    
    p = OptionParser(chain.__doc__)
    p.add_option("--dist", default=100000, help="Chaining distance")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    bedfile, = args
    cmd = "sort -k4,4 -k1,1 -k2,2n -k3,3n {0} -o {0}".format(bedfile)
    sh(cmd)
    bed = Bed(bedfile, sorted=False)
    newbed = Bed()
    for accn, bb in groupby(bed, key=lambda x: x.accn):
        bb = list(bb)
        g = Grouper()
        for a in bb:
            g.join(a)
        for a, b in pairwise(bb):
            if a.seqid == b.seqid and b.start - a.end < opts.dist:
                g.join(a, b)
        data = []
        for p in g:
            seqid = p[0].seqid
            start = min(x.start for x in p)
            end = max(x.end for x in p)
            score = sum(x.span for x in p)
            data.append((seqid, start - 1, end, accn, score))

        d = max(data, key=lambda x: x[-1])
        newbed.append(BedLine("\t".join(str(x) for x in d)))

    newbed.print_to_file(opts.outfile, sorted=True)