%prog bed blastfile

    Print out bed file based on coordinates in BLAST report. By default, write
    out subject positions. Use --swap to write query positions.
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
1164013
def _bed_to_bed6(orig_file, out_dir):
    
    bed6_file = os.path.join(out_dir, "%s-bed6%s" % os.path.splitext(os.path.basename(orig_file)))
    if not utils.file_exists(bed6_file):
        with open(bed6_file, "w") as out_handle:
            for i, region in enumerate(list(x) for x in pybedtools.BedTool(orig_file)):
                region = [x for x in list(region) if x]
                fillers = [str(i), "1.0", "+"]
                full = region + fillers[:6 - len(region)]
                out_handle.write("\t".join(full) + "\n")
    return bed6_file
59062
def save_bed(cls, query, filename=sys.stdout):
        
        out = _open(filename, 'w')
        for o in query:
            out.write(o.bed() + '\n')
472449
def pairs(args):
    
    p = OptionParser(pairs.__doc__)
    p.set_pairs()
    opts, args = p.parse_args(args)

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

    bedfile, = args

    basename = bedfile.split(".")[0]
    insertsfile = ".".join((basename, "inserts"))
    bedfile = sort([bedfile, "--accn"])

    fp = open(bedfile)
    data = [BedLine(row) for i, row in enumerate(fp) if i < opts.nrows]

    ascii = not opts.pdf
    return bedfile, report_pairs(data, opts.cutoff, opts.mateorientation,
           pairsfile=opts.pairsfile, insertsfile=insertsfile,
           rclip=opts.rclip, ascii=ascii, bins=opts.bins,
           distmode=opts.distmode)
509737
def blastp(seqs, blast_db="nr", e_value="1e-20", max_hits=200,
           working_dir=tempfile.gettempdir(), blast_mat_root=None,
           extra_params={}):
    

    # set up params to use with blastp
    params = {
        # matrix
        "-M":"BLOSUM62",

        # max procs
        "-a":"1",

        # expectation
        "-e":e_value,

        # max seqs to show
        "-b":max_hits,

        # max one line descriptions
        "-v":max_hits,

        # program
        "-p":"blastp"
    }
    params.update(extra_params)

    # blast
    blast_res =  blast_seqs(seqs,
        Blastall,
        blast_mat_root=blast_mat_root,
        blast_db=blast_db,
        params=params,
        add_seq_names=False,
        WorkingDir=working_dir
        )

    # get prot id map
    if blast_res['StdOut']:
        lines = [x for x in blast_res['StdOut']]
        return BlastResult(lines)

    return None
15571
def main(*args):
    '''
    
    '''
    # Parse the arguments...
    (args, opts) = _benson14_parser(*args)
    # help?
    if opts['help']:
        print(info, file=sys.stdout)
        return 1
    # verbose?
    if opts['verbose']:
        def note(s):
            print(s, file=sys.stdout)
            return True
    else:
        def note(s): return False
    # based on format, how do we export?
    sfmt = opts['surf_format'].lower()
    if sfmt in ['curv', 'auto', 'automatic', 'morph']:
        sfmt = 'freesurfer_morph'
        sext = ''
    elif sfmt == 'nifti':
        sext = '.nii.gz'
    elif sfmt in ['mgh', 'mgz', 'nii', 'nii.gz']:
        sext = '.' + sfmt
    else:
        raise ValueError('Unknown surface format: %s' % opts['surf_format'])
    vfmt = opts['vol_format'].lower()
    if vfmt == 'nifti':
        vext = '.nii.gz'
    elif vfmt in ['mgh', 'mgz', 'nii', 'nii.gz']:
        vext = '.' + vfmt
    else:
        raise ValueError('Unknown volume format: %s' % opts['vol_format'])
    # Add the subjects directory, if there is one
    if 'subjects_dir' in opts and opts['subjects_dir'] is not None:
        add_subject_path(opts['subjects_dir'])
    ow = not opts['no_overwrite']
    nse = opts['no_surf_export']
    nve = opts['no_vol_export']
    tr = {'angle': opts['angle_tag'],
          'eccen': opts['eccen_tag'],
          'varea': opts['label_tag'],
          'sigma': opts['sigma_tag']}
    # okay, now go through the subjects...
    for subnm in args:
        note('Processing subject %s:' % subnm)
        sub = subject(subnm)
        note('   - Interpolating template...')
        (lhdat, rhdat) = predict_retinotopy(sub,
                                            template=opts['template'],
                                            registration=opts['registration'])
        # Export surfaces
        if nse:
            note('   - Skipping surface export.')
        else:
            note('   - Exporting surfaces:')
            for (t,dat) in six.iteritems(lhdat):
                flnm = os.path.join(sub.path, 'surf', 'lh.' + tr[t] + sext)
                if ow or not os.path.exists(flnm):
                    note('    - Exporting LH prediction file: %s' % flnm)
                    nyio.save(flnm, dat, format=sfmt)
                else:
                    note('    - Not overwriting existing file: %s' % flnm)
            for (t,dat) in six.iteritems(rhdat):
                flnm = os.path.join(sub.path, 'surf', 'rh.' + tr[t] + sext)
                if ow or not os.path.exists(flnm):
                    note('    - Exporting RH prediction file: %s' % flnm)
                    nyio.save(flnm, dat, format=sfmt)
                else:
                    note('    - Not overwriting existing file: %s' % flnm)
        # Export volumes
        if nve:
            note('   - Skipping volume export.')
        else:
            note('   - Exporting Volumes:')
            for t in lhdat.keys():
                flnm = os.path.join(sub.path, 'mri', tr[t] + vext)
                if ow or not os.path.exists(flnm):
                    note('    - Preparing volume file: %s' % flnm)
                    dtyp = (np.int32 if t == 'varea' else np.float32)
                    vol = sub.cortex_to_image(
                        (lhdat[t], rhdat[t]),
                        method=('nearest' if t == 'varea' else 'linear'),
                        dtype=dtyp)
                    note('    - Exporting volume file: %s' % flnm)
                    nyio.save(flnm, vol, like=sub)
                else:
                    note('    - Not overwriting existing file: %s' % flnm)
        note('   Subject %s finished!' % sub.name)
    return 0
