1 | #!/usr/bin/env python |
---|
2 | |
---|
3 | """ |
---|
4 | Read a table dump in the UCSC gene table format and print a tab separated |
---|
5 | list of intervals corresponding to requested features of each gene. |
---|
6 | |
---|
7 | usage: ucsc_gene_table_to_intervals.py [options] |
---|
8 | |
---|
9 | options: |
---|
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 | |
---|
18 | import optparse, string, sys |
---|
19 | |
---|
20 | assert sys.version_info[:2] >= ( 2, 4 ) |
---|
21 | |
---|
22 | def 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 | for line in in_file: |
---|
60 | try: |
---|
61 | if line[0:1] == "#": |
---|
62 | continue |
---|
63 | # Parse fields from gene tabls |
---|
64 | fields = line.split( '\t' ) |
---|
65 | name = fields[0] |
---|
66 | chrom = fields[1] |
---|
67 | strand = fields[2].replace(" ","_") |
---|
68 | tx_start = int( fields[3] ) |
---|
69 | tx_end = int( fields[4] ) |
---|
70 | cds_start = int( fields[5] ) |
---|
71 | cds_end = int( fields[6] ) |
---|
72 | |
---|
73 | # Determine the subset of the transcribed region we are interested in |
---|
74 | if options.region == 'utr3': |
---|
75 | if strand == '-': region_start, region_end = tx_start, cds_start |
---|
76 | else: region_start, region_end = cds_end, tx_end |
---|
77 | elif options.region == 'utr5': |
---|
78 | if strand == '-': region_start, region_end = cds_end, tx_end |
---|
79 | else: region_start, region_end = tx_start, cds_start |
---|
80 | elif options.region == 'coding': |
---|
81 | region_start, region_end = cds_start, cds_end |
---|
82 | else: |
---|
83 | region_start, region_end = tx_start, tx_end |
---|
84 | |
---|
85 | # If only interested in exons, print the portion of each exon overlapping |
---|
86 | # the region of interest, otherwise print the span of the region |
---|
87 | if options.exons: |
---|
88 | exon_starts = map( int, fields[8].rstrip( ',\n' ).split( ',' ) ) |
---|
89 | exon_ends = map( int, fields[9].rstrip( ',\n' ).split( ',' ) ) |
---|
90 | for start, end in zip( exon_starts, exon_ends ): |
---|
91 | start = max( start, region_start ) |
---|
92 | end = min( end, region_end ) |
---|
93 | if start < end: |
---|
94 | if strand: print_tab_sep(out_file, chrom, start, end, name, "0", strand ) |
---|
95 | else: print_tab_sep(out_file, chrom, start, end ) |
---|
96 | else: |
---|
97 | if strand: print_tab_sep(out_file, chrom, region_start, region_end, name, "0", strand ) |
---|
98 | else: print_tab_sep(out_file, chrom, region_start, region_end ) |
---|
99 | except: |
---|
100 | continue |
---|
101 | |
---|
102 | def print_tab_sep(out_file, *args ): |
---|
103 | """Print items in `l` to stdout separated by tabs""" |
---|
104 | print >>out_file, string.join( [ str( f ) for f in args ], '\t' ) |
---|
105 | |
---|
106 | if __name__ == "__main__": main() |
---|