%prog paste flanks.bed flanks_vs_assembly.blast backbone.fasta

    Paste in good sequences in the final assembly.
348343
def paste(args):
    
    from jcvi.formats.bed import uniq

    p = OptionParser(paste.__doc__)
    p.add_option("--maxsize", default=300000, type="int",
            help="Maximum size of patchers to be replaced [default: %default]")
    p.add_option("--prefix", help="Prefix of the new object [default: %default]")
    p.set_rclip(rclip=1)
    opts, args = p.parse_args(args)

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

    pbed, blastfile, bbfasta = args
    maxsize = opts.maxsize  # Max DNA size to replace gap
    order = Bed(pbed).order

    beforebed, afterbed = blast_to_twobeds(blastfile, order, log=True,
                                           rclip=opts.rclip, maxsize=maxsize,
                                           flipbeds=True)
    beforebed = uniq([beforebed])

    afbed = Bed(beforebed)
    bfbed = Bed(afterbed)

    shuffle_twobeds(afbed, bfbed, bbfasta, prefix=opts.prefix)
731919
def paste(**kwargs):
    """"""

    clip.OpenClipboard() 
    d = clip.GetClipboardData(win32con.CF_UNICODETEXT)
    clip.CloseClipboard() 
    return d
132597
def paste(self):
        '''
        
        '''
        try:
            t = pygame.scrap.get(SCRAP_TEXT)
            if t:
                self.insert(t)
                return True
        except:
            # pygame.scrap is experimental, allow for changes
            return False
1324873
def sample_batch_transitions(self, batch_size, forward_steps=1):
        
        batches = [
            self._sample_batch_prioritized(
                segment_tree, batch_size, self.deque.frame_history, forward_steps=forward_steps
            )
            for segment_tree in self.segment_trees
        ]

        probs, idxs, tree_idxs = zip(*batches)

        return np.stack(probs, axis=1).astype(float), np.stack(idxs, axis=1), np.stack(tree_idxs, axis=1)
679174
def blast(args):
    
    p = OptionParser(blast.__doc__)
    p.add_option("--dist", default=100, type="int",
            help="Merge adjacent HSPs separated by [default: %default]")
    p.add_option("--db",
            help="Use a different database rather than UniVec_Core")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    fastaprefix = fastafile.split(".", 1)[0]

    univec = opts.db or download("ftp://ftp.ncbi.nih.gov/pub/UniVec/UniVec_Core")
    uniprefix = univec.split(".", 1)[0]

    fastablast = fastaprefix + ".{0}.blast".format(uniprefix)

    prog = run_megablast if opts.db else run_vecscreen
    prog(infile=fastafile, outfile=fastablast, db=univec, pctid=95, hitlen=50)

    fp = open(fastablast)
    ranges = []
    for row in fp:
        b = BlastLine(row)
        ranges.append((b.query, b.qstart, b.qstop))

    merged_ranges = range_merge(ranges, dist=opts.dist)
    bedfile = fastaprefix + ".{0}.bed".format(uniprefix)
    fw = must_open(bedfile, "w")
    for seqid, start, end in merged_ranges:
        print("\t".join(str(x) for x in (seqid, start - 1, end, uniprefix)), file=fw)

    return bedfile
