Mercurial > repos > xuebing > sharplabtool

comparison tools/filters/ucsc_gene_bed_to_exon_bed.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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, 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()