%prog translate cdsfasta

    Translate CDS to proteins. The tricky thing is that sometimes the CDS
    represents a partial gene, therefore disrupting the frame of the protein.
    Check all three frames to get a valid translation.
12751
def translate(args):
    
    transl_tables = [str(x) for x in xrange(1,25)]
    p = OptionParser(translate.__doc__)
    p.add_option("--ids", default=False, action="store_true",
                 help="Create .ids file with the complete/partial/gaps "
                      "label [default: %default]")
    p.add_option("--longest", default=False, action="store_true",
                 help="Find the longest ORF from each input CDS [default: %default]")
    p.add_option("--table", default=1, choices=transl_tables,
            help="Specify translation table to use [default: %default]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    cdsfasta, = args
    if opts.longest:
        cdsfasta = longestorf([cdsfasta])

    f = Fasta(cdsfasta, lazy=True)
    outfile = opts.outfile
    fw = must_open(outfile, "w")

    if opts.ids:
        idsfile = cdsfasta.rsplit(".", 1)[0] + ".ids"
        ids = open(idsfile, "w")
    else:
        ids = None

    five_prime_missing = three_prime_missing = 0
    contain_ns = complete = cannot_translate = total = 0

    for name, rec in f.iteritems_ordered():
        cds = rec.seq
        cdslen = len(cds)
        peplen = cdslen / 3
        total += 1

        # Try all three frames
        pep = ""
        for i in xrange(3):
            newcds = cds[i: i + peplen * 3]
            newpep = newcds.translate(table=opts.table)
            if len(newpep.split("*")[0]) > len(pep.split("*")[0]):
                pep = newpep

        labels = []
        if "*" in pep.rstrip("*"):
            logging.error("{0} cannot translate".format(name))
            cannot_translate += 1
            labels.append("cannot_translate")

        contains_start = pep.startswith("M")
        contains_stop = pep.endswith("*")
        contains_ns = "X" in pep
        start_ns = pep.startswith("X")
        end_ns = pep.endswith("X")

        if not contains_start:
            five_prime_missing += 1
            labels.append("five_prime_missing")
        if not contains_stop:
            three_prime_missing += 1
            labels.append("three_prime_missing")
        if contains_ns:
            contain_ns += 1
            labels.append("contain_ns")
        if contains_start and contains_stop:
            complete += 1
            labels.append("complete")
        if start_ns:
            labels.append("start_ns")
        if end_ns:
            labels.append("end_ns")

        if ids:
            print("\t".join((name, ",".join(labels))), file=ids)

        peprec = SeqRecord(pep, id=name, description=rec.description)
        SeqIO.write([peprec], fw, "fasta")
        fw.flush()

    print("Complete gene models: {0}".\
                        format(percentage(complete, total)), file=sys.stderr)
    print("Missing 5`-end: {0}".\
                        format(percentage(five_prime_missing, total)), file=sys.stderr)
    print("Missing 3`-end: {0}".\
                        format(percentage(three_prime_missing, total)), file=sys.stderr)
    print("Contain Ns: {0}".\
                        format(percentage(contain_ns, total)), file=sys.stderr)

    if cannot_translate:
        print("Cannot translate: {0}".\
                        format(percentage(cannot_translate, total)), file=sys.stderr)

    fw.close()

    return cdsfasta, outfile
307452
def validate(args):
    
    import pyfasta

    p = OptionParser(validate.__doc__)
    p.add_option("--prefix", help="Add prefix to seqid")
    opts, args = p.parse_args(args)

    vcffile, fastafile = args
    pf = opts.prefix
    genome = pyfasta.Fasta(fastafile, record_class=pyfasta.MemoryRecord)
    fp = must_open(vcffile)
    match_ref = match_alt = total = 0
    for row in fp:
        if row[0] == '#':
            continue
        seqid, pos, id, ref, alt = row.split()[:5]
        total += 1
        if pf:
            seqid = pf + seqid
        pos = int(pos)
        if seqid not in genome:
            continue
        true_ref = genome[seqid][pos - 1]
        if total % 100000 == 0:
            print(total, "sites parsed", file=sys.stderr)
        if ref == true_ref:
            match_ref += 1
        elif alt == true_ref:
            match_alt += 1

    logging.debug("Match REF: {}".format(percentage(match_ref, total)))
    logging.debug("Match ALT: {}".format(percentage(match_alt, total)))
1281335
def localize(self, *positions, **kwargs):
        
        cdna = kwargs.get('cdna', False)
        # TODO: account for strand ?? add kwarg ??
        # if it's to the CDNA, then it's based on the cdsStart
        start, end = (self.cdsStart, self.cdsEnd) if cdna else \
                                        (self.start, self.end)
        introns = self.introns or None
        if cdna:
            if not self.is_coding:
                return ([None] * len(positions)) if len(positions) > 1 else None
            introns = self._introns(self.cds) or None

        if introns is None:
            local_ps = [p - start if (start <= p < end) else None for p in positions]
            return local_ps[0] if len(positions) == 1 else local_ps

        introns = [(s - start, e - start) for s, e in introns]
        positions = [p - start for p in positions]
        # now both introns and positions are local starts based on cds/tx-Start
        local_ps = []
        l = end - start
        for original_p in positions:
            subtract = 0
            p = original_p
            print(p, l, file=sys.stderr)
            if p < 0 or p >= l: # outside of transcript
                local_ps.append(None)
                continue
            for s, e in introns:
                # within intron
                if s <= p <= e:
                    subtract = None
                    break
                # otherwise, adjust for intron length.
                elif p >= e:
                    subtract += (e - s)

            local_ps.append(p - subtract if subtract is not None else None)

        assert all(p is None or p >=0 for p in local_ps), (local_ps)
        return local_ps[0] if len(positions) == 1 else local_ps
950346
def from_translation_key(
            cls,
            translation_key,
            translations,
            overlapping_reads,
            ref_reads,
            alt_reads,
            alt_reads_supporting_protein_sequence,
            transcripts_overlapping_variant,
            transcripts_supporting_protein_sequence,
            gene):
        
        return cls(
            amino_acids=translation_key.amino_acids,
            variant_aa_interval_start=translation_key.variant_aa_interval_start,
            variant_aa_interval_end=translation_key.variant_aa_interval_end,
            ends_with_stop_codon=translation_key.ends_with_stop_codon,
            frameshift=translation_key.frameshift,
            translations=translations,
            overlapping_reads=overlapping_reads,
            ref_reads=ref_reads,
            alt_reads=alt_reads,
            alt_reads_supporting_protein_sequence=(
                alt_reads_supporting_protein_sequence),
            transcripts_overlapping_variant=transcripts_overlapping_variant,
            transcripts_supporting_protein_sequence=(
                transcripts_supporting_protein_sequence),
            gene=gene)
526958
def some(args):
    
    p = OptionParser(some.__doc__)
    p.add_option("--exclude", default=False, action="store_true",
            help="Output sequences not in the list file [default: %default]")
    p.add_option("--uniprot", default=False, action="store_true",
            help="Header is from uniprot [default: %default]")

    opts, args = p.parse_args(args)

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

    fastafile, listfile, outfastafile = args
    outfastahandle = must_open(outfastafile, "w")
    qualfile = get_qual(fastafile)

    names = set(x.strip() for x in open(listfile))
    if qualfile:
        outqualfile = outfastafile + ".qual"
        outqualhandle = open(outqualfile, "w")
        parser = iter_fasta_qual(fastafile, qualfile)
    else:
        parser = SeqIO.parse(fastafile, "fasta")

    num_records = 0
    for rec in parser:
        name = rec.id
        if opts.uniprot:
            name = name.split("|")[-1]

        if opts.exclude:
            if name in names:
                continue
        else:
            if name not in names:
                continue

        SeqIO.write([rec], outfastahandle, "fasta")
        if qualfile:
            SeqIO.write([rec], outqualhandle, "qual")

        num_records += 1

    logging.debug("A total of %d records written to `%s`" % \
            (num_records, outfastafile))
70351
def fromagp(args):
    
    from jcvi.formats.agp import AGP
    from jcvi.formats.sizes import Sizes

    p = OptionParser(fromagp.__doc__)
    p.add_option("--novalidate", default=False, action="store_true",
                 help="Do not validate AGP")
    opts, args = p.parse_args(args)

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

    agpfile, componentfasta, objectfasta = args
    chainfile = agpfile.rsplit(".", 1)[0] + ".chain"
    fw = open(chainfile, "w")
    agp = AGP(agpfile, validate=(not opts.novalidate))
    componentsizes = Sizes(componentfasta).mapping
    objectsizes = Sizes(objectfasta).mapping
    chain = "chain"
    score = 1000
    tStrand = "+"
    id = 0
    for a in agp:
        if a.is_gap:
            continue

        tName = a.component_id
        tSize = componentsizes[tName]
        tStart = a.component_beg
        tEnd = a.component_end
        tStart -= 1

        qName = a.object
        qSize = objectsizes[qName]
        qStrand = "-" if a.orientation == "-" else "+"
        qStart = a.object_beg
        qEnd = a.object_end
        if qStrand == '-':
            _qStart = qSize - qEnd + 1
            _qEnd = qSize - qStart + 1
            qStart, qEnd = _qStart, _qEnd
        qStart -= 1

        id += 1
        size = a.object_span
        headerline = "\t".join(str(x) for x in (
             chain, score, tName, tSize, tStrand, tStart,
             tEnd, qName, qSize, qStrand, qStart, qEnd, id
        ))
        alignmentline = size
        print(headerline, file=fw)
        print(alignmentline, file=fw)
        print(file=fw)

    fw.close()
    logging.debug("File written to `{0}`.".format(chainfile))
1003355
def translations_generator_to_dataframe(translations_generator):
    
    return dataframe_from_generator(
        element_class=Translation,
        variant_and_elements_generator=translations_generator,
        exclude=[],
        converters={
            "untrimmed_variant_sequence": lambda vs: vs.sequence,
            "variant_sequence_in_reading_frame": (
                lambda vs: vs.in_frame_cdna_sequence),
            "reference_context": (
                lambda rc: ";".join([
                    transcript.name for
                    transcript in rc.transcripts]))
        },
        extra_column_fns={
            "untrimmed_variant_sequence_read_count": (
                lambda _, t: len(t.untrimmed_variant_sequence.reads)),
        })
367067
function translate(cond) {

    let result = {};

    // If condition is not an array...
    if (!Array.isArray(cond)) {
        //...then it is a condition cluster
        if (cond.and) {
            result = translateAndList(cond.and)
        }
        else if (cond.or) {
            result = translateOrList(cond.or)
        }
        else {
            //This is an illegal subject.
            throw new Error('Invalid query Object keys. Objects must be either "and" or "or"');
        }
    }
    else {
        //...then it's a singular condition
        assert(cond.length === 3, 'Invalid condition item. Condition items must be an array of 3 items representing field, comparator, and value')

        let fieldName = cond[0];
        fieldName = fieldName.replace('.*', '');
        let comp = cond[1];
        let value = cond[2];

        if (fieldName === 'id') {
            fieldName = '_id';
            value = Array.isArray(value) ? value : new ObjectID(value)
        }

        // construct mini block
        if (comp === '=') {

            if (value === null || value === undefined) {
                value = {};
                value['$exists'] = false;
            }

            result[fieldName] = value;
        } else if (comp === '!=' && (value === null || value === undefined)) {
            let value = {};
            value['$exists'] = true;

            result[fieldName] = value;
        } else if (comp === 'like') {

            result[fieldName] = new RegExp(value);

        } else if (comp === 'nlike') {

            result[fieldName] = {
                $not: new RegExp(value)
            };
        }
        else {
            const subCondBlock = {};

            subCondBlock[convertComp(comp)] = value;

            result[fieldName] = subCondBlock;
        }
    }

    return result;

}
573154
public ExtendedResult<TranslationResult> translate(CdsFeature cds, Entry entry) throws RepositoryException {
    
        ExtendedResult<TranslationResult> validationResult = new ExtendedResult<TranslationResult>();
        if(cds == null || cds.getLocations() == null)
            return validationResult;
        
    	if((Entry.CON_DATACLASS.equals(entry.getDataClass()) && (entry.getSequence() == null || entry.getSequence().getSequenceByte() == null)) )
    		return validationResult;

        // Check translation amino acids.
		if (cds.getTranslation() != null) {
			CdsFeatureAminoAcidCheck validationCheck = new CdsFeatureAminoAcidCheck();
			validationResult.append(validationCheck.check(cds));
		}

        if(entry.getSequence() == null || entry.getSequence().getSequenceByte() == null){
            validationResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.ERROR, "CDSTranslator-5"));
        }

        if (validationResult.count(Severity.ERROR) > 0) {
			return validationResult;
		}

		// Check exceptional translation.
		if (cds.isException() && cds.getTranslation() == null) {
			// Missing exceptional translation.
            validationResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.ERROR, "CDSTranslator-1"));
        }
		if (validationResult.count(Severity.ERROR) > 0) {
			return validationResult;
		}

        createTranslator(cds, entry, validationResult);

        if(!checkLocations(entry.getSequence().getLength(), cds.getLocations())){
            validationResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.ERROR, "CDSTranslator-4"));
        }

		if (validationResult.count(Severity.ERROR) > 0) {
			return validationResult;
		}

		try
		{
		SegmentFactory factory=new SegmentFactory(planProperty.enproConnection.get()); 
		Segment segment = factory.createSegment(entry.getSequence(), cds.getLocations());
		if(segment==null)
		{
			return validationResult;
		}
		ExtendedResult<TranslationResult> extendedTranslatorResult = translator.translate(segment.getSequenceByte(), entry.getOrigin());
		
		for (ValidationMessage<Origin> message : extendedTranslatorResult.getMessages()) {
			message.append(cds.getOrigin());
// 			message.append(EntryOrigin.origin(entry));
		}
		extendedTranslatorResult.append(validationResult);

		TranslationResult translationResult = extendedTranslatorResult.getExtension();
		String expectedTranslation = cds.getTranslation();
		String conceptualTranslation = translationResult.getConceptualTranslation();


        if (extendedTranslatorResult.count(Severity.ERROR) > 0) {
            return extendedTranslatorResult;
        }

		// Note that if the compound location is a global complement
		// then the left and right partiality are reversed between CdsFeature
		// and Translator.

        if (translationResult.isFixedLeftPartial()) {
            if (!cds.getLocations().isComplement()) {
                cds.getLocations().setLeftPartial(translator.isLeftPartial());
            }
            else {
                cds.getLocations().setRightPartial(translator.isLeftPartial());
            }
        }

        if (translationResult.isFixedRightPartial()) {
            if (!cds.getLocations().isComplement()) {
                cds.getLocations().setRightPartial(translator.isRightPartial());
            }
            else {
                cds.getLocations().setLeftPartial(translator.isRightPartial());
            }
        }

		if(translationResult.isFixedPseudo()) {
		    cds.addQualifier(new QualifierFactory().createQualifier(Qualifier.PSEUDO_QUALIFIER_NAME));
		    cds.setTranslation(null);
        }
        if (expectedTranslation == null ||
			expectedTranslation.length() == 0 ) {
        	if(!(cds instanceof PeptideFeature))
			cds.setTranslation(conceptualTranslation); // Accept the translation.
		}
		else if (acceptTranslation 
				&& !translator.equalsTranslation(expectedTranslation, conceptualTranslation)) {
			cds.setTranslation(conceptualTranslation); // Accept the translation.
			// Add warning
			extendedTranslatorResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.WARNING,
                    "CDSTranslator-2"));
		}
		else if (!cds.isException() && !cds.isPseudo() &&
                !translator.equalsTranslation(expectedTranslation, conceptualTranslation)) {
			// Reject the translation.
			extendedTranslatorResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.ERROR,
                    "CDSTranslator-2"));
		}
		else if (cds.isException() && translator.equalsTranslation(expectedTranslation, conceptualTranslation)) {
			extendedTranslatorResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.WARNING,
                    "CDSTranslator-3"));
		}

		Set<String> fixes = translator.getFixes();
		if(fixes != null && !fixes.isEmpty()) {
		    fixes.forEach(f -> extendedTranslatorResult.append(EntryValidations.createMessage(cds.getOrigin(), Severity.FIX, f)));
        }
		
		if (severity != null) {
			Collection<ValidationMessage<Origin>> messages = 
				extendedTranslatorResult.getMessages();
			for (ValidationMessage<Origin> message : messages) {
				message.setSeverity(severity);
			}
		}
		
		return extendedTranslatorResult;
		
		}
		catch(Exception e)
		{
			throw new RepositoryException(e);
		}
	}