324821
def _fix_contig_orientation(contigs_fa, ref_fa, outfile, min_id=90, min_length=20, breaklen=200):
        ''''''
        if not os.path.exists(contigs_fa):
            raise Error('Cannot fix orientation of assembly contigs because file not found: ' + contigs_fa)

        tmp_coords = os.path.join(outfile + '.tmp.rename.coords')
        pymummer.nucmer.Runner(
            ref_fa,
            contigs_fa,
            tmp_coords,
            min_id=min_id,
            min_length=min_length,
            breaklen=breaklen,
            maxmatch=True,
        ).run()

        to_revcomp = set()
        not_revcomp = set()
        file_reader = pymummer.coords_file.reader(tmp_coords)
        for hit in file_reader:
            if hit.on_same_strand():
                not_revcomp.add(hit.qry_name)
            else:
                to_revcomp.add(hit.qry_name)

        os.unlink(tmp_coords)
        in_both = to_revcomp.intersection(not_revcomp)

        f = pyfastaq.utils.open_file_write(outfile)
        seq_reader = pyfastaq.sequences.file_reader(contigs_fa)
        for seq in seq_reader:
            if seq.id in to_revcomp and seq.id not in in_both:
                seq.revcomp()
            print(seq, file=f)
        pyfastaq.utils.close(f)

        return in_both
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)))
423050
def overlap(args):
    
    from jcvi.formats.blast import chain_HSPs

    p = OptionParser(overlap.__doc__)
    p.add_option("--dir", default=os.getcwd(),
            help="Download sequences to dir [default: %default]")
    p.add_option("--suffix", default="fasta",
            help="Suffix of the sequence file in dir [default: %default]")
    p.add_option("--qreverse", default=False, action="store_true",
            help="Reverse seq a [default: %default]")
    p.add_option("--nochain", default=False, action="store_true",
            help="Do not chain adjacent HSPs [default: chain HSPs]")
    p.set_align(pctid=GoodPct, hitlen=GoodOverlap, evalue=.01)
    p.set_outfile(outfile=None)
    opts, args = p.parse_args(args)

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

    afasta, bfasta = args
    dir = opts.dir
    chain = not opts.nochain
    suffix = opts.suffix
    evalue = opts.evalue
    pctid = opts.pctid
    hitlen = opts.hitlen
    cutoff = Cutoff(pctid, hitlen)

    # Check first whether it is file or accession name
    if not op.exists(afasta):
        af = op.join(dir, ".".join((afasta, suffix)))
        if not op.exists(af):  # Check to avoid redownload
            entrez([afasta, "--skipcheck", "--outdir=" + dir])
        afasta = af

    if not op.exists(bfasta):
        bf = op.join(dir, ".".join((bfasta, suffix)))
        if not op.exists(bf):
            entrez([bfasta, "--skipcheck", "--outdir=" + dir])
        bfasta = bf

    assert op.exists(afasta) and op.exists(bfasta)

    cmd = "blastn -dust no"
    cmd += " -query {0} -subject {1}".format(afasta, bfasta)
    cmd += " -evalue {0} -outfmt 6 -perc_identity {1}".format(evalue, pctid)

    fp = popen(cmd)
    hsps = fp.readlines()

    hsps = [BlastLine(x) for x in hsps]
    hsps = [x for x in hsps if x.hitlen >= hitlen]
    if chain:
        logging.debug("Chain HSPs in the Blast output.")
        dist = 2 * hitlen  # Distance to chain the HSPs
        hsps = chain_HSPs(hsps, xdist=dist, ydist=dist)

    if len(hsps) == 0:
        print("No match found.", file=sys.stderr)
        return None

    besthsp = hsps[0]

    aid, asize = next(Fasta(afasta).itersizes())
    bid, bsize = next(Fasta(bfasta).itersizes())
    o = Overlap(besthsp, asize, bsize, cutoff, qreverse=opts.qreverse)
    o.print_graphic()

    if opts.outfile:
        fw = must_open(opts.outfile, "w")
        print(str(o), file=fw)
        fw.close()

    return o
