def extract(args):
    """
    %prog extract fasta query

    extract query out of fasta file, query needs to be in the form of
    "seqname", or "seqname:start-stop", or "seqname:start-stop:-"
    """
    p = OptionParser(extract.__doc__)
    p.add_option('--include', default=False, action="store_true",
            help="search description line for match [default: %default]")
    p.add_option('--exclude', default=False, action="store_true",
            help="exclude description that matches [default: %default]")
    p.add_option('--idonly', default=False, action="store_true",
            help="Only search identifier")
    p.add_option('--bed', default=None,
            help="path to bed file to guide extraction by matching seqname")
    p.set_outfile()

    opts, args = p.parse_args(args)

    if len(args) == 2:
        fastafile, query = args
    elif len(args) == 1 and opts.bed:
        fastafile, = args
        bedaccns = Bed(opts.bed).accns
    else:
        sys.exit(p.print_help())

    if opts.bed:
        fw = must_open(opts.outfile, "w")
        f = Fasta(fastafile)
        for accn in bedaccns:
            try:
                rec = f[accn]
            except:
                logging.error("{0} not found in {1}".format(accn, fastafile))
                continue
            SeqIO.write([rec], fw, "fasta")
        return fw.name

    atoms = query.split(":")
    key = atoms[0]

    assert len(atoms) <= 3, "cannot have more than two ':' in your query"

    pos = ""
    if len(atoms) in (2, 3):
        pos = atoms[1]

    strand = "+"
    if len(atoms) == 3:
        strand = atoms[2]

    assert strand in ('+', '-'), "strand must be either '+' or '-'"

    feature = dict(chr=key)

    if "-" in pos:
        start, stop = pos.split("-")
        try:
            start, stop = int(start), int(stop)
        except ValueError as e:
            logging.error(e)
            sys.exit(p.print_help())

        feature["start"] = start
        feature["stop"] = stop
    else:
        start, stop = None, None

    assert start < stop or None in (start, stop), \
            "start must be < stop, you have ({0}, {1})".format(start, stop)
    feature["strand"] = strand

    include, exclude = opts.include, opts.exclude
    # conflicting options, cannot be true at the same time
    assert not (include and exclude), "--include and --exclude cannot be "\
            "on at the same time"
    fw = must_open(opts.outfile, "w")

    if include or exclude:
        f = Fasta(fastafile, lazy=True)
        fi = f.iteritems_ordered if opts.idonly else f.iterdescriptions_ordered
        for k, rec in fi():
            if include and key not in k:
                continue
            if exclude and key in k:
                continue

            seq = Fasta.subseq(rec, start, stop, strand)
            newid = rec.id
            if start is not None:
                newid += ":{0}-{1}:{2}".format(start, stop, strand)

            rec = SeqRecord(seq, id=newid, description=k)
            SeqIO.write([rec], fw, "fasta")
    else:
        f = Fasta(fastafile)
        try:
            seq = f.sequence(feature, asstring=False)
        except AssertionError as e:
            logging.error(e)
            return

        rec = SeqRecord(seq, id=query, description="")
        SeqIO.write([rec], fw, "fasta")

    return fw.name
dataset/ETHPy150Open tanghaibao/jcvi/formats/fasta.py/extract/original