Mercurial > repos > vipints > fml_gff3togtf
changeset 10:c42c69aa81f8
fixed manually the upload of version 2.1.0 - deleted accidentally added files to the repo
| author | vipints <vipin@cbio.mskcc.org> |
|---|---|
| date | Thu, 23 Apr 2015 18:01:45 -0400 |
| parents | 7d67331368f3 |
| children | 5c6f33e20fcc |
| files | GFFParser.py README.md bed_to_gff.py bed_to_gff.xml gbk_to_gff.py gbk_to_gff.xml gff_to_bed.py gff_to_bed.xml gff_to_gtf.py gff_to_gtf.xml gtf_to_gff.py gtf_to_gff.xml helper.py tool_conf.xml.sample tool_dependencies.xml |
| diffstat | 15 files changed, 1944 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/GFFParser.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,496 @@ +#!/usr/bin/env python +""" +Extract genome annotation from a GFF (a tab delimited format for storing sequence features and annotations) file. + +Requirements: + Numpy :- http://numpy.org/ + +Copyright (C) + +2009-2012 Friedrich Miescher Laboratory of the Max Planck Society, Tubingen, Germany. +2012-2015 Memorial Sloan Kettering Cancer Center, New York City, USA. +""" + +import re +import os +import sys +import urllib +import numpy as np +import helper as utils +from collections import defaultdict + +def attribute_tags(col9): + """ + Split the key-value tags from the attribute column, it takes column number 9 from GTF/GFF file + + @args col9: attribute column from GFF file + @type col9: str + """ + info = defaultdict(list) + is_gff = False + + if not col9: + return is_gff, info + + # trim the line ending semi-colon ucsc may have some white-space + col9 = col9.rstrip(';| ') + # attributes from 9th column + atbs = col9.split(" ; ") + if len(atbs) == 1: + atbs = col9.split("; ") + if len(atbs) == 1: + atbs = col9.split(";") + # check the GFF3 pattern which has key value pairs like: + gff3_pat = re.compile("\w+=") + # sometime GTF have: gene_id uc002zkg.1; + gtf_pat = re.compile("\s?\w+\s") + + key_vals = [] + + if gff3_pat.match(atbs[0]): # gff3 pattern + is_gff = True + key_vals = [at.split('=') for at in atbs] + elif gtf_pat.match(atbs[0]): # gtf pattern + for at in atbs: + key_vals.append(at.strip().split(" ",1)) + else: + # to handle attribute column has only single value + key_vals.append(['ID', atbs[0]]) + # get key, val items + for item in key_vals: + key, val = item + # replace the double qoutes from feature identifier + val = re.sub('"', '', val) + # replace the web formating place holders to plain text format + info[key].extend([urllib.unquote(v) for v in val.split(',') if v]) + + return is_gff, info + +def spec_features_keywd(gff_parts): + """ + Specify the feature key word according to the GFF specifications + + @args gff_parts: attribute field key + @type gff_parts: str + """ + for t_id in ["transcript_id", "transcriptId", "proteinId"]: + try: + gff_parts["info"]["Parent"] = gff_parts["info"][t_id] + break + except KeyError: + pass + for g_id in ["gene_id", "geneid", "geneId", "name", "gene_name", "genename"]: + try: + gff_parts["info"]["GParent"] = gff_parts["info"][g_id] + break + except KeyError: + pass + ## TODO key words + for flat_name in ["Transcript", "CDS"]: + if gff_parts["info"].has_key(flat_name): + # parents + if gff_parts['type'] in [flat_name] or re.search(r'transcript', gff_parts['type'], re.IGNORECASE): + if not gff_parts['id']: + gff_parts['id'] = gff_parts['info'][flat_name][0] + #gff_parts["info"]["ID"] = [gff_parts["id"]] + # children + elif gff_parts["type"] in ["intron", "exon", "three_prime_UTR", + "coding_exon", "five_prime_UTR", "CDS", "stop_codon", + "start_codon"]: + gff_parts["info"]["Parent"] = gff_parts["info"][flat_name] + break + return gff_parts + +def Parse(ga_file): + """ + Parsing GFF/GTF file based on feature relationship, it takes the input file. + + @args ga_file: input file name + @type ga_file: str + """ + child_map = defaultdict(list) + parent_map = dict() + + ga_handle = utils.open_file(ga_file) + + for rec in ga_handle: + rec = rec.strip('\n\r') + + # skip empty line fasta identifier and commented line + if not rec or rec[0] in ['#', '>']: + continue + # skip the genome sequence + if not re.search('\t', rec): + continue + + parts = rec.split('\t') + assert len(parts) >= 8, rec + + # process the attribute column (9th column) + ftype, tags = attribute_tags(parts[-1]) + if not tags: # skip the line if no attribute column. + continue + + # extract fields + if parts[1]: + tags["source"] = parts[1] + if parts[7]: + tags["phase"] = parts[7] + + gff_info = dict() + gff_info['info'] = dict(tags) + gff_info["is_gff3"] = ftype + gff_info['chr'] = parts[0] + gff_info['score'] = parts[5] + + if parts[3] and parts[4]: + gff_info['location'] = [int(parts[3]) , + int(parts[4])] + gff_info['type'] = parts[2] + gff_info['id'] = tags.get('ID', [''])[0] + if parts[6] in ['?', '.']: + parts[6] = None + gff_info['strand'] = parts[6] + + # key word according to the GFF spec. + # is_gff3 flag is false check this condition and get the attribute fields + if not ftype: + gff_info = spec_features_keywd(gff_info) + + # link the feature relationships + if gff_info['info'].has_key('Parent'): + for p in gff_info['info']['Parent']: + if p == gff_info['id']: + gff_info['id'] = '' + break + rec_category = 'child' + elif gff_info['id']: + rec_category = 'parent' + else: + rec_category = 'record' + + # depends on the record category organize the features + if rec_category == 'child': + for p in gff_info['info']['Parent']: + # create the data structure based on source and feature id + child_map[(gff_info['chr'], gff_info['info']['source'], p)].append( + dict( type = gff_info['type'], + location = gff_info['location'], + strand = gff_info['strand'], + score = gff_info['score'], + ID = gff_info['id'], + gene_id = gff_info['info'].get('GParent', '') + )) + elif rec_category == 'parent': + parent_map[(gff_info['chr'], gff_info['info']['source'], gff_info['id'])] = dict( + type = gff_info['type'], + location = gff_info['location'], + strand = gff_info['strand'], + score = gff_info['score'], + name = tags.get('Name', [''])[0]) + elif rec_category == 'record': + #TODO how to handle plain records? + c = 1 + ga_handle.close() + + # depends on file type create parent feature + if not ftype: + parent_map, child_map = create_missing_feature_type(parent_map, child_map) + + # connecting parent child relations + # essentially the parent child features are here from any type of GTF/GFF2/GFF3 file + gene_mat = format_gene_models(parent_map, child_map) + + return gene_mat + +def format_gene_models(parent_nf_map, child_nf_map): + """ + Genarate GeneObject based on the parsed file contents + + @args parent_nf_map: parent features with source and chromosome information + @type parent_nf_map: collections defaultdict + @args child_nf_map: transctipt and exon information are encoded + @type child_nf_map: collections defaultdict + """ + + g_cnt = 0 + gene = np.zeros((len(parent_nf_map),), dtype = utils.init_gene()) + + for pkey, pdet in parent_nf_map.items(): + # considering only gene features + #if not re.search(r'gene', pdet.get('type', '')): + # continue + + # infer the gene start and stop if not there in the + if not pdet.get('location', []): + GNS, GNE = [], [] + # multiple number of transcripts + for L1 in child_nf_map[pkey]: + GNS.append(L1.get('location', [])[0]) + GNE.append(L1.get('location', [])[1]) + GNS.sort() + GNE.sort() + pdet['location'] = [GNS[0], GNE[-1]] + + orient = pdet.get('strand', '') + gene[g_cnt]['id'] = g_cnt +1 + gene[g_cnt]['chr'] = pkey[0] + gene[g_cnt]['source'] = pkey[1] + gene[g_cnt]['name'] = pkey[-1] + gene[g_cnt]['start'] = pdet.get('location', [])[0] + gene[g_cnt]['stop'] = pdet.get('location', [])[1] + gene[g_cnt]['strand'] = orient + gene[g_cnt]['score'] = pdet.get('score','') + + # default value + gene[g_cnt]['is_alt_spliced'] = gene[g_cnt]['is_alt'] = 0 + if len(child_nf_map[pkey]) > 1: + gene[g_cnt]['is_alt_spliced'] = gene[g_cnt]['is_alt'] = 1 + + # complete sub-feature for all transcripts + dim = len(child_nf_map[pkey]) + TRS = np.zeros((dim,), dtype=np.object) + TR_TYP = np.zeros((dim,), dtype=np.object) + EXON = np.zeros((dim,), dtype=np.object) + UTR5 = np.zeros((dim,), dtype=np.object) + UTR3 = np.zeros((dim,), dtype=np.object) + CDS = np.zeros((dim,), dtype=np.object) + TISc = np.zeros((dim,), dtype=np.object) + TSSc = np.zeros((dim,), dtype=np.object) + CLV = np.zeros((dim,), dtype=np.object) + CSTOP = np.zeros((dim,), dtype=np.object) + TSTAT = np.zeros((dim,), dtype=np.object) + TSCORE = np.zeros((dim,), dtype=np.object) + + # fetching corresponding transcripts + for xq, Lv1 in enumerate(child_nf_map[pkey]): + + TID = Lv1.get('ID', '') + TRS[xq]= np.array([TID]) + + TYPE = Lv1.get('type', '') + TR_TYP[xq] = np.array('') + TR_TYP[xq] = np.array(TYPE) if TYPE else TR_TYP[xq] + + orient = Lv1.get('strand', '') + tr_score = Lv1.get('score', '') + + # fetching different sub-features + child_feat = defaultdict(list) + for Lv2 in child_nf_map[(pkey[0], pkey[1], TID)]: + E_TYP = Lv2.get('type', '') + child_feat[E_TYP].append(Lv2.get('location')) + + # make general ascending order of coordinates + if orient == '-': + for etype, excod in child_feat.items(): + if len(excod) > 1: + if excod[0][0] > excod[-1][0]: + excod.reverse() + child_feat[etype] = excod + + # make exon coordinate from cds and utr regions + if not child_feat.get('exon'): + if child_feat.get('CDS'): + exon_cod = utils.make_Exon_cod( orient, + NonetoemptyList(child_feat.get('five_prime_UTR')), + NonetoemptyList(child_feat.get('CDS')), + NonetoemptyList(child_feat.get('three_prime_UTR'))) + child_feat['exon'] = exon_cod + else: + # TODO only UTR's + # searching through keys to find a pattern describing exon feature + ex_key_pattern = [k for k in child_feat if k.endswith("exon")] + if ex_key_pattern: + child_feat['exon'] = child_feat[ex_key_pattern[0]] + + # stop_codon are seperated from CDS, add the coordinates based on strand + if child_feat.get('stop_codon'): + if orient == '+': + if child_feat.get('stop_codon')[0][0] - child_feat.get('CDS')[-1][1] == 1: + child_feat['CDS'][-1] = [child_feat.get('CDS')[-1][0], child_feat.get('stop_codon')[0][1]] + else: + child_feat['CDS'].append(child_feat.get('stop_codon')[0]) + elif orient == '-': + if child_feat.get('CDS')[0][0] - child_feat.get('stop_codon')[0][1] == 1: + child_feat['CDS'][0] = [child_feat.get('stop_codon')[0][0], child_feat.get('CDS')[0][1]] + else: + child_feat['CDS'].insert(0, child_feat.get('stop_codon')[0]) + + # transcript signal sites + TIS, cdsStop, TSS, cleave = [], [], [], [] + cds_status, exon_status, utr_status = 0, 0, 0 + + if child_feat.get('exon'): + TSS = [child_feat.get('exon')[-1][1]] + TSS = [child_feat.get('exon')[0][0]] if orient == '+' else TSS + cleave = [child_feat.get('exon')[0][0]] + cleave = [child_feat.get('exon')[-1][1]] if orient == '+' else cleave + exon_status = 1 + + if child_feat.get('CDS'): + if orient == '+': + TIS = [child_feat.get('CDS')[0][0]] + cdsStop = [child_feat.get('CDS')[-1][1]-3] + else: + TIS = [child_feat.get('CDS')[-1][1]] + cdsStop = [child_feat.get('CDS')[0][0]+3] + cds_status = 1 + # cds phase calculation + child_feat['CDS'] = utils.add_CDS_phase(orient, child_feat.get('CDS')) + + # sub-feature status + if child_feat.get('three_prime_UTR') or child_feat.get('five_prime_UTR'): + utr_status =1 + + if utr_status == cds_status == exon_status == 1: + t_status = 1 + else: + t_status = 0 + + # add sub-feature # make array for export to different out + TSTAT[xq] = t_status + EXON[xq] = np.array(child_feat.get('exon'), np.float64) + UTR5[xq] = np.array(NonetoemptyList(child_feat.get('five_prime_UTR'))) + UTR3[xq] = np.array(NonetoemptyList(child_feat.get('three_prime_UTR'))) + CDS[xq] = np.array(NonetoemptyList(child_feat.get('CDS'))) + TISc[xq] = np.array(TIS) + CSTOP[xq] = np.array(cdsStop) + TSSc[xq] = np.array(TSS) + CLV[xq] = np.array(cleave) + TSCORE[xq] = tr_score + + # add sub-features to the parent gene feature + gene[g_cnt]['transcript_status'] = TSTAT + gene[g_cnt]['transcripts'] = TRS + gene[g_cnt]['exons'] = EXON + gene[g_cnt]['utr5_exons'] = UTR5 + gene[g_cnt]['cds_exons'] = CDS + gene[g_cnt]['utr3_exons'] = UTR3 + gene[g_cnt]['transcript_type'] = TR_TYP + gene[g_cnt]['tis'] = TISc + gene[g_cnt]['cdsStop'] = CSTOP + gene[g_cnt]['tss'] = TSSc + gene[g_cnt]['cleave'] = CLV + gene[g_cnt]['transcript_score'] = TSCORE + + gene[g_cnt]['gene_info'] = dict( ID = pkey[-1], + Name = pdet.get('name'), + Source = pkey[1]) + # few empty fields // TODO fill this: + gene[g_cnt]['anno_id'] = [] + gene[g_cnt]['confgenes_id'] = [] + gene[g_cnt]['alias'] = '' + gene[g_cnt]['name2'] = [] + gene[g_cnt]['chr_num'] = [] + gene[g_cnt]['paralogs'] = [] + gene[g_cnt]['transcript_valid'] = [] + gene[g_cnt]['exons_confirmed'] = [] + gene[g_cnt]['tis_conf'] = [] + gene[g_cnt]['tis_info'] = [] + gene[g_cnt]['cdsStop_conf'] = [] + gene[g_cnt]['cdsStop_info'] = [] + gene[g_cnt]['tss_info'] = [] + gene[g_cnt]['tss_conf'] = [] + gene[g_cnt]['cleave_info'] = [] + gene[g_cnt]['cleave_conf'] = [] + gene[g_cnt]['polya_info'] = [] + gene[g_cnt]['polya_conf'] = [] + gene[g_cnt]['is_valid'] = [] + gene[g_cnt]['transcript_complete'] = [] + gene[g_cnt]['is_complete'] = [] + gene[g_cnt]['is_correctly_gff3_referenced'] = '' + gene[g_cnt]['splicegraph'] = [] + g_cnt += 1 + + ## deleting empty gene records from the main array + XPFLG=0 + for XP, ens in enumerate(gene): + if ens[0]==0: + XPFLG=1 + break + + if XPFLG==1: + XQC = range(XP, len(gene)+1) + gene = np.delete(gene, XQC) + + return gene + +def NonetoemptyList(XS): + """ + Convert a None type to empty list + + @args XS: None type + @type XS: str + """ + return [] if XS is None else XS + +def create_missing_feature_type(p_feat, c_feat): + """ + GFF/GTF file defines only child features. This function tries to create + the parent feature from the information provided in the attribute column. + + example: + chr21 hg19_knownGene exon 9690071 9690100 0.000000 + . gene_id "uc002zkg.1"; transcript_id "uc002zkg.1"; + chr21 hg19_knownGene exon 9692178 9692207 0.000000 + . gene_id "uc021wgt.1"; transcript_id "uc021wgt.1"; + chr21 hg19_knownGene exon 9711935 9712038 0.000000 + . gene_id "uc011abu.2"; transcript_id "uc011abu.2"; + + This function gets the parsed feature annotations. + + @args p_feat: Parent feature map + @type p_feat: collections defaultdict + @args c_feat: Child feature map + @type c_feat: collections defaultdict + """ + + child_n_map = defaultdict(list) + for fid, det in c_feat.items(): + # get the details from grand child + GID = STRD = SCR = None + SPOS, EPOS = [], [] + TYP = dict() + for gchild in det: + GID = gchild.get('gene_id', [''])[0] + SPOS.append(gchild.get('location', [])[0]) + EPOS.append(gchild.get('location', [])[1]) + STRD = gchild.get('strand', '') + SCR = gchild.get('score', '') + if gchild.get('type', '') == "gene": ## gencode GTF file has this problem + continue + TYP[gchild.get('type', '')] = 1 + SPOS.sort() + EPOS.sort() + + # infer transcript type + transcript_type = 'transcript' + transcript_type = 'mRNA' if TYP.get('CDS', '') or TYP.get('cds', '') else transcript_type + + # gene id and transcript id are same + transcript_id = fid[-1] + if GID == transcript_id: + transcript_id = 'Transcript:' + str(GID) + + # level -1 feature type + p_feat[(fid[0], fid[1], GID)] = dict( type = 'gene', + location = [], ## infer location based on multiple transcripts + strand = STRD, + name = GID ) + # level -2 feature type + child_n_map[(fid[0], fid[1], GID)].append( + dict( type = transcript_type, + location = [SPOS[0], EPOS[-1]], + strand = STRD, + score = SCR, + ID = transcript_id, + gene_id = '' )) + # reorganizing the grand child + for gchild in det: + child_n_map[(fid[0], fid[1], transcript_id)].append( + dict( type = gchild.get('type', ''), + location = gchild.get('location'), + strand = gchild.get('strand'), + ID = gchild.get('ID'), + score = gchild.get('score'), + gene_id = '' )) + return p_feat, child_n_map +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.md Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,71 @@ +GFFtools-GX +=========== + +A collection of tools for converting genome annotation between [GTF](https://genome.ucsc.edu/FAQ/FAQformat.html#format4), [BED](https://genome.ucsc.edu/FAQ/FAQformat.html#format1), [GenBank](http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html) and [GFF](https://genome.ucsc.edu/FAQ/FAQformat.html#format3). + +##### INTRODUCTION + +Several genome annotation centers provide their data in GTF, BED, GFF and GenBank format. I have few programs, they mainly deals with converting between GTF, BED GenBank and GFF formats. They are extensively tested with files from different centers like [ENSEMBL](http://www.ensembl.org), [UCSC](https://genome.ucsc.edu/), [JGI](http://genome.jgi.doe.gov/) and [NCBI AceView](http://www.ncbi.nlm.nih.gov/IEB/Research/Acembly/HelpJan.html). These programs can be easily integrated into your galaxy instance. + +##### CONTENTS + +Included utilities are: + + BED-to-GFF: convert data from a 12 column UCSC wiggle BED format to GFF + GBK-to-GFF: convert data from genbank format to GFF + GFF-to-BED: convert data from GFF to 12 column BED format + GFF-to-GTF: convert data from GFF to GTF + GTF-to-GFF: convert data from GTF to valid GFF + +test-data: Test data set. (move to your galaxy-root-folder/test-data/) + + You may need to move the test files into your test-data directory so galaxy can find them. + If you want to run the functional tests eg as: + + exmaple: + sh run_functional_tests.sh -id fml_gtf2gff + +##### REQUIREMENTS + + python2.6 or 2.7 and biopython + + Galaxy should be able to automatically install biopython via Galaxy toolshed. + +##### COMMENTS/QUESTIONS + +I can be reached at vipin [at] cbio.mskcc.org + +##### LICENSE + +Copyright (c) 2009-2012, Friedrich Miescher Laboratory of the Max Planck Society + + 2013-2015, Memorial Sloan Kettering Cancer Center + Vipin T Sreedharan <vipin@cbio.mskcc.org> +All rights reserved. + +Licensed under the BSD 2-Clause License: <http://opensource.org/licenses/BSD-2-Clause> + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +##### COURTESY + +To the Galaxy Team.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bed_to_gff.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,71 @@ +#!/usr/bin/env python +""" +Convert genome annotation data in a 12 column BED format to GFF3. + +Usage: + python bed_to_gff.py in.bed > out.gff + +Requirement: + helper.py : https://github.com/vipints/GFFtools-GX/blob/master/helper.py + +Copyright (C) + 2009-2012 Friedrich Miescher Laboratory of the Max Planck Society, Tubingen, Germany. + 2012-2015 Memorial Sloan Kettering Cancer Center New York City, USA. +""" + +import re +import sys +import helper + +def __main__(): + """ + main function + """ + + try: + bed_fname = sys.argv[1] + except: + print __doc__ + sys.exit(-1) + + bed_fh = helper.open_file(bed_fname) + + for line in bed_fh: + line = line.strip( '\n\r' ) + + if not line or line[0] in ['#']: + continue + + parts = line.split('\t') + assert len(parts) >= 12, line + + rstarts = parts[-1].split(',') + rstarts.pop() if rstarts[-1] == '' else rstarts + + exon_lens = parts[-2].split(',') + exon_lens.pop() if exon_lens[-1] == '' else exon_lens + + if len(rstarts) != len(exon_lens): + continue # checking the consistency col 11 and col 12 + + if len(rstarts) != int(parts[-3]): + continue # checking the number of exons and block count are same + + if not parts[5] in ['+', '-']: + parts[5] = '.' # replace the unknown strand with '.' + + # bed2gff result line + sys.stdout.write('%s\tbed2gff\tgene\t%d\t%s\t%s\t%s\t.\tID=Gene:%s;Name=Gene:%s\n' % (parts[0], int(parts[1])+1, parts[2], parts[4], parts[5], parts[3], parts[3])) + sys.stdout.write('%s\tbed2gff\ttranscript\t%d\t%s\t%s\t%s\t.\tID=%s;Name=%s;Parent=Gene:%s\n' % (parts[0], int(parts[1])+1, parts[2], parts[4], parts[5], parts[3], parts[3], parts[3])) + + st = int(parts[1]) + for ex_cnt in range(int(parts[-3])): + start = st + int(rstarts[ex_cnt]) + 1 + stop = start + int(exon_lens[ex_cnt]) - 1 + sys.stdout.write('%s\tbed2gff\texon\t%d\t%d\t%s\t%s\t.\tParent=%s\n' % (parts[0], start, stop, parts[4], parts[5], parts[3])) + + bed_fh.close() + + +if __name__ == "__main__": + __main__()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bed_to_gff.xml Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,95 @@ +<tool id="fml_bed2gff" name="BED-to-GFF" version="2.1.0"> + <description>converter</description> + <command interpreter="python">bed_to_gff.py $inf_bed > $gff_format + </command> + <inputs> + <param format="bed" name="inf_bed" type="data" label="Convert this query" help="Provide genome annotation in 12 column BED format."/> + </inputs> + <outputs> + <data format="gff" name="gff_format" label="${tool.name} on ${on_string}: Converted" /> + </outputs> + <tests> + <test> + <param name="inf_bed" value="CCDS30770.bed" /> + <output name="gff_format" file="CCDS30770.gff" /> + </test> + </tests> + <help> + +**What it does** + +This tool converts data from a 12 column UCSC wiggle BED format to GFF3 (scroll down for format description). + +-------- + +**Example** + +- The following data in UCSC Wiggle BED format:: + + chr1 11873 14409 uc001aaa.3 0 + 11873 11873 0 3 354,109,1189, 0,739,1347, + +- Will be converted to GFF3:: + + ##gff-version 3 + chr1 bed2gff gene 11874 14409 0 + . ID=Gene:uc001aaa.3;Name=Gene:uc001aaa.3 + chr1 bed2gff transcript 11874 14409 0 + . ID=uc001aaa.3;Name=uc001aaa.3;Parent=Gene:uc001aaa.3 + chr1 bed2gff exon 11874 12227 0 + . Parent=uc001aaa.3 + chr1 bed2gff exon 12613 12721 0 + . Parent=uc001aaa.3 + chr1 bed2gff exon 13221 14409 0 + . Parent=uc001aaa.3 + +-------- + +**Reference** + +**BED-to-GFF** is part of oqtans package and cited as [1]_. + +.. [1] Sreedharan VT, Schultheiss SJ, Jean G et.al., Oqtans: the RNA-seq workbench in the cloud for complete and reproducible quantitative transcriptome analysis. Bioinformatics (2014). `10.1093/bioinformatics/btt731`_ + +.. _10.1093/bioinformatics/btt731: http://goo.gl/I75poH + +-------- + +**About file formats** + +**BED format** Browser Extensible Data format was designed at UCSC for displaying data tracks in the Genome Browser. It has three required fields and several additional optional ones: + +The first three BED fields (required) are:: + + 1. chrom - The name of the chromosome (e.g. chr1, chrY_random). + 2. chromStart - The starting position in the chromosome. (The first base in a chromosome is numbered 0.) + 3. chromEnd - The ending position in the chromosome, plus 1 (i.e., a half-open interval). + +The additional BED fields (optional) are:: + + 4. name - The name of the BED line. + 5. score - A score between 0 and 1000. + 6. strand - Defines the strand - either '+' or '-'. + 7. thickStart - The starting position where the feature is drawn thickly at the Genome Browser. + 8. thickEnd - The ending position where the feature is drawn thickly at the Genome Browser. + 9. reserved - This should always be set to zero. + 10. blockCount - The number of blocks (exons) in the BED line. + 11. blockSizes - A comma-separated list of the block sizes. The number of items in this list should correspond to blockCount. + 12. blockStarts - A comma-separated list of block starts. All of the blockStart positions should be calculated relative to chromStart. The number of items in this list should correspond to blockCount. + +**GFF format** General Feature Format is a format for describing genes and other features associated with DNA, RNA and Protein sequences. GFF lines have nine tab-separated fields:: + + 1. seqid - Must be a chromosome or scaffold or contig. + 2. source - The program that generated this feature. + 3. type - The name of this type of feature. Some examples of standard feature types are "gene", "CDS", "protein", "mRNA", and "exon". + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. stop - The ending position of the feature (inclusive). + 6. score - A score between 0 and 1000. If there is no score value, enter ".". + 7. strand - Valid entries include '+', '-', or '.' (for don't know/care). + 8. phase - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'. + 9. attributes - All lines with the same group are linked together into a single item. + +-------- + +**Copyright** + +BED-to-GFF Wrapper Version 0.6 (Apr 2015) + +2009-2015 Max Planck Society, University of Tübingen & Memorial Sloan Kettering Cancer Center + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gbk_to_gff.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,216 @@ +#!/usr/bin/env python +""" +Convert data from Genbank format to GFF. + +Usage: +python gbk_to_gff.py in.gbk > out.gff + +Requirements: + BioPython:- http://biopython.org/ + helper.py:- https://github.com/vipints/GFFtools-GX/blob/master/helper.py + +Copyright (C) + 2009-2012 Friedrich Miescher Laboratory of the Max Planck Society, Tubingen, Germany. + 2012-2015 Memorial Sloan Kettering Cancer Center New York City, USA. +""" + +import os +import re +import sys +import helper +import collections +from Bio import SeqIO + +def feature_table(chr_id, source, orient, genes, transcripts, cds, exons, unk): + """ + Write the feature information + """ + for gname, ginfo in genes.items(): + line = [str(chr_id), + 'gbk2gff', + ginfo[3], + str(ginfo[0]), + str(ginfo[1]), + '.', + ginfo[2], + '.', + 'ID=%s;Name=%s' % (str(gname), str(gname))] + sys.stdout.write('\t'.join(line)+"\n") + ## construct the transcript line is not defined in the original file + t_line = [str(chr_id), 'gbk2gff', source, 0, 1, '.', ginfo[2], '.'] + + if not transcripts: + t_line.append('ID=Transcript:%s;Parent=%s' % (str(gname), str(gname))) + + if exons: ## get the entire transcript region from the defined feature + t_line[3] = str(exons[gname][0][0]) + t_line[4] = str(exons[gname][0][-1]) + elif cds: + t_line[3] = str(cds[gname][0][0]) + t_line[4] = str(cds[gname][0][-1]) + + if not cds: + t_line[2] = 'transcript' + else: + t_line[2] = 'mRNA' + sys.stdout.write('\t'.join(t_line)+"\n") + + if exons: + exon_line_print(t_line, exons[gname], 'Transcript:'+str(gname), 'exon') + + if cds: + exon_line_print(t_line, cds[gname], 'Transcript:'+str(gname), 'CDS') + if not exons: + exon_line_print(t_line, cds[gname], 'Transcript:'+str(gname), 'exon') + + else: ## transcript is defined + for idx in transcripts[gname]: + t_line[2] = idx[3] + t_line[3] = str(idx[0]) + t_line[4] = str(idx[1]) + t_line.append('ID='+str(idx[2])+';Parent='+str(gname)) + sys.stdout.write('\t'.join(t_line)+"\n") + + ## feature line print call + if exons: + exon_line_print(t_line, exons[gname], str(idx[2]), 'exon') + if cds: + exon_line_print(t_line, cds[gname], str(idx[2]), 'CDS') + if not exons: + exon_line_print(t_line, cds[gname], str(idx[2]), 'exon') + + if len(genes) == 0: ## feature entry with fragment information + + line = [str(chr_id), 'gbk2gff', source, 0, 1, '.', orient, '.'] + fStart = fStop = None + + for eid, ex in cds.items(): + fStart = ex[0][0] + fStop = ex[0][-1] + + for eid, ex in exons.items(): + fStart = ex[0][0] + fStop = ex[0][-1] + + if fStart or fStart: + + line[2] = 'gene' + line[3] = str(fStart) + line[4] = str(fStop) + line.append('ID=Unknown_Gene_' + str(unk) + ';Name=Unknown_Gene_' + str(unk)) + sys.stdout.write('\t'.join(line)+"\n") + + if not cds: + line[2] = 'transcript' + else: + line[2] = 'mRNA' + + line[8] = 'ID=Unknown_Transcript_' + str(unk) + ';Parent=Unknown_Gene_' + str(unk) + sys.stdout.write('\t'.join(line)+"\n") + + if exons: + exon_line_print(line, cds[None], 'Unknown_Transcript_' + str(unk), 'exon') + + if cds: + exon_line_print(line, cds[None], 'Unknown_Transcript_' + str(unk), 'CDS') + if not exons: + exon_line_print(line, cds[None], 'Unknown_Transcript_' + str(unk), 'exon') + + unk +=1 + + return unk + + +def exon_line_print(temp_line, trx_exons, parent, ftype): + """ + Print the EXON feature line + """ + for ex in trx_exons: + temp_line[2] = ftype + temp_line[3] = str(ex[0]) + temp_line[4] = str(ex[1]) + temp_line[8] = 'Parent=%s' % parent + sys.stdout.write('\t'.join(temp_line)+"\n") + + +def gbk_parse(fname): + """ + Extract genome annotation recods from genbank format + + @args fname: gbk file name + @type fname: str + """ + fhand = helper.open_file(gbkfname) + unk = 1 + + for record in SeqIO.parse(fhand, "genbank"): + gene_tags = dict() + tx_tags = collections.defaultdict(list) + exon = collections.defaultdict(list) + cds = collections.defaultdict(list) + mol_type, chr_id = None, None + + for rec in record.features: + + if rec.type == 'source': + try: + mol_type = rec.qualifiers['mol_type'][0] + except: + mol_type = '.' + pass + try: + chr_id = rec.qualifiers['chromosome'][0] + except: + chr_id = record.name + continue + + strand='-' + strand='+' if rec.strand>0 else strand + + fid = None + try: + fid = rec.qualifiers['gene'][0] + except: + pass + + transcript_id = None + try: + transcript_id = rec.qualifiers['transcript_id'][0] + except: + pass + + if re.search(r'gene', rec.type): + gene_tags[fid] = (rec.location._start.position+1, + rec.location._end.position, + strand, + rec.type + ) + elif rec.type == 'exon': + exon[fid].append((rec.location._start.position+1, + rec.location._end.position)) + elif rec.type=='CDS': + cds[fid].append((rec.location._start.position+1, + rec.location._end.position)) + else: + # get all transcripts + if transcript_id: + tx_tags[fid].append((rec.location._start.position+1, + rec.location._end.position, + transcript_id, + rec.type)) + # record extracted, generate feature table + unk = feature_table(chr_id, mol_type, strand, gene_tags, tx_tags, cds, exon, unk) + + fhand.close() + + +if __name__=='__main__': + + try: + gbkfname = sys.argv[1] + except: + print __doc__ + sys.exit(-1) + + ## extract gbk records + gbk_parse(gbkfname)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gbk_to_gff.xml Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,100 @@ +<tool id="fml_gbk2gff" name="GBK-to-GFF" version="2.1.0"> + <description>converter</description> + <command interpreter="python">gbk_to_gff.py $inf_gbk > $gff_format + </command> + <inputs> + <param format="gb,gbk,genbank" name="inf_gbk" type="data" label="Convert this query" help="GenBank flat file format consists of an annotation section and a sequence section."/> + </inputs> + <outputs> + <data format="gff" name="gff_format" label="${tool.name} on ${on_string}: Converted"/> + </outputs> + <tests> + <test> + <param name="inf_gbk" value="s_cerevisiae_SCU49845.gbk" /> + <output name="gff_format" file="s_cerevisiae_SCU49845.gff" /> + </test> + </tests> + <help> + +**What it does** + +This tool converts data from a GenBank_ flat file format to GFF (scroll down for format description). + +.. _GenBank: http://www.ncbi.nlm.nih.gov/genbank/ + +------ + +**Example** + +- The following data in GenBank format:: + + LOCUS NM_001202705 2406 bp mRNA linear PLN 28-MAY-2011 + DEFINITION Arabidopsis thaliana thiamine biosynthesis protein ThiC (THIC) + mRNA, complete cds. + ACCESSION NM_001202705 + VERSION NM_001202705.1 GI:334184566......... + FEATURES Location/Qualifiers + source 1..2406 + /organism="Arabidopsis thaliana" + /mol_type="mRNA" + /db_xref="taxon:3702"........ + gene 1..2406 + /gene="THIC" + /locus_tag="AT2G29630" + /gene_synonym="PY; PYRIMIDINE REQUIRING; T27A16.27;........ + ORIGIN + 1 aagcctttcg ctttaggctg cattgggccg tgacaatatt cagacgattc aggaggttcg + 61 ttcctttttt aaaggaccct aatcactctg agtaccactg actcactcag tgtgcgcgat + 121 tcatttcaaa aacgagccag cctcttcttc cttcgtctac tagatcagat ccaaagcttc + 181 ctcttccagc tatggctgct tcagtacact gtaccttgat gtccgtcgta tgcaacaaca + // + + +- Will be converted to GFF3:: + + NM_001202705 gbk2gff chromosome 1 2406 . + 1 ID=NM_001202705;Alias=2;Dbxref=taxon:3702;Name=NM_001202705 + NM_001202705 gbk2gff gene 1 2406 . + 1 ID=AT2G29630;Dbxref=GeneID:817513,TAIR:AT2G29630;Name=THIC + NM_001202705 gbk2gff mRNA 192 2126 . + 1 ID=AT2G29630.t01;Parent=AT2G29630 + NM_001202705 gbk2gff CDS 192 2126 . + 1 ID=AT2G29630.p01;Parent=AT2G29630.t01 + NM_001202705 gbk2gff exon 192 2126 . + 1 Parent=AT2G29630.t01 + +------ + +**Reference** + +**GBK-to-GFF** is part of oqtans package and cited as [1]_. + +.. [1] Sreedharan VT, Schultheiss SJ, Jean G et.al., Oqtans: the RNA-seq workbench in the cloud for complete and reproducible quantitative transcriptome analysis. Bioinformatics (2014). `10.1093/bioinformatics/btt731`_ + +.. _10.1093/bioinformatics/btt731: http://goo.gl/I75poH + +------ + +**About file formats** + +**GenBank format** An example of a GenBank record may be viewed here_ + +.. _here: http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html + +**GFF** Generic Feature Format is a format for describing genes and other features associated with DNA, RNA and Protein sequences. GFF lines have nine tab-separated fields:: + + 1. seqid - Must be a chromosome or scaffold or contig. + 2. source - The program that generated this feature. + 3. type - The name of this type of feature. Some examples of standard feature types are "gene", "CDS", "protein", "mRNA", and "exon". + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. stop - The ending position of the feature (inclusive). + 6. score - A score between 0 and 1000. If there is no score value, enter ".". + 7. strand - Valid entries include '+', '-', or '.' (for don't know/care). + 8. phase - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'. + 9. attributes - All lines with the same group are linked together into a single item. + +-------- + +**Copyright** + +GBK-to-GFF Wrapper Version 0.6 (Apr 2015) + +2009-2015 Max Planck Society, University of Tübingen & Memorial Sloan Kettering Cancer Center + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gff_to_bed.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,116 @@ +#!/usr/bin/env python +""" +Convert genome annotation data in GFF/GTF to a 12 column BED format. +BED format typically represents the transcript models. + +Usage: python gff_to_bed.py in.gff > out.bed + +Requirement: + GFFParser.py: https://github.com/vipints/GFFtools-GX/blob/master/GFFParser.py + +Copyright (C) + 2009-2012 Friedrich Miescher Laboratory of the Max Planck Society, Tubingen, Germany. + 2012-2015 Memorial Sloan Kettering Cancer Center New York City, USA. +""" + +import re +import sys +import GFFParser + +def limitBEDWrite(tinfo): + """ + Write a three column BED file + + @args tinfo: list of genes + @type tinfo: numpy object + """ + + for contig_id, feature in tinfo.items(): + uns_line = dict() + for tid, tloc in feature.items(): + uns_line[(int(tloc[0])-1, int(tloc[1]))]=1 + for ele in sorted(uns_line): + pline = [contig_id, + str(ele[0]-1), + str(ele[1])] + + sys.stdout.write('\t'.join(pline)+"\n") + + +def writeBED(tinfo): + """ + writing result files in bed format + + @args tinfo: list of genes + @type tinfo: numpy object + """ + + for ent1 in tinfo: + child_flag = False + + for idx, tid in enumerate(ent1['transcripts']): + child_flag = True + exon_cnt = len(ent1['exons'][idx]) + exon_len = '' + exon_cod = '' + rel_start = None + rel_stop = None + for idz, ex_cod in enumerate(ent1['exons'][idx]):#check for exons of corresponding transcript + exon_len += '%d,' % (ex_cod[1]-ex_cod[0]+1) + if idz == 0: #calculate the relative start position + exon_cod += '0,' + rel_start = int(ex_cod[0])-1 + rel_stop = int(ex_cod[1]) + else: + exon_cod += '%d,' % (ex_cod[0]-1-rel_start) ## shifting the coordinates to zero + rel_stop = int(ex_cod[1]) + + if exon_len: + score = 0 + score = ent1['transcript_score'][idx] if ent1['transcript_score'].any() else score ## getting the transcript score + out_print = [ent1['chr'], + str(rel_start), + str(rel_stop), + tid[0], + str(score), + ent1['strand'], + str(rel_start), + str(rel_stop), + '0', + str(exon_cnt), + exon_len, + exon_cod] + sys.stdout.write('\t'.join(out_print)+"\n") + + if not child_flag: # file just contains only a single parent type i.e, gff3 defines only one feature type + score = 0 + score = ent1['transcript_score'][0] if ent1['transcript_score'].any() else score + + out_print = [ent1['chr'], + '%d' % int(ent1['start'])-1, + '%d' % int(ent1['stop']), + ent1['name'], + str(score), + ent1['strand'], + '%d' % int(ent1['start']), + '%d' % int(ent1['stop']), + '0', + '1', + '%d,' % (int(ent1['stop'])-int(ent1['start'])+1), + '0,'] + + sys.stdout.write('\t'.join(out_print)+"\n") + + +def __main__(): + try: + query_file = sys.argv[1] + except: + print __doc__ + sys.exit(-1) + + Transcriptdb = GFFParser.Parse(query_file) + writeBED(Transcriptdb) + +if __name__ == "__main__": + __main__()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gff_to_bed.xml Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,96 @@ +<tool id="fml_gff2bed" name="GFF-to-BED" version="2.1.0"> + <description>converter</description> + <command interpreter="python">gff_to_bed.py $inf_gff > $bed_format + </command> + <inputs> + <param format="gtf,gff,gff3" name="inf_gff" type="data" label="Convert this query" help="Provide genome annotation file in GFF, GTF, GFF3."/> + </inputs> + <outputs> + <data format="bed" name="bed_format" label="${tool.name} on ${on_string}: Converted" /> + </outputs> + <tests> + <test> + <param name="inf_gff" value="MB7_3R.gff3" /> + <output name="bed_format" file="MB7_3R.bed" /> + </test> + </tests> + <help> + +**What it does** + +This tool converts gene transcript annotation from GTF or GFF or GFF3 to UCSC wiggle 12 column BED format. + +-------- + +**Example** + +- The following data in GFF3:: + + ##gff-version 3 + chr1 protein_coding gene 11874 14409 0 + . ID=Gene:uc001aaa.3;Name=Gene:uc001aaa.3 + chr1 protein_coding transcript 11874 14409 0 + . ID=uc001aaa.3;Name=uc001aaa.3;Parent=Gene:uc001aaa.3 + chr1 protein_coding exon 11874 12227 0 + . Parent=uc001aaa.3 + chr1 protein_coding exon 12613 12721 0 + . Parent=uc001aaa.3 + chr1 protein_coding exon 13221 14409 0 + . Parent=uc001aaa.3 + +- Will be converted to UCSC Wiggle BED format:: + + chr1 11874 14409 uc001aaa.3 0 + 11874 14409 0 3 354,109,1189, 0,739,1347, + +-------- + +**Reference** + +**GFF-to-BED** is part of oqtans package and cited as [1]_. + +.. [1] Sreedharan VT, Schultheiss SJ, Jean G et.al., Oqtans: the RNA-seq workbench in the cloud for complete and reproducible quantitative transcriptome analysis. Bioinformatics (2014). `10.1093/bioinformatics/btt731`_ + +.. _10.1093/bioinformatics/btt731: http://goo.gl/I75poH + +-------- + +**About file formats** + +**GFF format** General Feature Format is a format for describing genes and other features associated with DNA, RNA and Protein sequences. GFF lines have nine tab-separated fields:: + + + 1. seqid - Must be a chromosome or scaffold or contig. + 2. source - The program that generated this feature. + 3. type - The name of this type of feature. Some examples of standard feature types are "gene", "CDS", "protein", "mRNA", and "exon". + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. stop - The ending position of the feature (inclusive). + 6. score - A score between 0 and 1000. If there is no score value, enter ".". + 7. strand - Valid entries include '+', '-', or '.' (for don't know/care). + 8. phase - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'. + 9. attributes - All lines with the same group are linked together into a single item. + +**BED format** Browser Extensible Data format was designed at UCSC for displaying data tracks in the Genome Browser. It has three required fields and several additional optional ones: + +The first three BED fields (required) are:: + + 1. chrom - The name of the chromosome (e.g. chr1, chrY_random). + 2. chromStart - The starting position in the chromosome. (The first base in a chromosome is numbered 0.) + 3. chromEnd - The ending position in the chromosome, plus 1 (i.e., a half-open interval). + +The additional BED fields (optional) are:: + + 4. name - The name of the BED line. + 5. score - A score between 0 and 1000. + 6. strand - Defines the strand - either '+' or '-'. + 7. thickStart - The starting position where the feature is drawn thickly at the Genome Browser. + 8. thickEnd - The ending position where the feature is drawn thickly at the Genome Browser. + 9. reserved - This should always be set to zero. + 10. blockCount - The number of blocks (exons) in the BED line. + 11. blockSizes - A comma-separated list of the block sizes. The number of items in this list should correspond to blockCount. + 12. blockStarts - A comma-separated list of block starts. All of the blockStart positions should be calculated relative to chromStart. The number of items in this list should correspond to blockCount. + +-------- + +**Copyright** + +GFF-to-BED Wrapper Version 0.6 (Apr 2015) + +2009-2015 Max Planck Society, University of Tübingen & Memorial Sloan Kettering Cancer Center + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gff_to_gtf.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,76 @@ +#!/usr/bin/env python +""" +Program to convert data from GFF to GTF + +Usage: python gff_to_gtf.py in.gff > out.gtf + +Requirement: + GFFParser.py: https://github.com/vipints/GFFtools-GX/blob/master/GFFParser.py + +Copyright (C) + 2009-2012 Friedrich Miescher Laboratory of the Max Planck Society, Tubingen, Germany. + 2012-2015 Memorial Sloan Kettering Cancer Center New York City, USA. +""" + +import re +import sys +import GFFParser + +def printGTF(tinfo): + """ + writing result file in GTF format + + @args tinfo: parsed object from gff file + @type tinfo: numpy array + """ + + for ent1 in tinfo: + for idx, tid in enumerate(ent1['transcripts']): + + exons = ent1['exons'][idx] + cds_exons = ent1['cds_exons'][idx] + + stop_codon = start_codon = () + + if ent1['strand'] == '+': + if cds_exons.any(): + start_codon = (cds_exons[0][0], cds_exons[0][0]+2) + stop_codon = (cds_exons[-1][1]-2, cds_exons[-1][1]) + elif ent1['strand'] == '-': + if cds_exons.any(): + start_codon = (cds_exons[-1][1]-2, cds_exons[-1][1]) + stop_codon = (cds_exons[0][0], cds_exons[0][0]+2) + else: + sys.stdout.write('STRAND information missing - %s, skip the transcript - %s\n' % (ent1['strand'], tid[0])) + pass + + last_cds_cod = 0 + for idz, ex_cod in enumerate(exons): + + sys.stdout.write('%s\t%s\texon\t%d\t%d\t.\t%s\t.\tgene_id "%s"; transcript_id "%s"; exon_number "%d"; gene_name "%s"; \n' % (ent1['chr'], ent1['source'], ex_cod[0], ex_cod[1], ent1['strand'], ent1['name'], tid[0], idz+1, ent1['gene_info']['Name'])) + + if cds_exons.any(): + try: + sys.stdout.write('%s\t%s\tCDS\t%d\t%d\t.\t%s\t%d\tgene_id "%s"; transcript_id "%s"; exon_number "%d"; gene_name "%s"; \n' % (ent1['chr'], ent1['source'], cds_exons[idz][0], cds_exons[idz][1], ent1['strand'], cds_exons[idz][2], ent1['name'], tid[0], idz+1, ent1['gene_info']['Name'])) + last_cds_cod = idz + except: + pass + + if idz == 0: + sys.stdout.write('%s\t%s\tstart_codon\t%d\t%d\t.\t%s\t%d\tgene_id "%s"; transcript_id "%s"; exon_number "%d"; gene_name "%s"; \n' % (ent1['chr'], ent1['source'], start_codon[0], start_codon[1], ent1['strand'], cds_exons[idz][2], ent1['name'], tid[0], idz+1, ent1['gene_info']['Name'])) + + if stop_codon: + sys.stdout.write('%s\t%s\tstop_codon\t%d\t%d\t.\t%s\t%d\tgene_id "%s"; transcript_id "%s"; exon_number "%d"; gene_name "%s"; \n' % (ent1['chr'], ent1['source'], stop_codon[0], stop_codon[1], ent1['strand'], cds_exons[last_cds_cod][2], ent1['name'], tid[0], idz+1, ent1['gene_info']['Name'])) + + +if __name__ == "__main__": + + try: + gff_fname = sys.argv[1] + except: + print __doc__ + sys.exit(-1) + + Transcriptdb = GFFParser.Parse(gff_fname) + + printGTF(Transcriptdb)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gff_to_gtf.xml Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,92 @@ +<tool id="fml_gff2gtf" name="GFF-to-GTF" version="2.1.0"> + <description>converter</description> + <command interpreter="python">gff_to_gtf.py $inf_gff3 > $gtf_format + </command> + <inputs> + <param format="gff3,gff" name="inf_gff3" type="data" label="Convert this query" help="Provide genome annotation file in GFF or GFF3."/> + </inputs> + <outputs> + <data format="gtf" name="gtf_format" label="${tool.name} on ${on_string}: Converted" /> + </outputs> + <tests> + <test> + <param name="inf_gff3" value="ens_mm9_chr18.gff3" /> + <output name="gtf_format" file="ens_mm9_chr18.gtf" /> + </test> + </tests> + <help> + +**What it does** + +This tool converts data from GFF to GTF file format (scroll down for format description). + +-------- + +**Example** + +- The following data in GFF3:: + + ##gff-version 3 + 17 protein_coding gene 7255208 7258258 . + . ID=ENSG00000213859;Name=KCTD11 + 17 protein_coding mRNA 7255208 7258258 . + . ID=ENST00000333751;Name=KCTD11-001;Parent=ENSG00000213859 + 17 protein_coding protein 7256262 7256960 . + . ID=ENSP00000328352;Name=KCTD11-001;Parent=ENST00000333751 + 17 protein_coding five_prime_UTR 7255208 7256261 . + . Parent=ENST00000333751 + 17 protein_coding CDS 7256262 7256960 . + 0 Name=CDS:KCTD11;Parent=ENST00000333751,ENSP00000328352 + 17 protein_coding three_prime_UTR 7256961 7258258 . + . Parent=ENST00000333751 + 17 protein_coding exon 7255208 7258258 . + . Parent=ENST00000333751 + +- Will be converted to GTF:: + + 17 protein_coding exon 7255208 7258258 . + . gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; + 17 protein_coding CDS 7256262 7256957 . + 0 gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; protein_id "ENSP00000328352"; + 17 protein_coding start_codon 7256262 7256264 . + 0 gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; + 17 protein_coding stop_codon 7256958 7256960 . + 0 gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; + +-------- + +**Reference** + +**GFF-to-GTF** is part of oqtans package and cited as [1]_. + +.. [1] Sreedharan VT, Schultheiss SJ, Jean G et.al., Oqtans: the RNA-seq workbench in the cloud for complete and reproducible quantitative transcriptome analysis. Bioinformatics (2014). `10.1093/bioinformatics/btt731`_ + +.. _10.1093/bioinformatics/btt731: http://goo.gl/I75poH + +-------- + +**About formats** + +**GFF format** General Feature Format is a format for describing genes and other features associated with DNA, RNA and Protein sequences. GFF lines have nine tab-separated fields:: + + 1. seqid - Must be a chromosome or scaffold. + 2. source - The program that generated this feature. + 3. type - The name of this type of feature. Some examples of standard feature types are "gene", "CDS", "protein", "mRNA", and "exon". + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. stop - The ending position of the feature (inclusive). + 6. score - A score between 0 and 1000. If there is no score value, enter ".". + 7. strand - Valid entries include '+', '-', or '.' (for don't know/care). + 8. phase - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'. + 9. attributes - All lines with the same group are linked together into a single item. + +**GTF format** Gene Transfer Format, it borrows from GFF, but has additional structure that warrants a separate definition and format name. GTF lines have nine tab-seaparated fields:: + + 1. seqname - The name of the sequence. + 2. source - This indicating where the annotation came from. + 3. feature - The name of the feature types. The following feature types are required: 'CDS', 'start_codon' and 'stop_codon' + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. end - The ending position of the feature (inclusive). + 6. score - The score field indicates a degree of confidence in the feature's existence and coordinates. + 7. strand - Valid entries include '+', '-', or '.' + 8. frame - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. + 9. attributes - These attributes are designed for handling multiple transcripts from the same genomic region. + +-------- + +**Copyright** + +GFF-to-GTF Wrapper Version 0.6 (Apr 2015) + +2009-2015 Max Planck Society, University of Tübingen & Memorial Sloan Kettering Cancer Center + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gtf_to_gff.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,78 @@ +#!/usr/bin/env python +""" +Convert Gene Transfer Format [GTF] to Generic Feature Format Version 3 [GFF3]. + +Usage: python gtf_to_gff.py in.gtf > out.gff3 + +Requirement: + GFFParser.py: https://github.com/vipints/GFFtools-GX/blob/master/GFFParser.py + helper.py: https://github.com/vipints/GFFtools-GX/blob/master/helper.py + +Copyright (C) + 2009-2012 Friedrich Miescher Laboratory of the Max Planck Society, Tubingen, Germany. + 2012-2015 Memorial Sloan Kettering Cancer Center New York City, USA. +""" + +import re +import sys +import helper +import GFFParser + +def GFFWriter(gtf_content): + """ + write the feature information to GFF format + + @args gtf_content: Parsed object from gtf file + @type gtf_content: numpy array + """ + + sys.stdout.write('##gff-version 3\n') + for ent1 in gtf_content: + chr_name = ent1['chr'] + strand = ent1['strand'] + start = ent1['start'] + stop = ent1['stop'] + source = ent1['source'] + ID = ent1['name'] + Name = ent1['gene_info']['Name'] + Name = ID if not Name else Name + + sys.stdout.write('%s\t%s\tgene\t%d\t%d\t.\t%s\t.\tID=%s;Name=%s\n' % (chr_name, source, start, stop, strand, ID, Name)) + for idx, tid in enumerate(ent1['transcripts']): + + t_start = ent1['exons'][idx][0][0] + t_stop = ent1['exons'][idx][-1][-1] + t_type = ent1['transcript_type'][idx] + + utr5_exons, utr3_exons = [], [] + if ent1['exons'][idx].any() and ent1['cds_exons'][idx].any(): + utr5_exons, utr3_exons = helper.buildUTR(ent1['cds_exons'][idx], ent1['exons'][idx], strand) + + sys.stdout.write('%s\t%s\t%s\t%d\t%d\t.\t%s\t.\tID=%s;Parent=%s\n' % (chr_name, source, t_type, t_start, t_stop, strand, tid[0], ID)) + for ex_cod in utr5_exons: + sys.stdout.write('%s\t%s\tfive_prime_UTR\t%d\t%d\t.\t%s\t.\tParent=%s\n' % (chr_name, source, ex_cod[0], ex_cod[1], strand, tid[0])) + + for ex_cod in ent1['cds_exons'][idx]: + sys.stdout.write('%s\t%s\tCDS\t%d\t%d\t.\t%s\t%d\tParent=%s\n' % (chr_name, source, ex_cod[0], ex_cod[1], strand, ex_cod[2], tid[0])) + + for ex_cod in utr3_exons: + sys.stdout.write('%s\t%s\tthree_prime_UTR\t%d\t%d\t.\t%s\t.\tParent=%s\n' % (chr_name, source, ex_cod[0], ex_cod[1], strand, tid[0])) + + for ex_cod in ent1['exons'][idx]: + sys.stdout.write('%s\t%s\texon\t%d\t%d\t.\t%s\t.\tParent=%s\n' % (chr_name, source, ex_cod[0], ex_cod[1], strand, tid[0])) + + +def __main__(): + + try: + gtf_fname = sys.argv[1] + except: + print __doc__ + sys.exit(-1) + + gtf_file_content = GFFParser.Parse(gtf_fname) + + GFFWriter(gtf_file_content) + +if __name__ == "__main__": + __main__()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gtf_to_gff.xml Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,92 @@ +<tool id="fml_gtf2gff" name="GTF-to-GFF" version="2.1.0"> + <description>converter</description> + <command interpreter="python">gtf_to_gff.py $inf_gtf > $gff3_format + </command> + <inputs> + <param format="gtf" name="inf_gtf" type="data" label="Convert this query" help="Provide genome annotation file in GTF."/> + </inputs> + <outputs> + <data format="gff" name="gff3_format" label="${tool.name} on ${on_string}: Converted" /> + </outputs> + <tests> + <test> + <param name="inf_gtf" value="aceview_hs_37.gtf" /> + <output name="gff3_format" file="aceview_hs_37.gff3" /> + </test> + </tests> + <help> + +**What it does** + +This tool converts data from GTF to a valid GFF file (scroll down for format description). + +-------- + +**Example** + +- The following data in GTF:: + + 17 protein_coding exon 7255208 7258258 . + . gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; + 17 protein_coding CDS 7256262 7256957 . + 0 gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; protein_id "ENSP00000328352"; + 17 protein_coding start_codon 7256262 7256264 . + 0 gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; + 17 protein_coding stop_codon 7256958 7256960 . + 0 gene_id "ENSG00000213859"; transcript_id "ENST00000333751"; exon_number "1"; gene_name "KCTD11"; transcript_name "KCTD11-001"; + +- Will be converted to GFF3:: + + ##gff-version 3 + 17 protein_coding gene 7255208 7258258 . + . ID=ENSG00000213859;Name=KCTD11 + 17 protein_coding mRNA 7255208 7258258 . + . ID=ENST00000333751;Name=KCTD11-001;Parent=ENSG00000213859 + 17 protein_coding protein 7256262 7256960 . + . ID=ENSP00000328352;Name=KCTD11-001;Parent=ENST00000333751 + 17 protein_coding five_prime_UTR 7255208 7256261 . + . Parent=ENST00000333751 + 17 protein_coding CDS 7256262 7256960 . + 0 Name=CDS:KCTD11;Parent=ENST00000333751,ENSP00000328352 + 17 protein_coding three_prime_UTR 7256961 7258258 . + . Parent=ENST00000333751 + 17 protein_coding exon 7255208 7258258 . + . Parent=ENST00000333751 + +-------- + +**Reference** + +**GTF-to-GFF** is part of oqtans package and cited as [1]_. + +.. [1] Sreedharan VT, Schultheiss SJ, Jean G et.al., Oqtans: the RNA-seq workbench in the cloud for complete and reproducible quantitative transcriptome analysis. Bioinformatics (2014). `10.1093/bioinformatics/btt731`_ + +.. _10.1093/bioinformatics/btt731: http://goo.gl/I75poH + +------ + +**About formats** + +**GTF format** Gene Transfer Format, it borrows from GFF, but has additional structure that warrants a separate definition and format name. GTF lines have nine tab-seaparated fields:: + + 1. seqname - The name of the sequence. + 2. source - This indicating where the annotation came from. + 3. feature - The name of the feature types. The following feature types are required: 'CDS', 'start_codon' and 'stop_codon' + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. end - The ending position of the feature (inclusive). + 6. score - The score field indicates a degree of confidence in the feature's existence and coordinates. + 7. strand - Valid entries include '+', '-', or '.' + 8. frame - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. + 9. attributes - These attributes are designed for handling multiple transcripts from the same genomic region. + +**GFF format** General Feature Format is a format for describing genes and other features associated with DNA, RNA and Protein sequences. GFF lines have nine tab-separated fields:: + + 1. seqid - Must be a chromosome or scaffold. + 2. source - The program that generated this feature. + 3. type - The name of this type of feature. Some examples of standard feature types are "gene", "CDS", "protein", "mRNA", and "exon". + 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + 5. stop - The ending position of the feature (inclusive). + 6. score - A score between 0 and 1000. If there is no score value, enter ".". + 7. strand - Valid entries include '+', '-', or '.' (for don't know/care). + 8. phase - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'. + 9. attributes - All lines with the same group are linked together into a single item. + +-------- + +**Copyright** + +GTF-to-GFF Wrapper Version 0.6 (Apr 2015) + +2009-2015 Max Planck Society, University of Tübingen & Memorial Sloan Kettering Cancer Center + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/helper.py Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,333 @@ +#!/usr/bin/env python +""" +Common utility functions +""" + +import os +import re +import sys +import gzip +import bz2 +import numpy + +def init_gene(): + """ + Initializing the gene structure + """ + + gene_det = [('id', 'f8'), + ('anno_id', numpy.dtype), + ('confgenes_id', numpy.dtype), + ('name', 'S25'), + ('source', 'S25'), + ('gene_info', numpy.dtype), + ('alias', 'S15'), + ('name2', numpy.dtype), + ('strand', 'S2'), + ('score', 'S15'), + ('chr', 'S15'), + ('chr_num', numpy.dtype), + ('paralogs', numpy.dtype), + ('start', 'f8'), + ('stop', 'f8'), + ('transcripts', numpy.dtype), + ('transcript_type', numpy.dtype), + ('transcript_info', numpy.dtype), + ('transcript_score', numpy.dtype), + ('transcript_status', numpy.dtype), + ('transcript_valid', numpy.dtype), + ('exons', numpy.dtype), + ('exons_confirmed', numpy.dtype), + ('cds_exons', numpy.dtype), + ('utr5_exons', numpy.dtype), + ('utr3_exons', numpy.dtype), + ('tis', numpy.dtype), + ('tis_conf', numpy.dtype), + ('tis_info', numpy.dtype), + ('cdsStop', numpy.dtype), + ('cdsStop_conf', numpy.dtype), + ('cdsStop_info', numpy.dtype), + ('tss', numpy.dtype), + ('tss_info', numpy.dtype), + ('tss_conf', numpy.dtype), + ('cleave', numpy.dtype), + ('cleave_info', numpy.dtype), + ('cleave_conf', numpy.dtype), + ('polya', numpy.dtype), + ('polya_info', numpy.dtype), + ('polya_conf', numpy.dtype), + ('is_alt', 'f8'), + ('is_alt_spliced', 'f8'), + ('is_valid', numpy.dtype), + ('transcript_complete', numpy.dtype), + ('is_complete', numpy.dtype), + ('is_correctly_gff3_referenced', 'S5'), + ('splicegraph', numpy.dtype) ] + + return gene_det + +def open_file(fname): + """ + Open the file (supports .gz .bz2) and returns the handler + + @args fname: input file name for reading + @type fname: str + """ + + try: + if os.path.splitext(fname)[1] == ".gz": + FH = gzip.open(fname, 'rb') + elif os.path.splitext(fname)[1] == ".bz2": + FH = bz2.BZ2File(fname, 'rb') + else: + FH = open(fname, 'rU') + except Exception as error: + sys.exit(error) + + return FH + +def add_CDS_phase(strand, cds): + """ + Calculate CDS phase and add to the CDS exons + + @args strand: feature strand information + @type strand: +/- + @args cds: coding exon coordinates + @type cds: numpy array [[int, int, int]] + """ + + cds_region, cds_flag = [], 0 + if strand == '+': + for cdspos in cds: + if cds_flag == 0: + cdspos = (cdspos[0], cdspos[1], 0) + diff = (cdspos[1]-(cdspos[0]-1))%3 + else: + xy = 0 + if diff == 0: + cdspos = (cdspos[0], cdspos[1], 0) + elif diff == 1: + cdspos = (cdspos[0], cdspos[1], 2) + xy = 2 + elif diff == 2: + cdspos = (cdspos[0], cdspos[1], 1) + xy = 1 + diff = ((cdspos[1]-(cdspos[0]-1))-xy)%3 + cds_region.append(cdspos) + cds_flag = 1 + elif strand == '-': + cds.reverse() + for cdspos in cds: + if cds_flag == 0: + cdspos = (cdspos[0], cdspos[1], 0) + diff = (cdspos[1]-(cdspos[0]-1))%3 + else: + xy = 0 + if diff == 0: + cdspos = (cdspos[0], cdspos[1], 0) + elif diff == 1: + cdspos = (cdspos[0], cdspos[1], 2) + xy = 2 + elif diff == 2: + cdspos = (cdspos[0], cdspos[1], 1) + xy = 1 + diff = ((cdspos[1]-(cdspos[0]-1))-xy)%3 + cds_region.append(cdspos) + cds_flag = 1 + cds_region.reverse() + return cds_region + +def buildUTR(cc, ec, strand): + """ + Build UTR regions from a given set of CDS and exon coordiantes of a gene + + @args cc: coding exon coordinates + @type cc: numpy array [[int, int, int]] + @args ec: exon coordinates + @type ec: numpy array [[int, int]] + @args strand: feature strand information + @type strand: +/- + """ + + utr5 = [] + utr3 = [] + if strand == '+': + cds_s = cc[0][0] + for ex in ec: + if ex[0] <= cds_s and cds_s <= ex[1]: + if ex[0] != cds_s:utr5.append((ex[0], cds_s-1)) + break + else: + utr5.append(ex) + cds_e = cc[-1][1] + for i in range(len(ec)): + i += 1 + if ec[-i][0] <= cds_e and cds_e <= ec[-i][1]: + if ec[-i][1] != cds_e:utr3.append((cds_e +1, ec[-i][1])) + break + else: + utr3.append(ec[-i]) + utr3.reverse() + elif strand == '-': + cds_s = cc[-1][1] + for i in range(len(ec)): + i += 1 + if ec[-i][0] <= cds_s and cds_s <= ec[-i][1]: + if ec[-i][1] != cds_s:utr5.append((cds_s+1, ec[-i][1])) + break + else: + utr5.append(ec[-i]) + utr5.reverse() + cds_e = cc[0][0] + for ex in ec: + if ex[0] <= cds_e and cds_e <= ex[1]: + if ex[0] != cds_e:utr3.append((ex[0], cds_e-1)) + break + else: + utr3.append(ex) + return utr5, utr3 + +def make_Exon_cod(strand_p, five_p_utr, cds_cod, three_p_utr): + """ + Create exon cordinates from UTR's and CDS region + + @args strand_p: feature strand information + @type strand_p: +/- + @args five_p_utr: five prime utr exon coordinates + @type five_p_utr: numpy array [[int, int]] + @args cds_cod: coding exon coordinates + @type cds_cod: numpy array [[int, int, int]] + @args three_p_utr: three prime utr exon coordinates + @type three_p_utr: numpy array [[int, int]] + """ + + exon_pos = [] + if strand_p == '+': + utr5_start, utr5_end = 0, 0 + if five_p_utr != []: + utr5_start, utr5_end = five_p_utr[-1][0], five_p_utr[-1][1] + cds_5start, cds_5end = cds_cod[0][0], cds_cod[0][1] + jun_exon = [] + if cds_5start-utr5_end == 0 or cds_5start-utr5_end == 1: + jun_exon = [utr5_start, cds_5end] + if len(cds_cod) == 1: + five_prime_flag = 0 + if jun_exon != []: + five_p_utr = five_p_utr[:-1] + five_prime_flag = 1 + for utr5 in five_p_utr: + exon_pos.append(utr5) + jun_exon = [] + utr3_start, utr3_end = 0, 0 + if three_p_utr != []: + utr3_start = three_p_utr[0][0] + utr3_end = three_p_utr[0][1] + if utr3_start-cds_5end == 0 or utr3_start-cds_5end == 1: + jun_exon = [cds_5start, utr3_end] + three_prime_flag = 0 + if jun_exon != []: + cds_cod = cds_cod[:-1] + three_p_utr = three_p_utr[1:] + three_prime_flag = 1 + if five_prime_flag == 1 and three_prime_flag == 1: + exon_pos.append([utr5_start, utr3_end]) + if five_prime_flag == 1 and three_prime_flag == 0: + exon_pos.append([utr5_start, cds_5end]) + cds_cod = cds_cod[:-1] + if five_prime_flag == 0 and three_prime_flag == 1: + exon_pos.append([cds_5start, utr3_end]) + for cds in cds_cod: + exon_pos.append(cds) + for utr3 in three_p_utr: + exon_pos.append(utr3) + else: + if jun_exon != []: + five_p_utr = five_p_utr[:-1] + cds_cod = cds_cod[1:] + for utr5 in five_p_utr: + exon_pos.append(utr5) + exon_pos.append(jun_exon) if jun_exon != [] else '' + jun_exon = [] + utr3_start, utr3_end = 0, 0 + if three_p_utr != []: + utr3_start = three_p_utr[0][0] + utr3_end = three_p_utr[0][1] + cds_3start = cds_cod[-1][0] + cds_3end = cds_cod[-1][1] + if utr3_start-cds_3end == 0 or utr3_start-cds_3end == 1: + jun_exon = [cds_3start, utr3_end] + if jun_exon != []: + cds_cod = cds_cod[:-1] + three_p_utr = three_p_utr[1:] + for cds in cds_cod: + exon_pos.append(cds) + exon_pos.append(jun_exon) if jun_exon != [] else '' + for utr3 in three_p_utr: + exon_pos.append(utr3) + elif strand_p == '-': + utr3_start, utr3_end = 0, 0 + if three_p_utr != []: + utr3_start = three_p_utr[-1][0] + utr3_end = three_p_utr[-1][1] + cds_3start = cds_cod[0][0] + cds_3end = cds_cod[0][1] + jun_exon = [] + if cds_3start-utr3_end == 0 or cds_3start-utr3_end == 1: + jun_exon = [utr3_start, cds_3end] + if len(cds_cod) == 1: + three_prime_flag = 0 + if jun_exon != []: + three_p_utr = three_p_utr[:-1] + three_prime_flag = 1 + for utr3 in three_p_utr: + exon_pos.append(utr3) + jun_exon = [] + (utr5_start, utr5_end) = (0, 0) + if five_p_utr != []: + utr5_start = five_p_utr[0][0] + utr5_end = five_p_utr[0][1] + if utr5_start-cds_3end == 0 or utr5_start-cds_3end == 1: + jun_exon = [cds_3start, utr5_end] + five_prime_flag = 0 + if jun_exon != []: + cds_cod = cds_cod[:-1] + five_p_utr = five_p_utr[1:] + five_prime_flag = 1 + if three_prime_flag == 1 and five_prime_flag == 1: + exon_pos.append([utr3_start, utr5_end]) + if three_prime_flag == 1 and five_prime_flag == 0: + exon_pos.append([utr3_start, cds_3end]) + cds_cod = cds_cod[:-1] + if three_prime_flag == 0 and five_prime_flag == 1: + exon_pos.append([cds_3start, utr5_end]) + for cds in cds_cod: + exon_pos.append(cds) + for utr5 in five_p_utr: + exon_pos.append(utr5) + else: + if jun_exon != []: + three_p_utr = three_p_utr[:-1] + cds_cod = cds_cod[1:] + for utr3 in three_p_utr: + exon_pos.append(utr3) + if jun_exon != []: + exon_pos.append(jun_exon) + jun_exon = [] + (utr5_start, utr5_end) = (0, 0) + if five_p_utr != []: + utr5_start = five_p_utr[0][0] + utr5_end = five_p_utr[0][1] + cds_5start = cds_cod[-1][0] + cds_5end = cds_cod[-1][1] + if utr5_start-cds_5end == 0 or utr5_start-cds_5end == 1: + jun_exon = [cds_5start, utr5_end] + if jun_exon != []: + cds_cod = cds_cod[:-1] + five_p_utr = five_p_utr[1:] + for cds in cds_cod: + exon_pos.append(cds) + if jun_exon != []: + exon_pos.append(jun_exon) + for utr5 in five_p_utr: + exon_pos.append(utr5) + return exon_pos
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_conf.xml.sample Thu Apr 23 18:01:45 2015 -0400 @@ -0,0 +1,7 @@ +<section name="GFFtools" id="gfftools.web"> + <tool file="GFFtools-GX/gff_to_bed.xml"/> + <tool file="GFFtools-GX/bed_to_gff.xml"/> + <tool file="GFFtools-GX/gbk_to_gff.xml"/> + <tool file="GFFtools-GX/gff_to_gtf.xml"/> + <tool file="GFFtools-GX/gtf_to_gff.xml"/> +</section>