given a GTF file, produce a new GTF file with only the longest transcript
    for each gene
    function lifted from:
    https://pythonhosted.org/gffutils/_modules/gffutils/helpers.html
347685
def canonical_transcripts(gtf, out_file):
    
    if file_exists(out_file):
        return out_file
    db = get_gtf_db(gtf)
    with file_transaction(out_file) as tx_out_file:
        with open(tx_out_file, "w") as out_handle:
            for gene in db.features_of_type('gene'):
                exon_list = []
                for ti, transcript in enumerate(db.children(gene, level=1)):
                    cds_len = 0
                    total_len = 0
                    exons = list(db.children(transcript, level=1))
                    for exon in exons:
                        exon_length = len(exon)
                        if exon.featuretype == 'CDS':
                            cds_len += exon_length
                        total_len += exon_length

                    exon_list.append((cds_len, total_len, transcript, exons))

                # If we have CDS, then use the longest coding transcript
                if max(i[0] for i in exon_list) > 0:
                    best = sorted(exon_list)[0]
                # Otherwise, just choose the longest
                else:
                    best = sorted(exon_list, key=lambda x: x[1])[0]
                for exon in db.children(best[2], level=1):
                    out_handle.write(str(exon) + "\n")
    return out_file
34558
def transcript_to_gene(gtf):
    
    gene_lookup = {}
    for feature in complete_features(get_gtf_db(gtf)):
        gene_id = feature.attributes.get('gene_id', [None])[0]
        transcript_id = feature.attributes.get('transcript_id', [None])[0]
        gene_lookup[transcript_id] = gene_id
    return gene_lookup
