%prog alignextend bedpefile ref.fasta

    Similar idea to alignextend, using mates from BEDPE and FASTA ref. See AMOS
    script here:

    https://github.com/nathanhaigh/amos/blob/master/src/Experimental/alignextend.pl
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)
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)))
864079
def bed(args):
    
    p = OptionParser(bed.__doc__)
    p.add_option("-o", dest="output", default="stdout",
            help="Output file name [default: %default]")
    p.add_option("--cutoff", dest="cutoff", default=10, type="int",
            help="Minimum read depth to report intervals [default: %default]")
    opts, args = p.parse_args(args)

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

    binfile, fastafile = args
    fw = must_open(opts.output, "w")
    cutoff = opts.cutoff
    assert cutoff >= 0, "Need non-negative cutoff"

    b = BinFile(binfile)
    ar = b.array

    fastasize, sizes, offsets = get_offsets(fastafile)
    s = Sizes(fastafile)
    for ctg, ctglen in s.iter_sizes():
        offset = offsets[ctg]
        subarray = ar[offset:offset + ctglen]
        key = lambda x: x[1] >= cutoff
        for tf, array_elements in groupby(enumerate(subarray), key=key):
            array_elements = list(array_elements)
            if not tf:
                continue

            # 0-based system => 1-based system
            start = array_elements[0][0] + 1
            end = array_elements[-1][0] + 1

            mean_depth = sum([x[1] for x in array_elements]) / \
                len(array_elements)
            mean_depth = int(mean_depth)

            name = "na"
            print("\t".join(str(x) for x in (ctg, \
                    start - 1, end, name, mean_depth)), file=fw)
834742
def bed(args):
    
    p = OptionParser(bed.__doc__)
    opts, args = p.parse_args(args)

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

    frgscffile, = args
    bedfile = frgscffile.rsplit(".", 1)[0] + ".bed"
    fw = open(bedfile, "w")

    fp = open(frgscffile)
    for row in fp:
        f = FrgScfLine(row)
        print(f.bedline, file=fw)

    logging.debug("File written to `{0}`.".format(bedfile))

    return bedfile
281132
def lobstrindex(args):
    
    p = OptionParser(lobstrindex.__doc__)
    p.add_option("--notreds", default=False, action="store_true",
                 help="Remove TREDs from the bed file")
    p.set_home("lobstr")
    opts, args = p.parse_args(args)

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

    trfbed, fastafile = args
    pf = fastafile.split(".")[0]
    lhome = opts.lobstr_home
    mkdir(pf)

    if opts.notreds:
        newbedfile = trfbed + ".new"
        newbed = open(newbedfile, "w")
        fp = open(trfbed)
        retained = total = 0
        seen = set()
        for row in fp:
            r = STRLine(row)
            total += 1
            name = r.longname
            if name in seen:
                continue
            seen.add(name)
            print(r, file=newbed)
            retained += 1
        newbed.close()
        logging.debug("Retained: {0}".format(percentage(retained, total)))
    else:
        newbedfile = trfbed

    mm = MakeManager()
    cmd = "python {0}/scripts/lobstr_index.py".format(lhome)
    cmd += " --str {0} --ref {1} --out {2}".format(newbedfile, fastafile, pf)
    mm.add((newbedfile, fastafile), op.join(pf, "lobSTR_ref.fasta.rsa"), cmd)

    tabfile = "{0}/index.tab".format(pf)
    cmd = "python {0}/scripts/GetSTRInfo.py".format(lhome)
    cmd += " {0} {1} > {2}".format(newbedfile, fastafile, tabfile)
    mm.add((newbedfile, fastafile), tabfile, cmd)

    infofile = "{0}/index.info".format(pf)
    cmd = "cp {0} {1}".format(newbedfile, infofile)
    mm.add(trfbed, infofile, cmd)
    mm.write()
