annotate tools/filters/ucsc_gene_bed_to_intron_bed.py @ 1:cdcb0ce84a1b

Uploaded

#!/usr/bin/env python

"""
Read a table dump in the UCSC gene table format and print a tab separated
list of intervals corresponding to requested features of each gene.

usage: ucsc_gene_table_to_intervals.py [options]

options:
  -h, --help                  show this help message and exit
  -rREGION, --region=REGION
                              Limit to region: one of coding, utr3, utr5, transcribed [default]
  -e, --exons                 Only print intervals overlapping an exon
  -i, --input=inputfile       input file
  -o, --output=outputfile     output file
"""

import optparse, string, sys

assert sys.version_info[:2] >= ( 2, 4 )

def main():

    # Parse command line    
    parser = optparse.OptionParser( usage="%prog [options] " )
    #parser.add_option( "-r", "--region", dest="region", default="transcribed",
    #                   help="Limit to region: one of coding, utr3, utr5, transcribed [default]" )
    #parser.add_option( "-e", "--exons",  action="store_true", dest="exons",
    #                   help="Only print intervals overlapping an exon" )
    parser.add_option( "-s", "--strand",  action="store_true", dest="strand",
                       help="Print strand after interval" )
    parser.add_option( "-i", "--input",  dest="input",  default=None,
                       help="Input file" )
    parser.add_option( "-o", "--output", dest="output", default=None,
                       help="Output file" )
    options, args = parser.parse_args()
    #assert options.region in ( 'coding', 'utr3', 'utr5', 'transcribed' ), "Invalid region argument"
    
    try:
        out_file = open (options.output,"w")
    except:
        print >> sys.stderr, "Bad output file."
        sys.exit(0)
    
    try:
        in_file = open (options.input)
    except:
        print >> sys.stderr, "Bad input file."
        sys.exit(0)
    
    #print "Region:", options.region+";"
    #print "Only overlap with Exons:",
    #if options.exons:
    #    print "Yes"
    #else:
    #    print "No"
    
    # Read table and handle each gene
    
    for line in in_file:
        try:
	    #print ("len: %d", len(line))
            if line[0:1] == "#":
                continue
	   
            # Parse fields from gene tabls
            fields = line.split( '\t' )
            chrom     = fields[0]
            tx_start  = int( fields[1] )
            tx_end    = int( fields[2] )
            name      = fields[3]
            strand    = fields[5].replace(" ","_")
            cds_start = int( fields[6] )
            cds_end   = int( fields[7] )
	    	
	    exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
            exon_starts = map((lambda x: x + tx_start ), exon_starts)
            exon_ends = map( int, fields[10].rstrip( ',\n' ).split( ',' ) )
            exon_ends = map((lambda x, y: x + y ), exon_starts, exon_ends);
	    
	    i=0
	    while i < len(exon_starts)-1:
            	intron_starts = exon_ends[i] + 1
		intron_ends = exon_starts[i+1] - 1
		if strand: print_tab_sep(out_file, chrom, intron_starts, intron_ends, name, "0", strand )
                else: print_tab_sep(out_file, chrom, intron_starts, intron_ends )
		i+=1
            # If only interested in exons, print the portion of each exon overlapping
            # the region of interest, otherwise print the span of the region
            
        except:
            continue

def print_tab_sep(out_file, *args ):
    """Print items in `l` to stdout separated by tabs"""
    print >>out_file, string.join( [ str( f ) for f in args ], '\t' )

if __name__ == "__main__": main()