comparison utils/gff_util.py @ 0:90100b587723 draft

Imported from capsule None
author devteam
date Tue, 01 Apr 2014 10:52:59 -0400
parents
children 8307665c4b6c
comparison
equal deleted inserted replaced
-1:000000000000 0:90100b587723
1 """
2 Provides utilities for working with GFF files.
3 """
4
5 import copy
6 from bx.intervals.io import *
7 from bx.tabular.io import Header, Comment
8 from utils.odict import odict
9
10 class GFFInterval( GenomicInterval ):
11 """
12 A GFF interval, including attributes. If file is strictly a GFF file,
13 only attribute is 'group.'
14 """
15 def __init__( self, reader, fields, chrom_col=0, feature_col=2, start_col=3, end_col=4, \
16 strand_col=6, score_col=5, default_strand='.', fix_strand=False ):
17 # HACK: GFF format allows '.' for strand but GenomicInterval does not. To get around this,
18 # temporarily set strand and then unset after initing GenomicInterval.
19 unknown_strand = False
20 if not fix_strand and fields[ strand_col ] == '.':
21 unknown_strand = True
22 fields[ strand_col ] = '+'
23 GenomicInterval.__init__( self, reader, fields, chrom_col, start_col, end_col, strand_col, \
24 default_strand, fix_strand=fix_strand )
25 if unknown_strand:
26 self.strand = '.'
27 self.fields[ strand_col ] = '.'
28
29 # Handle feature, score column.
30 self.feature_col = feature_col
31 if self.feature_col >= self.nfields:
32 raise MissingFieldError( "No field for feature_col (%d)" % feature_col )
33 self.feature = self.fields[ self.feature_col ]
34 self.score_col = score_col
35 if self.score_col >= self.nfields:
36 raise MissingFieldError( "No field for score_col (%d)" % score_col )
37 self.score = self.fields[ self.score_col ]
38
39 # GFF attributes.
40 self.attributes = parse_gff_attributes( fields[8] )
41
42 def copy( self ):
43 return GFFInterval(self.reader, list( self.fields ), self.chrom_col, self.feature_col, self.start_col,
44 self.end_col, self.strand_col, self.score_col, self.strand)
45
46 class GFFFeature( GFFInterval ):
47 """
48 A GFF feature, which can include multiple intervals.
49 """
50 def __init__( self, reader, chrom_col=0, feature_col=2, start_col=3, end_col=4, \
51 strand_col=6, score_col=5, default_strand='.', fix_strand=False, intervals=[], \
52 raw_size=0 ):
53 # Use copy so that first interval and feature do not share fields.
54 GFFInterval.__init__( self, reader, copy.deepcopy( intervals[0].fields ), chrom_col, feature_col, \
55 start_col, end_col, strand_col, score_col, default_strand, \
56 fix_strand=fix_strand )
57 self.intervals = intervals
58 self.raw_size = raw_size
59 # Use intervals to set feature attributes.
60 for interval in self.intervals:
61 # Error checking. NOTE: intervals need not share the same strand.
62 if interval.chrom != self.chrom:
63 raise ValueError( "interval chrom does not match self chrom: %s != %s" % \
64 ( interval.chrom, self.chrom ) )
65 # Set start, end of interval.
66 if interval.start < self.start:
67 self.start = interval.start
68 if interval.end > self.end:
69 self.end = interval.end
70
71 def name( self ):
72 """ Returns feature's name. """
73 name = None
74 # Preference for name: GTF, GFF3, GFF.
75 for attr_name in [
76 # GTF:
77 'gene_id', 'transcript_id',
78 # GFF3:
79 'ID', 'id',
80 # GFF (TODO):
81 'group' ]:
82 name = self.attributes.get( attr_name, None )
83 if name is not None:
84 break
85 return name
86
87 def copy( self ):
88 intervals_copy = []
89 for interval in self.intervals:
90 intervals_copy.append( interval.copy() )
91 return GFFFeature(self.reader, self.chrom_col, self.feature_col, self.start_col, self.end_col, self.strand_col,
92 self.score_col, self.strand, intervals=intervals_copy )
93
94 def lines( self ):
95 lines = []
96 for interval in self.intervals:
97 lines.append( '\t'.join( interval.fields ) )
98 return lines
99
100
101 class GFFIntervalToBEDReaderWrapper( NiceReaderWrapper ):
102 """
103 Reader wrapper that reads GFF intervals/lines and automatically converts
104 them to BED format.
105 """
106
107 def parse_row( self, line ):
108 # HACK: this should return a GFF interval, but bx-python operations
109 # require GenomicInterval objects and subclasses will not work.
110 interval = GenomicInterval( self, line.split( "\t" ), self.chrom_col, self.start_col, \
111 self.end_col, self.strand_col, self.default_strand, \
112 fix_strand=self.fix_strand )
113 interval = convert_gff_coords_to_bed( interval )
114 return interval
115
116 class GFFReaderWrapper( NiceReaderWrapper ):
117 """
118 Reader wrapper for GFF files.
119
120 Wrapper has two major functions:
121
122 1. group entries for GFF file (via group column), GFF3 (via id attribute),
123 or GTF (via gene_id/transcript id);
124 2. convert coordinates from GFF format--starting and ending coordinates
125 are 1-based, closed--to the 'traditional'/BED interval format--0 based,
126 half-open. This is useful when using GFF files as inputs to tools that
127 expect traditional interval format.
128 """
129
130 def __init__( self, reader, chrom_col=0, feature_col=2, start_col=3, \
131 end_col=4, strand_col=6, score_col=5, fix_strand=False, convert_to_bed_coord=False, **kwargs ):
132 NiceReaderWrapper.__init__( self, reader, chrom_col=chrom_col, start_col=start_col, end_col=end_col, \
133 strand_col=strand_col, fix_strand=fix_strand, **kwargs )
134 self.feature_col = feature_col
135 self.score_col = score_col
136 self.convert_to_bed_coord = convert_to_bed_coord
137 self.last_line = None
138 self.cur_offset = 0
139 self.seed_interval = None
140 self.seed_interval_line_len = 0
141
142 def parse_row( self, line ):
143 interval = GFFInterval( self, line.split( "\t" ), self.chrom_col, self.feature_col, \
144 self.start_col, self.end_col, self.strand_col, self.score_col, \
145 self.default_strand, fix_strand=self.fix_strand )
146 return interval
147
148 def next( self ):
149 """ Returns next GFFFeature. """
150
151 #
152 # Helper function.
153 #
154
155 def handle_parse_error( parse_error ):
156 """ Actions to take when ParseError found. """
157 if self.outstream:
158 if self.print_delegate and hasattr(self.print_delegate,"__call__"):
159 self.print_delegate( self.outstream, e, self )
160 self.skipped += 1
161 # no reason to stuff an entire bad file into memmory
162 if self.skipped < 10:
163 self.skipped_lines.append( ( self.linenum, self.current_line, str( e ) ) )
164
165 # For debugging, uncomment this to propogate parsing exceptions up.
166 # I.e. the underlying reason for an unexpected StopIteration exception
167 # can be found by uncommenting this.
168 # raise e
169
170 #
171 # Get next GFFFeature
172 #
173 raw_size = self.seed_interval_line_len
174
175 # If there is no seed interval, set one. Also, if there are no more
176 # intervals to read, this is where iterator dies.
177 if not self.seed_interval:
178 while not self.seed_interval:
179 try:
180 self.seed_interval = GenomicIntervalReader.next( self )
181 except ParseError, e:
182 handle_parse_error( e )
183 # TODO: When no longer supporting python 2.4 use finally:
184 #finally:
185 raw_size += len( self.current_line )
186
187 # If header or comment, clear seed interval and return it with its size.
188 if isinstance( self.seed_interval, ( Header, Comment ) ):
189 return_val = self.seed_interval
190 return_val.raw_size = len( self.current_line )
191 self.seed_interval = None
192 self.seed_interval_line_len = 0
193 return return_val
194
195 # Initialize feature identifier from seed.
196 feature_group = self.seed_interval.attributes.get( 'group', None ) # For GFF
197 # For GFF3
198 feature_id = self.seed_interval.attributes.get( 'ID', None )
199 feature_parent_id = self.seed_interval.attributes.get( 'Parent', None )
200 # For GTF.
201 feature_gene_id = self.seed_interval.attributes.get( 'gene_id', None )
202 feature_transcript_id = self.seed_interval.attributes.get( 'transcript_id', None )
203
204 # Read all intervals associated with seed.
205 feature_intervals = []
206 feature_intervals.append( self.seed_interval )
207 while True:
208 try:
209 interval = GenomicIntervalReader.next( self )
210 raw_size += len( self.current_line )
211 except StopIteration, e:
212 # No more intervals to read, but last feature needs to be
213 # returned.
214 interval = None
215 raw_size += len( self.current_line )
216 break
217 except ParseError, e:
218 handle_parse_error( e )
219 raw_size += len( self.current_line )
220 continue
221 # TODO: When no longer supporting python 2.4 use finally:
222 #finally:
223 #raw_size += len( self.current_line )
224
225 # Ignore comments.
226 if isinstance( interval, Comment ):
227 continue
228
229 # Determine if interval is part of feature.
230 part_of = False
231 group = interval.attributes.get( 'group', None )
232 # GFF test:
233 if group and feature_group == group:
234 part_of = True
235 # GFF3 test:
236 parent_id = interval.attributes.get( 'Parent', None )
237 cur_id = interval.attributes.get( 'ID', None )
238 if ( cur_id and cur_id == feature_id ) or ( parent_id and parent_id == feature_id ):
239 part_of = True
240 # GTF test:
241 transcript_id = interval.attributes.get( 'transcript_id', None )
242 if transcript_id and transcript_id == feature_transcript_id:
243 part_of = True
244
245 # If interval is not part of feature, clean up and break.
246 if not part_of:
247 # Adjust raw size because current line is not part of feature.
248 raw_size -= len( self.current_line )
249 break
250
251 # Interval associated with feature.
252 feature_intervals.append( interval )
253
254 # Last interval read is the seed for the next interval.
255 self.seed_interval = interval
256 self.seed_interval_line_len = len( self.current_line )
257
258 # Return feature.
259 feature = GFFFeature( self, self.chrom_col, self.feature_col, self.start_col, \
260 self.end_col, self.strand_col, self.score_col, \
261 self.default_strand, fix_strand=self.fix_strand, \
262 intervals=feature_intervals, raw_size=raw_size )
263
264 # Convert to BED coords?
265 if self.convert_to_bed_coord:
266 convert_gff_coords_to_bed( feature )
267
268 return feature
269
270 def convert_bed_coords_to_gff( interval ):
271 """
272 Converts an interval object's coordinates from BED format to GFF format.
273 Accepted object types include GenomicInterval and list (where the first
274 element in the list is the interval's start, and the second element is
275 the interval's end).
276 """
277 if isinstance( interval, GenomicInterval ):
278 interval.start += 1
279 if isinstance( interval, GFFFeature ):
280 for subinterval in interval.intervals:
281 convert_bed_coords_to_gff( subinterval )
282 elif type ( interval ) is list:
283 interval[ 0 ] += 1
284 return interval
285
286 def convert_gff_coords_to_bed( interval ):
287 """
288 Converts an interval object's coordinates from GFF format to BED format.
289 Accepted object types include GFFFeature, GenomicInterval, and list (where
290 the first element in the list is the interval's start, and the second
291 element is the interval's end).
292 """
293 if isinstance( interval, GenomicInterval ):
294 interval.start -= 1
295 if isinstance( interval, GFFFeature ):
296 for subinterval in interval.intervals:
297 convert_gff_coords_to_bed( subinterval )
298 elif type ( interval ) is list:
299 interval[ 0 ] -= 1
300 return interval
301
302 def parse_gff_attributes( attr_str ):
303 """
304 Parses a GFF/GTF attribute string and returns a dictionary of name-value
305 pairs. The general format for a GFF3 attributes string is
306
307 name1=value1;name2=value2
308
309 The general format for a GTF attribute string is
310
311 name1 "value1" ; name2 "value2"
312
313 The general format for a GFF attribute string is a single string that
314 denotes the interval's group; in this case, method returns a dictionary
315 with a single key-value pair, and key name is 'group'
316 """
317 attributes_list = attr_str.split(";")
318 attributes = {}
319 for name_value_pair in attributes_list:
320 # Try splitting by '=' (GFF3) first because spaces are allowed in GFF3
321 # attribute; next, try double quotes for GTF.
322 pair = name_value_pair.strip().split("=")
323 if len( pair ) == 1:
324 pair = name_value_pair.strip().split("\"")
325 if len( pair ) == 1:
326 # Could not split for some reason -- raise exception?
327 continue
328 if pair == '':
329 continue
330 name = pair[0].strip()
331 if name == '':
332 continue
333 # Need to strip double quote from values
334 value = pair[1].strip(" \"")
335 attributes[ name ] = value
336
337 if len( attributes ) == 0:
338 # Could not split attributes string, so entire string must be
339 # 'group' attribute. This is the case for strictly GFF files.
340 attributes['group'] = attr_str
341 return attributes
342
343 def gff_attributes_to_str( attrs, gff_format ):
344 """
345 Convert GFF attributes to string. Supported formats are GFF3, GTF.
346 """
347 if gff_format == 'GTF':
348 format_string = '%s "%s"'
349 # Convert group (GFF) and ID, parent (GFF3) attributes to transcript_id, gene_id
350 id_attr = None
351 if 'group' in attrs:
352 id_attr = 'group'
353 elif 'ID' in attrs:
354 id_attr = 'ID'
355 elif 'Parent' in attrs:
356 id_attr = 'Parent'
357 if id_attr:
358 attrs['transcript_id'] = attrs['gene_id'] = attrs[id_attr]
359 elif gff_format == 'GFF3':
360 format_string = '%s=%s'
361 attrs_strs = []
362 for name, value in attrs.items():
363 attrs_strs.append( format_string % ( name, value ) )
364 return " ; ".join( attrs_strs )
365
366 def read_unordered_gtf( iterator, strict=False ):
367 """
368 Returns GTF features found in an iterator. GTF lines need not be ordered
369 or clustered for reader to work. Reader returns GFFFeature objects sorted
370 by transcript_id, chrom, and start position.
371 """
372
373 # -- Get function that generates line/feature key. --
374
375 get_transcript_id = lambda fields: parse_gff_attributes( fields[8] )[ 'transcript_id' ]
376 if strict:
377 # Strict GTF parsing uses transcript_id only to group lines into feature.
378 key_fn = get_transcript_id
379 else:
380 # Use lenient parsing where chromosome + transcript_id is the key. This allows
381 # transcripts with same ID on different chromosomes; this occurs in some popular
382 # datasources, such as RefGenes in UCSC.
383 key_fn = lambda fields: fields[0] + '_' + get_transcript_id( fields )
384
385
386 # Aggregate intervals by transcript_id and collect comments.
387 feature_intervals = odict()
388 comments = []
389 for count, line in enumerate( iterator ):
390 if line.startswith( '#' ):
391 comments.append( Comment( line ) )
392 continue
393
394 line_key = key_fn( line.split('\t') )
395 if line_key in feature_intervals:
396 feature = feature_intervals[ line_key ]
397 else:
398 feature = []
399 feature_intervals[ line_key ] = feature
400 feature.append( GFFInterval( None, line.split( '\t' ) ) )
401
402 # Create features.
403 chroms_features = {}
404 for count, intervals in enumerate( feature_intervals.values() ):
405 # Sort intervals by start position.
406 intervals.sort( lambda a,b: cmp( a.start, b.start ) )
407 feature = GFFFeature( None, intervals=intervals )
408 if feature.chrom not in chroms_features:
409 chroms_features[ feature.chrom ] = []
410 chroms_features[ feature.chrom ].append( feature )
411
412 # Sort features by chrom, start position.
413 chroms_features_sorted = []
414 for chrom_features in chroms_features.values():
415 chroms_features_sorted.append( chrom_features )
416 chroms_features_sorted.sort( lambda a,b: cmp( a[0].chrom, b[0].chrom ) )
417 for features in chroms_features_sorted:
418 features.sort( lambda a,b: cmp( a.start, b.start ) )
419
420 # Yield comments first, then features.
421 # FIXME: comments can appear anywhere in file, not just the beginning.
422 # Ideally, then comments would be associated with features and output
423 # just before feature/line.
424 for comment in comments:
425 yield comment
426
427 for chrom_features in chroms_features_sorted:
428 for feature in chrom_features:
429 yield feature
430