1250682
def _align_abi_input(arg_abi, arg):
    
    tuple_parts = get_tuple_type_str_parts(arg_abi['type'])

    if tuple_parts is None:
        # Arg is non-tuple.  Just return value.
        return arg

    tuple_prefix, tuple_dims = tuple_parts
    if tuple_dims is None:
        # Arg is non-list tuple.  Each sub arg in `arg` will be aligned
        # according to its corresponding abi.
        sub_abis = arg_abi['components']
    else:
        # Arg is list tuple.  A non-list version of its abi will be used to
        # align each element in `arg`.
        new_abi = copy.copy(arg_abi)
        new_abi['type'] = tuple_prefix

        sub_abis = itertools.repeat(new_abi)

    if isinstance(arg, abc.Mapping):
        # Arg is mapping.  Align values according to abi order.
        aligned_arg = tuple(arg[abi['name']] for abi in sub_abis)
    else:
        aligned_arg = arg

    if not is_list_like(aligned_arg):
        raise TypeError(
            'Expected non-string sequence for "{}" component type: got {}'.format(
                arg_abi['type'],
                aligned_arg,
            ),
        )

    return type(aligned_arg)(
        _align_abi_input(sub_abi, sub_arg)
        for sub_abi, sub_arg in zip(sub_abis, aligned_arg)
    )
231663
def bed(args):
    '''
    
    '''
    from jcvi.utils.cbook import gene_name

    p = OptionParser(bed.__doc__)
    p.add_option("--type", dest="type", default="gene",
            help="Feature type to extract, use comma for multiple [default: %default]")
    p.add_option("--key", default="ID", help="Key in the attributes to extract")
    p.add_option("--source",
            help="Source to extract from, use comma for multiple [default: %default]")
    p.add_option("--span", default=False, action="store_true",
            help="Use feature span in the score column")
    p.add_option("--score_attrib", dest="score_attrib", default=False,
            help="Attribute whose value is to be used as score in `bedline` [default: %default]")
    p.add_option("--append_source", default=False, action="store_true",
            help="Append GFF source name to extracted key value")
    p.add_option("--append_ftype", default=False, action="store_true",
            help="Append GFF feature type to extracted key value")
    p.add_option("--nosort", default=False, action="store_true",
            help="Do not sort the output bed file [default: %default]")
    p.set_stripnames(default=False)
    p.set_outfile()

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

    gffile, = args
    key = opts.key or None
    type = opts.type or set()
    source = opts.source or set()
    strip_names = opts.strip_names
    span = opts.span

    if opts.type:
        type = set(x.strip() for x in opts.type.split(","))
    if opts.source:
        source = set(x.strip() for x in opts.source.split(","))

    gff = Gff(gffile, key=key, append_source=opts.append_source, \
        append_ftype=opts.append_ftype, score_attrib=opts.score_attrib)
    b = Bed()

    for g in gff:
        if type and g.type not in type:
            continue
        if source and g.source not in source:
            continue

        bl = g.bedline
        if strip_names:
            bl.accn = gene_name(bl.accn)
        if span:
            bl.score = bl.span
        b.append(bl)

    sorted = not opts.nosort
    b.print_to_file(opts.outfile, sorted=sorted)
    logging.debug("Extracted {0} features (type={1} id={2})".\
                    format(len(b), ",".join(type), key))
160249
def fill(args):
    
    p = OptionParser(fill.__doc__)
    p.add_option("--extend", default=2000, type="int",
                 help="Extend seq flanking the gaps [default: %default]")
    opts, args = p.parse_args(args)

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

    gapsbed, badfasta = args
    Ext = opts.extend

    gapdist = 2 * Ext + 1  # This is to prevent to replacement ranges intersect
    gapsbed = mergeBed(gapsbed, d=gapdist, nms=True)

    bed = Bed(gapsbed)
    sizes = Sizes(badfasta).mapping
    pf = gapsbed.rsplit(".", 1)[0]
    extbed = pf + ".ext.bed"
    fw = open(extbed, "w")
    for b in bed:
        gapname = b.accn
        start, end = max(0, b.start - Ext - 1), b.start - 1
        print("\t".join(str(x) for x in \
                             (b.seqid, start, end, gapname + "L")), file=fw)
        start, end = b.end, min(sizes[b.seqid], b.end + Ext)
        print("\t".join(str(x) for x in \
                             (b.seqid, start, end, gapname + "R")), file=fw)
    fw.close()

    fastaFromBed(extbed, badfasta, name=True)
