%prog format infasta outfasta

    Reformat FASTA file and also clean up names.
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()
213106
def clean(inputline, strip=True, rem_indents=True, encoding_errors="replace"):
    """"""
    stdout_encoding = get_encoding(sys.stdout)
    inputline = str(inputline)
    if rem_indents:
        inputline = inputline.replace(openindent, "").replace(closeindent, "")
    if strip:
        inputline = inputline.strip()
    return inputline.encode(stdout_encoding, encoding_errors).decode(stdout_encoding)
378253
def run_RNAplfold(input_fasta, tmpdir, W=240, L=160, U=1):
    

    profiles = RNAplfold_PROFILES_EXECUTE
    for i, P in enumerate(profiles):
        print("running {P}_RNAplfold... ({i}/{N})".format(P=P, i=i + 1, N=len(profiles)))

        command = "{bin}/{P}_RNAplfold".format(bin=RNAplfold_BIN_DIR, P=P)
        file_out = "{tmp}/{P}_profile.fa".format(tmp=tmpdir, P=P)
        args = " -W {W} -L {L} -u {U} < {fa} > {file_out}".format(W=W, L=L, U=U, fa=input_fasta, file_out=file_out)

        os.system(command + args)

        # check if the file is empty
        if os.path.getsize(file_out) == 0:
            raise Exception("command wrote an empty file: {0}".format(file_out))
    print("done!")
574342
def _convert_fastq(srafn, outdir, single=False):
    ""
    cmd = "fastq-dump --split-files --gzip {srafn}"
    cmd = "%s %s" % (utils.local_path_export(), cmd)
    sraid = os.path.basename(utils.splitext_plus(srafn)[0])
    if not srafn:
        return None
    if not single:
        out_file = [os.path.join(outdir, "%s_1.fastq.gz" % sraid),
                    os.path.join(outdir, "%s_2.fastq.gz" % sraid)]
        if not utils.file_exists(out_file[0]):
            with utils.chdir(outdir):
                do.run(cmd.format(**locals()), "Covert to fastq %s" % sraid)
        if not utils.file_exists(out_file[0]):
            raise IOError("SRA %s didn't convert, something happened." % srafn)
        return [out for out in out_file if utils.file_exists(out)]
    else:
        raise ValueError("Not supported single-end sra samples for now.")
