Mercurial > repos > xuebing > sharplabtool
comparison tools/filters/ucsc_gene_table_to_intervals.py @ 0:9071e359b9a3
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:37:19 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:9071e359b9a3 |
---|---|
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() |