%prog bed fastafile kmer.dump.txt

    Map kmers on FASTA.
766645
def bed(args):
    
    from jcvi.formats.fasta import rc, parse_fasta

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

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

    fastafile, dumpfile = args
    fp = open(dumpfile)
    KMERS = set()
    for row in fp:
        kmer = row.split()[0]
        kmer_rc = rc(kmer)
        KMERS.add(kmer)
        KMERS.add(kmer_rc)

    K = len(kmer)
    logging.debug("Imported {} {}-mers".format(len(KMERS), K))

    for name, seq in parse_fasta(fastafile):
        name = name.split()[0]
        for i in range(len(seq) - K):
            if i % 5000000 == 0:
                print("{}:{}".format(name, i), file=sys.stderr)
            kmer = seq[i: i + K]
            if kmer in KMERS:
                print("\t".join(str(x) for x in (name, i, i + K, kmer)))
649387
def frombed(args):
    
    p = OptionParser(frombed.__doc__)
    p.add_option("--gapsize", default=100, type="int",
                 help="Insert gaps of size [default: %default]")
    opts, args = p.parse_args(args)

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

    bedfile, = args
    gapsize = opts.gapsize
    agpfile = bedfile.replace(".bed", ".agp")
    fw = open(agpfile, "w")

    bed = Bed(bedfile, sorted=False)
    for object, beds in groupby(bed, key=lambda x: x.accn):
        beds = list(beds)
        for i, b in enumerate(beds):
            if gapsize and i != 0:
                print("\t".join(str(x) for x in \
                    (object, 0, 0, 0, "U", \
                     gapsize, "scaffold", "yes", "map")), file=fw)

            print("\t".join(str(x) for x in \
                    (object, 0, 0, 0, "W", \
                     b.seqid, b.start, b.end, b.strand)), file=fw)

    fw.close()

    # Reindex
    return reindex([agpfile, "--inplace"])
296192
def clean(args):
    
    p = OptionParser(clean.__doc__)
    p.add_option("--fancy", default=False, action="store_true",
                 help="Pretty print the sequence [default: %default]")
    p.add_option("--canonical", default=False, action="store_true",
                 help="Use only acgtnACGTN [default: %default]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    fastafile, = args
    fw = must_open(opts.outfile, "w")
    if opts.fancy:
        for header, seq in iter_clean_fasta(fastafile):
            print(">" + header, file=fw)
            fancyprint(fw, seq)

        return 0

    iterator = iter_canonical_fasta if opts.canonical else iter_clean_fasta

    for header, seq in iterator(fastafile):
        seq = Seq(seq)
        s = SeqRecord(seq, id=header, description="")
        SeqIO.write([s], fw, "fasta")
532517
def filer(args):
        
        # List to store all the metadata objects
        samples = list()
        # Find all the sequence files in the path
        fastas = sorted(glob(os.path.join(args.sequencepath, '*.fa*')))
        for fasta in fastas:
            # Create a metadata object for each sample
            metadata = MetadataObject()
            # Populate the metadata object with the required attributes
            metadata.name = os.path.splitext(os.path.basename(fasta))[0]
            metadata.general = GenObject()
            metadata.commands = GenObject()
            metadata.general.bestassemblyfile = fasta
            metadata.general.outputdirectory = os.path.join(args.sequencepath, metadata.name)
            samples.append(metadata)
        return samples
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)))
345429
def write_prepareref_fasta_file(outfile, gene_coords, genes_need_upstream, genes_non_upstream, upstream_before=100, upstream_after=100):
    ''''''
    tmp_dict = {}
    fasta_in = os.path.join(data_dir, 'NC_000962.3.fa.gz')
    pyfastaq.tasks.file_to_dict(fasta_in, tmp_dict)
    ref_seq = tmp_dict['NC_000962.3']

    with open(outfile, 'w') as f:
        for gene in genes_non_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start:end+1])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[end:start+1])
                gene_fa.revcomp()

            print(gene_fa, file=f)

        for gene in genes_need_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_before:start + upstream_after])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_after + 1:start + upstream_before + 1])
                gene_fa.revcomp()

            gene_fa.id += '_upstream'
            print(gene_fa, file=f)
