%prog blat old.fasta new.fasta

    Generate psl file using blat.
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)
1283883
def _psed(text,
          before,
          after,
          limit,
          flags):
    '''
    
    '''
    atext = text
    if limit:
        limit = re.compile(limit)
        comps = text.split(limit)
        atext = ''.join(comps[1:])

    count = 1
    if 'g' in flags:
        count = 0
        flags = flags.replace('g', '')

    aflags = 0
    for flag in flags:
        aflags |= RE_FLAG_TABLE[flag]

    before = re.compile(before, flags=aflags)
    text = re.sub(before, after, atext, count=count)

    return text
58534
def generate_datafile_old(number_items=1000):
    
    from utils import get_names, generate_dataset
    from pprint import pprint
    filename = "samples.py"
    dataset = generate_dataset(number_items)
    fo = open(filename, "wb")
    fo.write("#!/usr/bin/env python\n")
    fo.write("# -*- coding: utf-8 -*-\n")
    fo.write("#Brainaetic: http://www.thenetplanet.com\n\n")
    fo.write("SAMPLES = ")
    pprint(dataset, fo)
    fo.close()
    print "%s generated with %d samples" % (filename, number_items)
58612
def writeFasta(sequence, sequence_name, output_file):
    
    i=0
    f=open(output_file,'w')
    f.write(">"+str(sequence_name)+"\n")
    while i <= len(sequence):
        f.write(sequence[i:i+60]+"\n")
        i=i+60
    f.close()
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
175312
def pasa(args):
    
    p = OptionParser(pasa.__doc__)
    opts, args = p.parse_args(args)

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

    pasa_db, fastafile = args

    termexons = "pasa.terminal_exons.gff3"
    if need_update(fastafile, termexons):
        cmd = "$ANNOT_DEVEL/PASA2/scripts/pasa_asmbls_to_training_set.dbi"
        cmd += ' -M "{0}:mysql.tigr.org" -p "access:access"'.format(pasa_db)
        cmd += ' -g {0}'.format(fastafile)
        sh(cmd)

        cmd = "$EVM/PasaUtils/retrieve_terminal_CDS_exons.pl"
        cmd += " trainingSetCandidates.fasta trainingSetCandidates.gff"
        sh(cmd, outfile=termexons)

    return termexons
593965
def run_bcl2fastq(run_folder, ss_csv, config):
    
    bc_dir = os.path.join(run_folder, "Data", "Intensities", "BaseCalls")
    output_dir = os.path.join(run_folder, "fastq")

    if not os.path.exists(os.path.join(output_dir, "Makefile")):
        subprocess.check_call(["configureBclToFastq.pl", "--no-eamss",
                               "--input-dir", bc_dir, "--output-dir", output_dir,
                               "--sample-sheet", ss_csv])
    with utils.chdir(output_dir):
        cores = str(utils.get_in(config, ("algorithm", "num_cores"), 1))
        cmd = ["make", "-j", cores]
        if "submit_cmd" in config["process"] and "bcl2fastq_batch" in config["process"]:
            _submit_and_wait(cmd, cores, config, output_dir)
        else:
            subprocess.check_call(cmd)
    return output_dir
1099262
def fastq_convert_pipe_cl(in_file, data):
    
    cmd = _seqtk_fastq_prep_cl(data, in_file)
    if not cmd:
        cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
        cmd = cat_cmd + " " + in_file
    return "<(%s)" % cmd
