/tools/filters/ucsc_gene_bed_to_intron_bed.py

 1#!/usr/bin/env python
 2
 3"""
 4Read a table dump in the UCSC gene table format and print a tab separated
 5list of intervals corresponding to requested features of each gene.
 6
 7usage: ucsc_gene_table_to_intervals.py [options]
 8
 9options:
10  -h, --help                  show this help message and exit
11  -rREGION, --region=REGION
12                              Limit to region: one of coding, utr3, utr5, transcribed [default]
13  -e, --exons                 Only print intervals overlapping an exon
14  -i, --input=inputfile       input file
15  -o, --output=outputfile     output file
16"""
17
18import optparse, string, sys
19
20assert sys.version_info[:2] >= ( 2, 4 )
21
22def main():
23
24    # Parse command line    
25    parser = optparse.OptionParser( usage="%prog [options] " )
26    #parser.add_option( "-r", "--region", dest="region", default="transcribed",
27    #                   help="Limit to region: one of coding, utr3, utr5, transcribed [default]" )
28    #parser.add_option( "-e", "--exons",  action="store_true", dest="exons",
29    #                   help="Only print intervals overlapping an exon" )
30    parser.add_option( "-s", "--strand",  action="store_true", dest="strand",
31                       help="Print strand after interval" )
32    parser.add_option( "-i", "--input",  dest="input",  default=None,
33                       help="Input file" )
34    parser.add_option( "-o", "--output", dest="output", default=None,
35                       help="Output file" )
36    options, args = parser.parse_args()
37    #assert options.region in ( 'coding', 'utr3', 'utr5', 'transcribed' ), "Invalid region argument"
38    
39    try:
40        out_file = open (options.output,"w")
41    except:
42        print >> sys.stderr, "Bad output file."
43        sys.exit(0)
44    
45    try:
46        in_file = open (options.input)
47    except:
48        print >> sys.stderr, "Bad input file."
49        sys.exit(0)
50    
51    #print "Region:", options.region+";"
52    #print "Only overlap with Exons:",
53    #if options.exons:
54    #    print "Yes"
55    #else:
56    #    print "No"
57    
58    # Read table and handle each gene
59    
60    for line in in_file:
61        try:
62	    #print ("len: %d", len(line))
63            if line[0:1] == "#":
64                continue
65	   
66            # Parse fields from gene tabls
67            fields = line.split( '\t' )
68            chrom     = fields[0]
69            tx_start  = int( fields[1] )
70            tx_end    = int( fields[2] )
71            name      = fields[3]
72            strand    = fields[5].replace(" ","_")
73            cds_start = int( fields[6] )
74            cds_end   = int( fields[7] )
75	    	
76	    exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
77            exon_starts = map((lambda x: x + tx_start ), exon_starts)
78            exon_ends = map( int, fields[10].rstrip( ',\n' ).split( ',' ) )
79            exon_ends = map((lambda x, y: x + y ), exon_starts, exon_ends);
80	    
81	    i=0
82	    while i < len(exon_starts)-1:
83            	intron_starts = exon_ends[i] + 1
84		intron_ends = exon_starts[i+1] - 1
85		if strand: print_tab_sep(out_file, chrom, intron_starts, intron_ends, name, "0", strand )
86                else: print_tab_sep(out_file, chrom, intron_starts, intron_ends )
87		i+=1
88            # If only interested in exons, print the portion of each exon overlapping
89            # the region of interest, otherwise print the span of the region
90            
91        except:
92            continue
93
94def print_tab_sep(out_file, *args ):
95    """Print items in `l` to stdout separated by tabs"""
96    print >>out_file, string.join( [ str( f ) for f in args ], '\t' )
97
98if __name__ == "__main__": main()