changeset 0:887e111c0919 draft

planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/proteogenomics/retrieve_ensembl_bed commit 3fd7be931712e7fa5b281bc8c48104c8583ef7f0
author galaxyp
date Sun, 14 Jan 2018 14:11:53 -0500
parents
children c3d600729b6f
files bedutil.py ensembl_rest.py macros.xml retrieve_ensembl_bed.py retrieve_ensembl_bed.xml
diffstat 5 files changed, 1011 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bedutil.py	Sun Jan 14 14:11:53 2018 -0500
@@ -0,0 +1,503 @@
+#!/usr/bin/env python
+"""
+#
+#------------------------------------------------------------------------------
+#                         University of Minnesota
+#         Copyright 2016, Regents of the University of Minnesota
+#------------------------------------------------------------------------------
+# Author:
+#
+#  James E Johnson
+#
+#------------------------------------------------------------------------------
+"""
+
+import sys
+
+from Bio.Seq import reverse_complement, translate
+
+
+def bed_from_line(line, ensembl=False, seq_column=None):
+    fields = line.rstrip('\r\n').split('\t')
+    if len(fields) < 12:
+        return None
+    (chrom, chromStart, chromEnd, name, score, strand,
+     thickStart, thickEnd, itemRgb,
+     blockCount, blockSizes, blockStarts) = fields[0:12]
+    bed_entry = BedEntry(chrom=chrom, chromStart=chromStart, chromEnd=chromEnd,
+                         name=name, score=score, strand=strand,
+                         thickStart=thickStart, thickEnd=thickEnd,
+                         itemRgb=itemRgb,
+                         blockCount=blockCount,
+                         blockSizes=blockSizes.rstrip(','),
+                         blockStarts=blockStarts.rstrip(','))
+    if seq_column is not None and -len(fields) <= seq_column < len(fields):
+        bed_entry.seq = fields[seq_column]
+    if ensembl and len(fields) >= 20:
+        bed_entry.second_name = fields[12]
+        bed_entry.cds_start_status = fields[13]
+        bed_entry.cds_end_status = fields[14]
+        bed_entry.exon_frames = fields[15].rstrip(',')
+        bed_entry.biotype = fields[16]
+        bed_entry.gene_name = fields[17]
+        bed_entry.second_gene_name = fields[18]
+        bed_entry.gene_type = fields[19]
+    return bed_entry
+
+
+class BedEntry(object):
+    def __init__(self, chrom=None, chromStart=None, chromEnd=None,
+                 name=None, score=None, strand=None,
+                 thickStart=None, thickEnd=None, itemRgb=None,
+                 blockCount=None, blockSizes=None, blockStarts=None):
+        self.chrom = chrom
+        self.chromStart = int(chromStart)
+        self.chromEnd = int(chromEnd)
+        self.name = name
+        self.score = int(score) if score is not None else 0
+        self.strand = '-' if str(strand).startswith('-') else '+'
+        self.thickStart = int(thickStart) if thickStart else self.chromStart
+        self.thickEnd = int(thickEnd) if thickEnd else self.chromEnd
+        self.itemRgb = str(itemRgb) if itemRgb is not None else r'100,100,100'
+        self.blockCount = int(blockCount)
+        if isinstance(blockSizes, str) or isinstance(blockSizes, unicode):
+            self.blockSizes = [int(x) for x in blockSizes.split(',')]
+        elif isinstance(blockSizes, list):
+            self.blockSizes = [int(x) for x in blockSizes]
+        else:
+            self.blockSizes = blockSizes
+        if isinstance(blockStarts, str) or isinstance(blockSizes, unicode):
+            self.blockStarts = [int(x) for x in blockStarts.split(',')]
+        elif isinstance(blockStarts, list):
+            self.blockStarts = [int(x) for x in blockStarts]
+        else:
+            self.blockStarts = blockStarts
+        self.second_name = None
+        self.cds_start_status = None
+        self.cds_end_status = None
+        self.exon_frames = None
+        self.biotype = None
+        self.gene_name = None
+        self.second_gene_name = None
+        self.gene_type = None
+        self.seq = None
+        self.cdna = None
+        self.pep = None
+        # T26C
+        self.aa_change = []
+        # p.Trp26Cys g.<pos><ref>><alt> # g.1304573A>G
+        self.variants = []
+
+    def __str__(self):
+        return '%s\t%d\t%d\t%s\t%d\t%s\t%d\t%d\t%s\t%d\t%s\t%s' % (
+            self.chrom, self.chromStart, self.chromEnd,
+            self.name, self.score, self.strand,
+            self.thickStart, self.thickEnd, str(self.itemRgb), self.blockCount,
+            ','.join([str(x) for x in self.blockSizes]),
+            ','.join([str(x) for x in self.blockStarts]))
+
+    def get_splice_junctions(self):
+        splice_juncs = []
+        for i in range(self.blockCount - 1):
+            splice_junc = "%s:%d_%d"\
+                % (self.chrom,
+                   self.chromStart + self.blockSizes[i],
+                   self.chromStart + self.blockStarts[i+1])
+            splice_juncs.append(splice_junc)
+        return splice_juncs
+
+    def get_exon_seqs(self):
+        if not self.seq:
+            return None
+        exons = []
+        for i in range(self.blockCount):
+            exons.append(self.seq[self.blockStarts[i]:self.blockStarts[i]
+                         + self.blockSizes[i]])
+        if self.strand == '-':  # reverse complement
+            exons.reverse()
+            for i, s in enumerate(exons):
+                exons[i] = reverse_complement(s)
+        return exons
+
+    def get_spliced_seq(self, strand=None):
+        if not self.seq:
+            return None
+        seq = ''.join(self.get_exon_seqs())
+        if strand and self.strand != strand:
+            seq = reverse_complement(seq)
+        return seq
+
+    def get_cdna(self):
+        if not self.cdna:
+            self.cdna = self.get_spliced_seq()
+        return self.cdna
+
+    def get_cds(self):
+        cdna = self.get_cdna()
+        if cdna:
+            if self.chromStart == self.thickStart\
+               and self.chromEnd == self.thickEnd:
+                return cdna
+            pos = [self.cdna_offset_of_pos(self.thickStart),
+                   self.cdna_offset_of_pos(self.thickEnd)]
+            if 0 <= min(pos) <= max(pos) <= len(cdna):
+                return cdna[min(pos):max(pos)]
+        return None
+
+    def set_cds(self, cdna_start, cdna_end):
+        cdna_len = sum(self.blockSizes)
+        if 0 <= cdna_start < cdna_end <= cdna_len:
+            cds_pos = [self.pos_of_cdna_offet(cdna_start),
+                       self.pos_of_cdna_offet(cdna_end)]
+            if all(cds_pos):
+                self.thickStart = min(cds_pos)
+                self.thickEnd = max(cds_pos)
+            return self
+        return None
+
+    def trim_cds(self, basepairs):
+        if self.chromStart <= self.thickStart < self.thickEnd <= self.chromEnd:
+            cds_pos = [self.cdna_offset_of_pos(self.thickStart),
+                       self.cdna_offset_of_pos(self.thickEnd)]
+            if basepairs > 0:
+                return self.set_cds(min(cds_pos) + basepairs, max(cds_pos))
+            else:
+                return self.set_cds(min(cds_pos), max(cds_pos) + basepairs)
+        return None
+
+    def get_cigar(self):
+        cigar = ''
+        r = range(self.blockCount)
+        xl = None
+        for x in r:
+            if xl is not None:
+                intronSize = abs(self.blockStarts[x] - self.blockSizes[xl]
+                                 - self.blockStarts[xl])
+                cigar += '%dN' % intronSize
+            cigar += '%dM' % self.blockSizes[x]
+            xl = x
+        return cigar
+
+    def get_cigar_md(self):
+        cigar = ''
+        md = ''
+        r = range(self.blockCount)
+        xl = None
+        for x in r:
+            if xl is not None:
+                intronSize = abs(self.blockStarts[x] - self.blockSizes[xl]
+                                 - self.blockStarts[xl])
+                cigar += '%dN' % intronSize
+            cigar += '%dM' % self.blockSizes[x]
+            xl = x
+        md = '%d' % sum(self.blockSizes)
+        return (cigar, md)
+
+    def get_translation(self, sequence=None):
+        translation = None
+        seq = sequence if sequence else self.get_spliced_seq()
+        if seq:
+            seqlen = len(seq) / 3 * 3
+            if seqlen >= 3:
+                translation = translate(seq[:seqlen])
+        return translation
+
+    def get_translations(self):
+        translations = []
+        seq = self.get_spliced_seq()
+        if seq:
+            for i in range(3):
+                translation = self.get_translation(sequence=seq[i:])
+                if translation:
+                    translations.append(translation)
+        return translations
+
+    def pos_of_cdna_offet(self, offset):
+        if offset is not None and 0 <= offset < sum(self.blockSizes):
+            r = range(self.blockCount)
+            rev = self.strand == '-'
+            if rev:
+                r.reverse()
+            nlen = 0
+            for x in r:
+                if offset < nlen + self.blockSizes[x]:
+                    if rev:
+                        return self.chromStart + self.blockStarts[x]\
+                               + self.blockSizes[x] - (offset - nlen)
+                    else:
+                        return self.chromStart + self.blockStarts[x]\
+                               + (offset - nlen)
+                nlen += self.blockSizes[x]
+        return None
+
+    def cdna_offset_of_pos(self, pos):
+        if not self.chromStart <= pos < self.chromEnd:
+            return -1
+        r = range(self.blockCount)
+        rev = self.strand == '-'
+        if rev:
+            r.reverse()
+        nlen = 0
+        for x in r:
+            bStart = self.chromStart + self.blockStarts[x]
+            bEnd = bStart + self.blockSizes[x]
+            if bStart <= pos < bEnd:
+                return nlen + (bEnd - pos if rev else pos - bStart)
+            nlen += self.blockSizes[x]
+
+    def apply_variant(self, pos, ref, alt):
+        pos = int(pos)
+        if not ref or not alt:
+            print >> sys.stderr, "variant requires ref and alt sequences"
+            return
+        if not self.chromStart <= pos <= self.chromEnd:
+            print >> sys.stderr, "variant not in entry %s: %s %d < %d < %d"\
+                % (self.name, self.strand, self.chromStart, pos, self.chromEnd)
+            print >> sys.stderr, "%s" % str(self)
+            return
+        if len(ref) != len(alt):
+            print >> sys.stderr, "variant only works for snp: %s  %s"\
+                % (ref, alt)
+            return
+        if not self.seq:
+            print >> sys.stderr, "variant entry %s has no seq" % self.name
+            return
+        """
+        if self.strand  == '-':
+            ref = reverse_complement(ref)
+            alt = reverse_complement(alt)
+        """
+        bases = list(self.seq)
+        offset = pos - self.chromStart
+        for i in range(len(ref)):
+            # offset = self.cdna_offset_of_pos(pos+i)
+            if offset is not None:
+                bases[offset+i] = alt[i]
+            else:
+                print >> sys.stderr,\
+                    "variant offset %s: %s %d < %d < %d"\
+                    % (self.name, self.strand, self.chromStart,
+                       pos+1, self.chromEnd)
+                print >> sys.stderr, "%s" % str(self)
+        self.seq = ''.join(bases)
+        self.variants.append("g.%d%s>%s" % (pos+1, ref, alt))
+
+    def get_variant_bed(self, pos, ref, alt):
+        pos = int(pos)
+        if not ref or not alt:
+            print >> sys.stderr, "variant requires ref and alt sequences"
+            return None
+        if not self.chromStart <= pos <= self.chromEnd:
+            print >> sys.stderr,\
+                "variant not in entry %s: %s %d < %d < %d"\
+                % (self.name, self.strand, self.chromStart, pos, self.chromEnd)
+            print >> sys.stderr, "%s" % str(self)
+            return None
+        if not self.seq:
+            print >> sys.stderr, "variant entry %s has no seq" % self.name
+            return None
+        tbed = BedEntry(chrom=self.chrom,
+                        chromStart=self.chromStart, chromEnd=self.chromEnd,
+                        name=self.name, score=self.score, strand=self.strand,
+                        thickStart=self.chromStart, thickEnd=self.chromEnd,
+                        itemRgb=self.itemRgb,
+                        blockCount=self.blockCount,
+                        blockSizes=self.blockSizes,
+                        blockStarts=self.blockStarts)
+        bases = list(self.seq)
+        offset = pos - self.chromStart
+        tbed.seq = ''.join(bases[:offset] + list(alt)
+                           + bases[offset+len(ref):])
+        if len(ref) != len(alt):
+            diff = len(alt) - len(ref)
+            rEnd = pos + len(ref)
+            # need to adjust blocks
+            #  change spans blocks,
+            for x in range(tbed.blockCount):
+                bStart = tbed.chromStart + tbed.blockStarts[x]
+                bEnd = bStart + tbed.blockSizes[x]
+                # change within a block or extends (last block)
+                #  adjust blocksize
+                #  seq:            GGGcatGGG
+                #  ref c alt tag:  GGGtagatGGG
+                #  ref cat alt a:  GGGaGGG
+                if bStart <= pos < rEnd < bEnd:
+                    tbed.blockSizes[x] += diff
+        return tbed
+
+    # (start, end)
+    def get_subrange(self, tstart, tstop, debug=False):
+        chromStart = self.chromStart
+        chromEnd = self.chromEnd
+        if debug:
+            print >> sys.stderr, "%s" % (str(self))
+        r = range(self.blockCount)
+        if self.strand == '-':
+            r.reverse()
+        bStart = 0
+        bEnd = 0
+        for x in r:
+            bEnd = bStart + self.blockSizes[x]
+            if bStart <= tstart < bEnd:
+                if self.strand == '+':
+                    chromStart = self.chromStart + self.blockStarts[x] +\
+                        (tstart - bStart)
+                else:
+                    chromEnd = self.chromStart + self.blockStarts[x] +\
+                        self.blockSizes[x] - (tstart - bStart)
+            if bStart <= tstop < bEnd:
+                if self.strand == '+':
+                    chromEnd = self.chromStart + self.blockStarts[x] +\
+                        (tstop - bStart)
+                else:
+                    chromStart = self.chromStart + self.blockStarts[x] +\
+                        self.blockSizes[x] - (tstop - bStart)
+            if debug:
+                print >> sys.stderr,\
+                    "%3d %s\t%d\t%d\t%d\t%d\t%d\t%d"\
+                    % (x, self.strand, bStart, bEnd,
+                       tstart, tstop, chromStart, chromEnd)
+            bStart += self.blockSizes[x]
+        return(chromStart, chromEnd)
+
+    # get the blocks for sub range
+    def get_blocks(self, chromStart, chromEnd):
+        tblockCount = 0
+        tblockSizes = []
+        tblockStarts = []
+        for x in range(self.blockCount):
+            bStart = self.chromStart + self.blockStarts[x]
+            bEnd = bStart + self.blockSizes[x]
+            if bStart > chromEnd:
+                break
+            if bEnd < chromStart:
+                continue
+            cStart = max(chromStart, bStart)
+            tblockStarts.append(cStart - chromStart)
+            tblockSizes.append(min(chromEnd, bEnd) - cStart)
+            tblockCount += 1
+        return (tblockCount, tblockSizes, tblockStarts)
+
+    def trim(self, tstart, tstop, debug=False):
+        (tchromStart, tchromEnd) =\
+            self.get_subrange(tstart, tstop, debug=debug)
+        (tblockCount, tblockSizes, tblockStarts) =\
+            self.get_blocks(tchromStart, tchromEnd)
+        tbed = BedEntry(
+            chrom=self.chrom, chromStart=tchromStart, chromEnd=tchromEnd,
+            name=self.name, score=self.score, strand=self.strand,
+            thickStart=tchromStart, thickEnd=tchromEnd, itemRgb=self.itemRgb,
+            blockCount=tblockCount,
+            blockSizes=tblockSizes, blockStarts=tblockStarts)
+        if self.seq:
+            ts = tchromStart-self.chromStart
+            te = tchromEnd - tchromStart + ts
+            tbed.seq = self.seq[ts:te]
+        return tbed
+
+    def get_filtered_translations(self, untrimmed=False, filtering=True,
+                                  ignore_left_bp=0, ignore_right_bp=0,
+                                  debug=False):
+        translations = [None, None, None]
+        seq = self.get_spliced_seq()
+        ignore = (ignore_left_bp if self.strand == '+'
+                  else ignore_right_bp) / 3
+        block_sum = sum(self.blockSizes)
+        exon_sizes = [x for x in self.blockSizes]
+        if self.strand == '-':
+            exon_sizes.reverse()
+        splice_sites = [sum(exon_sizes[:x]) / 3
+                        for x in range(1, len(exon_sizes))]
+        if debug:
+            print >> sys.stderr, "splice_sites: %s" % splice_sites
+        junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0]
+        if seq:
+            for i in range(3):
+                translation = self.get_translation(sequence=seq[i:])
+                if translation:
+                    tstart = 0
+                    tstop = len(translation)
+                    offset = (block_sum - i) % 3
+                    if debug:
+                        print >> sys.stderr,\
+                            "frame: %d\ttstart: %d  tstop: %d  offset: %d\t%s"\
+                            % (i, tstart, tstop, offset, translation)
+                    if not untrimmed:
+                        tstart = translation.rfind('*', 0, junc) + 1
+                        stop = translation.find('*', junc)
+                        tstop = stop if stop >= 0 else len(translation)
+                    offset = (block_sum - i) % 3
+                    trimmed = translation[tstart:tstop]
+                    if debug:
+                        print >> sys.stderr,\
+                            "frame: %d\ttstart: %d  tstop: %d  offset: %d\t%s"\
+                            % (i, tstart, tstop, offset, trimmed)
+                    if filtering and tstart > ignore:
+                        continue
+                    # get genomic locations for start and end
+                    if self.strand == '+':
+                        chromStart = self.chromStart + i + (tstart * 3)
+                        chromEnd = self.chromEnd - offset\
+                            - (len(translation) - tstop) * 3
+                    else:
+                        chromStart = self.chromStart + offset\
+                            + (len(translation) - tstop) * 3
+                        chromEnd = self.chromEnd - i - (tstart * 3)
+                    # get the blocks for this translation
+                    (tblockCount, tblockSizes, tblockStarts) =\
+                        self.get_blocks(chromStart, chromEnd)
+                    translations[i] = (chromStart, chromEnd, trimmed,
+                                       tblockCount, tblockSizes, tblockStarts)
+                    if debug:
+                        print >> sys.stderr,\
+                            "tblockCount: %d tblockStarts: %s tblockSizes: %s"\
+                            % (tblockCount, tblockStarts, tblockSizes)
+        return translations
+
+    def get_seq_id(self, seqtype='unk:unk', reference='', frame=None):
+        # Ensembl fasta ID format
+        # >ID SEQTYPE:STATUS LOCATION GENE TRANSCRIPT
+        # >ENSP00000328693 pep:splice chromosome:NCBI35:1:904515:910768:1\
+        #   gene:ENSG00000158815:transcript:ENST00000328693\
+        #    gene_biotype:protein_coding transcript_biotype:protein_coding
+        frame_name = ''
+        chromStart = self.chromStart
+        chromEnd = self.chromEnd
+        strand = 1 if self.strand == '+' else -1
+        if frame is not None:
+            block_sum = sum(self.blockSizes)
+            offset = (block_sum - frame) % 3
+            frame_name = '_' + str(frame + 1)
+            if self.strand == '+':
+                chromStart += frame
+                chromEnd -= offset
+            else:
+                chromStart += offset
+                chromEnd -= frame
+        location = "chromosome:%s:%s:%s:%s:%s"\
+            % (reference, self.chrom, chromStart, chromEnd, strand)
+        seq_id = "%s%s %s %s" % (self.name, frame_name, seqtype, location)
+        return seq_id
+
+    def get_line(self, start_offset=0, end_offset=0):
+        if start_offset or end_offset:
+            s_offset = start_offset if start_offset else 0
+            e_offset = end_offset if end_offset else 0
+            if s_offset > self.chromStart:
+                s_offset = self.chromStart
+            chrStart = self.chromStart - s_offset
+            chrEnd = self.chromEnd + e_offset
+            blkSizes = self.blockSizes
+            blkSizes[0] += s_offset
+            blkSizes[-1] += e_offset
+            blkStarts = self.blockStarts
+            for i in range(1, self.blockCount):
+                blkStarts[i] += s_offset
+            items = [str(x) for x in [self.chrom, chrStart, chrEnd, self.name,
+                                      self.score, self.strand, self.thickStart,
+                                      self.thickEnd, self.itemRgb,
+                                      self.blockCount,
+                                      ','.join([str(x) for x in blkSizes]),
+                                      ','.join([str(x) for x in blkStarts])]]
+            return '\t'.join(items) + '\n'
+        return self.line
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ensembl_rest.py	Sun Jan 14 14:11:53 2018 -0500
@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+"""
+#
+#------------------------------------------------------------------------------
+#                         University of Minnesota
+#         Copyright 2017, Regents of the University of Minnesota
+#------------------------------------------------------------------------------
+# Author:
+#
+#  James E Johnson
+#
+#------------------------------------------------------------------------------
+"""
+
+
+import sys
+
+from time import sleep
+
+import requests
+
+
+server = "https://rest.ensembl.org"
+ext = "/info/assembly/homo_sapiens?"
+max_region = 4000000
+debug = False
+
+
+def ensembl_rest(ext, headers):
+    if debug:
+        print >> sys.stderr, "%s" % ext
+    r = requests.get(server+ext, headers=headers)
+    if r.status_code == 429:
+        print >> sys.stderr, "response headers: %s\n" % r.headers
+        if 'Retry-After' in r.headers:
+            sleep(r.headers['Retry-After'])
+            r = requests.get(server+ext, headers=headers)
+    if not r.ok:
+        r.raise_for_status()
+    return r
+
+
+def get_species():
+    results = dict()
+    ext = "/info/species"
+    req_header = {"Content-Type": "application/json"}
+    r = ensembl_rest(ext, req_header)
+    for species in r.json()['species']:
+        results[species['name']] = species
+        print >> sys.stdout,\
+            "%s\t%s\t%s\t%s\t%s"\
+            % (species['name'], species['common_name'],
+               species['display_name'],
+               species['strain'],
+               species['taxon_id'])
+    return results
+
+
+def get_biotypes(species):
+    biotypes = []
+    ext = "/info/biotypes/%s?" % species
+    req_header = {"Content-Type": "application/json"}
+    r = ensembl_rest(ext, req_header)
+    for entry in r.json():
+        if 'biotype' in entry:
+            biotypes.append(entry['biotype'])
+    return biotypes
+
+
+def get_toplevel(species):
+    coord_systems = dict()
+    ext = "/info/assembly/%s?" % species
+    req_header = {"Content-Type": "application/json"}
+    r = ensembl_rest(ext, req_header)
+    toplevel = r.json()
+    for seq in toplevel['top_level_region']:
+        if seq['coord_system'] not in coord_systems:
+            coord_systems[seq['coord_system']] = dict()
+        coord_system = coord_systems[seq['coord_system']]
+        coord_system[seq['name']] = int(seq['length'])
+    return coord_systems
+
+
+def get_transcripts_bed(species, refseq, start, length, strand='',
+                        params=None):
+    bed = []
+    param = params if params else ''
+    req_header = {"Content-Type": "text/x-bed"}
+    regions = range(start, length, max_region)
+    if not regions or regions[-1] < length:
+        regions.append(length)
+    for end in regions[1:]:
+        ext = "/overlap/region/%s/%s:%d-%d%s?feature=transcript;%s"\
+            % (species, refseq, start, end, strand, param)
+        start = end + 1
+        r = ensembl_rest(ext, req_header)
+        if r.text:
+            bed += r.text.splitlines()
+    return bed
+
+
+def get_seq(id, seqtype, params=None):
+    param = params if params else ''
+    ext = "/sequence/id/%s?type=%s;%s" % (id, seqtype, param)
+    req_header = {"Content-Type": "text/plain"}
+    r = ensembl_rest(ext, req_header)
+    return r.text
+
+
+def get_cdna(id, params=None):
+    return get_seq(id, 'cdna', params=params)
+
+
+def get_cds(id, params=None):
+    return get_seq(id, 'cds', params=params)
+
+
+def get_genomic(id, params=None):
+    return get_seq(id, 'genomic', params=params)
+
+
+def get_transcript_haplotypes(species, transcript):
+    ext = "/transcript_haplotypes/%s/%s?aligned_sequences=1"\
+        % (species, transcript)
+    req_header = {"Content-Type": "application/json"}
+    r = ensembl_rest(ext, req_header)
+    decoded = r.json()
+    return decoded
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/macros.xml	Sun Jan 14 14:11:53 2018 -0500
@@ -0,0 +1,125 @@
+<macros>
+    <xml name="bedutil_requirements">
+        <requirement type="package" version="1.62">biopython</requirement>
+    </xml>
+    <xml name="ensembl_requirements">
+        <requirement type="package" version="0.4.10">requests-cache</requirement>
+    </xml>
+    <xml name="twobit_requirements">
+        <requirement type="package" version="3.1.4">twobitreader</requirement>
+    </xml>
+    <xml name="species_options">
+            <option value="homo_sapiens">homo_sapiens  (Human) taxon_id: 9606</option>
+            <option value="mus_musculus">mus_musculus  (Mouse) taxon_id: 10090</option>
+            <option value="ailuropoda_melanoleuca">ailuropoda_melanoleuca  (Panda) taxon_id: 9646</option>
+            <option value="anas_platyrhynchos">anas_platyrhynchos  (Duck) taxon_id: 8839</option>
+            <option value="anolis_carolinensis">anolis_carolinensis  (Anole lizard) taxon_id: 28377</option>
+            <option value="astyanax_mexicanus">astyanax_mexicanus  (Cave fish) taxon_id: 7994</option>
+            <option value="bos_taurus">bos_taurus  (Cow) taxon_id: 9913</option>
+            <option value="caenorhabditis_elegans">caenorhabditis_elegans  (Caenorhabditis elegans) taxon_id: 6239</option>
+            <option value="callithrix_jacchus">callithrix_jacchus  (Marmoset) taxon_id: 9483</option>
+            <option value="canis_familiaris">canis_familiaris  (Dog) taxon_id: 9615</option>
+            <option value="carlito_syrichta">carlito_syrichta  (Tarsier) taxon_id: 1868482</option>
+            <option value="cavia_aperea">cavia_aperea  (Brazilian guinea pig) taxon_id: 37548</option>
+            <option value="cavia_porcellus">cavia_porcellus  (Guinea Pig) taxon_id: 10141</option>
+            <option value="chinchilla_lanigera">chinchilla_lanigera  (Long-tailed chinchilla) taxon_id: 34839</option>
+            <option value="chlorocebus_sabaeus">chlorocebus_sabaeus  (Vervet-AGM) taxon_id: 60711</option>
+            <option value="choloepus_hoffmanni">choloepus_hoffmanni  (Sloth) taxon_id: 9358</option>
+            <option value="ciona_intestinalis">ciona_intestinalis  (C.intestinalis) taxon_id: 7719</option>
+            <option value="ciona_savignyi">ciona_savignyi  (C.savignyi) taxon_id: 51511</option>
+            <option value="cricetulus_griseus_chok1gshd">cricetulus_griseus_chok1gshd  (Chinese hamster CHOK1GS) taxon_id: 10029</option>
+            <option value="cricetulus_griseus_crigri">cricetulus_griseus_crigri  (Chinese hamster CriGri) taxon_id: 10029</option>
+            <option value="danio_rerio">danio_rerio  (Zebrafish) taxon_id: 7955</option>
+            <option value="dasypus_novemcinctus">dasypus_novemcinctus  (Armadillo) taxon_id: 9361</option>
+            <option value="dipodomys_ordii">dipodomys_ordii  (Kangaroo rat) taxon_id: 10020</option>
+            <option value="drosophila_melanogaster">drosophila_melanogaster  (Fruitfly) taxon_id: 7227</option>
+            <option value="echinops_telfairi">echinops_telfairi  (Lesser hedgehog tenrec) taxon_id: 9371</option>
+            <option value="equus_caballus">equus_caballus  (Horse) taxon_id: 9796</option>
+            <option value="erinaceus_europaeus">erinaceus_europaeus  (Hedgehog) taxon_id: 9365</option>
+            <option value="felis_catus">felis_catus  (Cat) taxon_id: 9685</option>
+            <option value="ficedula_albicollis">ficedula_albicollis  (Flycatcher) taxon_id: 59894</option>
+            <option value="fukomys_damarensis">fukomys_damarensis  (Damara mole rat) taxon_id: 885580</option>
+            <option value="gadus_morhua">gadus_morhua  (Cod) taxon_id: 8049</option>
+            <option value="gallus_gallus">gallus_gallus  (Chicken) taxon_id: 9031</option>
+            <option value="gasterosteus_aculeatus">gasterosteus_aculeatus  (Stickleback) taxon_id: 69293</option>
+            <option value="gorilla_gorilla">gorilla_gorilla  (Gorilla) taxon_id: 9595</option>
+            <option value="heterocephalus_glaber_female">heterocephalus_glaber_female  (Naked mole-rat female) taxon_id: 10181</option>
+            <option value="heterocephalus_glaber_male">heterocephalus_glaber_male  (Naked mole-rat male) taxon_id: 10181</option>
+            <option value="ictidomys_tridecemlineatus">ictidomys_tridecemlineatus  (Squirrel) taxon_id: 43179</option>
+            <option value="jaculus_jaculus">jaculus_jaculus  (Lesser Egyptian jerboa) taxon_id: 51337</option>
+            <option value="latimeria_chalumnae">latimeria_chalumnae  (Coelacanth) taxon_id: 7897</option>
+            <option value="lepisosteus_oculatus">lepisosteus_oculatus  (Spotted gar) taxon_id: 7918</option>
+            <option value="loxodonta_africana">loxodonta_africana  (Elephant) taxon_id: 9785</option>
+            <option value="macaca_mulatta">macaca_mulatta  (Macaque) taxon_id: 9544</option>
+            <option value="meleagris_gallopavo">meleagris_gallopavo  (Turkey) taxon_id: 9103</option>
+            <option value="mesocricetus_auratus">mesocricetus_auratus  (Golden Hamster) taxon_id: 10036</option>
+            <option value="microcebus_murinus">microcebus_murinus  (Mouse Lemur) taxon_id: 30608</option>
+            <option value="microtus_ochrogaster">microtus_ochrogaster  (Prairie vole) taxon_id: 79684</option>
+            <option value="monodelphis_domestica">monodelphis_domestica  (Opossum) taxon_id: 13616</option>
+            <option value="mus_caroli">mus_caroli  (Ryukyu mouse) taxon_id: 10089</option>
+            <option value="mus_musculus_129s1svimj">mus_musculus_129s1svimj  (Mouse 129S1/SvImJ) taxon_id: 10090</option>
+            <option value="mus_musculus_aj">mus_musculus_aj  (Mouse A/J) taxon_id: 10090</option>
+            <option value="mus_musculus_akrj">mus_musculus_akrj  (Mouse AKR/J) taxon_id: 10090</option>
+            <option value="mus_musculus_balbcj">mus_musculus_balbcj  (Mouse BALB/cJ) taxon_id: 10090</option>
+            <option value="mus_musculus_c3hhej">mus_musculus_c3hhej  (Mouse C3H/HeJ) taxon_id: 10090</option>
+            <option value="mus_musculus_c57bl6nj">mus_musculus_c57bl6nj  (Mouse C57BL/6NJ) taxon_id: 10090</option>
+            <option value="mus_musculus_casteij">mus_musculus_casteij  (Mouse CAST/EiJ) taxon_id: 10091</option>
+            <option value="mus_musculus_cbaj">mus_musculus_cbaj  (Mouse CBA/J) taxon_id: 10090</option>
+            <option value="mus_musculus_dba2j">mus_musculus_dba2j  (Mouse DBA/2J) taxon_id: 10090</option>
+            <option value="mus_musculus_fvbnj">mus_musculus_fvbnj  (Mouse FVB/NJ) taxon_id: 10090</option>
+            <option value="mus_musculus_lpj">mus_musculus_lpj  (Mouse LP/J) taxon_id: 10090</option>
+            <option value="mus_musculus_nodshiltj">mus_musculus_nodshiltj  (Mouse NOD/ShiLtJ) taxon_id: 10090</option>
+            <option value="mus_musculus_nzohlltj">mus_musculus_nzohlltj  (Mouse NZO/HlLtJ) taxon_id: 10090</option>
+            <option value="mus_musculus_pwkphj">mus_musculus_pwkphj  (Mouse PWK/PhJ) taxon_id: 39442</option>
+            <option value="mus_musculus_wsbeij">mus_musculus_wsbeij  (Mouse WSB/EiJ) taxon_id: 10092</option>
+            <option value="mus_pahari">mus_pahari  (Shrew mouse) taxon_id: 10093</option>
+            <option value="mus_spretus_spreteij">mus_spretus_spreteij  (Algerian mouse) taxon_id: 10096</option>
+            <option value="mustela_putorius_furo">mustela_putorius_furo  (Ferret) taxon_id: 9669</option>
+            <option value="myotis_lucifugus">myotis_lucifugus  (Microbat) taxon_id: 59463</option>
+            <option value="nannospalax_galili">nannospalax_galili  (Upper Galilee mountains blind mole rat) taxon_id: 1026970</option>
+            <option value="nomascus_leucogenys">nomascus_leucogenys  (Gibbon) taxon_id: 61853</option>
+            <option value="notamacropus_eugenii">notamacropus_eugenii  (Wallaby) taxon_id: 9315</option>
+            <option value="ochotona_princeps">ochotona_princeps  (Pika) taxon_id: 9978</option>
+            <option value="octodon_degus">octodon_degus  (Degu) taxon_id: 10160</option>
+            <option value="oreochromis_niloticus">oreochromis_niloticus  (Tilapia) taxon_id: 8128</option>
+            <option value="ornithorhynchus_anatinus">ornithorhynchus_anatinus  (Platypus) taxon_id: 9258</option>
+            <option value="oryctolagus_cuniculus">oryctolagus_cuniculus  (Rabbit) taxon_id: 9986</option>
+            <option value="oryzias_latipes">oryzias_latipes  (Medaka) taxon_id: 8090</option>
+            <option value="otolemur_garnettii">otolemur_garnettii  (Bushbaby) taxon_id: 30611</option>
+            <option value="ovis_aries">ovis_aries  (Sheep) taxon_id: 9940</option>
+            <option value="pan_troglodytes">pan_troglodytes  (Chimpanzee) taxon_id: 9598</option>
+            <option value="papio_anubis">papio_anubis  (Olive baboon) taxon_id: 9555</option>
+            <option value="pelodiscus_sinensis">pelodiscus_sinensis  (Chinese softshell turtle) taxon_id: 13735</option>
+            <option value="peromyscus_maniculatus_bairdii">peromyscus_maniculatus_bairdii  (Northern American deer mouse) taxon_id: 230844</option>
+            <option value="petromyzon_marinus">petromyzon_marinus  (Lamprey) taxon_id: 7757</option>
+            <option value="poecilia_formosa">poecilia_formosa  (Amazon molly) taxon_id: 48698</option>
+            <option value="pongo_abelii">pongo_abelii  (Orangutan) taxon_id: 9601</option>
+            <option value="procavia_capensis">procavia_capensis  (Hyrax) taxon_id: 9813</option>
+            <option value="pteropus_vampyrus">pteropus_vampyrus  (Megabat) taxon_id: 132908</option>
+            <option value="rattus_norvegicus">rattus_norvegicus  (Rat) taxon_id: 10116</option>
+            <option value="saccharomyces_cerevisiae">saccharomyces_cerevisiae  (Saccharomyces cerevisiae) taxon_id: 4932</option>
+            <option value="sarcophilus_harrisii">sarcophilus_harrisii  (Tasmanian devil) taxon_id: 9305</option>
+            <option value="sorex_araneus">sorex_araneus  (Shrew) taxon_id: 42254</option>
+            <option value="sus_scrofa">sus_scrofa  (Pig) taxon_id: 9823</option>
+            <option value="taeniopygia_guttata">taeniopygia_guttata  (Zebra Finch) taxon_id: 59729</option>
+            <option value="takifugu_rubripes">takifugu_rubripes  (Fugu) taxon_id: 31033</option>
+            <option value="tetraodon_nigroviridis">tetraodon_nigroviridis  (Tetraodon) taxon_id: 99883</option>
+            <option value="tupaia_belangeri">tupaia_belangeri  (Tree Shrew) taxon_id: 37347</option>
+            <option value="tursiops_truncatus">tursiops_truncatus  (Dolphin) taxon_id: 9739</option>
+            <option value="vicugna_pacos">vicugna_pacos  (Alpaca) taxon_id: 30538</option>
+            <option value="xenopus_tropicalis">xenopus_tropicalis  (Xenopus) taxon_id: 8364</option>
+            <option value="xiphophorus_maculatus">xiphophorus_maculatus  (Platyfish) taxon_id: 8083</option>
+    </xml>
+    <xml name="biotypes_help">
+            <help><![CDATA[
+Example biotypes: 
+protein_coding, non_coding, pseudogene, nonsense_mediated_decay, non_stop_decay, 
+translated_processed_pseudogene, transcribed_processed_pseudogene, transcribed_unitary_pseudogene, transcribed_unprocessed_pseudogene, 
+polymorphic_pseudogene, processed_pseudogene, unprocessed_pseudogene, unitary_pseudogene, processed_transcript, 
+retained_intron, ccds_gene, sense_overlapping, sense_intronic, cdna_update, antisense, 
+LRG_gene, IG_C_gene, IG_D_gene, IG_J_gene, IG_LV_gene IG_V_gene, TR_C_gene, TR_D_gene, TR_J_gene, TR_V_gene, 
+IG_pseudogene, IG_C_pseudogene, IG_D_pseudogene, IG_J_pseudogene, IG_V_pseudogene, TR_J_pseudogene, TR_V_pseudogene, TEC, 
+ribozyme, RNase_P_RNA, guide_RNA, macro_lncRNA, bidirectional_promoter_lncRNA, 3prime_overlapping_ncRNA, antisense_RNA, vaultRNA, Y_RNA, SRP_RNA, RNase_MRP_RNA, IG_C_pseudogene, lncRNA, lincRNA, miRNA, snRNA, sRNA, telomerase_RNA, Mt_tRNA, Mt_rRNA, scaRNA, misc_RNA, rRNA, tRNA, scRNA, snoRNA, other
+            ]]></help>
+    </xml>
+</macros>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/retrieve_ensembl_bed.py	Sun Jan 14 14:11:53 2018 -0500
@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+"""
+#
+#------------------------------------------------------------------------------
+#                         University of Minnesota
+#         Copyright 2017, Regents of the University of Minnesota
+#------------------------------------------------------------------------------
+# Author:
+#
+#  James E Johnson
+#
+#------------------------------------------------------------------------------
+"""
+
+import argparse
+import re
+import sys
+
+from bedutil import bed_from_line
+
+from ensembl_rest import get_toplevel, get_transcripts_bed, max_region
+
+
+def __main__():
+    parser = argparse.ArgumentParser(
+        description='Retrieve Ensembl cDNAs in BED format')
+    parser.add_argument(
+        '-s', '--species', default='human',
+        help='Ensembl Species to retrieve')
+    parser.add_argument(
+        '-R', '--regions', action='append', default=[],
+        help='Restrict Ensembl retrieval to regions e.g.:'
+             + ' X,2:20000-25000,3:100-500+')
+    parser.add_argument(
+        '-B', '--biotypes', action='append', default=[],
+        help='Restrict Ensembl biotypes to retrieve')
+    parser.add_argument(
+        '-X', '--extended_bed', action='store_true', default=False,
+        help='Include the extended columns returned from Ensembl')
+    parser.add_argument(
+        '-U', '--ucsc_chrom_names', action='store_true', default=False,
+        help='Use the UCSC names for Chromosomes')
+    parser.add_argument(
+        '-t', '--toplevel', action='store_true',
+        help='Print Ensembl toplevel for species')
+    parser.add_argument(
+        'output',
+        help='Output BED filepath, or for stdout: "-"')
+    parser.add_argument('-v', '--verbose', action='store_true', help='Verbose')
+    parser.add_argument('-d', '--debug', action='store_true', help='Debug')
+    args = parser.parse_args()
+    # print >> sys.stderr, "args: %s" % args
+    species = args.species
+    out_wtr = open(args.output, 'w') if args.output != '-' else sys.stdout
+    biotypes = ';'.join(['biotype=%s' % bt.strip()
+                         for biotype in args.biotypes
+                         for bt in biotype.split(',') if bt.strip()])
+
+    selected_regions = dict()  # chrom:(start, end)
+    region_pat = '^([^:]+)(?::(\d*)(?:-(\d+)([+-])?)?)?'
+    if args.regions:
+        for entry in args.regions:
+            if not entry:
+                continue
+            regs = [x.strip() for x in entry.split(',') if x.strip()]
+            for reg in regs:
+                m = re.match(region_pat, reg)
+                if m:
+                    (chrom, start, end, strand) = m.groups()
+                    if chrom:
+                        if chrom not in selected_regions:
+                            selected_regions[chrom] = []
+                        selected_regions[chrom].append([start, end, strand])
+        if args.debug:
+            print >> sys.stderr, "selected_regions: %s" % selected_regions
+
+    def retrieve_region(species, ref, start, stop, strand):
+        transcript_count = 0
+        regions = range(start, stop, max_region)
+        if not regions or regions[-1] < stop:
+            regions.append(stop)
+        for end in regions[1:]:
+            bedlines = get_transcripts_bed(species, ref, start, end,
+                                           strand=strand, params=biotypes)
+            if args.debug:
+                print >> sys.stderr,\
+                    "%s\t%s\tstart: %d\tend: %d\tcDNA transcripts:%d"\
+                    % (species, ref, start, end, len(bedlines))
+            # start, end, seq
+            for i, bedline in enumerate(bedlines):
+                if args.debug:
+                    print >> sys.stderr, "%s\n" % (bedline)
+                if not args.ucsc_chrom_names:
+                    bedline = re.sub('^[^\t]+', ref, bedline)
+                try:
+                    if out_wtr:
+                        out_wtr.write(bedline.replace(',\t', '\t')
+                                      if args.extended_bed
+                                      else str(bed_from_line(bedline)))
+                        out_wtr.write("\n")
+                        out_wtr.flush()
+                except Exception as e:
+                    print >> sys.stderr,\
+                        "BED error (%s) : %s\n" % (e, bedline)
+            start = end + 1
+        return transcript_count
+
+    coord_systems = get_toplevel(species)
+    if 'chromosome' in coord_systems:
+        ref_lengths = dict()
+        for ref in sorted(coord_systems['chromosome'].keys()):
+            length = coord_systems['chromosome'][ref]
+            ref_lengths[ref] = length
+            if args.toplevel:
+                print >> sys.stderr,\
+                    "%s\t%s\tlength: %d" % (species, ref, length)
+        if selected_regions:
+            transcript_count = 0
+            for ref in sorted(selected_regions.keys()):
+                if ref in ref_lengths:
+                    for reg in selected_regions[ref]:
+                        (_start, _stop, _strand) = reg
+                        start = int(_start) if _start else 0
+                        stop = int(_stop) if _stop else ref_lengths[ref]
+                        strand = '' if not _strand else ':1'\
+                            if _strand == '+' else ':-1'
+                        transcript_count += retrieve_region(species, ref,
+                                                            start, stop,
+                                                            strand)
+                        if args.debug or args.verbose:
+                            length = stop - start
+                            print >> sys.stderr,\
+                                "%s\t%s:%d-%d%s\tlength: %d\ttrancripts:%d"\
+                                % (species, ref, start, stop, strand,
+                                   length, transcript_count)
+        else:
+            strand = ''
+            start = 0
+            for ref in sorted(ref_lengths.keys()):
+                length = ref_lengths[ref]
+                transcript_count = 0
+                if args.debug:
+                    print >> sys.stderr,\
+                        "Retrieving transcripts: %s\t%s\tlength: %d"\
+                        % (species, ref, length)
+                transcript_count += retrieve_region(species, ref, start,
+                                                    length, strand)
+                if args.debug or args.verbose:
+                    print >> sys.stderr,\
+                        "%s\t%s\tlength: %d\ttrancripts:%d"\
+                        % (species, ref, length, transcript_count)
+
+
+if __name__ == "__main__":
+    __main__()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/retrieve_ensembl_bed.xml	Sun Jan 14 14:11:53 2018 -0500
@@ -0,0 +1,100 @@
+<tool id="retrieve_ensembl_bed" name="Retrieve Ensembl features in BED format" version="0.1.0">
+    <description>using Ensembl REST API</description>
+    <macros>
+        <import>macros.xml</import>
+    </macros>
+    <requirements>
+        <expand macro="ensembl_requirements" />
+        <expand macro="bedutil_requirements" />
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        python '$__tool_directory__/retrieve_ensembl_bed.py'  
+            --species '$species'
+            #if $extended_bed:
+                --extended_bed
+            #end if
+            $ucsc_chrom_names
+            #if $biotypes:
+                --biotypes '$biotypes'
+            #end if
+            #if $regions:
+                --regions '$regions'
+            #end if
+            '$transcript_bed'
+    ]]></command>
+    <inputs>
+        <param name="species" type="text" value="" label="Ensembl species" >
+            <help>
+            </help>
+            <expand macro="species_options" />
+            <validator type="regex" message="Enter an Ensembl organism">^\w+.*$</validator>
+        </param>
+        <param name="extended_bed" type="boolean" truevalue=",second_name,cds_start_status,cds_end_status,exon_frames,type,gene_name,second_gene_name,gene_type" falsevalue="" checked="true" 
+               label="Keep extra columns from ensembl BED"/>
+        <param name="ucsc_chrom_names" type="boolean" truevalue="--ucsc_chrom_names" falsevalue="" checked="false" 
+               label="Use the UCSC names for Chromosomes"/>
+        <param name="biotypes" type="text" value="" optional="true" label="Restrict Feature retrieval to these biotypes" >
+            <expand macro="biotypes_help" />
+        </param>
+        <param name="regions" type="text" value="" optional="true" label="Restrict Feature retrieval to comma-separated list of regions" >
+            <help>Each region is specifed as: chr or chr:pos or chr:from-to</help>
+            <validator type="regex" message="">^(\w+(:\d+(-\d+)?)?(,\w+(:\d+(-\d+)?)?)*)?$</validator>
+        </param>
+    </inputs>
+    <outputs>
+        <data name="transcript_bed" format="bed" label="Ensembl ${species} transcripts.bed">
+            <actions>
+                <action name="column_names" type="metadata" 
+                 default="chrom,chromStart,chromEnd,name,score,strand,thickStart,thickEnd,itemRgb,blockCount,blockSizes,blockStarts${extended_bed}"/>
+            </actions>
+        </data>
+    </outputs>
+    <tests>
+        <test>
+            <param name="species" value="human"/>
+            <param name="biotypes" value="protein_coding"/>
+            <param name="regions" value="1:51194990-51275150"/>
+            <output name="transcript_bed">
+                <assert_contents>
+                    <has_text_matching expression="(chr)?1\t\d+\t\d+\tENST" />
+                </assert_contents>
+            </output>
+        </test>
+    </tests>
+    <help><![CDATA[
+Retrieve Ensembl cDNAs in BED format
+
+usage: retrieve_ensembl_bed.py [-h] [-s SPECIES] [-R REGIONS] [-B BIOTYPES]
+                               [-X] [-U] [-t] [-v] [-d]
+                               output
+
+positional arguments:
+  output                Output BED filepath, or for stdout: "-"
+
+optional arguments:
+  -h, --help            show this help message and exit
+  -s SPECIES, --species SPECIES
+                        Ensembl Species to retrieve
+  -R REGIONS, --regions REGIONS
+                        Restrict Ensembl retrieval to regions e.g.:
+                        X,2:20000-25000,3:100-500+
+  -B BIOTYPES, --biotypes BIOTYPES
+                        Restrict Ensembl biotypes to retrieve
+  -X, --extended_bed    Include the extended columns returned from Ensembl
+  -U, --ucsc_chrom_names
+                        Use the UCSC names for Chromosomes
+  -t, --toplevel        Print Ensembl toplevel for species
+  -v, --verbose         Verbose
+  -d, --debug           Debug
+
+
+Ensembl REST API returns an extended BED format with these additional columns::
+
+  second_name, cds_start_status, cds_end_status, exon_frames, type, gene_name, second_gene_name, gene_type
+
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1093/bioinformatics/btu613</citation>
+        <citation type="doi">10.1093/nar/gku1010</citation>
+    </citations>
+</tool>