991339
def _prep_callable_bed(in_file, work_dir, stats, data):
    
    out_file = os.path.join(work_dir, "%s-merge.bed.gz" % utils.splitext_plus(os.path.basename(in_file))[0])
    gsort = config_utils.get_program("gsort", data)
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(data, out_file) as tx_out_file:
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            cmd = ("{gsort} {in_file} {fai_file} | bedtools merge -i - -d {stats[merge_size]} | "
                   "bgzip -c > {tx_out_file}")
            do.run(cmd.format(**locals()), "Prepare SV callable BED regions")
    return vcfutils.bgzip_and_index(out_file, data["config"])
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)
367255
def clr(args):
    
    p = OptionParser(clr.__doc__)
    p.set_bedpe()
    opts, args = p.parse_args(args)

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

    bedpe, ref = args
    if bedpe.endswith(".bam"):
        bedpefile = bedpe.replace(".bam", ".bedpe")
        if need_update(bedpe, bedpefile):
            cmd = "bamToBed -bedpe -i {0}".format(bedpe)
            sh(cmd, outfile=bedpefile)
        bedpe = bedpefile

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

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

    converted = rmdup + ".converted"
    if need_update(rmdup, converted):
        fp = open(rmdup)
        fw = open(converted, "w")
        for row in fp:
            r = BedpeLine(row)
            print(r.bedline, file=fw)
        fw.close()

    merged = converted + ".merge.bed"
    if need_update(converted, merged):
        mergeBed(converted)
572038
def mcscanx(args):
    
    p = OptionParser(mcscanx.__doc__)
    opts, args = p.parse_args(args)

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

    blastfile = args[0]
    bedfiles = args[1:]
    prefix = "_".join(op.basename(x)[:2] for x in bedfiles)
    symlink(blastfile, prefix + ".blast")
    allbedfile = prefix + ".gff"
    fw = open(allbedfile, "w")
    for i, bedfile in enumerate(bedfiles):
        prefix = chr(ord('A') + i)
        make_gff(bedfile, prefix, fw)
    fw.close()
1002819
def run_coral(clus_obj, out_dir, args):
    
    if not args.bed:
        raise ValueError("This module needs the bed file output from cluster subcmd.")
    workdir = op.abspath(op.join(args.out, 'coral'))
    safe_dirs(workdir)
    bam_in = op.abspath(args.bam)
    bed_in = op.abspath(args.bed)
    reference = op.abspath(args.ref)
    with chdir(workdir):
        bam_clean = coral.prepare_bam(bam_in, bed_in)
        out_dir = op.join(workdir, "regions")
        safe_dirs(out_dir)
        prefix = "seqcluster"
        loci_file = coral.detect_regions(bam_clean, bed_in, out_dir, prefix)
        coral.create_features(bam_clean, loci_file, reference, out_dir)
497162
def blat(args):
    
    p = OptionParser(blat.__doc__)
    p.add_option("--minscore", default=100, type="int",
                 help="Matches minus mismatches gap penalty [default: %default]")
    p.add_option("--minid", default=98, type="int",
                 help="Minimum sequence identity [default: %default]")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    oldfasta, newfasta = args
    twobitfiles = []
    for fastafile in args:
        tbfile = faToTwoBit(fastafile)
        twobitfiles.append(tbfile)

    oldtwobit, newtwobit = twobitfiles
    cmd = "pblat -threads={0}".format(opts.cpus) if which("pblat") else "blat"
    cmd += " {0} {1}".format(oldtwobit, newfasta)
    cmd += " -tileSize=12 -minScore={0} -minIdentity={1} ".\
                format(opts.minscore, opts.minid)
    pslfile = "{0}.{1}.psl".format(*(op.basename(x).split('.')[0] \
                for x in (newfasta, oldfasta)))
    cmd += pslfile
    sh(cmd)
796035
def summary(args):
    
    p = OptionParser(summary.__doc__)
    p.add_option("--sizes", default=False, action="store_true",
                 help="Write .sizes file")
    p.add_option("--all", default=False, action="store_true",
                 help="Write summary stats per seqid")
    opts, args = p.parse_args(args)

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

    bedfile, = args
    bed = Bed(bedfile)
    bs = BedSummary(bed)
    if opts.sizes:
        sizesfile = bedfile + ".sizes"
        fw = open(sizesfile, "w")
        for span, accn in bs.mspans:
            print(span, file=fw)
        fw.close()
        logging.debug("Spans written to `{0}`.".format(sizesfile))
        return bs

    if not opts.all:
        bs.report()
        return bs

    for seqid, subbeds in bed.sub_beds():
        bs = BedSummary(subbeds)
        print("\t".join((seqid, str(bs))))
319485
def tofasta(args):
    
    p = OptionParser(tofasta.__doc__)
    p.add_option("--prefix", default="gbfasta",
            help="prefix of output files [default: %default]")
    filenames, accessions, idfile, opts, args = preparegb(p, args)
    prefix = opts.prefix

    GenBank(filenames=filenames, accessions=accessions, idfile=idfile).\
        write_fasta(output=prefix, individual=opts.individual)

    if opts.individual:
        logging.debug("Output written dir {0}".format(prefix))
    else:
        logging.debug("Output written to {0}.fasta".format(prefix))
11993
def merge(args):
    
    p = OptionParser(merge.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    bedfiles = args
    fw = must_open(opts.outfile, "w")
    for bedfile in bedfiles:
        bed = Bed(bedfile)
        pf = op.basename(bedfile).split(".")[0]
        for b in bed:
            b.seqid = "_".join((pf, b.seqid))
            print(b, file=fw)
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))