[Bio] / FigKernelPackages / gjogff.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/gjogff.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download) (as text)

1 : golsen 1.1 package gjogff;
2 :    
3 :     use strict;
4 :     use gjoseqlib;
5 :     use Data::Dumper;
6 :    
7 :     #-------------------------------------------------------------------------------
8 :     # Parse the descriptions of the features into key-value pairs
9 :     #
10 : golsen 1.2 # ( \%features_by_contig, \@dna ) = read_gff( \*FH, \%options );
11 :     #
12 :     # Options:
13 :     #
14 :     # contig_ids => \%ids
15 :     # contig_ids => \@ids
16 :     # dna => \@dna
17 : golsen 1.1 #
18 :     # Features are:
19 :     #
20 :     # [ $type, $loc, $id, $name, \%keys_values ]
21 :     #
22 : golsen 1.2 # Location is Sapling style contig_begĀ±length
23 : golsen 1.1 #
24 :     #-------------------------------------------------------------------------------
25 :     sub read_gff
26 :     {
27 : golsen 1.2 my $opts = @_ && $_[ 0] && ref $_[ 0] eq 'HASH' ? shift :
28 :     @_ && $_[-1] && ref $_[-1] eq 'HASH' ? pop :
29 :     {};
30 :    
31 :     my $contig_ids = $opts->{ contig_ids } && ref $opts->{ contig_ids } eq 'HASH' ? $opts->{ contig_ids } :
32 :     $opts->{ contig_ids } && ref $opts->{ contig_ids } eq 'ARRAY' ? { map { $_ => 1 } @{ $opts->{ contig_ids } } } :
33 :     $opts->{ contig_len } && ref $opts->{ contig_len } eq 'HASH' ? $opts->{ contig_len } :
34 :     $opts->{ contig_len } && ref $opts->{ contig_len } eq 'ARRAY' ? { map { $_ => 1 } @{ $opts->{ contig_len } } } :
35 :     $opts->{ dna } && ref $opts->{ dna } eq 'ARRAY' ? { map { $_->[0] => length $_->[2] } @{ $opts->{ dna } } } :
36 :     undef;
37 :    
38 : golsen 1.1 my ( $fh, $close ) = input_filehandle( $_[0] );
39 :    
40 :     my $dna;
41 :     my %ftrs;
42 :     local $_;
43 :    
44 : golsen 1.2 my %no_contig;
45 : golsen 1.1 while ( defined( $_ = <$fh> ) )
46 :     {
47 :     if ( m/^##FASTA/i ) { $dna = 1; last } # GFF DNA introduction
48 :     next if m/^track/i; # JGI header
49 :     next if m/^#/; # GFF comment
50 :    
51 :     chomp;
52 :     my ( $contig, $auth, $type, $l, $r, $scr, $dir, $fr, $rest ) = split /\t/;
53 :     my ( $id, $name, $key_value ) = parse_description( $rest );
54 :     push @{ $ftrs{ $id }->{ $type } },
55 :     [ $id, $name, $contig, $l, $r, $dir, $fr, $key_value, $l+$r ];
56 :     }
57 :    
58 :     my @dna;
59 :     @dna = gjoseqlib::read_fasta( $fh ) if $dna;
60 : golsen 1.2 $contig_ids ||= { map { $_->[0] => 1 } @dna } if @dna;
61 : golsen 1.1
62 :     close( $fh ) if $close;
63 :    
64 :     # Merge segments of same feature. Resulting descriptions are:
65 :     #
66 :     # [ $type, $loc, $id, $name, $frame, \%keys_values ]
67 :     #
68 :    
69 : golsen 1.2 my $features_by_contig = process_gff_features( \%ftrs, $contig_ids );
70 : golsen 1.1
71 :     wantarray ? ( $features_by_contig, \@dna ) : $features_by_contig;
72 :     }
73 :    
74 :    
75 :     #-------------------------------------------------------------------------------
76 :     # Parse the descriptions of the features into key-value pairs
77 :     #
78 :     # ( $id, $name, \%key_values ) = parse_description( $gff_last_field );
79 :     #
80 :     #-------------------------------------------------------------------------------
81 :     sub parse_description
82 :     {
83 :     # Check for %xx escaped text:
84 :    
85 :     my @candidates = $_[0] =~ m/(\%..)/gi;
86 :     my $escaped = @candidates && ( @candidates == grep { m/\%[0-9a-f][0-9a-f]/i } @candidates );
87 :    
88 :     # Split into key=value pairs, remove surounding quotes, and unescape %..
89 :     # Some keys can have multiple values, so keep lists.
90 :    
91 :     my %key_values;
92 :     foreach ( split / *; */, $_[0] )
93 :     {
94 :     my ( $key, $value ) = /^(\S+)[ =](.+)$/ ? ( $1 => $2 ) : ( 'name' => $_ );
95 :     $value =~ s/""/"/g if $value =~ s/^"(.*)"$/$1/;
96 :     $value =~ s/''/'/g if $value =~ s/^'(.*)'$/$1/;
97 :     $value = unescape( $value ) if $escaped;
98 :     push @{$key_values{$key}}, $value;
99 :     }
100 :    
101 :     # Look for an ID and name
102 :    
103 :     my $key;
104 :    
105 :     ( $key ) = grep { lc $_ eq 'name' } keys %key_values;
106 :     ( $key ) = grep { m/name/i } keys %key_values if ! $key;
107 :     ( $key ) = grep { /I[dD]$/ } keys %key_values if ! $key;
108 :     ( $key ) = grep { /id$/ } keys %key_values if ! $key;
109 :     my $name = $key ? $key_values{ $key }->[0] : $_[0];
110 :    
111 :     ( $key ) = grep { m/^ID$/i } keys %key_values;
112 :     ( $key ) = grep { /I[dD]$/ } keys %key_values if ! $key;
113 :     ( $key ) = grep { /id$/ } keys %key_values if ! $key;
114 :     ( $key ) = grep { lc eq 'parent' } keys %key_values if ! $key;
115 :     my $id = $key ? $key_values{ $key }->[0] : $name;
116 :    
117 :     ( $id, $name, \%key_values );
118 :     }
119 :    
120 :    
121 :     #-------------------------------------------------------------------------------
122 :     # Integrate and sort feature information:
123 :     #
124 : golsen 1.2 # \%features_by_contig = process_gff_features( \%features, \%contig_ids );
125 : golsen 1.1 #
126 :     # Features are:
127 :     #
128 :     # [ $type, $loc, $id, $name, \%keys_values ]
129 :     #
130 : golsen 1.2 # Location is Sapling style contig_begĀ±length
131 : golsen 1.1 #
132 :     #-------------------------------------------------------------------------------
133 :     sub process_gff_features
134 :     {
135 : golsen 1.2 my $ftrs = shift || {};
136 :     my $contig_ids = $_[0] && ref $_[0] eq 'HASH' ? shift : {};
137 :     my $have_ids = keys %$contig_ids;
138 :     my $have_len = ( values %$contig_ids )[0] > 1;
139 :    
140 : golsen 1.1 my @ftrs;
141 : golsen 1.2 my %missing_contig; # Not currently reporting
142 : golsen 1.1 foreach my $id ( keys %$ftrs )
143 :     {
144 :     my $types = $ftrs->{ $id };
145 :     foreach my $type ( keys %$types )
146 :     {
147 :     my $okay = 1;
148 :     my $data = $types->{ $type };
149 :    
150 :     # 0 1 2 3 4 5 6 7 8
151 :     # $datum = [ $id, $name, $contig, $left, $right, $dir, $frame, \%keys_values, $mid ]
152 :    
153 : golsen 1.2 #
154 :     # Contig anlysis:
155 :     #
156 :     # Number of contigs contig_data ids_matching response
157 :     # -----------------------------------------------------------------
158 :     # 1 yes 0 contig missing
159 :     # 1 yes 1 okay
160 :     # 1 no okay
161 :     # >1 yes 0 contig missing
162 :     # >1 yes 1 keep the 1
163 :     # >1 yes >1 multiple contigs
164 :     # >1 no multiple contigs
165 :     # -----------------------------------------------------------------
166 :     #
167 :    
168 :     my %contigs = map { $_->[2] => 1 } @$data;
169 :     my @contigs = keys %contigs;
170 :     my @have_contig = grep { $contig_ids->{ $_ } } @contigs;
171 :     if ( @contigs == 1 )
172 :     {
173 :     $missing_contig{ $contigs[0] } = 1 if $have_ids && ! @have_contig;
174 :     }
175 :     elsif ( @have_contig == 1 )
176 :     {
177 :     # Filter raw data to the contig we have:
178 :     @$data = grep { $_->[2] eq $have_contig[0] } @$data;
179 :     }
180 :     elsif ( $have_ids && ! @have_contig )
181 :     {
182 :     $missing_contig{ $contigs[0] } = 1;
183 :     }
184 :     else
185 :     {
186 :     print STDERR "Feature '$id' of type '$type' has multiple contigs: ",
187 :     join( ', ', sort @contigs ), "\n";
188 :     $okay = 0;
189 :     }
190 :    
191 : golsen 1.1 my %names = map { $_->[1] => 1 } @$data;
192 :     if ( keys %names > 1 )
193 :     {
194 :     print STDERR "Feature '$id' of type '$type' has multiple names: ",
195 :     join( ', ', sort keys %names ), "\n";
196 :     $okay = 0;
197 :     }
198 :    
199 :     my %dirs = map { $_->[5] => 1 } @$data;
200 :     if ( keys %dirs > 1 )
201 :     {
202 :     print STDERR "Feature '$id' of type '$type' has multiple directions.\n";
203 :     $okay = 0;
204 :     }
205 :    
206 :     next if ! $okay;
207 :    
208 :     # Grab shared information from the first exon:
209 :    
210 :     my ( $name, $cont, $dir, $frame ) = @{ $data->[0] }[ 1, 2, 5, 6 ];
211 :    
212 :     # Sort my midpoint
213 :    
214 :     my $ndir = $dir =~ /^-/ ? -1 : +1;
215 :     my @segs = sort { $ndir * ( $a->[8] <=> $b->[8] ) } @$data;
216 :    
217 :     # Location in Sapling format
218 :    
219 :     my @locs = map { [ $cont, # contig
220 :     ( $ndir > 0 ? $_->[3] : $_->[4] ), # start
221 :     $_->[5], # dir
222 :     $_->[4] - $_->[3] + 1 # length
223 :     ] }
224 :     @segs;
225 :     my $loc = join( ',', map { "$_->[0]_$_->[1]$_->[2]$_->[3]" } @locs );
226 :    
227 :     my $mid = 0.5 * ( $segs[0]->[3] + $segs[-1]->[4] );
228 :    
229 :     # Merge the keys and values for the merged gff lines:
230 :    
231 :     my %keys_values;
232 :     my %n_seen;
233 :     foreach my $keys_values0 ( map { $_->[7] } @$data )
234 :     {
235 :     foreach my $key ( keys %$keys_values0 )
236 :     {
237 :     my $values = $keys_values0->{ $key };
238 :     foreach ( @$values )
239 :     {
240 :     push @{ $keys_values{ $key } }, $_ if ! $n_seen{ $key }->{ $_ }++;
241 :     }
242 :     }
243 :     }
244 :    
245 : golsen 1.2 if ( $have_len && $contig_ids->{ $cont } )
246 :     {
247 :     if ( max( map { $_->[4] } @$data ) > $contig_ids->{ $cont } )
248 :     {
249 :     print STDERR "Feature $type: $id $name extends beyond end of the contig DNA.\n";
250 :     next;
251 :     }
252 :     }
253 :    
254 :     # The JGI files list frame 0 in all forward genes, so kill this.
255 : golsen 1.1 # $keys_values{ codon_start } = [ $frame =~ /\d/ ? $frame+1 : 1 ];
256 :    
257 :     push @ftrs, [ [ $type, $loc, $id, $name, \%keys_values ], $cont, $mid ];
258 :     }
259 :     }
260 :    
261 :     my %features_by_contig;
262 :     foreach ( sort { $a->[1] cmp $b->[1] || $a->[2] <=> $b->[2] } @ftrs )
263 :     {
264 :     push @{ $features_by_contig{ $_->[1] } }, $_->[0];
265 :     }
266 :    
267 : golsen 1.2 if ( $have_ids && keys %missing_contig )
268 :     {
269 :     print STDERR "No data for contigs:\n";
270 :     print STDERR map { " $_\n" } sort keys %missing_contig;
271 :     }
272 :    
273 : golsen 1.1 \%features_by_contig;
274 :     }
275 :    
276 :    
277 : golsen 1.2 sub max { my $max = shift; foreach ( @_ ) { $max = $_ if $_ > $max } $max }
278 :    
279 :    
280 : golsen 1.1 #-------------------------------------------------------------------------------
281 :     # Deal with %.. escaped text
282 :     #
283 :     # $text = unescape( $text );
284 :     #
285 :     #-------------------------------------------------------------------------------
286 :     sub unescape
287 :     {
288 :     join '', map { /\%(..)/ ? chr( hex( $1 ) ) : $_ } # convert
289 :     $_[0] =~ m/(\%[0-9a-f][0-9a-f]|.)/gi; # split
290 :     }
291 :    
292 :    
293 :     #-------------------------------------------------------------------------------
294 :     # Helper function for defining an input filehandle:
295 :     #
296 :     # filehandle is passed through
297 :     # string is taken as file name to be openend
298 :     # undef or "" defaults to STDOUT
299 :     #
300 :     # \*FH = input_filehandle( $file );
301 :     # ( \*FH, $close ) = input_filehandle( $file );
302 :     #
303 :     #-------------------------------------------------------------------------------
304 :     sub input_filehandle
305 :     {
306 :     my $file = shift;
307 :    
308 :     # Null string or undef
309 :    
310 :     if ( ! defined( $file ) || ( $file eq '' ) )
311 :     {
312 :     return wantarray ? ( \*STDIN, 0 ) : \*STDIN;
313 :     }
314 :    
315 :     # FILEHANDLE?
316 :    
317 :     if ( ref( $file ) eq "GLOB" )
318 :     {
319 :     return wantarray ? ( $file, 0 ) : $file;
320 :     }
321 :    
322 :     # File name
323 :    
324 :     if ( ! ref( $file ) )
325 :     {
326 :     -f $file or die "Could not find input file \"$file\"\n";
327 :     my $fh;
328 :     open( $fh, "<$file" ) || die "Could not open \"$file\" for input\n";
329 :     return wantarray ? ( $fh, 1 ) : $fh;
330 :     }
331 :    
332 :     return wantarray ? ( \*STDIN, undef ) : \*STDIN;
333 :     }
334 :    
335 :    
336 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3