%prog condense OM.bed

    Merge split alignments in OM bed.
160499
def condense(args):
    
    from itertools import groupby
    from jcvi.assembly.patch import merge_ranges

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

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

    bedfile, = args
    bed = Bed(bedfile, sorted=False)
    key = lambda x: (x.seqid, x.start, x.end)
    for k, sb in groupby(bed, key=key):
        sb = list(sb)
        b = sb[0]
        chr, start, end, strand = merge_ranges(sb)

        id = "{0}:{1}-{2}".format(chr, start, end)
        b.accn = id
        print(b)
1332559
def mergebydepth(args):
    
    p = OptionParser(mergebydepth.__doc__)
    p.add_option("--mindepth", default=3, type="int",
                 help="Minimum depth required")
    opts, args = p.parse_args(args)

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

    bedfile, fastafile = args
    mindepth = opts.mindepth
    bedgraph = make_bedgraph(bedfile)

    bedgraphfiltered = bedgraph + ".d{0}".format(mindepth)
    if need_update(bedgraph, bedgraphfiltered):
        filter([bedgraph, "--minaccn={0}".format(mindepth),
                "--outfile={0}".format(bedgraphfiltered)])

    merged = bedgraphfiltered + ".merge.fasta"
    if need_update(bedgraphfiltered, merged):
        mergeBed(bedgraphfiltered, sorted=True)
1349689
def join(l, conj=CONJUNCTION, im_a_moron=MORON_MODE, separator=COMMA):
    """"""

    collector = []
    left = len(l)
    separator = separator + SPACE
    conj = conj + SPACE

    for _l in l[:]:

        left += -1

        collector.append(_l)
        if left == 1:
            if len(l) == 2 or im_a_moron:
                collector.append(SPACE)
            else:
                collector.append(separator)

            collector.append(conj)

        elif left is not 0:
            collector.append(separator)

    return unicode(str().join(collector))
35024
def union(*argv):
        
        res = OrderedSet()
        for ods in argv:
            res = res | ods
        return res
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))
1285643
def merge_coords(objs, compat='minimal', join='outer', priority_arg=None,
                 indexes=None):
    
    _assert_compat_valid(compat)
    coerced = coerce_pandas_values(objs)
    aligned = deep_align(coerced, join=join, copy=False, indexes=indexes)
    expanded = expand_variable_dicts(aligned)
    priority_vars = _get_priority_vars(aligned, priority_arg, compat=compat)
    variables = merge_variables(expanded, priority_vars, compat=compat)
    assert_unique_multiindex_level_names(variables)

    return variables
217588
def bed(args):
    
    p = OptionParser(bed.__doc__)
    p.add_option("--switch", default=False, action="store_true",
                 help="Switch reference and aligned map elements [default: %default]")
    opts, args = p.parse_args(args)

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

    mapout, = args
    pf = mapout.split(".")[0]
    mapbed = pf + ".bed"
    bm = BinMap(mapout)
    bm.print_to_bed(mapbed, switch=opts.switch)

    return mapbed
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"])
1134276
def split_diff(old, new):
        
        return map(lambda l: l.rstrip(),
                   icdiff.ConsoleDiff(cols=COLUMNS).make_table(old.splitlines(), new.splitlines()))
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)
221062
def bed(args):
    
    p = OptionParser(main.__doc__)
    opts, args = p.parse_args(args)

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

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

    for rec in c.iter_records():
        for r in rec.reads:
            print(r.bedline, file=fw)

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

    return bedfile
888641
def pasteprepare(args):
    
    p = OptionParser(pasteprepare.__doc__)
    p.add_option("--flank", default=5000, type="int",
                 help="Get the seq of size on two ends [default: %default]")
    opts, args = p.parse_args(args)

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

    goodfasta, = args
    flank = opts.flank
    pf = goodfasta.rsplit(".", 1)[0]
    extbed = pf + ".ext.bed"

    sizes = Sizes(goodfasta)
    fw = open(extbed, "w")
    for bac, size in sizes.iter_sizes():
        print("\t".join(str(x) for x in \
                               (bac, 0, min(flank, size), bac + "L")), file=fw)
        print("\t".join(str(x) for x in \
                               (bac, max(size - flank, 0), size, bac + "R")), file=fw)
    fw.close()

    fastaFromBed(extbed, goodfasta, name=True)