796096
def translate(self, exc):
        """"""
        from boto.exception import StorageResponseError

        if isinstance(exc, StorageResponseError):
            if exc.status == 404:
                return self.error_cls(str(exc))

        return None
228150
function translateUiTitles (elements) {
  if (elements.title && typeof elements.title == 'string') {
    elements.title = flattenTitles(elements.title)
  }
  if (elements.elements) {
    translateUiTitles(elements.elements)
  }
  for (section in elements) {
    if (elements[section].ui) {
      var sectionTitle = elements[section].ui.title
      elements[section].ui.title = flattenTitles(sectionTitle)
    }
    if (elements[section].components) {
      translateUiTitles(elements[section].components.elements)
    }
  }
  return elements
}
705163
def _retranslate(seq):
        '''
        
        '''
        if len(seq['vdj_nt']) % 3 != 0:
            trunc = len(seq['vdj_nt']) % 3
            seq['vdj_nt'] = seq['vdj_nt'][:-trunc]
        seq['vdj_aa'] = Seq(seq['vdj_nt'], generic_dna).translate()
255055
def handle_program_options():
    
    parser = argparse.ArgumentParser(description="Convert Sanger-sequencing \
                                     derived data files for use with the \
                                     metagenomics analysis program QIIME, by \
                                     extracting Sample ID information, adding\
                                     barcodes and primers to the sequence \
                                     data, and outputting a mapping file and\
                                     single FASTA-formatted sequence file \
                                     formed by concatenating all input data.")
    parser.add_argument('-i', '--input_dir', required=True,
                        help="The directory containing sequence data files. \
                              Assumes all data files are placed in this \
                              directory. For files organized within folders by\
                              sample, use -s in addition.")
    parser.add_argument('-m', '--map_file', default='map.txt',
                        help="QIIME-formatted mapping file linking Sample IDs \
                              with barcodes and primers.")
    parser.add_argument('-o', '--output', default='output.fasta',
                        metavar='OUTPUT_FILE',
                        help="Single file containing all sequence data found \
                              in input_dir, FASTA-formatted with barcode and \
                              primer preprended to sequence. If the -q option \
                              is passed, any quality data will also be output \
                              to a single file of the same name with a .qual \
                              extension.")
    parser.add_argument('-b', '--barcode_length', type=int, default=12,
                        help="Length of the generated barcode sequences. \
                              Default is 12 (QIIME default), minimum is 8.")

    parser.add_argument('-q', '--qual', action='store_true', default=False,
                        help="Instruct the program to look for quality \
                              input files")
    parser.add_argument('-u', '--utf16', action='store_true', default=False,
                        help="UTF-16 encoded input files")

    parser.add_argument('-t', '--treatment',
                        help="Inserts an additional column into the mapping \
                              file specifying some treatment or other variable\
                              that separates the current set of sequences \
                              from any other set of seqeunces. For example:\
                              -t DiseaseState=healthy")

    # data input options
    sidGroup = parser.add_mutually_exclusive_group(required=True)
    sidGroup.add_argument('-d', '--identifier_pattern',
                          action=ValidateIDPattern,
                          nargs=2, metavar=('SEPARATOR', 'FIELD_NUMBER'),
                          help="Indicates how to extract the Sample ID from \
                               the description line. Specify two things: \
                               1. Field separator, 2. Field number of Sample \
                               ID (1 or greater). If the separator is a space \
                               or tab, use \s or \\t respectively. \
                               Example: >ka-SampleID-2091, use -i - 2, \
                               indicating - is the separator and the Sample ID\
                               is field #2.")
    sidGroup.add_argument('-f', '--filename_sample_id', action='store_true',
                          default=False, help='Specify that the program should\
                          the name of each fasta file as the Sample ID for use\
                          in the mapping file. This is meant to be used when \
                          all sequence data for a sample is stored in a single\
                          file.')

    return parser.parse_args()