1082530
def write(self):
        

        rabs = [x.seqrecord for x in self.gproteins.values() if x.is_rab()]
        return SeqIO.write(rabs, self.tmpfname + '.phase2', 'fasta')
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)
153363
def build_one(self, object, lines, fasta, fw, newagp=None):
        
        components = []

        total_bp = 0
        for line in lines:

            if line.is_gap:
                seq = 'N' * line.gap_length
                if newagp:
                    print(line, file=newagp)
            else:
                seq = fasta.sequence(dict(chr=line.component_id,
                        start=line.component_beg,
                        stop=line.component_end,
                        strand=line.orientation))
                # Check for dangling N's
                if newagp:
                    trimNs(seq, line, newagp)

            components.append(seq)
            total_bp += len(seq)

            if self.validate:
                assert total_bp == line.object_end, \
                        "cumulative base pairs (%d) does not match (%d)" % \
                        (total_bp, line.object_end)

        if not newagp:
            rec = SeqRecord(Seq(''.join(components)), id=object, description="")
            SeqIO.write([rec], fw, "fasta")
            if len(rec) > 1000000:
                logging.debug("Write object %s to `%s`" % (object, fw.name))
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"])
416497
def tidy(args):
    
    p = OptionParser(tidy.__doc__)
    p.add_option("--gapsize", dest="gapsize", default=0, type="int",
            help="Set all gaps to the same size [default: %default]")
    p.add_option("--minlen", dest="minlen", default=100, type="int",
            help="Minimum component size [default: %default]")

    opts, args = p.parse_args(args)

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

    fastafile, = args
    gapsize = opts.gapsize
    minlen = opts.minlen

    tidyfastafile = fastafile.rsplit(".", 1)[0] + ".tidy.fasta"
    fw = must_open(tidyfastafile, "w")

    removed = normalized = 0
    fasta = Fasta(fastafile, lazy=True)
    for name, rec in fasta.iteritems_ordered():
        rec.seq = rec.seq.upper()
        if minlen:
            removed += remove_small_components(rec, minlen)
        trim_terminal_Ns(rec)
        if gapsize:
            normalized += normalize_gaps(rec, gapsize)

        if len(rec) == 0:
            logging.debug("Drop seq {0}".format(rec.id))
            continue
        SeqIO.write([rec], fw, "fasta")

    # Print statistics
    if removed:
        logging.debug("Total discarded bases: {0}".format(removed))
    if normalized:
        logging.debug("Gaps normalized: {0}".format(normalized))

    logging.debug("Tidy FASTA written to `{0}`.".format(tidyfastafile))
    fw.close()

    return tidyfastafile
794895
def executePlateBiasAnalysis(options):
    
    # The plink command
    plinkCommand = ["plink", "--noweb", "--bfile", options.bfile,
                    "--loop-assoc", options.loop_assoc, "--fisher",
                    "--pfilter", str(options.pfilter), "--out", options.out]
    runCommand(plinkCommand)
1013135
def assemble_amplicon_skesa(self):
        
        with progressbar(self.metadata) as bar:
            for sample in bar:
                # Initialise variables
                sample[self.analysistype].skesa_outdir = os.path.join(
                    sample[self.analysistype].outputdir, self.analysistype)
                make_path(sample[self.analysistype].skesa_outdir)
                sample[self.analysistype].assemblyfile = os.path.join(sample[self.analysistype].skesa_outdir,
                                                                      'contigs.fasta')
                # If there are two fastq files
                if len(sample.general.fastqfiles) == 2:
                    # Set the reverse fastq name https://github.com/ncbi/SKESA/issues/7
                    sample[self.analysistype].spadescommand = 'skesa --fastq {fastqfiles} --cores {threads} ' \
                                                              '--use_paired_ends --vector_percent 1 ' \
                                                              '--contigs_out {contigs}'\
                        .format(fastqfiles=','.join(sample.general.fastqfiles),
                                threads=self.cpus,
                                contigs=sample[self.analysistype].assemblyfile)
                # Same as above, but use single read settings for the assembler
                else:
                    sample[self.analysistype].spadescommand = 'skesa --fastq {fastqfiles} --cores {threads} ' \
                                                              '--vector_percent 1 --contigs_out {contigs}'\
                        .format(fastqfiles=','.join(sample.general.fastqfiles),
                                threads=self.cpus,
                                contigs=sample[self.analysistype].assemblyfile)
                if not os.path.isfile(sample[self.analysistype].assemblyfile):
                    # Create and run the sequence call as required
                    run_subprocess(sample[self.analysistype].spadescommand)
                    if not os.path.isfile(sample[self.analysistype].assemblyfile):
                        sample[self.analysistype].assemblyfile = 'NA'
1114815
def strip_inserts(fasta):
    
    for seq in parse_fasta(fasta):
        seq[1] = ''.join([b for b in seq[1] if b == '-' or b.isupper()])
        yield seq