<filename>taco/lib/transcript.py
'''
Created on Nov 30, 2011

@author: mkiyer

AssemblyLine: transcriptome meta-assembly from RNA-Seq

Copyright (C) 2012 <NAME>

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
'''
import collections
import bisect
import os
from gtf_comp import parse_loci, GTFFeature, GTFError

# attributes
TRANSCRIPT_ID = "transcript_id"
GENE_ID = "gene_id"

# strand constants
POS_STRAND = 0
NEG_STRAND = 1
NO_STRAND = 2

# strand format conversion
_strand_int_to_str = ["+", "-", "."]
_strand_str_to_int = {"+": POS_STRAND,
                      "-": NEG_STRAND,
                      ".": NO_STRAND}
def strand_str_to_int(strand):
    return _strand_str_to_int[strand]
def strand_int_to_str(strand):
    return _strand_int_to_str[strand]

def merge_strand(strand1, strand2):
    if strand1 == strand2:
        return strand1
    elif strand1 == NO_STRAND:
        return strand2
    elif strand2 == NO_STRAND:
        return strand1
    else:
        assert False

def cmp_strand(a, b):
    if (a == NO_STRAND) or (b == NO_STRAND):
        return True
    return a == b

def interval_overlap(a, b):
    return (a.start < b.end) and (b.start < a.end)

def interval_overlap_threshold(a, b, d=0):
    return (a.start < b.end + d) and (b.start < a.end + d)

class Exon(object):
    __slots__ = ('start', 'end')
    def __init__(self, start, end):
        self.start = start
        self.end = end
    def __str__(self):
        return ("%s-%s" % (self.start, self.end))
    def __repr__(self):
        return ("<%s(start=%d end=%d)>" %
                (self.__class__.__name__, self.start, self.end))
    def __eq__(self, other):
        return (self.start == other.start) and (self.end == other.end)
    def __ne__(self, other):
        return (self.start != other.start) or (self.end != other.end)
    def __lt__(self, other):
        return self.start < other.start
    def __le__(self, other):
        return self.start <= other.start
    def __gr__(self, other):
        return self.start > other.start
    def __ge__(self, other):
        return self.start >= other.start
    def __hash__(self):
        return (self.start << 16) | (self.end)
    def is_overlapping(self, other):
        return interval_overlap(self, other)

class Transcript(object):
    __slots__ = ('chrom', 'start', 'end', 'strand', 'score', 'exons', 'attrs')

    def __init__(self):
        self.chrom = None
        self.start = -1
        self.end = -1
        self.strand = NO_STRAND
        self.score = 0.0
        self.exons = None
        self.attrs = {}

    def __str__(self):
        return ("<%s(chrom='%s', start='%d', end='%d', strand='%s', "
                "score='%s' exons='%s', attrs='%s'" %
                (self.__class__.__name__, self.chrom, self.start, self.end,
                 strand_int_to_str(self.strand), str(self.score), self.exons,
                 self.attrs))

    @property
    def length(self):
        return sum((e.end - e.start) for e in self.exons)

    def iterintrons(self):
        #e1 = self.exons[0]
        #for e2 in self.exons[1:]:
        #    yield (e1.end,e2.start)
        #    e1 = e2
        e1 = self.exons[0]
        for j in xrange(1, len(self.exons)):
            e2 = self.exons[j]
            yield e1.end, e2.start
            e1 = e2

    def introns(self):
        return list(self.iterintrons())

    def to_bed12(self):
        block_sizes = []
        block_starts = []
        for e0, e1 in self.exons:
            block_starts.append(e0 - self.tx_start)
            block_sizes.append(e1 - e0)
        # write
        s = '\t'.join([self.chrom,
                       str(self.start),
                       str(self.end),
                       str(self.attrs["transcript_id"]),
                       '0',
                       strand_int_to_str(self.strand),
                       str(self.start),
                       str(self.start),
                       '0',
                       str(len(self.exons)),
                       ','.join(map(str,block_sizes)) + ',',
                       ','.join(map(str,block_starts)) + ','])
        return s

    def to_gtf_features(self, source=None, score=1000):
        if source is None:
            source = 'assemblyline'
        # transcript feature
        f = GTFFeature()
        f.seqid = self.chrom
        f.source = source
        f.feature_type = 'transcript'
        f.start = self.start
        f.end = self.end
        f.score = score
        f.strand = strand_int_to_str(self.strand)
        f.phase = '.'
        f.attrs = self.attrs
        features = [f]
        # exon features
        for i,e in enumerate(self.exons):
            f = GTFFeature()
            f.seqid = self.chrom
            f.source = source
            f.feature_type = 'exon'
            f.start = e.start
            f.end = e.end
            f.score = score
            f.strand = strand_int_to_str(self.strand)
            f.phase = '.'
            f.attrs = self.attrs.copy()
            f.attrs["exon_number"] = i
            features.append(f)
        return features