945580
def format(args):
    
    from jcvi.formats.obo import load_GODag, validate_term

    valid_multiparent_ops = ["split", "merge"]

    p = OptionParser(format.__doc__)

    g1 = OptionGroup(p, "Parameter(s) used to modify GFF attributes (9th column)")
    g1.add_option("--name", help="Add Name attribute from two-column file [default: %default]")
    g1.add_option("--note", help="Add Note attribute from two-column file [default: %default]")
    g1.add_option("--add_attribute", dest="attrib_files", help="Add new attribute(s) " + \
                "from two-column file(s); attribute name comes from filename; " + \
                "accepts comma-separated list of files [default: %default]")
    g1.add_option("--add_dbxref", dest="dbxref_files", help="Add new Dbxref value(s) (DBTAG:ID) " + \
                "from two-column file(s). DBTAG comes from filename, ID comes from 2nd column; " + \
                "accepts comma-separated list of files [default: %default]")
    g1.add_option("--nostrict", default=False, action="store_true",
                 help="Disable strict parsing of GFF file and/or mapping file [default: %default]")
    g1.add_option("--remove_attr", dest="remove_attrs", help="Specify attributes to remove; " + \
                "accepts comma-separated list of attribute names [default: %default]")
    g1.add_option("--copy_id_attr_to_name", default=False, action="store_true",
                 help="Copy `ID` attribute value to `Name`, when `Name` is not defined")
    g1.add_option("--invent_name_attr", default=False, action="store_true",
                 help="Invent `Name` attribute for 2nd level child features; " + \
                "Formatted like PARENT:FEAT_TYPE:FEAT_INDEX")
    g1.add_option("--no_keep_attr_order", default=False, action="store_true",
                 help="Do not maintain attribute order [default: %default]")
    p.add_option_group(g1)

    g2 = OptionGroup(p, "Parameter(s) used to modify content within columns 1-8")
    g2.add_option("--seqid", help="Switch seqid from two-column file. If not" + \
                " a file, value will globally replace GFF seqid [default: %default]")
    g2.add_option("--source", help="Switch GFF source from two-column file. If not" + \
                " a file, value will globally replace GFF source [default: %default]")
    g2.add_option("--type", help="Switch GFF feature type from two-column file. If not" + \
                " a file, value will globally replace GFF type [default: %default]")
    g2.add_option("--fixphase", default=False, action="store_true",
                 help="Change phase 1<->2, 2<->1 [default: %default]")
    p.add_option_group(g2)

    g3 = OptionGroup(p, "Other parameter(s) to perform manipulations to the GFF " + \
                 "file content")
    g3.add_option("--unique", default=False, action="store_true",
                 help="Make IDs unique [default: %default]")
    g3.add_option("--chaindup", default=None, dest="duptype",
                 help="Chain duplicate features of a particular GFF3 `type`," + \
                      " sharing the same ID attribute [default: %default]")
    g3.add_option("--multiparents", default=None, choices=valid_multiparent_ops,
                 help="Split/merge identical features (same `seqid`, `source`, `type` " + \
                 "`coord-range`, `strand`, `phase`) mapping to multiple parents " + \
                 "[default: %default]")
    g3.add_option("--remove_feats", help="Comma separated list of features to remove by type" + \
                " [default: %default]")
    g3.add_option("--remove_feats_by_ID", help="List of features to remove by ID;" + \
                " accepts comma-separated list or list file [default: %default]")
    g3.add_option("--gsac", default=False, action="store_true",
                 help="Fix GSAC GFF3 file attributes [default: %default]")
    g3.add_option("--invent_protein_feat", default=False, action="store_true",
                 help="Invent a protein feature span (chain CDS feats)")
    g3.add_option("--process_ftype", default=None, type="str",
                 help="Specify feature types to process; "
                 "accepts comma-separated list of feature types [default: %default]")
    g3.add_option("--gff3", default=False, action="store_true",
                 help="Print output in GFF3 format [default: %default]")
    g3.add_option("--make_gff_store", default=False, action="store_true",
                 help="Store entire GFF file in memory during first iteration [default: %default]")
    p.add_option_group(g3)

    p.set_outfile()
    p.set_SO_opts()

    opts, args = p.parse_args(args)

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

    gffile, = args
    mapfile = opts.seqid
    names = opts.name
    note = opts.note
    source = opts.source
    ftype = opts.type
    attrib_files = opts.attrib_files.split(",") if opts.attrib_files else None
    dbxref_files = opts.dbxref_files.split(",") if opts.dbxref_files else None
    remove_attrs = opts.remove_attrs.split(",") if opts.remove_attrs else None
    process_ftype = opts.process_ftype.split(",") if opts.process_ftype else None
    gsac = opts.gsac
    assert not (opts.unique and opts.duptype), \
        "Cannot use `--unique` and `--chaindup` together"
    assert not(opts.type and opts.duptype), \
        "Cannot use `--type` and `--chaindup` together"
    unique = opts.unique
    duptype = opts.duptype
    fixphase = opts.fixphase
    phaseT = {"1":"2", "2":"1"}
    remove_feats = opts.remove_feats.split(",") if opts.remove_feats else None
    remove_feats_by_ID = None
    if opts.remove_feats_by_ID:
        remove_feats_by_ID = LineFile(opts.remove_feats_by_ID, load=True).lines \
                if op.isfile(opts.remove_feats_by_ID) else \
                opts.remove_feats_by_ID.split(",")
    strict = False if opts.nostrict else True
    make_gff_store = True if gffile in ("-", "stdin") else opts.make_gff_store
    assert not (opts.copy_id_attr_to_name and opts.invent_name_attr), \
        "Cannot use `--copy_id_attr_to_name` and `--invent_name_attr` together"
    copy_id_attr_to_name = opts.copy_id_attr_to_name
    invent_name_attr = opts.invent_name_attr
    invent_protein_feat = opts.invent_protein_feat
    compute_signature = False

    outfile = opts.outfile

    mapping = None
    mod_attrs = set()
    if mapfile and op.isfile(mapfile):
        mapping = DictFile(mapfile, delimiter="\t", strict=strict)
        mod_attrs.add("ID")
    if note:
        note = DictFile(note, delimiter="\t", strict=strict)
        mod_attrs.add("Note")
    if source and op.isfile(source):
        source = DictFile(source, delimiter="\t", strict=strict)
    if ftype and op.isfile(ftype):
        ftype = DictFile(ftype, delimiter="\t", strict=strict)
    if names:
        names = DictFile(names, delimiter="\t", strict=strict)
        mod_attrs.add("Name")

    if attrib_files:
        attr_values = {}
        for fn in attrib_files:
            attr_name = op.basename(fn).rsplit(".", 1)[0]
            if attr_name not in reserved_gff_attributes:
                attr_name = attr_name.lower()
            attr_values[attr_name] = DictFile(fn, delimiter="\t", strict=strict)
            mod_attrs.add(attr_name)
    if dbxref_files:
        dbxref_values = {}
        for fn in dbxref_files:
            dbtag = op.basename(fn).rsplit(".", 1)[0]
            dbxref_values[dbtag] = DictFile(fn, delimiter="\t", strict=strict)
        mod_attrs.add("Dbxref")

    if remove_attrs:
        mod_remove_attrs = []
        for remove_attr in remove_attrs:
            if remove_attr in mod_attrs:
                mod_remove_attrs.append(remove_attr)

        if mod_remove_attrs:
            logging.error("Attributes `{0}` cannot be removed and modified".format( \
                    ",".join(mod_remove_attrs)))
            sys.exit()

    if gsac:  # setting gsac will force IDs to be unique
        unique = True
        notes = {}

    remove = set()
    if unique or duptype or remove_feats or remove_feats_by_ID \
            or opts.multiparents == "merge" or invent_name_attr or make_gff_store \
            or invent_protein_feat:
        if unique:
            dupcounts = defaultdict(int)
            seen = defaultdict(int)
            newparentid = {}
        elif duptype:
            dupranges = AutoVivification()
            skip = defaultdict(int)
        if opts.multiparents == "merge":
            merge_feats = AutoVivification()
        if invent_name_attr:
            ft = GffFeatureTracker()
        elif copy_id_attr_to_name:
            pass
        if invent_protein_feat:
            cds_track = {}
        if opts.multiparents == "merge" or invent_name_attr:
            make_gff_store = compute_signature = True
        gff = Gff(gffile, keep_attr_order=(not opts.no_keep_attr_order), \
                make_gff_store=make_gff_store, compute_signature=compute_signature, \
                strict=strict)
        for g in gff:
            if process_ftype and g.type not in process_ftype:
                continue
            id = g.accn
            if remove_feats and g.type in remove_feats:
                remove.add(id)
            if remove_feats_by_ID and id in remove_feats_by_ID:
                remove.add(id)
            if unique:
                dupcounts[id] += 1
            elif duptype and g.type == duptype:
                dupranges[g.seqid][id][g.idx] = (g.start, g.end)
            if opts.multiparents == "merge" and g.type != "CDS": #don't merge CDS
                pp = g.get_attr("Parent", first=False)
                if pp and len(pp) > 0:
                    for parent in pp:
                        if parent not in remove:
                            sig = g.sign
                            if sig not in merge_feats:
                                merge_feats[sig]['parents'] = []
                            if parent not in merge_feats[sig]['parents']:
                                merge_feats[sig]['parents'].append(parent)
            if invent_name_attr:
                parent, iso = atg_name(g.get_attr("Parent"), retval="locus,iso")
                if not parent:
                    parent = g.get_attr("Parent")
                ft.track(parent, g)
            if invent_protein_feat:
                if g.type == 'CDS':
                    cds_parent = g.get_attr("Parent")
                    if cds_parent not in cds_track:
                        cds_track[cds_parent] = []
                    cds_track[cds_parent].append((g.start, g.end))

    if opts.verifySO:
        so = load_GODag()
        valid_soterm = {}

    fw = must_open(outfile, "w")
    if not make_gff_store:
        gff = Gff(gffile, keep_attr_order=(not opts.no_keep_attr_order), \
                strict=strict)
    for g in gff:
        if process_ftype and g.type not in process_ftype:
            print(g, file=fw)
            continue

        id = g.accn

        if opts.multiparents == "merge" and g.type != "CDS": #don't merge CDS
            sig = g.sign
            if len(merge_feats[sig]['parents']) > 1:
                if 'candidate' not in merge_feats[sig]:
                    merge_feats[sig]['candidate'] = id
                    g.set_attr("Parent", merge_feats[sig]['parents'])
                else:
                    continue

        if remove_feats or remove_feats_by_ID:
            if id in remove:
                continue
            else:
                if "Parent" in g.attributes:
                    keep, parent = [], g.get_attr("Parent", first=False)
                    for i, pid in enumerate(parent):
                        if pid not in remove:
                            keep.append(parent[i])
                        else:
                            remove.add(id)
                    if len(keep) == 0:
                        continue
                    parent = g.set_attr("Parent", keep)

        if remove_attrs:
            for remove_attr in remove_attrs:
                if remove_attr in g.attributes:
                    g.set_attr(remove_attr, None)

        if opts.verifySO:
            if g.type not in valid_soterm:
                valid_soterm[g.type] = validate_term(g.type, so=so, method=opts.verifySO)
            ntype = valid_soterm[g.type]
            if ntype and g.type != ntype:
                g.type = ntype

        origid = g.seqid
        if fixphase:
            phase = g.phase
            g.phase = phaseT.get(phase, phase)

        if mapfile:
            if isinstance(mapping, dict):
                if origid in mapping:
                    g.seqid = mapping[origid]
                else:
                    logging.error("{0} not found in `{1}`. ID unchanged.".\
                            format(origid, mapfile))
            else:
                g.seqid = mapfile

        if source:
            if isinstance(source, dict) and g.source in source:
                g.source = source[g.source]
            else:
                g.source = source

        if names:
            if id in names:
                g.set_attr("Name", names[id])

        if note:
            name = g.get_attr("Name")
            tag = None
            if id in note:
                tag = note[id]
            elif name and name in note:
                tag = note[name]

            if tag:
                g.set_attr("Note", tag, update=False)

        if attrib_files:
            for attr_name in attr_values:
                name = g.get_attr("Name")
                if id in attr_values[attr_name]:
                    g.set_attr(attr_name, attr_values[attr_name][id])
                elif name and name in attr_values[attr_name]:
                    g.set_attr(attr_name, attr_values[attr_name][name])

        if dbxref_files:
            for dbtag in dbxref_values:
                if id in dbxref_values[dbtag]:
                    g.set_attr("Dbxref", dbxref_values[dbtag][id], dbtag=dbtag, append=True)

        if unique:
            if dupcounts[id] > 1:
                seen[id] += 1
                old_id = id
                id = "{0}-{1}".format(old_id, seen[old_id])
                newparentid[old_id] = id
                g.set_attr("ID", id)

            if "Parent" in g.attributes:
                parent = g.attributes["Parent"][0]
                if dupcounts[parent] > 1:
                    g.set_attr("Parent", newparentid[parent])

        if duptype:
            if duptype == g.type and len(dupranges[g.seqid][id]) > 1:
                p = sorted(dupranges[g.seqid][id])
                s, e = dupranges[g.seqid][id][p[0]][0:2]  # get coords of first encountered feature
                if g.start == s and g.end == e and p[0] == g.idx:
                    r = [dupranges[g.seqid][id][x] for x in dupranges[g.seqid][id]]
                    g.start, g.end = range_minmax(r)
                else:
                    skip[(g.seqid, g.idx, id, g.start, g.end)] = 1

        if gsac and g.type == "gene":
            notes[id] = g.attributes["Name"]

        if ftype:
            if isinstance(ftype, dict) and g.type in ftype:
                g.type = ftype[g.type]
            else:
                g.type = ftype

        if invent_name_attr:
            ft.store_symbol(g)
            if re.search(ft.ftype, g.type):
                parent, iso = atg_name(g.get_attr("Parent"), retval="locus,iso")
                if not parent:
                    parent = g.get_attr("Parent")
                if parent in ft.tracker:
                    fidx = ft.feat_index(parent, g.type, g.strand, (g.start, g.end, g.sign))
                    symbol = ft.get_symbol(parent)
                    attr = "ID" if symbol == parent else "Name"
                    g.set_attr(attr, "{0}:{1}:{2}".format(symbol, g.type, fidx + 1))
                    if opts.multiparents == "merge" and attr == "Name":
                        g.set_attr("ID", "{0}:{1}:{2}".format(parent, g.type, fidx + 1))
        elif copy_id_attr_to_name:
            if "Name" not in g.attributes.keys():
                g.set_attr("Name", g.get_attr("ID"))

        protein_feat = None
        if invent_protein_feat:
            if g.type == 'mRNA':
                if id in cds_track:
                    pstart, pstop = range_minmax(cds_track[id])
                    protein_feat = GffLine("\t".join(str(x) for x in [g.seqid, g.source, "protein", pstart, pstop, \
                            ".", g.strand, ".", "ID={0}-Protein;Name={0};Derives_from={0}".format(id)]))
            elif g.type == 'CDS':
                parent = g.get_attr("Parent")
                if parent in cds_track:
                    _parent = [parent, "{0}-Protein".format(parent)]
                    g.set_attr("Parent", _parent)

        pp = g.get_attr("Parent", first=False)
        if opts.multiparents == "split" and (pp and len(pp) > 1) and g.type != "CDS":  # separate features with multiple parents
            id = g.get_attr("ID")
            for i, parent in enumerate(pp):
                if id: g.set_attr("ID", "{0}-{1}".format(id, i + 1))
                g.set_attr("Parent", parent, update=True, urlquote=True)
                if gsac:
                    fix_gsac(g, notes)
                print(g, file=fw)
        else:
            if g.gff3 and not opts.gff3:
                opts.gff3 = True
            g.update_attributes(gff3=opts.gff3)
            if gsac:
                fix_gsac(g, notes)
            if duptype == g.type and skip[(g.seqid, g.idx, id, g.start, g.end)] == 1:
                continue
            print(g, file=fw)
            if g.type == 'mRNA' and invent_protein_feat:
                print(protein_feat, file=fw)

    fw.close()
