Mercurial > repos > galaxyp > retrieve_ensembl_bed
changeset 1:c3d600729b6f draft
planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/proteogenomics/retrieve_ensembl_bed commit 88cf1e923a8c9e5bc6953ad412d15a7c70f054d1
| author | galaxyp |
|---|---|
| date | Mon, 22 Jan 2018 13:13:26 -0500 |
| parents | 887e111c0919 |
| children | e385fe93df68 |
| files | bedutil.py ensembl_rest.py retrieve_ensembl_bed.py |
| diffstat | 3 files changed, 89 insertions(+), 76 deletions(-) [+] |
line wrap: on
line diff
--- a/bedutil.py Sun Jan 14 14:11:53 2018 -0500 +++ b/bedutil.py Mon Jan 22 13:13:26 2018 -0500 @@ -12,6 +12,8 @@ #------------------------------------------------------------------------------ """ +from __future__ import print_function + import sys from Bio.Seq import reverse_complement, translate @@ -45,6 +47,17 @@ return bed_entry +def as_int_list(obj): + if obj is None: + return None + if isinstance(obj, list): + return [int(x) for x in obj] + elif isinstance(obj, str): + return [int(x) for x in obj.split(',')] + else: # python2 unicode? + return [int(x) for x in str(obj).split(',')] + + class BedEntry(object): def __init__(self, chrom=None, chromStart=None, chromEnd=None, name=None, score=None, strand=None, @@ -60,18 +73,8 @@ 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.blockSizes = as_int_list(blockSizes) + self.blockStarts = as_int_list(blockStarts) self.second_name = None self.cds_start_status = None self.cds_end_status = None @@ -165,6 +168,11 @@ return self.set_cds(min(cds_pos), max(cds_pos) + basepairs) return None + def get_cds_bed(self): + cds_pos = [self.cdna_offset_of_pos(self.thickStart), + self.cdna_offset_of_pos(self.thickEnd)] + return self.trim(min(cds_pos), max(cds_pos)) + def get_cigar(self): cigar = '' r = range(self.blockCount) @@ -214,7 +222,7 @@ def pos_of_cdna_offet(self, offset): if offset is not None and 0 <= offset < sum(self.blockSizes): - r = range(self.blockCount) + r = list(range(self.blockCount)) rev = self.strand == '-' if rev: r.reverse() @@ -233,7 +241,7 @@ def cdna_offset_of_pos(self, pos): if not self.chromStart <= pos < self.chromEnd: return -1 - r = range(self.blockCount) + r = list(range(self.blockCount)) rev = self.strand == '-' if rev: r.reverse() @@ -248,19 +256,21 @@ 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" + print("variant requires ref and alt sequences", file=sys.stderr) 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) + print("variant not in entry %s: %s %d < %d < %d" % + (self.name, self.strand, + self.chromStart, pos, self.chromEnd), + file=sys.stderr) + print("%s" % str(self), file=sys.stderr) return if len(ref) != len(alt): - print >> sys.stderr, "variant only works for snp: %s %s"\ - % (ref, alt) + print("variant only works for snp: %s %s" % (ref, alt), + file=sys.stderr) return if not self.seq: - print >> sys.stderr, "variant entry %s has no seq" % self.name + print("variant entry %s has no seq" % self.name, file=sys.stderr) return """ if self.strand == '-': @@ -274,27 +284,27 @@ 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) + print("variant offset %s: %s %d < %d < %d" % + (self.name, self.strand, self.chromStart, + pos+1, self.chromEnd), file=sys.stderr) + print("%s" % str(self), file=sys.stderr) 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" + print("variant requires ref and alt sequences", file=sys.stderr) 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) + print("variant not in entry %s: %s %d < %d < %d" % + (self.name, self.strand, + self.chromStart, pos, self.chromEnd), + file=sys.stderr) + print("%s" % str(self), file=sys.stderr) return None if not self.seq: - print >> sys.stderr, "variant entry %s has no seq" % self.name + print("variant entry %s has no seq" % self.name, file=sys.stderr) return None tbed = BedEntry(chrom=self.chrom, chromStart=self.chromStart, chromEnd=self.chromEnd, @@ -330,8 +340,8 @@ chromStart = self.chromStart chromEnd = self.chromEnd if debug: - print >> sys.stderr, "%s" % (str(self)) - r = range(self.blockCount) + print("%s" % (str(self)), file=sys.stderr) + r = list(range(self.blockCount)) if self.strand == '-': r.reverse() bStart = 0 @@ -353,10 +363,9 @@ 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) + print("%3d %s\t%d\t%d\t%d\t%d\t%d\t%d" % + (x, self.strand, bStart, bEnd, + tstart, tstop, chromStart, chromEnd), file=sys.stderr) bStart += self.blockSizes[x] return(chromStart, chromEnd) @@ -409,7 +418,7 @@ 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 + print("splice_sites: %s" % splice_sites, file=sys.stderr) junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0] if seq: for i in range(3): @@ -419,9 +428,10 @@ 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) + print("frame: %d\ttstart: %d tstop: %d " + + "offset: %d\t%s" % + (i, tstart, tstop, offset, translation), + file=sys.stderr) if not untrimmed: tstart = translation.rfind('*', 0, junc) + 1 stop = translation.find('*', junc) @@ -429,9 +439,10 @@ 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) + print("frame: %d\ttstart: %d tstop: %d " + + "offset: %d\t%s" % + (i, tstart, tstop, offset, trimmed), + file=sys.stderr) if filtering and tstart > ignore: continue # get genomic locations for start and end @@ -449,9 +460,10 @@ translations[i] = (chromStart, chromEnd, trimmed, tblockCount, tblockSizes, tblockStarts) if debug: - print >> sys.stderr,\ - "tblockCount: %d tblockStarts: %s tblockSizes: %s"\ - % (tblockCount, tblockStarts, tblockSizes) + print("tblockCount: %d tblockStarts: %s " + + "tblockSizes: %s" % + (tblockCount, tblockStarts, tblockSizes), + file=sys.stderr) return translations def get_seq_id(self, seqtype='unk:unk', reference='', frame=None):
--- a/ensembl_rest.py Sun Jan 14 14:11:53 2018 -0500 +++ b/ensembl_rest.py Mon Jan 22 13:13:26 2018 -0500 @@ -12,6 +12,8 @@ #------------------------------------------------------------------------------ """ +from __future__ import print_function +from __future__ import unicode_literals import sys @@ -28,10 +30,10 @@ def ensembl_rest(ext, headers): if debug: - print >> sys.stderr, "%s" % ext + print("%s" % ext, file=sys.stderr) r = requests.get(server+ext, headers=headers) if r.status_code == 429: - print >> sys.stderr, "response headers: %s\n" % r.headers + print("response headers: %s\n" % r.headers, file=sys.stderr) if 'Retry-After' in r.headers: sleep(r.headers['Retry-After']) r = requests.get(server+ext, headers=headers) @@ -47,12 +49,11 @@ 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'], + print("%s\t%s\t%s\t%s\t%s" % + (species['name'], species['common_name'], species['display_name'], species['strain'], - species['taxon_id']) + species['taxon_id']), file=sys.stdout) return results @@ -86,7 +87,7 @@ bed = [] param = params if params else '' req_header = {"Content-Type": "text/x-bed"} - regions = range(start, length, max_region) + regions = list(range(start, length, max_region)) if not regions or regions[-1] < length: regions.append(length) for end in regions[1:]:
--- a/retrieve_ensembl_bed.py Sun Jan 14 14:11:53 2018 -0500 +++ b/retrieve_ensembl_bed.py Mon Jan 22 13:13:26 2018 -0500 @@ -12,6 +12,8 @@ #------------------------------------------------------------------------------ """ +from __future__ import print_function + import argparse import re import sys @@ -49,7 +51,6 @@ 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() @@ -72,24 +73,24 @@ selected_regions[chrom] = [] selected_regions[chrom].append([start, end, strand]) if args.debug: - print >> sys.stderr, "selected_regions: %s" % selected_regions + print("selected_regions: %s" % selected_regions, file=sys.stderr) def retrieve_region(species, ref, start, stop, strand): transcript_count = 0 - regions = range(start, stop, max_region) + regions = list(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)) + print("%s\t%s\tstart: %d\tend: %d\tcDNA transcripts:%d" % + (species, ref, start, end, len(bedlines)), + file=sys.stderr) # start, end, seq for i, bedline in enumerate(bedlines): if args.debug: - print >> sys.stderr, "%s\n" % (bedline) + print("%s\n" % (bedline), file=sys.stderr) if not args.ucsc_chrom_names: bedline = re.sub('^[^\t]+', ref, bedline) try: @@ -100,8 +101,8 @@ out_wtr.write("\n") out_wtr.flush() except Exception as e: - print >> sys.stderr,\ - "BED error (%s) : %s\n" % (e, bedline) + print("BED error (%s) : %s\n" % (e, bedline), + file=sys.stderr) start = end + 1 return transcript_count @@ -112,8 +113,8 @@ length = coord_systems['chromosome'][ref] ref_lengths[ref] = length if args.toplevel: - print >> sys.stderr,\ - "%s\t%s\tlength: %d" % (species, ref, length) + print("%s\t%s\tlength: %d" % (species, ref, length), + file=sys.stderr) if selected_regions: transcript_count = 0 for ref in sorted(selected_regions.keys()): @@ -129,10 +130,10 @@ 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) + print("%s\t%s:%d-%d%s\tlength: %d\ttrancripts:%d" % + (species, ref, start, stop, strand, + length, transcript_count), + file=sys.stderr) else: strand = '' start = 0 @@ -140,15 +141,14 @@ length = ref_lengths[ref] transcript_count = 0 if args.debug: - print >> sys.stderr,\ - "Retrieving transcripts: %s\t%s\tlength: %d"\ - % (species, ref, length) + print("Retrieving transcripts: %s\t%s\tlength: %d" % + (species, ref, length), file=sys.stderr) 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) + print("%s\t%s\tlength: %d\ttrancripts:%d" % + (species, ref, length, transcript_count), + file=sys.stderr) if __name__ == "__main__":