1243475
def location(args):
    
    from jcvi.formats.bed import BedLine
    from jcvi.graphics.histogram import stem_leaf_plot

    p = OptionParser(location.__doc__)
    p.add_option("--dist", default=100, type="int",
                 help="Distance cutoff to call 5` and 3` [default: %default]")
    opts, args = p.parse_args(args)

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

    bedfile, fastafile = args
    dist = opts.dist
    sizes = Sizes(fastafile).mapping
    fp = open(bedfile)
    fiveprime = threeprime = total = 0
    percentages = []
    for row in fp:
        b = BedLine(row)
        pos = b.start
        size = sizes[b.seqid]
        if pos < dist:
            fiveprime += 1
        if size - pos < dist:
            threeprime += 1
        total += 1
        percentages.append(100 * pos / size)

    m = "Five prime (within {0}bp of start codon): {1}\n".format(dist, fiveprime)
    m += "Three prime (within {0}bp of stop codon): {1}\n".format(dist, threeprime)
    m += "Total: {0}".format(total)
    print(m, file=sys.stderr)

    bins = 10
    title = "Locations within the gene [0=Five-prime, 100=Three-prime]"
    stem_leaf_plot(percentages, 0, 100, bins, title=title)
1307243
def blast(args):
    
    p = OptionParser(blast.__doc__)
    opts, args = p.parse_args(args)

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

    deltafile, = args
    blastfile = deltafile.rsplit(".", 1)[0] + ".blast"

    if need_update(deltafile, blastfile):
        coords = Coords(deltafile)
        fw = open(blastfile, "w")
        for c in coords:
            print(c.blastline, file=fw)
354243
def phase(args):
    
    p = OptionParser(phase.__doc__)
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    fw = must_open(opts.outfile, "w")
    for gbfile in args:
        for rec in SeqIO.parse(gbfile, "gb"):
            bac_phase, keywords = get_phase(rec)
            chr, clone = get_clone(rec)
            keyword_field = ";".join(keywords)
            print("\t".join((rec.id, str(bac_phase), keyword_field,
                    chr, clone)), file=fw)
79535
def breakpoint(args):
    
    from jcvi.formats.blast import bed
    from jcvi.utils.range import range_interleave

    p = OptionParser(breakpoint.__doc__)
    p.add_option("--xdist", type="int", default=20,
                 help="xdist (in related genome) cutoff [default: %default]")
    p.add_option("--ydist", type="int", default=200000,
                 help="ydist (in current genome) cutoff [default: %default]")
    p.add_option("-n", type="int", default=5,
                 help="number of markers in a block [default: %default]")
    opts, args = p.parse_args(args)

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

    blastfile, bedfile = args
    order = Bed(bedfile).order
    blastbedfile = bed([blastfile])
    bbed = Bed(blastbedfile)
    for scaffold, bs in bbed.sub_beds():
        blocks = get_blocks(scaffold, bs, order,
                            xdist=opts.xdist, ydist=opts.ydist, N=opts.n)
        sblocks = []
        for block in blocks:
            xx, yy = zip(*block)
            sblocks.append((scaffold, min(yy), max(yy)))
        iblocks = range_interleave(sblocks)
        for ib in iblocks:
            ch, start, end = ib
            print("{0}\t{1}\t{2}".format(ch, start - 1, end))