1291224
def fixpartials(args):
    
    p = OptionParser(fixpartials.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    gffile, gfasta, partials, = args

    gff = make_index(gffile)
    genome = Fasta(gfasta, index=True)
    partials = LineFile(partials, load=True).lines

    # all_transcripts = [f.id for f in gff.features_of_type("mRNA", \
    #     order_by=("seqid", "start"))]
    seen = set()
    fw = must_open(opts.outfile, "w")
    for gene in gff.features_of_type('gene',  order_by=('seqid', 'start')):
        children = AutoVivification()
        cflag = False
        transcripts = list(gff.children(gene, level=1, order_by=('start')))
        for transcript in transcripts:
            trid, seqid, strand = transcript.id, transcript.seqid, transcript.strand

            for child in gff.children(transcript, order_by=('start')):
                ftype = child.featuretype
                if ftype not in children[trid]: children[trid][ftype] = []
                children[trid][ftype].append(child)

            five_prime, three_prime = True, True
            nstart, nstop = (None, None), (None, None)
            cds_span = [children[trid]['CDS'][0].start, \
                children[trid]['CDS'][-1].stop]
            new_cds_span = [x for x in cds_span]

            start_codon = (cds_span[0], cds_span[0] + 2)
            stop_codon = (cds_span[1] - 2, cds_span[1])
            if strand == '-': start_codon, stop_codon = stop_codon, start_codon

            if trid in partials:
                seen.add(trid)
                start_codon_fasta = _fasta_slice(genome, seqid, \
                    start_codon[0], start_codon[1], strand)
                stop_codon_fasta = _fasta_slice(genome, seqid, \
                    stop_codon[0], stop_codon[1], strand)

                if not is_valid_codon(start_codon_fasta, type='start'):
                    five_prime = False
                    nstart = scan_for_valid_codon(start_codon, strand, \
                        seqid, genome, type='start')

                if not is_valid_codon(stop_codon_fasta, type='stop'):
                    three_prime = False
                    nstop = scan_for_valid_codon(stop_codon, strand, \
                        seqid, genome, type='stop')

                logging.debug("feature={0} ({1})".format(trid, strand) + \
                ", 5'={0}, 3'={1}".format(five_prime, three_prime) + \
                ", {0} <== {1} ==> {2}".format(nstart if strand == '+' else nstop, \
                    cds_span, nstop if strand == '+' else nstart))

            if not five_prime or not three_prime:
                if nstart != (None, None) and (start_codon != nstart):
                    i = 0 if strand == '+' else 1
                    new_cds_span[i] = nstart[i]
                if nstop != (None, None) and (stop_codon != nstop):
                    i = 1 if strand == '+' else 0
                    new_cds_span[i] = nstop[i]
                new_cds_span.sort()

            if set(cds_span) != set(new_cds_span):
                cflag = True
                # if CDS has been extended, appropriately adjust all relevent
                # child feature (CDS, exon, UTR) coordinates
                for ftype in children[trid]:
                    for idx in range(len(children[trid][ftype])):
                        child_span = (children[trid][ftype][idx].start, \
                            children[trid][ftype][idx].stop)
                        if ftype in ('exon', 'CDS'):
                            # if exons/CDSs, adjust start and stop according to
                            # new CDS start and stop, respectively
                            if child_span[0] == cds_span[0]:
                                children[trid][ftype][idx].start = new_cds_span[0]
                            if child_span[1] == cds_span[1]:
                                children[trid][ftype][idx].stop = new_cds_span[1]
                        elif ftype.endswith('UTR'):
                            # if *_prime_UTR, adjust stop according to new CDS start and
                            #                 adjust start according to new CDS stop
                            if child_span[1] == cds_span[0]:
                                children[trid][ftype][idx].stop = new_cds_span[0]
                            if child_span[0] == cds_span[1]:
                                children[trid][ftype][idx].start = new_cds_span[1]

                transcript.start, transcript.stop = \
                    children[trid]['exon'][0].start, children[trid]['exon'][-1].stop

        if cflag:
            _gene_span = range_minmax([(tr.start, tr.stop) for tr in transcripts])
            gene.start, gene.stop = _gene_span[0], _gene_span[1]

        # print gff file
        print(gene, file=fw)
        for transcript in transcripts:
            trid = transcript.id
            print(transcript, file=fw)
            for cftype in children[trid]:
                for child in children[trid][cftype]:
                    print(child, file=fw)

    fw.close()
703843
def sanitize_gff_file(gff_fname,
                      in_memory=True,
                      in_place=False):
    
    db = None
    if is_gff_db(gff_fname):
        # It's a database filename, so load it
        db = gffutils.FeatureDB(gff_fname)
    else:
        # Need to create a database for file
        if in_memory:
            db = gffutils.create_db(gff_fname, ":memory:",
                                    verbose=False)
        else:
            db = get_gff_db(gff_fname)
    if in_place:
        gff_out = gffwriter.GFFWriter(gff_fname,
                                      in_place=in_place)
    else:
        gff_out = gffwriter.GFFWriter(sys.stdout)
    sanitized_db = sanitize_gff_db(db)
    for gene_rec in sanitized_db.all_features(featuretype="gene"):
        gff_out.write_gene_recs(sanitized_db, gene_rec.id)
    gff_out.close()
596422
def get_gtf_db(gtf, in_memory=False):
    
    db_file = gtf + '.db'
    if gtf.endswith('.gz'):
        db_file = gtf[:-3] + '.db'
    if file_exists(db_file):
        return gffutils.FeatureDB(db_file)
    db_file = ':memory:' if in_memory else db_file
    if in_memory or not file_exists(db_file):
        debug('GTF database does not exist, creating...')
        infer_extent = guess_infer_extent(gtf)
        db = gffutils.create_db(gtf, dbfn=db_file,
                                infer_gene_extent=infer_extent)
        return db
    else:
        return gffutils.FeatureDB(db_file)
142945
private List<String> calculateTranscriptHgvs(Variant variant, Transcript transcript, String geneId) {

        String mutationType = ">";

        // Populate HGVSName parse tree.
        HgvsStringBuilder hgvsStringBuilder = new HgvsStringBuilder();

        // Populate coordinates.
        // Use cDNA coordinates.
        hgvsStringBuilder.setKind(isCoding(transcript) ? HgvsStringBuilder.Kind.CODING : HgvsStringBuilder.Kind.NON_CODING);

        hgvsStringBuilder.setCdnaStart(genomicToCdnaCoord(transcript, variant.getStart()));
        hgvsStringBuilder.setCdnaEnd(hgvsStringBuilder.getCdnaStart());

        // Populate prefix.
        hgvsStringBuilder.setTranscriptId(transcript.getId());
        hgvsStringBuilder.setGeneId(geneId);

        String reference;
        String alternate;
        // Convert alleles to transcript strand.
        if (transcript.getStrand().equals("-")) {
            reference = String.valueOf(VariantAnnotationUtils.COMPLEMENTARY_NT.get(variant.getReference().charAt(0)));
            alternate = String.valueOf(VariantAnnotationUtils.COMPLEMENTARY_NT.get(variant.getAlternate().charAt(0)));
        } else {
            reference = variant.getReference();
            alternate = variant.getAlternate();
        }

        // Populate alleles.
        hgvsStringBuilder.setMutationType(mutationType);
        hgvsStringBuilder.setReference(reference);
        hgvsStringBuilder.setAlternate(alternate);

        return Collections.singletonList(hgvsStringBuilder.format());
    }
172990
def GeneReader( fh, format='gff' ):
    

    known_formats = ( 'gff', 'gtf', 'bed')
    if format not in known_formats: 
        print('%s format not in %s' % (format, ",".join( known_formats )), file=sys.stderr)
        raise Exception('?')

    if format == 'bed':
        for line in fh:    
            f = line.strip().split()
            chrom = f[0]
            chrom_start = int(f[1])
            name = f[4]
            strand = f[5]
            cdsStart = int(f[6])
            cdsEnd = int(f[7])
            blockCount = int(f[9])
            blockSizes = [ int(i) for i in f[10].strip(',').split(',') ]
            blockStarts = [ chrom_start + int(i) for i in f[11].strip(',').split(',') ]

            # grab cdsStart - cdsEnd
            gene_exons = []
            for base,offset in zip( blockStarts, blockSizes ):
                exon_start = base
                exon_end = base+offset
                gene_exons.append( (exon_start, exon_end) )
            
    genelist = {}
    grouplist = []
    if format == 'gff' or format == 'gtf':
        for line in fh:
            if line.startswith('#'): continue
            fields = line.strip().split('\t')
            if len( fields ) < 9: continue

            # fields

            chrom = fields[0]
            ex_st = int( fields[3] ) - 1 # make zero-centered
            ex_end = int( fields[4] ) #+ 1 # make exclusive
            strand = fields[6]

            if format == 'gtf':
                group = fields[8].split(';')[0]
            else:
                group = fields[8]

            if group not in grouplist: grouplist.append( group )
            if group not in genelist:
                genelist[group] = (chrom, strand, [])
            exons_i = 2
            genelist[group][exons_i].append( ( ex_st, ex_end ) )

        sp = lambda a,b: cmp( a[0], b[0] )

        #for gene in genelist.values():
        for gene in grouplist:
            chrom, strand, gene_exons = genelist[ gene ]
            gene_exons = bitset_union( gene_exons )
            yield chrom, strand, gene_exons, gene
1242986
def _add_genes_to_bed(in_file, gene_file, fai_file, out_file, data, max_distance=10000):
    
    try:
        input_rec = next(iter(pybedtools.BedTool(in_file)))
    except StopIteration:  # empty file
        utils.copy_plus(in_file, out_file)
        return
    # keep everything after standard chrom/start/end, 1-based
    extra_fields = list(range(4, len(input_rec.fields) + 1))
    # keep the new gene annotation
    gene_index = len(input_rec.fields) + 4
    extra_fields.append(gene_index)
    columns = ",".join([str(x) for x in extra_fields])
    max_column = max(extra_fields) + 1
    ops = ",".join(["distinct"] * len(extra_fields))
    # swap over gene name to '.' if beyond maximum distance
    # cut removes the last distance column which can cause issues
    # with bedtools merge: 'ERROR: illegal character '.' found in integer conversion of string'
    distance_filter = (r"""awk -F$'\t' -v OFS='\t' '{if ($NF > %s || $NF < -%s) $%s = "."} {print}'""" %
                       (max_distance, max_distance, gene_index))
    sort_cmd = bedutils.get_sort_cmd(os.path.dirname(out_file))
    cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
    # Ensure gene transcripts match reference genome
    ready_gene_file = os.path.join(os.path.dirname(out_file), "%s-genomeonly.bed" %
                                   (utils.splitext_plus(os.path.basename(gene_file))[0]))
    ready_gene_file = bedutils.subset_to_genome(gene_file, ready_gene_file, data)
    exports = "export TMPDIR=%s && %s" % (os.path.dirname(out_file), utils.local_path_export())
    bcbio_py = sys.executable
    gsort = config_utils.get_program("gsort", data)
    cmd = ("{exports}{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^# | "
           "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
           "{gsort} - {fai_file} | "
            "bedtools closest -g {fai_file} "
            "-D ref -t first -a - -b <({gsort} {ready_gene_file} {fai_file}) | "
            "{distance_filter} | cut -f 1-{max_column} | "
            "bedtools merge -i - -c {columns} -o {ops} -delim ',' -d -10 > {out_file}")
    do.run(cmd.format(**locals()), "Annotate BED file with gene info")
677722
def gb(args):
    
    from Bio.Alphabet import generic_dna
    try:
        from BCBio import GFF
    except ImportError:
        print("You need to install dep first: $ easy_install bcbio-gff", file=sys.stderr)

    p = OptionParser(gb.__doc__)
    opts, args = p.parse_args(args)

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

    gff_file, fasta_file = args
    pf = op.splitext(gff_file)[0]
    out_file = pf + ".gb"
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")
31071
def from_seqfeature(s, **kwargs):
    
    source = s.qualifiers.get('source', '.')[0]
    score = s.qualifiers.get('score', '.')[0]
    seqid = s.qualifiers.get('seqid', '.')[0]
    frame = s.qualifiers.get('frame', '.')[0]
    strand = _feature_strand[s.strand]

    # BioPython parses 1-based GenBank positions into 0-based for use within
    # Python.  We need to convert back to 1-based GFF format here.
    start = s.location.start.position + 1
    stop = s.location.end.position
    featuretype = s.type
    id = s.id
    attributes = dict(s.qualifiers)
    attributes.pop('source', '.')
    attributes.pop('score', '.')
    attributes.pop('seqid', '.')
    attributes.pop('frame', '.')
    return Feature(seqid, source, featuretype, start, stop, score, strand,
                   frame, attributes, id=id, **kwargs)
303752
def fromgtf(args):
    
    p = OptionParser(fromgtf.__doc__)
    p.add_option("--transcript_id", default="transcript_id",
                 help="Field name for transcript [default: %default]")
    p.add_option("--gene_id", default="gene_id",
                 help="Field name for gene [default: %default]")
    p.add_option("--augustus", default=False, action="store_true",
                 help="Input is AUGUSTUS gtf [default: %default]")
    p.set_home("augustus")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    gtffile, = args
    outfile = opts.outfile
    if opts.augustus:
        ahome = opts.augustus_home
        s = op.join(ahome, "scripts/gtf2gff.pl")
        cmd = "{0} --gff3 < {1} --out={2}".format(s, gtffile, outfile)
        sh(cmd)
        return

    gff = Gff(gtffile)
    fw = must_open(outfile, "w")
    transcript_id = opts.transcript_id
    gene_id = opts.gene_id
    nfeats = 0
    for g in gff:
        if g.type in ("transcript", "mRNA"):
            g.type = "mRNA"
            g.update_tag(transcript_id, "ID")
            g.update_tag("mRNA", "ID")
            g.update_tag(gene_id, "Parent")
            g.update_tag("Gene", "Parent")
        elif g.type in ("exon", "CDS") or "UTR" in g.type:
            g.update_tag("transcript_id", "Parent")
            g.update_tag(g.type, "Parent")
        elif g.type == "gene":
            g.update_tag(gene_id, "ID")
            g.update_tag("Gene", "ID")
        else:
            assert 0, "Don't know how to deal with {0}".format(g.type)

        g.update_attributes()
        print(g, file=fw)
        nfeats += 1

    logging.debug("A total of {0} features written.".format(nfeats))
1231139
def add_genes(in_file, data, max_distance=10000, work_dir=None):
    
    gene_file = regions.get_sv_bed(data, "exons", out_dir=os.path.dirname(in_file))
    if gene_file and utils.file_exists(in_file):
        out_file = "%s-annotated.bed" % utils.splitext_plus(in_file)[0]
        if work_dir:
            out_file = os.path.join(work_dir, os.path.basename(out_file))
        if not utils.file_uptodate(out_file, in_file):
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            with file_transaction(data, out_file) as tx_out_file:
                _add_genes_to_bed(in_file, gene_file, fai_file, tx_out_file, data, max_distance)
        return out_file
    else:
        return in_file
144001
def load_transcript_fpkm_dict_from_gtf(
        gtf_path,
        transcript_id_column_name="reference_id",
        fpkm_column_name="FPKM",
        feature_column_name="feature"):
    
    df = gtfparse.read_gtf(
        gtf_path,
        column_converters={fpkm_column_name: float})
    transcript_ids = _get_gtf_column(transcript_id_column_name, gtf_path, df)
    fpkm_values = _get_gtf_column(fpkm_column_name, gtf_path, df)
    features = _get_gtf_column(feature_column_name, gtf_path, df)
    logging.info("Loaded %d rows from %s" % (len(transcript_ids), gtf_path))
    logging.info("Found %s transcript entries" % sum(
        feature == "transcript" for feature in features))
    result = {
        transcript_id: float(fpkm)
        for (transcript_id, fpkm, feature)
        in zip(transcript_ids, fpkm_values, features)
        if (
            (transcript_id is not None) and
            (len(transcript_id) > 0) and
            (feature == "transcript")
        )
    }
    logging.info("Keeping %d transcript rows with reference IDs" % (
        len(result),))
    return result
298556
def CDSReader( fh, format='gff' ):
    

    known_formats = ( 'gff', 'gtf', 'bed')
    if format not in known_formats: 
        print('%s format not in %s' % (format, ",".join( known_formats )), file=sys.stderr)
        raise Exception('?')

    if format == 'bed':
        for line in fh:    
            f = line.strip().split()
            chrom = f[0]
            chrom_start = int(f[1])
            name = f[4]
            strand = f[5]
            cdsStart = int(f[6])
            cdsEnd = int(f[7])
            blockCount = int(f[9])
            blockSizes = [ int(i) for i in f[10].strip(',').split(',') ]
            blockStarts = [ chrom_start + int(i) for i in f[11].strip(',').split(',') ]

            # grab cdsStart - cdsEnd
            cds_exons = []
            cds_seq = ''
            genome_seq_index = []
            for base,offset in zip( blockStarts, blockSizes ):
                if (base + offset) < cdsStart: continue
                if base > cdsEnd: continue
                exon_start = max( base, cdsStart )
                exon_end = min( base+offset, cdsEnd ) 
                cds_exons.append( (exon_start, exon_end) )
            

    genelist = {}
    grouplist = []
    if format == 'gff' or format == 'gtf':
        for line in fh:
            if line.startswith('#'): continue
            fields = line.strip().split('\t')
            if len( fields ) < 9: continue
            if fields[2] not in ('CDS', 'stop_codon', 'start_codon'): continue

            # fields

            chrom = fields[0]
            ex_st = int( fields[3] ) - 1 # make zero-centered
            ex_end = int( fields[4] ) #+ 1 # make exclusive
            strand = fields[6]

            if format == 'gtf':
                group = fields[8].split(';')[0]
            else:
                group = fields[8]

            if group not in grouplist: grouplist.append( group )
            if group not in genelist:
                genelist[group] = (chrom, strand, [])
            
            genelist[group][2].append( ( ex_st, ex_end ) )

        sp = lambda a,b: cmp( a[0], b[0] )

        #for gene in genelist.values():
        for gene in grouplist:
            chrom, strand, cds_exons = genelist[ gene ]
            seqlen = sum([ a[1]-a[0] for a in cds_exons ])
            overhang = seqlen % 3
            if overhang > 0:
                #print >>sys.stderr, "adjusting ", gene  
                if strand == '+': 
                    cds_exons[-1] = ( cds_exons[-1][0], cds_exons[-1][1] - overhang )
                else:
                    cds_exons[0] = ( cds_exons[0][0] + overhang, cds_exons[0][1] )
            cds_exons = bitset_union( cds_exons )
            yield chrom, strand, cds_exons, gene
1230699
def g_to_n(self, g_interval):
        """"""

        # This code is extremely convoluted. To begin with, it
        # confuses interbase intervals for a *single* hgvs position
        # with intervals of *two* hgvs positions.  This is the origin
        # of bug #205 (now closed). Although the bug is itself fixed,
        # the mappers really need an overhaul. See #208.

        def _hgvs_offset(g_position, grs, gre, strand):
            """ Calculates the HGVS coordinate offset from a given genomic position"""
            if g_position == grs or g_position == gre:
                return 0
            mid = (grs + gre) / 2
            if g_position < mid or (g_position == mid and strand == 1):
                offset = g_position - grs
            else:
                offset = g_position - gre
            return offset

        def map_g_to_n_pos(pos):
            g_ci = _hgvs_coord_to_ci(pos, pos)
            # frs, fre = (f)orward (r)na (s)tart & (e)nd; forward w.r.t. genome
            frs, fre = self.im.map_ref_to_tgt(
                g_ci[0] - self.gc_offset, g_ci[1] - self.gc_offset, max_extent=False)
            if self.strand == -1:
                frs, fre = self.tgt_len - fre, self.tgt_len - frs
            # get BaseOffsetPosition information for n.
            if self.strand == 1:
                grs, gre = self.im.map_tgt_to_ref(frs, fre, max_extent=False)
            elif self.strand == -1:
                grs, gre = self.im.map_tgt_to_ref((self.tgt_len - fre), (self.tgt_len - frs),
                                                  max_extent=False)
            grs, gre = grs + self.gc_offset, gre + self.gc_offset
            # 0-width interval indicates an intron.  Need to calculate offsets but we"re are in ci coordinates
            # requires adding 1 strategically to get the HGVS position
            # (shift coordinates to 3' end of the ref nucleotide)

            if frs == fre:
                start_offset = _hgvs_offset(g_ci[0] + 1, grs, gre + 1, self.strand)
                end_offset = _hgvs_offset(g_ci[1], grs, gre + 1, self.strand)
            else:
                start_offset = _hgvs_offset(g_ci[0], grs, gre, self.strand)
                end_offset = _hgvs_offset(g_ci[1], grs, gre, self.strand)
            if self.strand == -1:
                start_offset, end_offset = self.strand * end_offset, self.strand * start_offset
            if start_offset > 0:
                frs -= 1
            if end_offset < 0:
                fre += 1
            start_base = _ci_to_hgvs_coord(frs, fre)[0]
            end_base = _ci_to_hgvs_coord(frs, fre)[1]
            return start_base, start_offset, end_base, end_offset

        start_bo = map_g_to_n_pos(g_interval.start.base)[0:2]
        end_bo = start_bo if g_interval.start.base == g_interval.end.base else map_g_to_n_pos(
            g_interval.end.base)[2:4]
        if self.strand == -1:
            start_bo, end_bo = end_bo, start_bo

        return hgvs.location.BaseOffsetInterval(
            start=hgvs.location.BaseOffsetPosition(
                base=start_bo[0], offset=start_bo[1], datum=Datum.SEQ_START),
            end=hgvs.location.BaseOffsetPosition(
                base=end_bo[0], offset=end_bo[1], datum=Datum.SEQ_START),
            uncertain=g_interval.uncertain)
90783
def guess_infer_extent(gtf_file):
    
    _, ext = os.path.splitext(gtf_file)
    tmp_out = tempfile.NamedTemporaryFile(suffix=".gtf", delete=False).name
    with open(tmp_out, "w") as out_handle:
        count = 0
        in_handle = utils.open_gzipsafe(gtf_file)
        for line in in_handle:
            if count > 1000:
                break
            out_handle.write(line)
            count += 1
        in_handle.close()
    db = gffutils.create_db(tmp_out, dbfn=":memory:", infer_gene_extent=False)
    os.remove(tmp_out)
    features = [x for x in db.featuretypes()]
    if "gene" in features and "transcript" in features:
        return False
    else:
        return True
142964
def index_transcriptome(gtf_file, ref_file, data):
    
    gtf_fasta = gtf.gtf_to_fasta(gtf_file, ref_file)
    bowtie2_index = os.path.splitext(gtf_fasta)[0]
    bowtie2_build = config_utils.get_program("bowtie2", data["config"]) + "-build"
    cmd = "{bowtie2_build} --offrate 1 {gtf_fasta} {bowtie2_index}".format(**locals())
    message = "Creating transcriptome index of %s with bowtie2." % (gtf_fasta)
    do.run(cmd, message)
    return bowtie2_index