1236091
function alignLeft(str, len, ch) {
    return str + new Array(len - str.length + 1).join(ch || " ");
}
726691
def align(args):
    
    valid_modes = ("bwasw", "aln", "mem")
    p = OptionParser(align.__doc__)
    p.add_option("--mode", default="mem", choices=valid_modes, help="BWA mode")
    set_align_options(p)
    p.set_sam_options()

    opts, args = p.parse_args(args)
    mode = opts.mode
    nargs = len(args)

    if nargs not in (2, 3):
        sys.exit(not p.print_help())

    tag = "bwa-{0}: ".format(mode)
    c = mem
    if nargs == 2:
        tag += "Single-end alignment"
        if mode == "bwasw":
            c = bwasw
        elif mode == "aln":
            c = samse
    else:
        assert mode != "bwasw", "Cannot use --bwasw with paired-end mode"
        tag += "Paired-end alignment"
        if mode == "aln":
            c = sampe

    logging.debug(tag)
    cmd, samfile = c(args, opts)
    if cmd:
        cmd = output_bam(cmd, samfile)

    bam = opts.bam
    unmapped = opts.unmapped

    sh(cmd)
    if unmapped:
        dbfile, readfile = args[:2]
        mopts = [samfile, "--unmapped"]
        if not bam:
            mopts += ["--sam"]
        mapped(mopts)
        FileShredder([samfile])

    return samfile, None
174085
def tophat_align(fastq_file, pair_file, ref_file, out_base, align_dir, data,
                 names=None):
    
    config = data["config"]
    options = get_in(config, ("resources", "tophat", "options"), {})
    options = _set_fusion_mode(options, data)
    options = _set_quality_flag(options, data)
    options = _set_transcriptome_option(options, data, ref_file)
    options = _set_cores(options, config)
    options = _set_rg_options(options, names)
    options = _set_stranded_flag(options, config)

    ref_file, runner = _determine_aligner_and_reference(ref_file, data)

    # fusion search does not work properly with Bowtie2
    if options.get("fusion-search", False):
        ref_file = ref_file.replace("/bowtie2", "/bowtie")

    if _tophat_major_version(config) == 1:
        raise NotImplementedError("Tophat versions < 2.0 are not supported, please "
                                  "download the newest version of Tophat here: "
                                  "http://tophat.cbcb.umd.edu")

    if _ref_version(ref_file) == 1 or options.get("fusion-search", False):
        options["bowtie1"] = True

    out_dir = os.path.join(align_dir, "%s_tophat" % out_base)
    final_out = os.path.join(out_dir, "{0}.bam".format(names["sample"]))
    if file_exists(final_out):
        return final_out

    out_file = os.path.join(out_dir, "accepted_hits.bam")
    unmapped = os.path.join(out_dir, "unmapped.bam")
    files = [ref_file, fastq_file]
    if not file_exists(out_file):
        with file_transaction(config, out_dir) as tx_out_dir:
            safe_makedir(tx_out_dir)
            if pair_file and not options.get("mate-inner-dist", None):
                d, d_stdev = _estimate_paired_innerdist(fastq_file, pair_file,
                                                        ref_file, out_base,
                                                        tx_out_dir, data)
                options["mate-inner-dist"] = d
                options["mate-std-dev"] = d_stdev
                files.append(pair_file)
            options["output-dir"] = tx_out_dir
            options["no-coverage-search"] = True
            options["no-mixed"] = True
            cmd = [utils.get_program_python("tophat"), config_utils.get_program("tophat", config)]
            for k, v in options.items():
                if v is True:
                    cmd.append("--%s" % k)
                else:
                    assert not isinstance(v, bool)
                    cmd.append("--%s=%s" % (k, v))
            # tophat requires options before arguments, otherwise it silently ignores them
            cmd += files
            do.run(cmd, "Running Tophat on %s and %s." % (fastq_file, pair_file))
    if pair_file and _has_alignments(out_file):
        fixed = _fix_mates(out_file, os.path.join(out_dir, "%s-align.bam" % out_base),
                           ref_file, config)
    else:
        fixed = out_file
    fixed_unmapped = _fix_unmapped(fixed, unmapped, data)
    fixed = merge_unmapped(fixed, fixed_unmapped, config)
    fixed = _add_rg(fixed, config, names)
    fixed = bam.sort(fixed, config)
    picard = broad.runner_from_path("picard", config)
    # set the contig order to match the reference file so GATK works
    fixed = picard.run_fn("picard_reorder", fixed, data["sam_ref"],
                          os.path.splitext(fixed)[0] + ".picard.bam")
    fixed = fix_insert_size(fixed, config)
    if not file_exists(final_out):
        symlink_plus(fixed, final_out)
    return final_out