667111
def run(self):
        
        # thanks to deluge guys ;)
        po_dir = os.path.join(os.path.dirname(__file__), 'webant', 'translations')
        print('Compiling po files from "{}"...'.format(po_dir))
        for lang in os.listdir(po_dir):
            sys.stdout.write("\tCompiling {}... ".format(lang))
            sys.stdout.flush()
            curr_lang_path = os.path.join(po_dir, lang)
            for path, dirs, filenames in os.walk(curr_lang_path):
                for f in filenames:
                    if f.endswith('.po'):
                        src = os.path.join(path, f)
                        dst = os.path.join(path, f[:-3] + ".mo")
                        if not os.path.exists(dst) or self.force:
                            msgfmt.make(src, dst)
                            print("ok.")
                        else:
                            src_mtime = os.stat(src)[8]
                            dst_mtime = os.stat(dst)[8]
                            if src_mtime > dst_mtime:
                                msgfmt.make(src, dst)
                                print("ok.")
                            else:
                                print("already up to date.")
        print('Finished compiling translation files.')
648438
def validate(args):
    
    p = OptionParser(validate.__doc__)
    opts, args = p.parse_args(args)

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

    imputed, withheld = args
    register = {}
    fp = open(withheld)
    for row in fp:
        if row[0] == "#":
            continue
        v = VcfLine(row)
        register[(v.seqid, v.pos)] = v.genotype

    logging.debug("Imported {0} records from `{1}`".\
                    format(len(register), withheld))

    fp = must_open(imputed)
    hit = concordant = 0
    seen = set()
    for row in fp:
        if row[0] == "#":
            continue
        v = VcfLine(row)
        chr, pos, genotype = v.seqid, v.pos, v.genotype
        if (chr, pos) in seen:
            continue
        seen.add((chr, pos))
        if (chr, pos) not in register:
            continue
        truth = register[(chr, pos)]
        imputed = genotype.split(":")[0]
        if "|" in imputed:
            imputed = "/".join(sorted(genotype.split(":")[0].split("|")))
            #probs = [float(x) for x in genotype.split(":")[-1].split(",")]
            #imputed = max(zip(probs, ["0/0", "0/1", "1/1"]))[-1]
        hit += 1
        if truth == imputed:
            concordant += 1
        else:
            print(row.strip(), "truth={0}".format(truth), file=sys.stderr)

    logging.debug("Total concordant: {0}".\
            format(percentage(concordant, hit)))