1105895
def _clean_file(self, filename):
        """"""
        if not os.path.exists(filename):
            return
        self.announce("removing '{}'".format(filename))
        if not self.dry_run:
            os.remove(filename)
451199
def process_fastq_plain(fastq, **kwargs):
    """"""
    logging.info("Nanoget: Starting to collect statistics from plain fastq file.")
    inputfastq = handle_compressed_input(fastq)
    return ut.reduce_memory_usage(pd.DataFrame(
        data=[res for res in extract_from_fastq(inputfastq) if res],
        columns=["quals", "lengths"]
    ).dropna())
691805
def parse_fasta(infile, upper=False):
    '''
    
    '''
    try:
        fp = must_open(infile)
    except:
        fp = infile
    # keep header
    fa_iter = (x[1] for x in groupby(fp, lambda row: row[0] == '>'))
    for header in fa_iter:
        header = next(header)
        if header[0] != '>':
            continue
        # drop '>'
        header = header.strip()[1:]
        # stitch the sequence lines together and make into upper case
        seq = "".join(s.strip() for s in next(fa_iter))
        if upper:
            seq = seq.upper()
        yield header, seq
22211
def reformat(self, snip, index=None):
        """"""
        if index is not None:
            return self.reformat(snip), len(self.reformat(snip[:index]))
        else:
            return self.repl_proc(snip, reformatting=True, log=False)