def transcripts_from_gtf_lines(lines, attr_defs=None, log=None):
    transcripts = collections.OrderedDict()
    for line in lines:
        feature = GTFFeature.from_string(line, attr_defs, log)
        if feature.feature_type == 'gene':
            continue
        t_id = feature.attrs["transcript_id"]
        ref = feature.attrs["ref"]
        if (t_id, ref) not in transcripts.keys():
            if feature.feature_type != "transcript":
                raise GTFError("Feature type '%s' found before 'transcript' record: %s" %
                               (feature.feature_type, str(feature)))
            t = Transcript()
            t.chrom = feature.seqid
            t.start = feature.start
            t.end = feature.end
            # convert from string strand notation ("+", "-", ".")
            # to integer (0, 1)
            t.strand = strand_str_to_int(feature.strand)
            t.exons = []
            t.attrs = feature.attrs
            transcripts[(t_id, ref)] = t
        else:
            t = transcripts[(t_id, ref)]
        if feature.feature_type == "exon":
            t.exons.append(Exon(feature.start, feature.end))
    # sort transcript exons by genomic position
    for t in transcripts.itervalues():
        t.exons.sort()
        t.end = t.exons[-1].end
    return transcripts.values()

def parse_gtf(fileh, attr_defs=None):
    for locus_features in parse_loci(fileh):
        yield transcripts_from_gtf_lines(locus_features, attr_defs, fileh.name)

def genes_from_gtf_lines(lines, attr_defs=None):
    transcripts = collections.OrderedDict()
    for line in lines:
        feature = GTFFeature.from_string(line, attr_defs)
        if feature.feature_type != 'gene':
            continue
        t_id = feature.attrs["gene_id"]

        if t_id not in transcripts:
            if feature.feature_type != "gene":
                raise GTFError("Feature type '%s' found before 'gene' record: %s" %
                               (feature.feature_type, str(feature)))
            t = Transcript()
            t.chrom = feature.seqid
            t.start = feature.start
            t.end = feature.end
            # convert from string strand notation ("+", "-", ".")
            # to integer (0, 1)
            t.strand = strand_str_to_int(feature.strand)
            t.attrs = feature.attrs
            transcripts[t_id] = t
        else:
            t = transcripts[t_id]
    return transcripts.values()

def parse_gtf_gene(fileh, attr_defs=None):
    for locus_features in parse_loci(fileh):
        yield genes_from_gtf_lines(locus_features, attr_defs)

def find_exon_boundaries(transcripts):
    '''
    input: a list of transcripts (not Node objects, these are transcripts)
    parsed directly from an input file and not added to an isoform graph
    yet.

    output: sorted list of exon boundaries
    '''
    exon_boundaries = set()
    # first add introns to the graph and keep track of
    # all intron boundaries
    for transcript in transcripts:
        # add transcript exon boundaries
        for exon in transcript.exons:
            # keep track of positions where introns can be joined to exons
            exon_boundaries.add(exon.start)
            exon_boundaries.add(exon.end)
    # sort the intron boundary positions and add them to interval trees
    return sorted(exon_boundaries)


def split_exon(exon, boundaries):
    """
    partition the exon given list of node boundaries

    generator yields (start,end) intervals for exon
    """
    if exon.start == exon.end:
        return
    # find the indexes into the intron boundaries list that
    # border the exon.  all the indexes in between these two
    # are overlapping the exon and we must use them to break
    # the exon into pieces
    start_ind = bisect.bisect_right(boundaries, exon.start)
    end_ind = bisect.bisect_left(boundaries, exon.end)
    if start_ind == end_ind:
        yield exon.start, exon.end
    else:
        yield exon.start, boundaries[start_ind]
        for j in xrange(start_ind, end_ind-1):
            yield boundaries[j], boundaries[j+1]
        yield boundaries[end_ind-1], exon.end

def split_exons(t, boundaries):
    # split exons that cross boundaries and to get the
    # nodes in the transcript path
    for exon in t.exons:
        for start, end in split_exon(exon, boundaries):
            yield start, end

<filename>taco/lib/transcript.py
'''
Created on Nov 30, 2011

@author: mkiyer

AssemblyLine: transcriptome meta-assembly from RNA-Seq

Copyright (C) 2012 <NAME>

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
'''
import collections
import bisect
import os
from gtf_comp import parse_loci, GTFFeature, GTFError

# attributes
TRANSCRIPT_ID = "transcript_id"
GENE_ID = "gene_id"

# strand constants
POS_STRAND = 0
NEG_STRAND = 1
NO_STRAND = 2

