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

Annotation of /FigKernelPackages/MapIDs.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : wilke 1.1 #
2 :     # Package to retrieve accession numbers
3 :     #
4 :    
5 :    
6 :     package MapIDs;
7 :    
8 :     use Carp;
9 :     use Data::Dumper;
10 :     use strict;
11 :     use warnings;
12 :    
13 :     use FIG;
14 :    
15 :    
16 :    
17 :    
18 :     =head3 get_gi
19 :    
20 :     get a gi number for a given peg id
21 :     returns -1 or -2,-3,-4 if no gi number is found or results are amibgious
22 :    
23 :     =cut
24 :    
25 :    
26 :     sub get_gi_extended{
27 :     my ($fig,$peg,$org) = @_;
28 :    
29 :     my $error = 0;
30 :     my $error_code = 0;
31 :     my $msg = "";
32 :     my $gi_list = [];
33 :     my @return_values = ( $error , $gi_list, $msg , $error_code );
34 :    
35 :     # initialize gi number, returns
36 :     # -1 no gi number found
37 :     # -2 too many matched gis
38 :     # -3 no gi in corresponding organism
39 :     # -4 gi found but length is unequal
40 :    
41 :    
42 :     # get all synonyms for a given id
43 :     my $pid = $fig->maps_to_id($peg);
44 :     my @mapped = $fig->mapped_prot_ids($pid);
45 :    
46 :     # get organism for peg
47 :     my ($genome_id, $peg_number) = $fig->genome_and_peg_of($peg);
48 :    
49 :    
50 :     # get all gi numbers from synonym list and map to organism
51 :     my @gis; # list of all gi synonyms for a given peg
52 :     my @pegs; # list of all peg synonyms for a given peg
53 :     my $synonyms = {};
54 :     foreach my $map ( @mapped ){
55 :     push ( @gis, $map->[0] ) if ( $map->[0] =~ m/^gi\|/ );
56 :     push ( @pegs, $map->[0] ) if ( $map->[0] =~ m/^fig\|$genome_id/ );
57 :     $synonyms->{ $map->[0] } = $map->[1];
58 :    
59 :     }
60 :    
61 :    
62 :     if ( scalar(@gis) > 0){
63 :    
64 :    
65 :     # sort gis by organism
66 :     my $org_name = {};
67 :     foreach my $gi (@gis){
68 :    
69 :     if ( $org_name->{ $fig->org_of($gi) } or $org_name->{ unknown }){
70 :    
71 :     if ( $fig->org_of($gi) ){
72 :     push @{$org_name->{ $fig->org_of($gi) }} ,$gi;
73 :     }
74 :     else{
75 :     push @{$org_name->{ unknown }} ,$gi;
76 :     }
77 :    
78 :    
79 :     }
80 :     else{
81 :    
82 :     # in different genomes
83 :     if ( $fig->org_of($gi) ) {
84 :     $org_name->{ $fig->org_of($gi) } = [$gi];
85 :     }
86 :     else{
87 :     $org_name->{ unknown } = [$gi];
88 :     }
89 :     }
90 :     }
91 :    
92 :    
93 :     # check genome name
94 :     # hope that fig genome name match genome name from refseq/genbank
95 :    
96 :     my $hits = [ ] ;
97 :     $hits = $org_name->{ $fig->genus_species( $genome_id) } if ( defined $fig->genus_species( $genome_id) );
98 :    
99 :     if ( ref $hits and scalar @$hits ){
100 :    
101 :     # check for hit counts
102 :     # resolve if more than one hit per organism
103 :    
104 :     #print scalar @$hits," hit(s): ", @$hits, "\n" ;
105 :    
106 :     if ( scalar @$hits > 1 ){
107 :    
108 :     # paralog within the same genome
109 :     # do some magic stuff to resolve
110 :    
111 :     # 1. compare by length
112 :     # 2. compare by sequence
113 :     # 3. compare by start and stop / position in genome
114 :    
115 :    
116 :    
117 :     my $pid = $fig->maps_to_id($peg);
118 :     my @mapped = $fig->mapped_prot_ids($pid);
119 :    
120 :     my ($genome_id, $peg_number) = $fig->genome_and_peg_of($peg);
121 :    
122 :    
123 :     my @gis; # list of all gi synonyms for a given peg
124 :     my @pegs; # list of all peg synonyms for a given peg
125 :    
126 :    
127 :    
128 :     # compare length of peg with length of hits
129 :     my $same_length = [];
130 :     foreach my $hit ( @$hits ){
131 :    
132 :     if ( $synonyms-> {$hit} ){
133 :     push ( @$same_length, $hit) if ( $synonyms-> {$hit} == $synonyms-> {$peg} );
134 :     }
135 :     else{
136 :     print STDERR "ERROR, shouldn't be here\n";
137 :     die;
138 :     }
139 :    
140 :     }
141 :    
142 :     if ( scalar @$same_length == 1){
143 :     #$gi_number = $same_length->[0];
144 :     $error = 1;
145 :     @return_values = ( $error ,
146 :     $same_length ,
147 :     "gi found for $peg in $genome_id", 1 );
148 :     }
149 :     elsif ( scalar @$same_length > 1){
150 :     $error = 0;
151 :     @return_values = ( $error , $same_length , "multiple gi's with same length found for $peg" , -2 );
152 :    
153 :     # need some magic to resolve
154 :    
155 :     }
156 :     else{
157 :     # hits but length is not matching
158 :     $error = 0;
159 :    
160 :     @return_values = ( $error , $hits , "no gi's length is matching $peg" , -4 );
161 :     }
162 :    
163 :    
164 :     }
165 :     elsif( scalar @$hits == 1){
166 :     $gi_list = $hits;
167 :     $return_values[1] = $gi_list;
168 :    
169 :     if ( _compare_length($fig, $peg, $hits->[0], $synonyms) ){
170 :     $return_values[3] = 1;
171 :     }
172 :     else{
173 :     $error = 0;
174 :     @return_values = ( $error , $hits , "gi found for $peg and $genome_id, but length does not match" , -4 ) ;
175 :     # return -4;
176 :     }
177 :     }
178 :     else{
179 :     print STDERR "ERROR, shoudn't be here \n";
180 :     exit;
181 :     }
182 :    
183 :     }
184 :     elsif ( ref $hits ){
185 :     print STDERR scalar @$hits," hit(s): ",$hits, "\n"
186 :     }
187 :     else{
188 :     print STDERR "no gi for $peg in ".$fig->genus_species( $genome_id)." \n";
189 :     $return_values[2] = "no gi found for $peg in $genome_id";
190 :     $return_values[3] = -3;
191 :     #$gi_number = -3;
192 :     }
193 :    
194 :     }
195 :     else{
196 :     $return_values[2] = "no gi's found for $peg at all";
197 :     $return_values[3] = -1;
198 :     #print STDERR "no gis for $peg \n";
199 :     }
200 :    
201 :     return @return_values;
202 :     }
203 :    
204 :    
205 :    
206 :     sub get_gi{
207 :     my ($fig,$peg,$org) = @_;
208 :    
209 :     # initialize gi number, returns
210 :     # -1 no gi number found
211 :     # -2 too many matched gis
212 :     # -3 no gi in corresponding organism
213 :     # -4 gi found but length is unequal
214 :    
215 :     my $gi_number = -1;
216 :    
217 :     # get all synonyms for a given id
218 :     my $pid = $fig->maps_to_id($peg);
219 :     my @mapped = $fig->mapped_prot_ids($pid);
220 :    
221 :     # get organism for peg
222 :     my ($genome_id, $peg_number) = $fig->genome_and_peg_of($peg);
223 :    
224 :    
225 :     # get all gi numbers from synonym list and map to organism
226 :     my @gis; # list of all gi synonyms for a given peg
227 :     my @pegs; # list of all peg synonyms for a given peg
228 :     my $synonyms = {};
229 :     foreach my $map ( @mapped ){
230 :     push ( @gis, $map->[0] ) if ( $map->[0] =~ m/^gi\|/ );
231 :     push ( @pegs, $map->[0] ) if ( $map->[0] =~ m/^fig\|$genome_id/ );
232 :     $synonyms->{ $map->[0] } = $map->[1];
233 :    
234 :     }
235 :    
236 :    
237 :     if ( scalar(@gis) > 0){
238 :    
239 :    
240 :     # sort gis by organism
241 :     my $org_name = {};
242 :     foreach my $gi (@gis){
243 :    
244 :     if ( $org_name->{ $fig->org_of($gi) } or $org_name->{ unknown }){
245 :    
246 :     if ( $fig->org_of($gi) ){
247 :     push @{$org_name->{ $fig->org_of($gi) }} ,$gi;
248 :     }
249 :     else{
250 :     push @{$org_name->{ unknown }} ,$gi;
251 :     }
252 :    
253 :    
254 :     }
255 :     else{
256 :    
257 :     # in different genomes
258 :     if ( $fig->org_of($gi) ) {
259 :     $org_name->{ $fig->org_of($gi) } = [$gi];
260 :     }
261 :     else{
262 :     $org_name->{ unknown } = [$gi];
263 :     }
264 :     }
265 :     }
266 :    
267 :    
268 :     # check genome name
269 :     # hope that fig genome name match genome name from refseq/genbank
270 :    
271 :     my $hits = [ ] ;
272 :     $hits = $org_name->{ $fig->genus_species( $genome_id) } if ( defined $fig->genus_species( $genome_id) );
273 :    
274 :     if ( ref $hits and scalar @$hits ){
275 :    
276 :     # check for hit counts
277 :     # resolve if more than one hit per organism
278 :    
279 :     #print scalar @$hits," hit(s): ", @$hits, "\n" ;
280 :    
281 :     if ( scalar @$hits > 1 ){
282 :    
283 :     # paralog within the same genome
284 :     # do some magic stuff to resolve
285 :    
286 :     # 1. compare by length
287 :     # 2. compare by sequence
288 :     # 3. compare by start and stop / position in genome
289 :    
290 :    
291 :    
292 :     my $pid = $fig->maps_to_id($peg);
293 :     my @mapped = $fig->mapped_prot_ids($pid);
294 :    
295 :     my ($genome_id, $peg_number) = $fig->genome_and_peg_of($peg);
296 :    
297 :    
298 :     my @gis; # list of all gi synonyms for a given peg
299 :     my @pegs; # list of all peg synonyms for a given peg
300 :    
301 :    
302 :    
303 :     # compare length of peg with length of hits
304 :     my $same_length = [];
305 :     foreach my $hit ( @$hits ){
306 :    
307 :     if ( $synonyms-> {$hit} ){
308 :     push ( @$same_length, $hit) if ( $synonyms-> {$hit} == $synonyms-> {$peg} );
309 :     }
310 :     else{
311 :     print STDERR "ERROR, shouldn't be here\n";
312 :     die;
313 :     }
314 :    
315 :     }
316 :    
317 :     if ( scalar @$same_length == 1){
318 :     $gi_number = $same_length->[0];
319 :     }
320 :     else{
321 :     # more magic stuff
322 :     return -2 ;
323 :     }
324 :    
325 :    
326 :     }
327 :     elsif( scalar @$hits == 1){
328 :     if ( _compare_length($fig, $peg, $hits->[0], $synonyms) ){
329 :     $gi_number = $hits->[0];
330 :     }
331 :     else{
332 :     return -4;
333 :     }
334 :     }
335 :     else{
336 :     print STDERR "ERROR, shoudn't be here \n";
337 :     exit;
338 :     }
339 :    
340 :     }
341 :     elsif ( ref $hits ){
342 :     print STDERR scalar @$hits," hit(s): ",$hits, "\n"
343 :     }
344 :     else{
345 :     print STDERR "no gi for $peg in ".$fig->genus_species( $genome_id)." \n";
346 :     $gi_number = -3;
347 :     }
348 :    
349 :     }
350 :     else{
351 :     print STDERR "no gis for $peg \n";
352 :     }
353 :    
354 :     return $gi_number;
355 :     }
356 :    
357 :    
358 :    
359 :    
360 :    
361 :     sub _compare_length{
362 :     my ($fig, $peg, $gi, $synonyms) = @_;
363 :    
364 :     if (defined $synonyms and ref $synonyms) {
365 :     return $synonyms->{ $peg } if ( $synonyms->{ $peg } == $synonyms->{ $gi } );
366 :     }
367 :     else{
368 :    
369 :     my $pid = $fig->maps_to_id($peg);
370 :     my @mapped = $fig->mapped_prot_ids($pid);
371 :    
372 :     my ($genome_id, $peg_number) = $fig->genome_and_peg_of($peg);
373 :    
374 :     my $synonyms = {};
375 :    
376 :     foreach my $map ( @mapped ){
377 :     $synonyms->{ $map->[0] } = $map->[1];
378 :     }
379 :    
380 :     return $synonyms->{ $peg } if ( $synonyms->{ $peg } == $synonyms->{ $gi } );
381 :    
382 :     }
383 :     return -1
384 :     }
385 :    
386 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3