122892
def dump(args):
    
    p = OptionParser(dump.__doc__)
    p.add_option("--dir",
                help="Working directory [default: %default]")
    p.add_option("--nosim", default=False, action="store_true",
                 help="Do not simulate qual to 50 [default: %default]")
    opts, args = p.parse_args(args)

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

    fastbfile, = args
    d = opts.dir
    if d:
        from jcvi.assembly.preprocess import export_fastq

        rc = "jump" in fastbfile
        export_fastq(d, fastbfile, rc=rc)
        return

    sim = not opts.nosim
    pf = "j" if "jump" in fastbfile else "f"

    statsfile = "{0}.lib_stats".format(pf)
    if op.exists(statsfile):
        os.remove(statsfile)

    cmd = "SplitReadsByLibrary READS_IN={0}".format(fastbfile)
    cmd += " READS_OUT={0} QUALS=True".format(pf)
    sh(cmd)

    libs = []
    fp = open(statsfile)
    next(fp); next(fp)  # skip two rows
    for row in fp:
        if row.strip() == "":
            continue

        libname = row.split()[0]
        if libname == "Unpaired":
            continue

        libs.append(libname)

    logging.debug("Found libraries: {0}".format(",".join(libs)))

    cmds = []
    for libname in libs:
        cmd = "FastbQualbToFastq"
        cmd += " HEAD_IN={0}.{1}.AB HEAD_OUT={1}".format(pf, libname)
        cmd += " PAIRED=True PHRED_OFFSET=33"
        if sim:
            cmd += " SIMULATE_QUALS=True"
        if pf == 'j':
            cmd += " FLIP=True"

        cmds.append((cmd, ))

    m = Jobs(target=sh, args=cmds)
    m.run()

    for libname in libs:
        cmd = "mv {0}.A.fastq {0}.1.fastq".format(libname)
        sh(cmd)
        cmd = "mv {0}.B.fastq {0}.2.fastq".format(libname)
        sh(cmd)