13945
def write_equalwidth_bedfile(bedfile, width, outfile):
    """"""

    BUFSIZE = 10000
    f = open(bedfile)
    out = open(outfile, "w")
    lines = f.readlines(BUFSIZE)
    line_count = 0
    while lines:
        for line in lines:
            line_count += 1
            if not line.startswith("#") and not line.startswith("track") and not line.startswith("browser"):
                vals = line.strip().split("\t")
                try:
                    start, end = int(vals[1]), int(vals[2])
                except ValueError:
                    print("Error on line %s while reading %s. Is the file in BED or WIG format?" % (line_count, bedfile))
                    sys.exit(1)

                start = (start + end) // 2 - (width // 2)
                # This shifts the center, but ensures the width is identical... maybe not ideal
                if start < 0:
                    start = 0
                end = start + width
                # Keep all the other information in the bedfile if it's there
                if len(vals) > 3:
                    out.write("%s\t%s\t%s\t%s\n" % (vals[0], start, end, "\t".join(vals[3:])))
                else:
                    out.write("%s\t%s\t%s\n" % (vals[0], start, end))
        lines = f.readlines(BUFSIZE)
    
    out.close()
    f.close()
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()
516396
def cyntenator(args):
    
    p = OptionParser(cyntenator.__doc__)
    opts, args = p.parse_args(args)

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

    lastfile = args[0]
    fp = open(lastfile)
    filteredlastfile = lastfile + ".blast"
    fw = open(filteredlastfile, "w")
    for row in fp:
        b = BlastLine(row)
        if b.query == b.subject:
            continue
        print("\t".join((b.query, b.subject, str(b.score))), file=fw)
    fw.close()

    bedfiles = args[1:]
    fp = open(lastfile)
    b = BlastLine(next(fp))
    subject = b.subject
    txtfiles = []
    for bedfile in bedfiles:
        order = Bed(bedfile).order
        if subject in order:
            db = op.basename(bedfile).split(".")[0][:20]
            logging.debug("Found db: {0}".format(db))
        txtfile = write_txt(bedfile)
        txtfiles.append(txtfile)

    db += ".txt"
    mm = MakeManager()
    for txtfile in txtfiles:
        outfile = txtfile + ".alignment"
        cmd = 'cyntenator -t "({0} {1})" -h blast {2} > {3}'\
                .format(txtfile, db, filteredlastfile, outfile)
        mm.add((txtfile, db, filteredlastfile), outfile, cmd)
    mm.write()
527735
def align(self, out_path=None):
        """"""
        if out_path is None: out_path = self.prefix_path + '.aln'
        sh.muscle38("-in", self.path, "-out", out_path)
        return AlignedFASTA(out_path)
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)))
15114
def bed(args):
    
    from jcvi.formats.blast import BlastLine
    p = OptionParser(bed.__doc__)

    opts, args = p.parse_args(args)

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

    btabfile, = args
    btab = Btab(btabfile)
    for b in btab:
        Bline = BlastLine(b.blastline)
        print(Bline.bedline)