# strand format conversion
_strand_int_to_str = ["+", "-", "."]
_strand_str_to_int = {"+": POS_STRAND,
                      "-": NEG_STRAND,
                      ".": NO_STRAND}
def strand_str_to_int(strand):
    return _strand_str_to_int[strand]
def strand_int_to_str(strand):
    return _strand_int_to_str[strand]

def merge_strand(strand1, strand2):
    if strand1 == strand2:
        return strand1
    elif strand1 == NO_STRAND:
        return strand2
    elif strand2 == NO_STRAND:
        return strand1
    else:
        assert False

def cmp_strand(a, b):
    if (a == NO_STRAND) or (b == NO_STRAND):
        return True
    return a == b

def interval_overlap(a, b):
    return (a.start < b.end) and (b.start < a.end)

def interval_overlap_threshold(a, b, d=0):
    return (a.start < b.end + d) and (b.start < a.end + d)

class Exon(object):
    __slots__ = ('start', 'end')
    def __init__(self, start, end):
        self.start = start
        self.end = end
    def __str__(self):
        return ("%s-%s" % (self.start, self.end))
    def __repr__(self):
        return ("<%s(start=%d end=%d)>" %
                (self.__class__.__name__, self.start, self.end))
    def __eq__(self, other):
        return (self.start == other.start) and (self.end == other.end)
    def __ne__(self, other):
        return (self.start != other.start) or (self.end != other.end)
    def __lt__(self, other):
        return self.start < other.start
    def __le__(self, other):
        return self.start <= other.start
    def __gr__(self, other):
        return self.start > other.start
    def __ge__(self, other):
        return self.start >= other.start
    def __hash__(self):
        return (self.start << 16) | (self.end)
    def is_overlapping(self, other):
        return interval_overlap(self, other)

class Transcript(object):
    __slots__ = ('chrom', 'start', 'end', 'strand', 'score', 'exons', 'attrs')

    def __init__(self):
        self.chrom = None
        self.start = -1
        self.end = -1
        self.strand = NO_STRAND
        self.score = 0.0
        self.exons = None
        self.attrs = {}

    def __str__(self):
        return ("<%s(chrom='%s', start='%d', end='%d', strand='%s', "
                "score='%s' exons='%s', attrs='%s'" %
                (self.__class__.__name__, self.chrom, self.start, self.end,
                 strand_int_to_str(self.strand), str(self.score), self.exons,
                 self.attrs))

    @property
    def length(self):
        return sum((e.end - e.start) for e in self.exons)

    def iterintrons(self):
        #e1 = self.exons[0]
        #for e2 in self.exons[1:]:
        #    yield (e1.end,e2.start)
        #    e1 = e2
        e1 = self.exons[0]
        for j in xrange(1, len(self.exons)):
            e2 = self.exons[j]
            yield e1.end, e2.start
            e1 = e2

    def introns(self):
        return list(self.iterintrons())

    def to_bed12(self):
        block_sizes = []
        block_starts = []
        for e0, e1 in self.exons:
            block_starts.append(e0 - self.tx_start)
            block_sizes.append(e1 - e0)
        # write
        s = '\t'.join([self.chrom,
                       str(self.start),
                       str(self.end),
                       str(self.attrs["transcript_id"]),
                       '0',
                       strand_int_to_str(self.strand),
                       str(self.start),
                       str(self.start),
                       '0',
                       str(len(self.exons)),
                       ','.join(map(str,block_sizes)) + ',',
                       ','.join(map(str,block_starts)) + ','])
        return s

    def to_gtf_features(self, source=None, score=1000):
        if source is None:
            source = 'assemblyline'
        # transcript feature
        f = GTFFeature()
        f.seqid = self.chrom
        f.source = source
        f.feature_type = 'transcript'
        f.start = self.start
        f.end = self.end
        f.score = score
        f.strand = strand_int_to_str(self.strand)
        f.phase = '.'
        f.attrs = self.attrs
        features = [f]
        # exon features
        for i,e in enumerate(self.exons):
            f = GTFFeature()
            f.seqid = self.chrom
            f.source = source
            f.feature_type = 'exon'
            f.start = e.start
            f.end = e.end
            f.score = score
            f.strand = strand_int_to_str(self.strand)
            f.phase = '.'
            f.attrs = self.attrs.copy()
            f.attrs["exon_number"] = i
            features.append(f)
        return features