666909
def alignextend(args):
    
    choices = "prepare,align,filter,rmdup,genreads".split(",")
    p = OptionParser(alignextend.__doc__)
    p.add_option("--nosuffix", default=False, action="store_true",
                 help="Do not add /1/2 suffix to the read [default: %default]")
    p.add_option("--rc", default=False, action="store_true",
                 help="Reverse complement the reads before alignment")
    p.add_option("--len", default=100, type="int",
                 help="Extend to this length")
    p.add_option("--stage", default="prepare", choices=choices,
                 help="Start from certain stage")
    p.add_option("--dup", default=10, type="int",
                 help="Filter duplicates with coordinates within this distance")
    p.add_option("--maxdiff", default=1, type="int",
                 help="Maximum number of differences")
    p.set_home("amos")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    ref, r1, r2 = args
    pf = op.basename(r1).split(".")[0]
    cmd = op.join(opts.amos_home, "src/Experimental/alignextend.pl")
    if not opts.nosuffix:
        cmd += " -suffix"
    bwa_idx = "{0}.ref.fa.sa".format(pf)
    if not need_update(ref, bwa_idx):
        cmd += " -noindex"
    cmd += " -threads {0}".format(opts.cpus)
    offset = guessoffset([r1])
    if offset == 64:
        cmd += " -I"
    if opts.rc:
        cmd += " -rc"
    cmd += " -allow -len {0} -dup {1}".format(opts.len, opts.dup)
    cmd += " -min {0} -max {1}".format(2 * opts.len, 20 * opts.len)
    cmd += " -maxdiff {0}".format(opts.maxdiff)
    cmd += " -stage {0}".format(opts.stage)
    cmd += " ".join(("", pf, ref, r1, r2))
    sh(cmd)
549849
def format(args):
    
    sequential_choices = ("replace", "prefix", "suffix")
    p = OptionParser(format.__doc__)
    p.add_option("--pairs", default=False, action="store_true",
            help="Add trailing /1 and /2 for interleaved pairs [default: %default]")
    p.add_option("--sequential", default=None, choices=sequential_choices,
            help="Add sequential IDs [default: %default]")
    p.add_option("--sequentialoffset", default=0, type="int",
            help="Sequential IDs start at [default: %default]")
    p.add_option("--pad0", default=0, type="int",
            help="Pad a few zeros in front of sequential [default: %default]")
    p.add_option("--gb", default=False, action="store_true",
            help="For Genbank ID, get the accession [default: %default]")
    p.add_option("--sep", default=None,
            help="Split description by certain symbol [default: %default]")
    p.add_option("--index", default=0, type="int",
            help="Extract i-th field after split with --sep [default: %default]")
    p.add_option("--noversion", default=False, action="store_true",
            help="Remove the gb trailing version [default: %default]")
    p.add_option("--prefix", help="Prepend prefix to sequence ID")
    p.add_option("--suffix", help="Append suffix to sequence ID")
    p.add_option("--template", default=False, action="store_true",
            help="Extract `template=aaa dir=x library=m` to `m-aaa/x` [default: %default]")
    p.add_option("--switch", help="Switch ID from two-column file [default: %default]")
    p.add_option("--annotation", help="Add functional annotation from "
                        "two-column file ('ID <--> Annotation') [default: %default]")
    p.add_option("--ids", help="Generate ID conversion table [default: %default]")
    p.add_option("--upper", default=False, action="store_true",
            help="Convert sequence to upper case [default: %default]")
    p.add_option("--nodesc", default=False, action="store_true",
            help="Remove description after identifier")
    opts, args = p.parse_args(args)

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

    infasta, outfasta = args
    gb = opts.gb
    pairs = opts.pairs
    prefix = opts.prefix
    suffix = opts.suffix
    noversion = opts.noversion
    sequential = opts.sequential
    sequentialoffset = opts.sequentialoffset
    sep = opts.sep
    idx = opts.index
    mapfile = opts.switch
    annotfile = opts.annotation
    desc = not opts.nodesc
    idsfile = opts.ids
    idsfile = open(idsfile, "w") if idsfile else None
    upper = opts.upper

    if mapfile:
        mapping = DictFile(mapfile, delimiter="\t")
    if annotfile:
        annotation = DictFile(annotfile, delimiter="\t")

    fp = SeqIO.parse(must_open(infasta), "fasta")
    fw = must_open(outfasta, "w")
    for i, rec in enumerate(fp):
        origid = rec.id
        description = rec.description.replace(origid, "").strip()
        if sep:
            rec.id = rec.description.split(sep)[idx].strip()
        if gb:
            # gi|262233616|gb|GU123895.1| Coffea arabica clone BAC
            atoms = rec.id.split("|")
            if len(atoms) >= 3:
                rec.id = atoms[3]
            elif len(atoms) == 2:
                rec.id = atoms[1]
        if pairs:
            id = "/1" if (i % 2 == 0) else "/2"
            rec.id += id
        if noversion:
            rec.id = rec.id.rsplit(".", 1)[0]
        if sequential:
            rec.id = "{0:0{1}d}".format(sequentialoffset, opts.pad0)
            if sequential == "prefix":
                rec.id = "{0}-{1}".format(rec.id, origid)
            elif sequential == "suffix":
                rec.id = "{0}-{1}".format(origid, rec.id)
            sequentialoffset += 1
        if opts.template:
            template, dir, lib = [x.split("=")[-1] for x in
                    rec.description.split()[1:4]]
            rec.id = "{0}-{1}/{2}".format(lib, template, dir)
        if mapfile:
            if origid in mapping:
                rec.id = mapping[origid]
            else:
                logging.error("{0} not found in `{1}`. ID unchanged.".\
                        format(origid, mapfile))
        if prefix:
            rec.id = prefix + rec.id
        if suffix:
            rec.id += suffix
        if annotfile:
            rec.description = annotation.get(origid, "") if not mapfile \
                    else annotation.get(rec.id, "")
        else:
            rec.description = description if desc else ""
        if idsfile:
            print("\t".join((origid, rec.id)), file=idsfile)
        if upper:
            rec.seq = rec.seq.upper()

        SeqIO.write(rec, fw, "fasta")

    if idsfile:
        logging.debug("Conversion table written to `{0}`.".\
                      format(idsfile.name))
        idsfile.close()
