| Next changeset 1:9c4a48f5d4e7 (2019-10-07) |
|
Commit message:
planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/proteogenomics/retrieve_ensembl_bed commit 88cf1e923a8c9e5bc6953ad412d15a7c70f054d1 |
|
added:
bedutil.py ensembl_rest.py macros.xml retrieve_ensembl_bed.py retrieve_ensembl_bed.xml |
| b |
| diff -r 000000000000 -r da1b538b87e5 bedutil.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bedutil.py Mon Jan 22 13:13:47 2018 -0500 |
| [ |
| b'@@ -0,0 +1,515 @@\n+#!/usr/bin/env python\n+"""\n+#\n+#------------------------------------------------------------------------------\n+# University of Minnesota\n+# Copyright 2016, Regents of the University of Minnesota\n+#------------------------------------------------------------------------------\n+# Author:\n+#\n+# James E Johnson\n+#\n+#------------------------------------------------------------------------------\n+"""\n+\n+from __future__ import print_function\n+\n+import sys\n+\n+from Bio.Seq import reverse_complement, translate\n+\n+\n+def bed_from_line(line, ensembl=False, seq_column=None):\n+ fields = line.rstrip(\'\\r\\n\').split(\'\\t\')\n+ if len(fields) < 12:\n+ return None\n+ (chrom, chromStart, chromEnd, name, score, strand,\n+ thickStart, thickEnd, itemRgb,\n+ blockCount, blockSizes, blockStarts) = fields[0:12]\n+ bed_entry = BedEntry(chrom=chrom, chromStart=chromStart, chromEnd=chromEnd,\n+ name=name, score=score, strand=strand,\n+ thickStart=thickStart, thickEnd=thickEnd,\n+ itemRgb=itemRgb,\n+ blockCount=blockCount,\n+ blockSizes=blockSizes.rstrip(\',\'),\n+ blockStarts=blockStarts.rstrip(\',\'))\n+ if seq_column is not None and -len(fields) <= seq_column < len(fields):\n+ bed_entry.seq = fields[seq_column]\n+ if ensembl and len(fields) >= 20:\n+ bed_entry.second_name = fields[12]\n+ bed_entry.cds_start_status = fields[13]\n+ bed_entry.cds_end_status = fields[14]\n+ bed_entry.exon_frames = fields[15].rstrip(\',\')\n+ bed_entry.biotype = fields[16]\n+ bed_entry.gene_name = fields[17]\n+ bed_entry.second_gene_name = fields[18]\n+ bed_entry.gene_type = fields[19]\n+ return bed_entry\n+\n+\n+def as_int_list(obj):\n+ if obj is None:\n+ return None\n+ if isinstance(obj, list):\n+ return [int(x) for x in obj]\n+ elif isinstance(obj, str):\n+ return [int(x) for x in obj.split(\',\')]\n+ else: # python2 unicode?\n+ return [int(x) for x in str(obj).split(\',\')]\n+\n+\n+class BedEntry(object):\n+ def __init__(self, chrom=None, chromStart=None, chromEnd=None,\n+ name=None, score=None, strand=None,\n+ thickStart=None, thickEnd=None, itemRgb=None,\n+ blockCount=None, blockSizes=None, blockStarts=None):\n+ self.chrom = chrom\n+ self.chromStart = int(chromStart)\n+ self.chromEnd = int(chromEnd)\n+ self.name = name\n+ self.score = int(score) if score is not None else 0\n+ self.strand = \'-\' if str(strand).startswith(\'-\') else \'+\'\n+ self.thickStart = int(thickStart) if thickStart else self.chromStart\n+ self.thickEnd = int(thickEnd) if thickEnd else self.chromEnd\n+ self.itemRgb = str(itemRgb) if itemRgb is not None else r\'100,100,100\'\n+ self.blockCount = int(blockCount)\n+ self.blockSizes = as_int_list(blockSizes)\n+ self.blockStarts = as_int_list(blockStarts)\n+ self.second_name = None\n+ self.cds_start_status = None\n+ self.cds_end_status = None\n+ self.exon_frames = None\n+ self.biotype = None\n+ self.gene_name = None\n+ self.second_gene_name = None\n+ self.gene_type = None\n+ self.seq = None\n+ self.cdna = None\n+ self.pep = None\n+ # T26C\n+ self.aa_change = []\n+ # p.Trp26Cys g.<pos><ref>><alt> # g.1304573A>G\n+ self.variants = []\n+\n+ def __str__(self):\n+ return \'%s\\t%d\\t%d\\t%s\\t%d\\t%s\\t%d\\t%d\\t%s\\t%d\\t%s\\t%s\' % (\n+ self.chrom, self.chromStart, self.chromEnd,\n+ self.name, self.score, self.strand,\n+ self.thickStart, self.thickEnd, str(self.itemRgb), self.blockCount,\n+ \',\'.join([str(x) for x in self.blockSizes]),\n+ \',\'.join([str(x) for x in self.blockStarts]))\n+\n+ def get_splice_junctions(self):\n+ spl'..b' translation.find(\'*\', junc)\n+ tstop = stop if stop >= 0 else len(translation)\n+ offset = (block_sum - i) % 3\n+ trimmed = translation[tstart:tstop]\n+ if debug:\n+ print("frame: %d\\ttstart: %d tstop: %d " +\n+ "offset: %d\\t%s" %\n+ (i, tstart, tstop, offset, trimmed),\n+ file=sys.stderr)\n+ if filtering and tstart > ignore:\n+ continue\n+ # get genomic locations for start and end\n+ if self.strand == \'+\':\n+ chromStart = self.chromStart + i + (tstart * 3)\n+ chromEnd = self.chromEnd - offset\\\n+ - (len(translation) - tstop) * 3\n+ else:\n+ chromStart = self.chromStart + offset\\\n+ + (len(translation) - tstop) * 3\n+ chromEnd = self.chromEnd - i - (tstart * 3)\n+ # get the blocks for this translation\n+ (tblockCount, tblockSizes, tblockStarts) =\\\n+ self.get_blocks(chromStart, chromEnd)\n+ translations[i] = (chromStart, chromEnd, trimmed,\n+ tblockCount, tblockSizes, tblockStarts)\n+ if debug:\n+ print("tblockCount: %d tblockStarts: %s " +\n+ "tblockSizes: %s" %\n+ (tblockCount, tblockStarts, tblockSizes),\n+ file=sys.stderr)\n+ return translations\n+\n+ def get_seq_id(self, seqtype=\'unk:unk\', reference=\'\', frame=None):\n+ # Ensembl fasta ID format\n+ # >ID SEQTYPE:STATUS LOCATION GENE TRANSCRIPT\n+ # >ENSP00000328693 pep:splice chromosome:NCBI35:1:904515:910768:1\\\n+ # gene:ENSG00000158815:transcript:ENST00000328693\\\n+ # gene_biotype:protein_coding transcript_biotype:protein_coding\n+ frame_name = \'\'\n+ chromStart = self.chromStart\n+ chromEnd = self.chromEnd\n+ strand = 1 if self.strand == \'+\' else -1\n+ if frame is not None:\n+ block_sum = sum(self.blockSizes)\n+ offset = (block_sum - frame) % 3\n+ frame_name = \'_\' + str(frame + 1)\n+ if self.strand == \'+\':\n+ chromStart += frame\n+ chromEnd -= offset\n+ else:\n+ chromStart += offset\n+ chromEnd -= frame\n+ location = "chromosome:%s:%s:%s:%s:%s"\\\n+ % (reference, self.chrom, chromStart, chromEnd, strand)\n+ seq_id = "%s%s %s %s" % (self.name, frame_name, seqtype, location)\n+ return seq_id\n+\n+ def get_line(self, start_offset=0, end_offset=0):\n+ if start_offset or end_offset:\n+ s_offset = start_offset if start_offset else 0\n+ e_offset = end_offset if end_offset else 0\n+ if s_offset > self.chromStart:\n+ s_offset = self.chromStart\n+ chrStart = self.chromStart - s_offset\n+ chrEnd = self.chromEnd + e_offset\n+ blkSizes = self.blockSizes\n+ blkSizes[0] += s_offset\n+ blkSizes[-1] += e_offset\n+ blkStarts = self.blockStarts\n+ for i in range(1, self.blockCount):\n+ blkStarts[i] += s_offset\n+ items = [str(x) for x in [self.chrom, chrStart, chrEnd, self.name,\n+ self.score, self.strand, self.thickStart,\n+ self.thickEnd, self.itemRgb,\n+ self.blockCount,\n+ \',\'.join([str(x) for x in blkSizes]),\n+ \',\'.join([str(x) for x in blkStarts])]]\n+ return \'\\t\'.join(items) + \'\\n\'\n+ return self.line\n' |
| b |
| diff -r 000000000000 -r da1b538b87e5 ensembl_rest.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ensembl_rest.py Mon Jan 22 13:13:47 2018 -0500 |
| [ |
| @@ -0,0 +1,129 @@ +#!/usr/bin/env python +""" +# +#------------------------------------------------------------------------------ +# University of Minnesota +# Copyright 2017, Regents of the University of Minnesota +#------------------------------------------------------------------------------ +# Author: +# +# James E Johnson +# +#------------------------------------------------------------------------------ +""" + +from __future__ import print_function +from __future__ import unicode_literals + +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("%s" % ext, file=sys.stderr) + r = requests.get(server+ext, headers=headers) + if r.status_code == 429: + 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) + 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("%s\t%s\t%s\t%s\t%s" % + (species['name'], species['common_name'], + species['display_name'], + species['strain'], + species['taxon_id']), file=sys.stdout) + 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 = list(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 |
| b |
| diff -r 000000000000 -r da1b538b87e5 macros.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros.xml Mon Jan 22 13:13:47 2018 -0500 |
| [ |
| b'@@ -0,0 +1,125 @@\n+<macros>\n+ <xml name="bedutil_requirements">\n+ <requirement type="package" version="1.62">biopython</requirement>\n+ </xml>\n+ <xml name="ensembl_requirements">\n+ <requirement type="package" version="0.4.10">requests-cache</requirement>\n+ </xml>\n+ <xml name="twobit_requirements">\n+ <requirement type="package" version="3.1.4">twobitreader</requirement>\n+ </xml>\n+ <xml name="species_options">\n+ <option value="homo_sapiens">homo_sapiens (Human) taxon_id: 9606</option>\n+ <option value="mus_musculus">mus_musculus (Mouse) taxon_id: 10090</option>\n+ <option value="ailuropoda_melanoleuca">ailuropoda_melanoleuca (Panda) taxon_id: 9646</option>\n+ <option value="anas_platyrhynchos">anas_platyrhynchos (Duck) taxon_id: 8839</option>\n+ <option value="anolis_carolinensis">anolis_carolinensis (Anole lizard) taxon_id: 28377</option>\n+ <option value="astyanax_mexicanus">astyanax_mexicanus (Cave fish) taxon_id: 7994</option>\n+ <option value="bos_taurus">bos_taurus (Cow) taxon_id: 9913</option>\n+ <option value="caenorhabditis_elegans">caenorhabditis_elegans (Caenorhabditis elegans) taxon_id: 6239</option>\n+ <option value="callithrix_jacchus">callithrix_jacchus (Marmoset) taxon_id: 9483</option>\n+ <option value="canis_familiaris">canis_familiaris (Dog) taxon_id: 9615</option>\n+ <option value="carlito_syrichta">carlito_syrichta (Tarsier) taxon_id: 1868482</option>\n+ <option value="cavia_aperea">cavia_aperea (Brazilian guinea pig) taxon_id: 37548</option>\n+ <option value="cavia_porcellus">cavia_porcellus (Guinea Pig) taxon_id: 10141</option>\n+ <option value="chinchilla_lanigera">chinchilla_lanigera (Long-tailed chinchilla) taxon_id: 34839</option>\n+ <option value="chlorocebus_sabaeus">chlorocebus_sabaeus (Vervet-AGM) taxon_id: 60711</option>\n+ <option value="choloepus_hoffmanni">choloepus_hoffmanni (Sloth) taxon_id: 9358</option>\n+ <option value="ciona_intestinalis">ciona_intestinalis (C.intestinalis) taxon_id: 7719</option>\n+ <option value="ciona_savignyi">ciona_savignyi (C.savignyi) taxon_id: 51511</option>\n+ <option value="cricetulus_griseus_chok1gshd">cricetulus_griseus_chok1gshd (Chinese hamster CHOK1GS) taxon_id: 10029</option>\n+ <option value="cricetulus_griseus_crigri">cricetulus_griseus_crigri (Chinese hamster CriGri) taxon_id: 10029</option>\n+ <option value="danio_rerio">danio_rerio (Zebrafish) taxon_id: 7955</option>\n+ <option value="dasypus_novemcinctus">dasypus_novemcinctus (Armadillo) taxon_id: 9361</option>\n+ <option value="dipodomys_ordii">dipodomys_ordii (Kangaroo rat) taxon_id: 10020</option>\n+ <option value="drosophila_melanogaster">drosophila_melanogaster (Fruitfly) taxon_id: 7227</option>\n+ <option value="echinops_telfairi">echinops_telfairi (Lesser hedgehog tenrec) taxon_id: 9371</option>\n+ <option value="equus_caballus">equus_caballus (Horse) taxon_id: 9796</option>\n+ <option value="erinaceus_europaeus">erinaceus_europaeus (Hedgehog) taxon_id: 9365</option>\n+ <option value="felis_catus">felis_catus (Cat) taxon_id: 9685</option>\n+ <option value="ficedula_albicollis">ficedula_albicollis (Flycatcher) taxon_id: 59894</option>\n+ <option value="fukomys_damarensis">fukomys_damarensis (Damara mole rat) taxon_id: 885580</option>\n+ <option value="gadus_morhua">gadus_morhua (Cod) taxon_id: 8049</option>\n+ <option value="gallus_gallus">gallus_gallus (Chicken) taxon_id: 9031</option>\n+ <option value="gasterosteus_aculeatus">gasterosteus_aculeatus (Stickleback) taxon_id: 69293</option>\n+ <option value="gorilla_gorilla">gorilla_gorilla (Gorilla) taxon_id: 9595</option>\n+ '..b'n value="oreochromis_niloticus">oreochromis_niloticus (Tilapia) taxon_id: 8128</option>\n+ <option value="ornithorhynchus_anatinus">ornithorhynchus_anatinus (Platypus) taxon_id: 9258</option>\n+ <option value="oryctolagus_cuniculus">oryctolagus_cuniculus (Rabbit) taxon_id: 9986</option>\n+ <option value="oryzias_latipes">oryzias_latipes (Medaka) taxon_id: 8090</option>\n+ <option value="otolemur_garnettii">otolemur_garnettii (Bushbaby) taxon_id: 30611</option>\n+ <option value="ovis_aries">ovis_aries (Sheep) taxon_id: 9940</option>\n+ <option value="pan_troglodytes">pan_troglodytes (Chimpanzee) taxon_id: 9598</option>\n+ <option value="papio_anubis">papio_anubis (Olive baboon) taxon_id: 9555</option>\n+ <option value="pelodiscus_sinensis">pelodiscus_sinensis (Chinese softshell turtle) taxon_id: 13735</option>\n+ <option value="peromyscus_maniculatus_bairdii">peromyscus_maniculatus_bairdii (Northern American deer mouse) taxon_id: 230844</option>\n+ <option value="petromyzon_marinus">petromyzon_marinus (Lamprey) taxon_id: 7757</option>\n+ <option value="poecilia_formosa">poecilia_formosa (Amazon molly) taxon_id: 48698</option>\n+ <option value="pongo_abelii">pongo_abelii (Orangutan) taxon_id: 9601</option>\n+ <option value="procavia_capensis">procavia_capensis (Hyrax) taxon_id: 9813</option>\n+ <option value="pteropus_vampyrus">pteropus_vampyrus (Megabat) taxon_id: 132908</option>\n+ <option value="rattus_norvegicus">rattus_norvegicus (Rat) taxon_id: 10116</option>\n+ <option value="saccharomyces_cerevisiae">saccharomyces_cerevisiae (Saccharomyces cerevisiae) taxon_id: 4932</option>\n+ <option value="sarcophilus_harrisii">sarcophilus_harrisii (Tasmanian devil) taxon_id: 9305</option>\n+ <option value="sorex_araneus">sorex_araneus (Shrew) taxon_id: 42254</option>\n+ <option value="sus_scrofa">sus_scrofa (Pig) taxon_id: 9823</option>\n+ <option value="taeniopygia_guttata">taeniopygia_guttata (Zebra Finch) taxon_id: 59729</option>\n+ <option value="takifugu_rubripes">takifugu_rubripes (Fugu) taxon_id: 31033</option>\n+ <option value="tetraodon_nigroviridis">tetraodon_nigroviridis (Tetraodon) taxon_id: 99883</option>\n+ <option value="tupaia_belangeri">tupaia_belangeri (Tree Shrew) taxon_id: 37347</option>\n+ <option value="tursiops_truncatus">tursiops_truncatus (Dolphin) taxon_id: 9739</option>\n+ <option value="vicugna_pacos">vicugna_pacos (Alpaca) taxon_id: 30538</option>\n+ <option value="xenopus_tropicalis">xenopus_tropicalis (Xenopus) taxon_id: 8364</option>\n+ <option value="xiphophorus_maculatus">xiphophorus_maculatus (Platyfish) taxon_id: 8083</option>\n+ </xml>\n+ <xml name="biotypes_help">\n+ <help><![CDATA[\n+Example biotypes: \n+protein_coding, non_coding, pseudogene, nonsense_mediated_decay, non_stop_decay, \n+translated_processed_pseudogene, transcribed_processed_pseudogene, transcribed_unitary_pseudogene, transcribed_unprocessed_pseudogene, \n+polymorphic_pseudogene, processed_pseudogene, unprocessed_pseudogene, unitary_pseudogene, processed_transcript, \n+retained_intron, ccds_gene, sense_overlapping, sense_intronic, cdna_update, antisense, \n+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, \n+IG_pseudogene, IG_C_pseudogene, IG_D_pseudogene, IG_J_pseudogene, IG_V_pseudogene, TR_J_pseudogene, TR_V_pseudogene, TEC, \n+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\n+ ]]></help>\n+ </xml>\n+</macros>\n' |
| b |
| diff -r 000000000000 -r da1b538b87e5 retrieve_ensembl_bed.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/retrieve_ensembl_bed.py Mon Jan 22 13:13:47 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 +# +#------------------------------------------------------------------------------ +""" + +from __future__ import print_function + +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() + 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("selected_regions: %s" % selected_regions, file=sys.stderr) + + def retrieve_region(species, ref, start, stop, strand): + transcript_count = 0 + 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("%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("%s\n" % (bedline), file=sys.stderr) + 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("BED error (%s) : %s\n" % (e, bedline), + file=sys.stderr) + 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("%s\t%s\tlength: %d" % (species, ref, length), + file=sys.stderr) + 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("%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 + for ref in sorted(ref_lengths.keys()): + length = ref_lengths[ref] + transcript_count = 0 + if args.debug: + 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("%s\t%s\tlength: %d\ttrancripts:%d" % + (species, ref, length, transcript_count), + file=sys.stderr) + + +if __name__ == "__main__": + __main__() |
| b |
| diff -r 000000000000 -r da1b538b87e5 retrieve_ensembl_bed.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/retrieve_ensembl_bed.xml Mon Jan 22 13:13:47 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> |