563867
def best(args):
    
    p = OptionParser(best.__doc__)

    p.add_option("-n", default=1, type="int",
            help="get best N hits [default: %default]")
    p.add_option("--nosort", default=False, action="store_true",
            help="assume BLAST is already sorted [default: %default]")
    p.add_option("--hsps", default=False, action="store_true",
            help="get all HSPs for the best pair [default: %default]")
    p.add_option("--subject", default=False, action="store_true",
            help="get best hit(s) for subject genome instead [default: %default]")
    p.set_tmpdir()
    opts, args = p.parse_args(args)

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

    blastfile, = args
    n = opts.n
    hsps = opts.hsps
    tmpdir = opts.tmpdir
    ref = "query" if not opts.subject else "subject"

    if not opts.nosort:
        sargs = [blastfile]
        if tmpdir:
            sargs += ["-T {0}".format(tmpdir)]
        if ref != "query":
            sargs += ["--refscore"]
        sort(sargs)
    else:
        logging.debug("Assuming sorted BLAST")

    if not opts.subject:
        bestblastfile = blastfile + ".best"
    else:
        bestblastfile = blastfile + ".subject.best"
    fw = open(bestblastfile, "w")

    b = Blast(blastfile)
    for q, bline in b.iter_best_hit(N=n, hsps=hsps, ref=ref):
        print(bline, file=fw)

    return bestblastfile
1209385
def frombed(args):
    
    from jcvi.formats.fasta import Fasta
    from jcvi.formats.bed import Bed
    from jcvi.utils.cbook import fill

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

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

    bedfile, contigfasta, readfasta = args
    prefix = bedfile.rsplit(".", 1)[0]
    contigfile = prefix + ".contig"
    idsfile = prefix + ".ids"

    contigfasta = Fasta(contigfasta)
    readfasta = Fasta(readfasta)

    bed = Bed(bedfile)
    checksum = "00000000 checksum."
    fw_ids = open(idsfile, "w")
    fw = open(contigfile, "w")

    for ctg, reads in bed.sub_beds():
        ctgseq = contigfasta[ctg]
        ctgline = "##{0} {1} {2} bases, {3}".format(\
                ctg, len(reads), len(ctgseq), checksum)

        print(ctg, file=fw_ids)
        print(ctgline, file=fw)
        print(fill(ctgseq.seq), file=fw)

        for b in reads:
            read = b.accn
            strand = b.strand
            readseq = readfasta[read]
            rc = " [RC]" if strand == "-" else ""
            readlen = len(readseq)
            rstart, rend = 1, readlen
            if strand == "-":
                rstart, rend = rend, rstart

            readrange = "{{{0} {1}}}".format(rstart, rend)
            conrange = "<{0} {1}>".format(b.start, b.end)
            readline = "#{0}(0){1} {2} bases, {3} {4} {5}".format(\
                    read, rc, readlen, checksum, readrange, conrange)
            print(readline, file=fw)
            print(fill(readseq.seq), file=fw)

    logging.debug("Mapped contigs written to `{0}`.".format(contigfile))
    logging.debug("Contig IDs written to `{0}`.".format(idsfile))
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))
951276
def __run_blast(blast_command, input_file, *args, **kwargs):
    '''
    
    '''

    # XXX: Eventually, translate results on the fly as requested? Or
    #      just always use our parsed object?
    if 'outfmt' in kwargs:
        raise Exception('Use of the -outfmt option is not supported')

    num_processes = kwargs.get(
        'pb_num_processes', os.sysconf('SC_NPROCESSORS_ONLN'))
    fields = kwargs.get('pb_fields', DEFAULT_HIT_FIELDS)

    blast_args = [blast_command]
    blast_args += ['-outfmt', '7 {}'.format(' '.join(fields))]
    for a in args:
        blast_args += ['-' + a]
    for k, v in kwargs.iteritems():
        if not k.startswith('pb_'):
            blast_args += ['-' + k, str(v)]

    popens = []
    for _ in range(num_processes):
        popens.append(
            subprocess.Popen(
                args=blast_args, stdin=subprocess.PIPE, stdout=subprocess.PIPE,
                stderr=None, close_fds=True))

    try:
        for r in __run_blast_select_loop(input_file, popens, fields):
            yield r
    finally:
        for p in popens:
            if p.poll() is None:
                p.terminate()
            p.wait()
1055863
def blast(args):
    '''
    
    '''
    p = OptionParser(blast.__doc__)
    opts, args = p.parse_args(args)

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

    flist = args

    for f in flist:
        maf_to_blast8(f)
553085
def dump(self, title, coordinates):
        
        print("% 8i" % len(self.symbols), file=self._f)
        print(str(title), file=self._f)
        for symbol, coordinate in zip(self.symbols, coordinates):
            print("% 2s % 12.9f % 12.9f % 12.9f" % ((symbol, ) + tuple(coordinate/self.file_unit)), file=self._f)
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"])