gtf2juncs

#!/usr/local/bin/python
"""
Takes a GTF file of transcripts with exon features and outputs a
list of all splice junctions in juncs format.
The file must not have non-exon features in it or else this will create
some false junctions.
This does not filter for canonical splice junctions, since it is
sequence-unaware.
This also can report all exon skipping events.
It also throws out all transcripts with less than 3 exons.
"""

from gtfUtils import parseGTFlineToDict, aggregateFeaturesByTranscript
from gtfUtils import aggregateFeaturesByGene, mergeOverlappedExons
import logging
import sys

JUNC_HEADER = ["seqname", "left", "right", "id", "score", "strand"]
IE_HEADER = ["exon_id", "inc", "exc"]

logging.basicConfig(level=logging.INFO)

def _read_gtf():
    logging.info("Processing GTF input.")
    gtflines = []
    numlines = 0
    for line in sys.stdin:
        numlines = numlines + 1
        gtflines.append(parseGTFlineToDict(line))
    logging.info("Read %d lines." %(numlines))

    return gtflines

def _make_junction(exon1, exon2):
    """
    makes a junction of the form junc
    example: chr1:5:200:+
    """
    return ":".join([exon1['seqname'], exon1['end'], exon2['start'],
                     exon1['strand']])

def _make_junction_line(exon1, exon2):
    id = _make_junction(exon1, exon2)
    # convert to 0 based
    end_of_one = str(int(exon1['end']) - 1)
    start_of_two = str(int(exon2['start']) - 1)
    
    return "\t".join([exon1['seqname'], end_of_one, start_of_two, id, ".",
                      exon1['strand']]) + "\n"    

def _make_tophat_junction_line(exon1, exon2):
    end_of_one = str(int(exon1['end']))
    start_of_two = str(int(exon2['start']) - 1)
    if 'exon_id' not in exon1:
        return "\t".join([exon1['seqname'], end_of_one, start_of_two, exon1['transcript_id'], ".", exon1['strand']])
    else:
        junc_id = ".".join([exon1['exon_id'], exon2['exon_id']])
        return "\t".join([exon1['seqname'], end_of_one, start_of_two, junc_id, ".", exon1['strand']])

def _normal_mode(transcripts):
    numjuncs = 0
    for transcript, exons in transcripts.items():
        last_pair = len(exons) - 1
        for i in range(0, last_pair):
            sys.stdout.write(_make_junction(exons[i], exons[i + 1]) + "\n")
            numjuncs = numjuncs + 1

    logging.info("Wrote %d junctions specified from the GTF." %(numjuncs))

def _tophat_mode(transcripts):
    numjuncs = 0
    for transcript, exons in transcripts.items():
        last_pair = len(exons) - 1
        for i in range(0, last_pair):
            sys.stdout.write(_make_tophat_junction_line(exons[i], exons[i + 1]) + "\n")
            numjuncs = numjuncs + 1

    logging.info("Wrote %d junctions specified form the GTF." %(numjuncs))

def _single_skip_mode(transcripts):
    """
    currently not used, calculates just the single skips
    """
    numjuncs = 0
    numskips = 0
    for transcript, exons in transcripts.items():
        last_pair = len(exons) - 1
        for i in range(0, last_pair):
            sys.stdout.write(_make_junction_line(exons[i], exons[i + 1]))
            numjuncs = numjuncs + 1
            if (i + 2) < len(exons):
                sys.stdout.write(_make_junction_line(exons[i], exons[i + 2]))
                numskips = numskips + 1

    logging.info("Wrote %d junctions specified from the GTF." %(numjuncs))
    logging.info("Wrote %d single exon skipping events." %(numskips))

def _skip_mode(transcripts):
    numjuncs = 0
    numskips = 0
    for transcript, exons in transcripts.items():
        last_pair = len(exons) - 1
        for i in range(0, last_pair):
            sys.stdout.write(_make_junction_line(exons[i], exons[i + 1]))
            numjuncs = numjuncs + 1
            for j in range(i + 2, len(exons)):
                sys.stdout.write(_make_junction_line(exons[i], exons[j]))
                numskips = numskips + 1

    logging.info("Wrote %d junctions specified from the GTF." %(numjuncs))
    logging.info("Wrote %d exon skipping events." %(numskips))

def _build_ie_id(exon):
    return exon['seqname'] + ":" + exon['start'] + "-" + exon['end']
                
def _ie_mode(transcripts):

    sys.stdout.write(",".join(IE_HEADER) + "\n")
    for transcript, exons in transcripts.items():
        # skip if the number of exons is less than 3, since skipping is
        # meaningless if that is true
        if len(exons) < 3:
            continue

        last_skip = len(exons) - 2
        
        for i in range(1, last_skip):
            """
            calculate inclusion junctions for exon i (junctions that
            connect to i)
            """
            ie_id = _build_ie_id(exons[i])
            inc = []
            for j in range(0, len(exons)):
                if j < i:
                    inc.append(_make_junction(exons[j], exons[i]))
                if i < j:
                    inc.append(_make_junction(exons[i], exons[j]))
            inc_str = ",".join(inc)

            """
            calculate exclusion junctions for exon i (junctions that
            skip over exon i
            """
            exc = []
            for j in range(0, i):
                for k in range(i + 1, len(exons)):
                    exc.append(_make_junction(exons[j], exons[k]))
            exc_str = ",".join(exc)

            sys.stdout.write("\t".join([ie_id, inc_str, exc_str]) + "\n")

def main():

    description = "Outputs all the splice junctions specified by a GTF file."
    usage = "gtf2juncs [--skip]:\n\n\t" \
            "options:\n\t" \
            "--tophat: output these in the tophat raw juncs format\n " \
            "--skip: enumerate all possible single exon skipping events for " \
            "each transcript\n\t" \
            "--ie: format as inclusion/exclusion events\n\t" \
            "--help: display help"

    if sys.argv.count("-h") or sys.argv.count("--help"):
        print description
        print usage
        exit(1)

    gtflines = _read_gtf()

    
    if sys.argv.count("--skip"):
        #transcripts = aggregateFeaturesByTranscript(gtflines)
        genes = aggregateFeaturesByGene(gtflines)
        genes = mergeOverlappedExons(genes)
        logging.info("Will calculate all exon skipping events.")
        _skip_mode(genes)
        exit(1)

    if sys.argv.count("--ie"):
        genes = aggregateFeaturesByGene(gtflines)
        genes = mergeOverlappedExons(genes)
        logging.info("Writing out inclusion/exclusion events file.")
        _ie_mode(genes)
        exit(1)

    transcripts = aggregateFeaturesByTranscript(gtflines)

    if sys.argv.count("--tophat"):
        _tophat_mode(transcripts)
        exit(1)
        
    _normal_mode(transcripts)

if __name__ == "__main__":
    main()