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

Annotation of /FigKernelPackages/gjogff.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3