0
|
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, codon, intron, 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', 'intron', 'codon' ), "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
|
|
59 # Read table and handle each gene
|
|
60 for line in in_file:
|
|
61 try:
|
|
62 if line[0:1] == "#":
|
|
63 continue
|
|
64 # Parse fields from gene tabls
|
|
65 fields = line.split( '\t' )
|
|
66 chrom = fields[0]
|
|
67 tx_start = int( fields[1] )
|
|
68 tx_end = int( fields[2] )
|
|
69 name = fields[3]
|
|
70 strand = fields[5].replace(" ","_")
|
|
71 cds_start = int( fields[6] )
|
|
72 cds_end = int( fields[7] )
|
|
73
|
|
74 # Determine the subset of the transcribed region we are interested in
|
|
75 if options.region == 'utr3':
|
|
76 if strand == '-': region_start, region_end = tx_start, cds_start
|
|
77 else: region_start, region_end = cds_end, tx_end
|
|
78 elif options.region == 'utr5':
|
|
79 if strand == '-': region_start, region_end = cds_end, tx_end
|
|
80 else: region_start, region_end = tx_start, cds_start
|
|
81 elif options.region == 'coding' or options.region == 'codon':
|
|
82 region_start, region_end = cds_start, cds_end
|
|
83 else:
|
|
84 region_start, region_end = tx_start, tx_end
|
|
85
|
|
86 # If only interested in exons, print the portion of each exon overlapping
|
|
87 # the region of interest, otherwise print the span of the region
|
|
88 # options.exons is always TRUE
|
|
89 if options.exons:
|
|
90 exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
|
|
91 exon_starts = map((lambda x: x + tx_start ), exon_starts)
|
|
92 exon_ends = map( int, fields[10].rstrip( ',\n' ).split( ',' ) )
|
|
93 exon_ends = map((lambda x, y: x + y ), exon_starts, exon_ends);
|
|
94
|
|
95 #for Intron regions:
|
|
96 if options.region == 'intron':
|
|
97 i=0
|
|
98 while i < len(exon_starts)-1:
|
|
99 intron_starts = exon_ends[i]
|
|
100 intron_ends = exon_starts[i+1]
|
|
101 if strand: print_tab_sep(out_file, chrom, intron_starts, intron_ends, name, "0", strand )
|
|
102 else: print_tab_sep(out_file, chrom, intron_starts, intron_ends )
|
|
103 i+=1
|
|
104 #for non-intron regions:
|
|
105 else:
|
|
106 for start, end in zip( exon_starts, exon_ends ):
|
|
107 start = max( start, region_start )
|
|
108 end = min( end, region_end )
|
|
109 if start < end:
|
|
110 if options.region == 'codon':
|
|
111 start += (3 - ((start-region_start)%3))%3
|
|
112 c_start = start
|
|
113 while c_start+3 <= end:
|
|
114 if strand:
|
|
115 print_tab_sep(out_file, chrom, c_start, c_start+3, name, "0", strand )
|
|
116 else:
|
|
117 print_tab_sep(out_file, chrom, c_start, c_start+3)
|
|
118 c_start += 3
|
|
119 else:
|
|
120 if strand:
|
|
121 print_tab_sep(out_file, chrom, start, end, name, "0", strand )
|
|
122 else:
|
|
123 print_tab_sep(out_file, chrom, start, end )
|
|
124 """
|
|
125 else:
|
|
126 if options.region == 'codon':
|
|
127 c_start = start
|
|
128 c_end = end
|
|
129 if c_start > c_end:
|
|
130 t = c_start
|
|
131 c_start = c_end
|
|
132 c_end = t
|
|
133 while c_start+3 <= c_end:
|
|
134 if strand:
|
|
135 print_tab_sep(out_file, chrom, c_start, c_start+3, name, "0", strand )
|
|
136 else:
|
|
137 print_tab_sep(out_file, chrom, c_start, c_start+3)
|
|
138 c_start += 3
|
|
139 else:
|
|
140 if strand:
|
|
141 print_tab_sep(out_file, chrom, region_start, region_end, name, "0", strand )
|
|
142 else:
|
|
143 print_tab_sep(out_file, chrom, region_start, region_end )
|
|
144 """
|
|
145 except:
|
|
146 continue
|
|
147
|
|
148 def print_tab_sep(out_file, *args ):
|
|
149 """Print items in `l` to stdout separated by tabs"""
|
|
150 print >>out_file, string.join( [ str( f ) for f in args ], '\t' )
|
|
151
|
|
152 if __name__ == "__main__": main()
|