490619
def get_splits(split_bed, gff_file, stype, key):
    
    bed_file = get_bed_file(gff_file, stype, key)
    cmd = "intersectBed -a {0} -b {1} -wao".format(split_bed, bed_file)
    cmd += " | cut -f4,10"
    p = popen(cmd)
    splits = defaultdict(set)
    for row in p:
        a, b = row.split()
        splits[a].add(b)

    return splits
217924
def merge_split_alignments(samples, run_parallel):
    
    ready = []
    file_key = "work_bam"
    to_merge = collections.defaultdict(list)
    for data in (xs[0] for xs in samples):
        if data.get("combine"):
            out_key = tz.get_in(["combine", file_key, "out"], data)
            if not out_key:
                out_key = data["rgnames"]["lane"]
            to_merge[out_key].append(data)
        else:
            ready.append([data])
    ready_merge = []
    hla_merges = []
    for mgroup in to_merge.values():
        cur_data = mgroup[0]
        del cur_data["align_split"]
        for x in mgroup[1:]:
            cur_data["combine"][file_key]["extras"].append(x[file_key])
        ready_merge.append([cur_data])
        cur_hla = None
        for d in mgroup:
            hla_files = tz.get_in(["hla", "fastq"], d)
            if hla_files:
                if not cur_hla:
                    cur_hla = {"rgnames": {"sample": dd.get_sample_name(cur_data)},
                               "config": cur_data["config"], "dirs": cur_data["dirs"],
                               "hla": {"fastq": []}}
                cur_hla["hla"]["fastq"].append(hla_files)
        if cur_hla:
            hla_merges.append([cur_hla])
    if not tz.get_in(["config", "algorithm", "kraken"], data):
        # kraken requires fasta filenames from data['files'] as input.
        # We don't want to remove those files if kraken qc is required.
        _save_fastq_space(samples)
    merged = run_parallel("delayed_bam_merge", ready_merge)
    hla_merge_raw = run_parallel("merge_split_alignments", hla_merges)
    hla_merges = {}
    for hla_merge in [x[0] for x in hla_merge_raw]:
        hla_merges[dd.get_sample_name(hla_merge)] = tz.get_in(["hla", "fastq"], hla_merge)

    # Add stable 'align_bam' target to use for retrieving raw alignment
    out = []
    for data in [x[0] for x in merged + ready]:
        if data.get("work_bam"):
            data["align_bam"] = data["work_bam"]
        if dd.get_sample_name(data) in hla_merges:
            data["hla"]["fastq"] = hla_merges[dd.get_sample_name(data)]
        else:
            hla_files = glob.glob(os.path.join(dd.get_work_dir(data), "align",
                                               dd.get_sample_name(data), "hla", "*.fq"))
            if hla_files:
                data["hla"]["fastq"] = hla_files
        out.append([data])
    return out
19449
def sort(args):
    
    p = OptionParser(sort.__doc__)
    p.add_option("-i", "--inplace", dest="inplace",
            default=False, action="store_true",
            help="Sort bed file in place [default: %default]")
    p.add_option("-u", dest="unique",
            default=False, action="store_true",
            help="Uniqify the bed file")
    p.add_option("--accn", default=False, action="store_true",
            help="Sort based on the accessions [default: %default]")
    p.set_outfile(outfile=None)
    p.set_tmpdir()
    opts, args = p.parse_args(args)

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

    bedfile, = args
    inplace = opts.inplace

    if not inplace and ".sorted." in bedfile:
        return bedfile

    sortedbed = opts.outfile
    if inplace:
        sortedbed = bedfile
    elif opts.outfile is None:
        pf, sf = op.basename(bedfile).rsplit(".", 1)
        sortedbed = pf + ".sorted." + sf

    sortopt = "-k1,1 -k2,2n -k3,3n -k4,4" if not opts.accn else \
              "-k4,4 -k1,1 -k2,2n -k3,3n"
    cmd = "sort"
    if opts.tmpdir:
        cmd += " -T {0}".format(opts.tmpdir)
    if opts.unique:
        cmd += " -u"
    cmd += " {0} {1} -o {2}".format(sortopt, bedfile, sortedbed)

    if inplace or need_update(bedfile, sortedbed):
        sh(cmd)

    return sortedbed