442345
def patch(args):
    
    from jcvi.formats.base import write_file
    from jcvi.formats.fasta import format

    p = OptionParser(patch.__doc__)
    p.add_option("--cleanfasta", default=False, action="store_true",
                 help="Clean FASTA to remove description [default: %default]")
    p.add_option("--highqual", default=False, action="store_true",
                 help="Reads are of high quality [default: %default]")
    p.set_home("pbjelly")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    ref, reads = args
    cpus = opts.cpus
    cmd = op.join(opts.pbjelly_home, "setup.sh")
    setup = "source {0}".format(cmd)
    if not which("fakeQuals.py"):
        sh(setup)

    pf = ref.rsplit(".", 1)[0]
    pr, px = reads.rsplit(".", 1)
    # Remove description line
    if opts.cleanfasta:
        oref = pf + ".f.fasta"
        oreads = pr + ".f.fasta"
        format([ref, oref])
        format([reads, oreads])
        ref, reads = oref, oreads

    # Check if the FASTA has qual
    ref, refq = fake_quals(ref)
    convert_reads = not px in ("fq", "fastq", "txt")
    if convert_reads:
        reads, readsq = fake_quals(reads)
        readsfiles = " ".join((reads, readsq))
    else:
        readsfiles = reads

    # Make directory structure
    dref, dreads = "data/reference", "data/reads"
    cwd = os.getcwd()
    reference = op.join(cwd, "{0}/{1}".format(dref, ref))
    reads = op.join(cwd, "{0}/{1}".format(dreads, reads))
    if not op.exists(reference):
        sh("mkdir -p {0}".format(dref))
        sh("cp {0} {1}/".format(" ".join((ref, refq)), dref))
    if not op.exists(reads):
        sh("mkdir -p {0}".format(dreads))
        sh("cp {0} {1}/".format(readsfiles, dreads))

    outputDir = cwd
    p = Protocol(outputDir, reference, reads, highqual=opts.highqual)
    p.write_xml()

    # Build the pipeline
    runsh = [setup]
    for action in "setup|mapping|support|extraction".split("|"):
        runsh.append("Jelly.py {0} Protocol.xml".format(action))

    runsh.append('Jelly.py assembly Protocol.xml -x "--nproc={0}"'.format(cpus))
    runsh.append("Jelly.py output Protocol.xml")

    runfile = "run.sh"
    contents = "\n".join(runsh)
    write_file(runfile, contents)
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)
1026979
def fastaWrite(fileHandleOrFile, name, seq, mode="w"):
    
    fileHandle = _getFileHandle(fileHandleOrFile, mode)
    valid_chars = {x for x in string.ascii_letters + "-"}
    try:
        assert any([isinstance(seq, unicode), isinstance(seq, str)])
    except AssertionError:
        raise RuntimeError("Sequence is not unicode or string")
    try:
        assert all(x in valid_chars for x in seq)
    except AssertionError:
        bad_chars = {x for x in seq if x not in valid_chars}
        raise RuntimeError("Invalid FASTA character(s) see in fasta sequence: {}".format(bad_chars))
    fileHandle.write(">%s\n" % name)
    chunkSize = 100
    for i in xrange(0, len(seq), chunkSize):
        fileHandle.write("%s\n" % seq[i:i+chunkSize])
    if isinstance(fileHandleOrFile, "".__class__):
        fileHandle.close()