719576
def postprocess(options):
    

    resdir = options.resdir
    out_file = options.outfile
    tol = options.tol

    print('.. load permutation results')
    file_name = os.path.join(resdir,'perm*','*.res')
    files = glob.glob(file_name)
    LLR0 = []
    for _file in files:
        print(_file)
        LLR0.append(NP.loadtxt(_file,usecols=[6]))
    LLR0 = NP.concatenate(LLR0)

    print('.. fit test statistics')
    t0 = time.time()
    c2m = C2M.Chi2mixture(tol=4e-3)
    c2m.estimate_chi2mixture(LLR0)
    pv0 = c2m.sf(LLR0)
    t1 = time.time()
    print(('finished in %s seconds'%(t1-t0)))

    print('.. export permutation results')
    perm_file = out_file+'.perm'
    RV = NP.array([LLR0,pv0]).T
    NP.savetxt(perm_file,RV,delimiter='\t',fmt='%.6f %.6e')

    print('.. load test results')
    file_name = os.path.join(resdir,'test','*.res')
    files = glob.glob(file_name)
    RV_test = []
    for _file in files:
        print(_file)
        RV_test.append(NP.loadtxt(_file))
    RV_test = NP.concatenate(RV_test)

    print('.. calc pvalues')
    pv = c2m.sf(RV_test[:,-1])[:,NP.newaxis]

    print('.. export test results')
    perm_file = out_file+'.test'
    RV_test = NP.hstack([RV_test,pv])
    NP.savetxt(perm_file,RV_test,delimiter='\t',fmt='%d %d %d %d %d %d %.6e %.6e')

    if options.manhattan:
        manhattan_file = out_file+'.manhattan.jpg'
        plot_manhattan(pv,manhattan_file)
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)))
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)
423875
def fastq(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(fastq.__doc__)
    p.add_option("--outtie", dest="outtie", default=False, action="store_true",
            help="Are these outie reads? [default: %default]")
    p.set_phred()
    p.set_size()

    opts, args = p.parse_args(args)

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

    fastqfiles = [get_abs_path(x) for x in args]
    size = opts.size
    outtie = opts.outtie
    if size > 1000 and (not outtie):
        logging.debug("[warn] long insert size {0} but not outtie".format(size))

    mated = (size != 0)
    libname = op.basename(args[0]).split(".")[0]
    libname = libname.replace("_1_sequence", "")

    frgfile = libname + ".frg"
    mean, sv = get_mean_sv(opts.size)

    cmd = "fastqToCA"
    cmd += " -libraryname {0} ".format(libname)
    fastqs = " ".join("-reads {0}".format(x) for x in fastqfiles)
    if mated:
        assert len(args) in (1, 2), "you need one or two fastq files for mated library"
        fastqs = "-mates {0}".format(",".join(fastqfiles))
        cmd += "-insertsize {0} {1} ".format(mean, sv)
    cmd += fastqs

    offset = int(opts.phred) if opts.phred else guessoffset([fastqfiles[0]])
    illumina = (offset == 64)
    if illumina:
        cmd += " -type illumina"
    if outtie:
        cmd += " -outtie"

    sh(cmd, outfile=frgfile)
1190268
def blast(args):
    
    from jcvi.apps.align import run_megablast

    p = OptionParser(blast.__doc__)
    p.add_option("-n", type="int", default=2,
            help="Take best N hits [default: %default]")
    opts, args = p.parse_args(args)

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

    allfasta, clonename = args
    fastadir = "fasta"
    infile = op.join(fastadir, clonename + ".fasta")
    if not op.exists(infile):
        entrez([clonename, "--skipcheck", "--outdir=" + fastadir])

    outfile = "{0}.{1}.blast".format(clonename, allfasta.split(".")[0])
    run_megablast(infile=infile, outfile=outfile, db=allfasta, \
            pctid=GoodPct, hitlen=GoodOverlap)

    blasts = [BlastLine(x) for x in open(outfile)]
    besthits = []
    for b in blasts:
        if b.query.count("|") >= 3:
            b.query = b.query.split("|")[3]

        if b.subject.count("|") >= 3:
            b.subject = b.subject.split("|")[3]

        b.query = b.query.rsplit(".", 1)[0]
        b.subject = b.subject.rsplit(".", 1)[0]

        if b.query == b.subject:
            continue

        if b.subject not in besthits:
            besthits.append(b.subject)
        if len(besthits) == opts.n:
            break

    for b in besthits:
        overlap([clonename, b, "--dir=" + fastadir])
1085785
def ispcr(args):
    
    from jcvi.utils.iter import grouper

    p = OptionParser(ispcr.__doc__)
    p.add_option("-r", dest="rclip", default=1, type="int",
            help="pair ID is derived from rstrip N chars [default: %default]")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    ispcrfile = fastafile + ".isPcr"
    fw = open(ispcrfile, "w")

    N = opts.rclip
    strip_name = lambda x: x[:-N] if N else str

    npairs = 0
    fastaiter = SeqIO.parse(fastafile, "fasta")
    for a, b in grouper(fastaiter, 2):

        aid, bid = [strip_name(x) for x in (a.id, b.id)]
        assert aid == bid, "Name mismatch {0}".format((aid, bid))

        print("\t".join((aid, str(a.seq), str(b.seq))), file=fw)
        npairs += 1

    fw.close()
    logging.debug("A total of {0} pairs written to `{1}`.".\
                  format(npairs, ispcrfile))
642252
def fasta2bed(fastafile):
    
    dustfasta = fastafile.rsplit(".", 1)[0] + ".dust.fasta"
    for name, seq in parse_fasta(dustfasta):
        for islower, ss in groupby(enumerate(seq), key=lambda x: x[-1].islower()):
            if not islower:
                continue
            ss = list(ss)
            ms, mn = min(ss)
            xs, xn = max(ss)
            print("\t".join(str(x) for x in (name, ms, xs)))
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)
167965
def needle(args):
    
    from jcvi.formats.fasta import Fasta, SeqIO

    p = OptionParser(needle.__doc__)

    opts, args = p.parse_args(args)

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

    manager = mp.Manager()
    results = manager.list()
    needle_pool = mp.Pool(processes=mp.cpu_count())

    pairsfile, apep, bpep = args
    afasta, bfasta = Fasta(apep), Fasta(bpep)
    fp = must_open(pairsfile)
    for i, row in enumerate(fp):
        a, b = row.split()
        a, b = afasta[a], bfasta[b]
        fa, fb = must_open("{0}_{1}_a.fasta".format(pairsfile, i), "w"), \
            must_open("{0}_{1}_b.fasta".format(pairsfile, i), "w")
        SeqIO.write([a], fa, "fasta")
        SeqIO.write([b], fb, "fasta")
        fa.close()
        fb.close()

        needlefile = "{0}_{1}_ab.needle".format(pairsfile, i)
        needle_pool.apply_async(_needle, \
            (fa.name, fb.name, needlefile, a.id, b.id, results))

    needle_pool.close()
    needle_pool.join()

    fp.close()

    scoresfile = "{0}.scores".format(pairsfile.rsplit(".")[0])
    fw = must_open(scoresfile, "w")
    for result in results:
        print(result, file=fw)
    fw.close()