1062372
def parse_fasta_annotations(fastas, annot_tables, trans_table):
    
    if annot_tables is not False:
        annots = {}
        for table in annot_tables:
            for cds in open(table):
                ID, start, end, strand = cds.strip().split()
                annots[ID] = [start, end, int(strand)]
    for fasta in fastas:
        for seq in parse_fasta(fasta):
            if ('# ;gc_cont' not in seq[0] and '# ID=' not in seq[0]) and annot_tables is False:
                print('# specify fasta from Prodigal or annotations table (-t)', file=sys.stderr)
                exit()
            if 'ID=' in seq[0]:
                ID = seq[0].rsplit('ID=', 1)[1].split(';', 1)[0]
                contig = seq[0].split()[0].split('>')[1].rsplit('_%s' % (ID), 1)[0]
            else:
                contig = seq[0].split()[0].split('>')[1].rsplit('_', 1)[0]
            locus = seq[0].split()[0].split('>')[1]
            # annotation info from Prodigal
            if ('# ;gc_cont' in seq[0] or '# ID=' in seq[0]):
                info = seq[0].split(' # ')
                start, end, strand = int(info[1]), int(info[2]), info[3]
                if strand == '1':
                    strand = 1
                else:
                    strand = -1
                product = [''.join(info[4].split()[1:])]
            # annotation info from table
            else:
                start, end, strand = annots[locus]
                product = seq[0].split(' ', 1)[1]
            info = {'transl_table':[trans_table], \
                    'translation':[seq[1]], \
                    'product':product}
            yield contig, [locus, [start, end, strand], info]
886330
def fasta(args):
    
    from jcvi.formats.sizes import Sizes

    p = OptionParser(fasta.__doc__)
    p.add_option("--extend", default=1000, type="int",
                 help="Extend seq flanking the gaps [default: %default]")
    opts, args = p.parse_args(args)

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

    mapout, sfasta = args
    Flank = opts.extend
    pf = mapout.split(".")[0]
    mapbed = pf + ".bed"
    bm = BinMap(mapout)
    bm.print_to_bed(mapbed)

    bed = Bed(mapbed, sorted=False)
    markersbed = pf + ".markers.bed"
    fw = open(markersbed, "w")
    sizes = Sizes(sfasta).mapping
    for b in bed:
        accn = b.accn
        scf, pos = accn.split(".")
        pos = int(pos)
        start = max(0, pos - Flank)
        end = min(pos + Flank, sizes[scf])
        print("\t".join(str(x) for x in \
                    (scf, start, end, accn)), file=fw)

    fw.close()

    fastaFromBed(markersbed, sfasta, name=True)