624838
def fix_fasta(fasta):
    
    for seq in parse_fasta(fasta):
        seq[0] = remove_char(seq[0])
        if len(seq[1]) > 0:
            yield seq
218877
def extract(args):
    
    p = OptionParser(extract.__doc__)
    p.add_option("--format", default=False, action="store_true",
            help="enable flag to reformat header into a symbol separated list of constituent reads "+ \
            "[default: %default]")
    p.add_option("--singlets", default=False, action="store_true",
            help="ask the program to look in the singlets file (should be in the same folder) for " +\
            "unused reads and put them in the resultant fasta file [default: %default]")
    p.set_sep(sep="|", help="Separator used to list the reads in the FASTA header")
    opts, args = p.parse_args(args)

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

    acefile, = args
    ace = Ace.read(must_open(acefile))
    logging.debug('Loaded ace file {0}'.format(acefile))

    fastafile = acefile.rsplit(".", 1)[0] + ".fasta"
    fw = open(fastafile, "w")
    for c in ace.contigs:
        id = c.name
        if opts.format:
           id = opts.sep.join([read.name for read in c.af])

        seqrec = SeqRecord(Seq(c.sequence), id=id, description="")
        SeqIO.write([seqrec], fw, "fasta")

    if opts.singlets:
        singletsfile = acefile.rsplit(".", 1)[0] + ".singlets"
        if os.path.getsize(singletsfile) > 0:
            fp = SeqIO.parse(must_open(singletsfile), "fasta")
            for rec in fp:
                SeqIO.write(rec, fw, "fasta")

    fw.close()
    logging.debug('Wrote contigs to fasta file {0}'.format(fastafile))
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
590978
def format_stack(skip=0, length=6, _sep=os.path.sep):
    
    return ' < '.join("%s:%s:%s" % (
        '/'.join(f.f_code.co_filename.split(_sep)[-2:]),
        f.f_lineno,
        f.f_code.co_name
    ) for f in islice(frame_iterator(sys._getframe(1 + skip)), length))