def transcripts_from_gtf_lines(lines, attr_defs=None, log=None):
    transcripts = collections.OrderedDict()
    for line in lines:
        feature = GTFFeature.from_string(line, attr_defs, log)
        if feature.feature_type == 'gene':
            continue
        t_id = feature.attrs["transcript_id"]
        ref = feature.attrs["ref"]
        if (t_id, ref) not in transcripts.keys():
            if feature.feature_type != "transcript":
                raise GTFError("Feature type '%s' found before 'transcript' record: %s" %
                               (feature.feature_type, str(feature)))
            t = Transcript()
            t.chrom = feature.seqid
            t.start = feature.start
            t.end = feature.end
            # convert from string strand notation ("+", "-", ".")
            # to integer (0, 1)
            t.strand = strand_str_to_int(feature.strand)
            t.exons = []
            t.attrs = feature.attrs
            transcripts[(t_id, ref)] = t
        else:
            t = transcripts[(t_id, ref)]
        if feature.feature_type == "exon":
            t.exons.append(Exon(feature.start, feature.end))
    # sort transcript exons by genomic position
    for t in transcripts.itervalues():
        t.exons.sort()
        t.end = t.exons[-1].end
    return transcripts.values()

def parse_gtf(fileh, attr_defs=None):
    for locus_features in parse_loci(fileh):
        yield transcripts_from_gtf_lines(locus_features, attr_defs, fileh.name)

def genes_from_gtf_lines(lines, attr_defs=None):
    transcripts = collections.OrderedDict()
    for line in lines:
        feature = GTFFeature.from_string(line, attr_defs)
        if feature.feature_type != 'gene':
            continue
        t_id = feature.attrs["gene_id"]

        if t_id not in transcripts:
            if feature.feature_type != "gene":
                raise GTFError("Feature type '%s' found before 'gene' record: %s" %
                               (feature.feature_type, str(feature)))
            t = Transcript()
            t.chrom = feature.seqid
            t.start = feature.start
            t.end = feature.end
            # convert from string strand notation ("+", "-", ".")
            # to integer (0, 1)
            t.strand = strand_str_to_int(feature.strand)
            t.attrs = feature.attrs
            transcripts[t_id] = t
        else:
            t = transcripts[t_id]
    return transcripts.values()

def parse_gtf_gene(fileh, attr_defs=None):
    for locus_features in parse_loci(fileh):
        yield genes_from_gtf_lines(locus_features, attr_defs)

def find_exon_boundaries(transcripts):
    '''
    input: a list of transcripts (not Node objects, these are transcripts)
    parsed directly from an input file and not added to an isoform graph
    yet.

    output: sorted list of exon boundaries
    '''
    exon_boundaries = set()
    # first add introns to the graph and keep track of
    # all intron boundaries
    for transcript in transcripts:
        # add transcript exon boundaries
        for exon in transcript.exons:
            # keep track of positions where introns can be joined to exons
            exon_boundaries.add(exon.start)
            exon_boundaries.add(exon.end)
    # sort the intron boundary positions and add them to interval trees
    return sorted(exon_boundaries)


def split_exon(exon, boundaries):
    """
    partition the exon given list of node boundaries

    generator yields (start,end) intervals for exon
    """
    if exon.start == exon.end:
        return
    # find the indexes into the intron boundaries list that
    # border the exon.  all the indexes in between these two
    # are overlapping the exon and we must use them to break
    # the exon into pieces
    start_ind = bisect.bisect_right(boundaries, exon.start)
    end_ind = bisect.bisect_left(boundaries, exon.end)
    if start_ind == end_ind:
        yield exon.start, exon.end
    else:
        yield exon.start, boundaries[start_ind]
        for j in xrange(start_ind, end_ind-1):
            yield boundaries[j], boundaries[j+1]
        yield boundaries[end_ind-1], exon.end

def split_exons(t, boundaries):
    # split exons that cross boundaries and to get the
    # nodes in the transcript path
    for exon in t.exons:
        for start, end in split_exon(exon, boundaries):
            yield start, end

Created on Nov 30, 2011  @author: mkiyer  AssemblyLine: transcriptome meta-assembly from RNA-Seq  Copyright (C) 2012 <NAME>  Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:  The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. # attributes # strand constants # strand format conversion #e1 = self.exons[0] #for e2 in self.exons[1:]: #    yield (e1.end,e2.start) #    e1 = e2 # write # transcript feature # exon features # convert from string strand notation ("+", "-", ".") # to integer (0, 1) # sort transcript exons by genomic position # convert from string strand notation ("+", "-", ".") # to integer (0, 1) input: a list of transcripts (not Node objects, these are transcripts)     parsed directly from an input file and not added to an isoform graph     yet.      output: sorted list of exon boundaries # first add introns to the graph and keep track of # all intron boundaries # add transcript exon boundaries # keep track of positions where introns can be joined to exons # sort the intron boundary positions and add them to interval trees partition the exon given list of node boundaries      generator yields (start,end) intervals for exon # find the indexes into the intron boundaries list that # border the exon.  all the indexes in between these two # are overlapping the exon and we must use them to break # the exon into pieces # split exons that cross boundaries and to get the # nodes in the transcript path