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

Annotation of /FigKernelPackages/FigGFF.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 #
2 :     # FIG GFF utilities.
3 :     #
4 :    
5 :     #
6 :     # A GFFWriter handles the generation of GFF files from SEED data structures.
7 :     #
8 :    
9 : olson 1.4 package FigGFF;
10 :    
11 :     use strict;
12 :    
13 :     use base qw(Exporter);
14 :     use vars qw(@EXPORT);
15 :     @EXPORT = qw(map_seed_alias_to_dbxref);
16 :    
17 :     #
18 :     # General GFF-related routines.
19 :     #
20 :    
21 :    
22 :     #
23 :     # Alias translation.
24 :     #
25 :     # These routines map between the SEED aliases and the standard
26 :     # dbxref names as defined here:
27 :     #
28 :     # ftp://ftp.geneontology.org/pub/go/doc/GO.xrf_abbs
29 :     #
30 :     # In particular:
31 :     #
32 :     # abbreviation: NCBI_gi
33 :     # database: NCBI databases.
34 :     # object: Identifier.
35 :     # example_id: NCBI_gi:10727410
36 :     # generic_url: http://www.ncbi.nlm.nih.gov/
37 :     # url_syntax:
38 :     # url_example: http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=nucleotide&list_uids=10727410&dopt=GenBank
39 :     #
40 :     # abbreviation: NCBI_NP
41 :     # database: NCBI RefSeq.
42 :     # object: Protein identifier.
43 :     # example_id: NCBI_NP:123456
44 :     # generic_url: http://www.ncbi.nlm.nih.gov/
45 :     # url_syntax:
46 :     #
47 :     # abbreviation: KEGG
48 :     # database: Kyoto Encyclopedia of Genes and Genomes.
49 :     # generic_url: http://www.genome.ad.jp/kegg/
50 :    
51 :     sub map_seed_alias_to_dbxref
52 :     {
53 :     my($alias) = @_;
54 :    
55 :     $_ = $alias;
56 :     if (/^NP_(\d+.*)/)
57 :     {
58 :     return "NCBI_NP:$1";
59 :     }
60 :     elsif (/^gi\|(\d+)/)
61 :     {
62 :     return "NCBI_gi:$1";
63 :     }
64 :     elsif (/^kegg\|(\S+):(\S+)/)
65 :     {
66 :     return "KEGG:$1 $2";
67 :     }
68 :     elsif (/^uni\|(\S+)/)
69 :     {
70 :     return "UniProt:$1";
71 :     }
72 :     elsif (/^sp\|(\S+)/)
73 :     {
74 :     return "Swiss-Prot:$1";
75 :     }
76 :    
77 :     return undef;
78 :     }
79 :    
80 :     #
81 :     # And map back again.
82 :     #
83 :    
84 :     sub map_dbxref_to_seed_alias
85 :     {
86 :     my($dbxref) = @_;
87 :    
88 :     my($type, $ref) = split(/:/, $dbxref, 2);
89 :    
90 :     if ($type eq "NCBI_NP")
91 :     {
92 :     return "NP_$ref";
93 :     }
94 :     elsif ($type eq "NCBI_gi")
95 :     {
96 :     return "gi|$ref";
97 :     }
98 :     elsif ($type eq "KEGG")
99 :     {
100 :     $ref =~ s/ /:/;
101 :     return "kegg|$ref";
102 :     }
103 :     elsif ($type eq "UniProt")
104 :     {
105 :     return "uni|$ref";
106 :     }
107 :     elsif ($type eq "Swiss-Prot")
108 :     {
109 :     return "sp|$ref";
110 :     }
111 :    
112 :     return undef;
113 :     }
114 : olson 1.1
115 :     package GFFWriter;
116 :    
117 :     use strict;
118 :     use FIG;
119 : olson 1.4 use FigGFF;
120 : olson 1.1
121 :     use Carp;
122 :     use URI::Escape;
123 :     use Data::Dumper;
124 :    
125 :     sub new
126 :     {
127 :     my($class, $fig, %options) = @_;
128 :    
129 :     my $default_options = {
130 :     escapespace => 0,
131 :     outputfasta => 1,
132 :     linelength => 60,
133 :     };
134 :    
135 :     map { $default_options->{$_} = $options{$_} } keys(%options);
136 :    
137 :     my $self = {
138 :     options => $default_options,
139 :     contig_length_cache => {},
140 :     fig => $fig,
141 :     };
142 :    
143 :     return bless $self, $class;
144 :     }
145 :    
146 :    
147 :     =head1 gff3_for_feature
148 :    
149 :     Returns the GFF3 information for a given feature.
150 :    
151 :     The return is a pair ($contig_data, $fasta_sequences) that can be passed
152 :     into write_gff3().
153 :    
154 :     $contig_data is a hashref mapping a contig name to a list of GFF3 file lines
155 :     for the sequences in that contig.
156 :    
157 :     =cut
158 :    
159 :     sub gff3_for_feature
160 :     {
161 : olson 1.4 my($self, $fid, $user, $user_note, $in_aliases, $in_loc) = @_;
162 : olson 1.1
163 :     #
164 :     # Options we need to figure out somehow.
165 :     #
166 :     my $options = $self->{options};
167 :    
168 :     my $escapespace = $options->{escapespace};
169 :     my $outputfasta = $options->{outputfasta};
170 : olson 1.2
171 :     my %outputtype;
172 :     map { $outputtype{$_} = 1 } @{$options->{outputtype}};
173 : olson 1.1
174 :     my $fastasequences = '';
175 :     my $contig_data;
176 :     my $linelength = $options->{linelength};
177 :    
178 : olson 1.3 my $beg = $self->{options}->{beg};
179 :     my $end = $self->{options}->{end};
180 :    
181 : olson 1.1 my $fig = $self->{fig};
182 :    
183 :     #
184 :     # Do this first to make sure that we really have a feature.
185 :     #
186 : olson 1.4 my @location = ref($in_loc) ? @$in_loc : $fig->feature_location($fid);
187 : olson 1.1 if (@location == 0 or !defined($location[0]))
188 :     {
189 :     warn "No location found for feature $fid\n";
190 :     return ({}, "");
191 :     }
192 :    
193 :     ###########
194 :     #
195 :     # Begin figuring out the column 9 information about notes and aliases and GO terms
196 :     # All the information is temporarily stored in @alias or @note, and at the end is joined
197 :     # into $allnote
198 :     #
199 :     ###########
200 :    
201 :     #
202 :     # the notes for the last column
203 :     #
204 :     my $note;
205 :     #
206 :     # all the aliases we are going to use
207 :     #
208 :     my @alias;
209 : olson 1.4 my @xref;
210 : olson 1.1
211 : olson 1.2 if ($options->{with_assignments})
212 : olson 1.1 {
213 : olson 1.2 my $func = $fig->function_of($fid, $user);
214 :     if ($func)
215 :     {
216 : olson 1.4 push @$note, ("Note=" . uri_escape($func));
217 : olson 1.2 }
218 : olson 1.1 }
219 : olson 1.2
220 :     if ($options->{with_aliases})
221 :     {
222 : olson 1.4 # now find aliases
223 :     my @feat_aliases = ref($in_aliases) ? @$in_aliases : $fig->feature_aliases($fid);
224 :     foreach my $alias (@feat_aliases)
225 : olson 1.1 {
226 : olson 1.4 my $mapped = FigGFF::map_seed_alias_to_dbxref($alias);
227 :     if ($mapped)
228 : olson 1.2 {
229 : olson 1.4 push(@xref, $mapped);
230 : olson 1.2 }
231 :     else
232 :     {
233 : olson 1.4 push(@alias, $alias);
234 : olson 1.2 }
235 : olson 1.1 }
236 : olson 1.2 }
237 : olson 1.1
238 :     # now just join all the aliases and put them into @$note so we can add it to the array
239 :     if (@alias)
240 :     {
241 : olson 1.4 push @$note, "Alias=". join (",", map { uri_escape($_) } @alias);
242 : olson 1.1 }
243 : olson 1.2
244 :     #
245 :     # If we have user note passed in, add it.
246 :     #
247 :    
248 :     if ($user_note)
249 :     {
250 :     push @$note, $user_note;
251 :     }
252 : olson 1.1
253 :     # the LAST thing I am going to add as a note is the FIG id so that I can grep it out easily
254 : olson 1.4 #
255 :     # for now, make SEED xref the first in the list so we can search for DBxref=SEEd
256 :     #
257 :    
258 :     unshift(@xref, "SEED:$fid");
259 :    
260 :     push(@$note, "Dbxref=" . join(",", map { uri_escape($_) } @xref));
261 : olson 1.1
262 :     # finally join all the notes into a long string that can be added as column 9
263 :     my $allnotes;
264 :     $allnotes = join ";", @$note;
265 :    
266 :     # do we want to convert '%20' to ' '
267 :     unless ($escapespace)
268 :     {
269 :     $allnotes =~ s/\%20/ /g;
270 :     }
271 :    
272 :     ###########
273 :     #
274 :     # End figuring out the column 9 information about notes and aliases and GO terms
275 :     #
276 :     ###########
277 :    
278 :     #
279 :     # Cache contig lengths.
280 :     #
281 :     my $len = $self->{contig_length_cache};
282 :    
283 :     my $genome = $fig->genome_of($fid);
284 :    
285 :     foreach my $loc (@location)
286 :     {
287 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
288 :     my ($contig, $start, $stop) = ($1, $2, $3);
289 :     my $original_contig=$contig;
290 :    
291 :     #
292 :     # the contig name must be escaped
293 :     #
294 :     $contig = uri_escape($contig);
295 :    
296 :     #my $contig_key = "$genome:$contig";
297 :     my $contig_key = $contig;
298 :    
299 :     unless (defined $len->{$contig})
300 :     {
301 :     $len->{$contig}=$fig->contig_ln($genome, $original_contig);
302 :     }
303 :     my $strand='+';
304 :    
305 :     #
306 :     # These were bounds-checking for dumping all of a genome.
307 :     #
308 :     #next if (defined $beg && ($start < $beg || $stop < $beg));
309 :     #next if (defined $end && ($start > $end || $stop > $end));
310 :    
311 :     if ($start > $stop)
312 :     {
313 :     ($start, $stop, $strand)=($stop, $start, '-');
314 :     }
315 :     elsif ($start == $stop)
316 :     {
317 :     $strand=".";
318 :     }
319 :    
320 :     my $type=$fig->ftype($fid);
321 :    
322 :     if ($type eq "peg")
323 :     {
324 :     # it is a protein coding gene
325 :     # create an artificial id that is just the fid.(\d+) information
326 :     # we will use this to create ids in the form cds.xxx; trn.xxxx; pro.xxx; gen.xxx;
327 :     $fid =~ /\.peg\.(\d+)/;
328 :     my $geneid=$1;
329 :    
330 :     ############## KLUDGE
331 :     #
332 :     # At the moment the outputs for transcript, gene, CDS, and pro are all the same.
333 :     # This is clearly a kludge and wrong, but it will work at the moment
334 :     #
335 :    
336 :     # defined some truncations
337 :     my %trunc=(
338 :     "transcript" => "trn",
339 :     "gene" => "gen",
340 :     "protein" => "pro",
341 :     "cds" => "cds",
342 :     );
343 :    
344 :     # SO terms:
345 :     # transcript: SO:0000673
346 :     # gene: SO:0000704
347 :     # cds: SO:0000316
348 :     # protein: NOT VALID: should be protein_coding_primary_transcript SO:0000120
349 : olson 1.3
350 :     #
351 :     # For now, we will only output CDS features, and include
352 :     # the translation as an attribute (attribuute name is
353 :     # translation_id, value is a key that corresponds to a FASTA
354 :     # section at the end of the file).
355 :     #
356 :    
357 :     my $type = "cds";
358 :    
359 :     my $protein_id = "pro.$geneid";
360 :     my $cds_id = "cds.$geneid";
361 :    
362 :     # we want to store some sequences to be output
363 :     if ($outputfasta)
364 : olson 1.1 {
365 : olson 1.3 my $addseq = $fig->get_translation($fid);
366 :    
367 :     #
368 :     # the chomp is so that we know for sure to add the line back
369 :     #
370 :     $addseq =~ s/(.{$linelength})/$1\n/g;
371 :     chomp($addseq);
372 :     $fastasequences .= ">$protein_id\n$addseq\n";
373 :    
374 : olson 1.4 my $addseq = uc($fig->dna_seq($genome, @location));
375 : olson 1.3 $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
376 :    
377 :     $fastasequences .= ">$cds_id\n$addseq\n";
378 :    
379 :     $allnotes .= ";translation_id=$protein_id";
380 :     }
381 :    
382 :     push (@{$contig_data->{$contig_key}},
383 :     (join "\t",
384 :     ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$cds_id;$allnotes")));
385 :     }
386 : olson 1.1 elsif ($type eq "rna")
387 :     {
388 :     $fid =~ /\.rna\.(\d+)/;
389 :     my $geneid=$1;
390 :     #
391 :     # tRNA is a valid SOFA term == SO:0000253
392 :     #
393 :     my ($id, $type)=("rna.$geneid", "tRNA");
394 :     if ($outputfasta)
395 :     {
396 : olson 1.4 my $addseq = $fig->dna_seq($genome, @location);
397 : olson 1.1 $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
398 :     $fastasequences .= ">$id\n$addseq\n";
399 :     }
400 : olson 1.2 push (@{$contig_data->{$contig_key}}, (join "\t", ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$id;$allnotes")));
401 : olson 1.1 } # end the if type == rna
402 :     else
403 :     {
404 :     die "Don't know what type: |$type| is";
405 :     }
406 :     }
407 :     return ($contig_data, $fastasequences);
408 :     }
409 :    
410 :     =head1 write_gff3
411 :    
412 :     Write a set of gff3 per-contig data and fasta sequence data to a file or filehandle.
413 :    
414 :     $genome is the genome these contigs are a part of.
415 :     $contig_list is a list of contig-data hashes as returned by gff_for_feature.
416 :     $fast_list is a list of fasta data strings.
417 :    
418 :     =cut
419 :    
420 :     sub write_gff3
421 :     {
422 :     my($self, $output, $genome, $contig_list, $fasta_list) = @_;
423 :    
424 :     my $fig = $self->{fig};
425 :    
426 :     my $len = $self->{contig_length_cache};
427 :     my $fh;
428 :    
429 :     my $beg = $self->{options}->{beg};
430 :     my $end = $self->{options}->{end};
431 :    
432 :     my $close_output;
433 :    
434 :     if (ref($output))
435 :     {
436 :     $fh = $output;
437 :     }
438 :     else
439 :     {
440 :     open($fh, ">$output") or confess "Cannot open output '$output': $!";
441 :     $close_output = 1;
442 :     }
443 :    
444 :     #
445 :     # Build a data structure from the list of contigs
446 :     # that has a list of lists of data per contig name.
447 :     # (Do this so we don't copy all of the contig data itself, as it
448 :     # could be quite large).
449 :     #
450 :     my %contigs;
451 :    
452 :     #
453 :     # iterate over the given list of contig hashes.
454 :     #
455 :     for my $chash (@$contig_list)
456 :     {
457 :     #
458 :     # Then for each contig in the individual contig hashes,
459 :     # add the data list to %contigs.
460 :     #
461 :     for my $contig (keys %$chash)
462 :     {
463 :     push(@{$contigs{$contig}}, $chash->{$contig});
464 :     }
465 :     }
466 :    
467 :     foreach my $contig (sort keys %contigs)
468 :     {
469 :     print $fh "##sequence-region\t$contig\t";
470 :     if (defined $beg) {
471 :     print $fh "$beg\t";
472 :     } else {
473 :     print $fh "1\t";
474 :     }
475 :     if (defined $end) {
476 :     print $fh "$end\n";
477 :     } else {
478 :     print $fh "$len->{$contig}\n";
479 :     }
480 :     for my $list (@{$contigs{$contig}})
481 :     {
482 : olson 1.2 print $fh join("\n", @$list), "\n";
483 : olson 1.1 }
484 :     }
485 :    
486 :     print $fh "##FASTA\n";
487 :     # print out the cds and pro if we need them
488 :    
489 :     if ($self->{options}->{outputfasta})
490 :     {
491 :     for my $fastasequences (@$fasta_list)
492 :     {
493 :     print $fh $fastasequences;
494 :     }
495 :     }
496 :    
497 :     my $ll = $self->{options}->{linelength};
498 :     foreach my $contig (sort keys %contigs)
499 :     {
500 :     my $len=$fig->contig_ln($genome, $contig);
501 :     my $dna_seq=$fig->dna_seq($genome, $contig . "_1_". $len);
502 :     if (defined $beg)
503 :     {
504 :     unless (defined $end) {
505 :     $end=$len;
506 :     }
507 :     $dna_seq = substr($dna_seq, $beg, $end);
508 :     }
509 :     elsif (defined $end)
510 :     {
511 :     $beg=1;
512 :     $dna_seq = substr($dna_seq, $beg, $end);
513 :     }
514 :    
515 :     my $contig=uri_escape($contig);
516 :    
517 :     $dna_seq =~ s/(.{$ll})/$1\n/g;
518 :     chomp($dna_seq); # just remove the last \n if there is one
519 :     print $fh ">$contig\n$dna_seq\n";
520 :     }
521 :    
522 :     close($fh) if $close_output;
523 :     }
524 :    
525 : olson 1.2 package GFFParser;
526 :    
527 :     use strict;
528 :     use URI::Escape;
529 :     use Carp;
530 :     use Data::Dumper;
531 :    
532 :     use base qw(Class::Accessor);
533 :    
534 :     __PACKAGE__->mk_accessors(qw(fig current_file));
535 :    
536 :     my $count;
537 :    
538 :    
539 :     #
540 :     # GFF file parser. Creates GFFFiles.
541 :     #
542 :    
543 :     sub new
544 :     {
545 :     my($class, $fig) = @_;
546 :    
547 :     my $self = {
548 :     fig => $fig,
549 :     };
550 :    
551 :     return bless($self, $class);
552 :     }
553 :    
554 :     sub parse
555 :     {
556 :     my($self, $file) = @_;
557 :    
558 :     my($fh, $close_handle);
559 :    
560 :     my $fobj = GFFFile->new($self->fig);
561 :     $self->current_file($fobj);
562 :    
563 :     if (ref($file) ? (ref($file) eq 'GLOB'
564 :     || UNIVERSAL::isa($file, 'GLOB')
565 :     || UNIVERSAL::isa($file, 'IO::Handle'))
566 :     : (ref(\$file) eq 'GLOB'))
567 :     {
568 :     $fh = $file;
569 :     }
570 :     else
571 :     {
572 :     open($fh, "<$file") or confess "Cannot open $file: $!";
573 :     $fobj->filename($file);
574 :     $close_handle = 1;
575 :     }
576 :    
577 :     #
578 :     # Start parsing by verifying this is a gff3 file.
579 :     #
580 :    
581 :     $_ = <$fh>;
582 :    
583 :     if (m,^\#gff-version\t(\S+),)
584 :     {
585 :     if ($1 != 3)
586 :     {
587 :     confess "Invalid GFF File: version is not 3";
588 :     }
589 :     }
590 :    
591 :     #
592 :     # Now parse.
593 :     #
594 :    
595 :     while (<$fh>)
596 :     {
597 :     chomp;
598 :     #
599 :     # Check first for the fasta directive so we can run off and parse that
600 :     # separately.
601 :     #
602 :    
603 :     if (/^>/)
604 :     {
605 :     $self->parse_fasta($fh, $_);
606 :     last;
607 :     }
608 :     elsif (/^\#\#FASTA/)
609 :     {
610 :     print "Got fasta directive\n";
611 :     $_ = <$fh>;
612 :     chomp;
613 :     $self->parse_fasta($fh, $_);
614 :     last;
615 :     }
616 :     elsif (/^\#\s/)
617 :     {
618 :     #
619 :     # comment.
620 :     #
621 :     next;
622 :     }
623 :     elsif (/^\#\#(\S+)(?:\t(.*))?/)
624 :     {
625 :     #
626 :     # GFF3 directive.
627 :     #
628 :    
629 :     $self->parse_gff3_directive($1, $2);
630 :    
631 :     }
632 :     elsif (/^\#(\S+)(?:\t(.*))?/)
633 :     {
634 :     #
635 :     # Directive.
636 :     #
637 :    
638 : redwards 1.5 if (lc($1) eq "seed")
639 : olson 1.2 {
640 :     $self->parse_seed_directive($2);
641 :     }
642 :     else
643 :     {
644 :     $self->parse_local_directive($1, $2);
645 :     }
646 :    
647 :     }
648 :     elsif (/^([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)$/)
649 :     {
650 :     $self->parse_feature($1, $2, $3, $4, $5, $6, $7, $8, $9);
651 :     }
652 :     else
653 :     {
654 :     die "bad line: '$_'\n";
655 :     }
656 :     }
657 :    
658 :     return $fobj;
659 :     }
660 :    
661 :     sub parse_gff3_directive
662 :     {
663 :     my($self, $directive, $rest) = @_;
664 :    
665 :     print "Have gff3 directive '$directive' rest='$rest'\n";
666 :     }
667 :    
668 :     sub parse_seed_directive
669 :     {
670 :     my($self, $rest) = @_;
671 :    
672 :     my($verb, @rest) = split(/\t/, $rest);
673 :    
674 :     if ($verb eq "anno_start")
675 :     {
676 :     $self->current_file->anno_start($rest[0]);
677 :     }
678 :     elsif ($verb eq "anno_end")
679 :     {
680 :     $self->current_file->anno_start($rest[0]);
681 :     }
682 :     elsif ($verb eq "genome_md5")
683 :     {
684 :     $self->current_file->set_genome_checksum(@rest[0,1]);
685 :     }
686 :     elsif ($verb eq "contig_md5")
687 :     {
688 :     $self->current_file->set_contig_checksum(@rest[0,1,2]);
689 :     }
690 :     }
691 :    
692 :     sub parse_local_directive
693 :     {
694 :     my($self, $directive, $rest) = @_;
695 :    
696 :     print "Have local directive '$directive' rest='$rest'\n";
697 :     }
698 :    
699 :     sub parse_feature
700 :     {
701 :     my($self, $seqid, $source, $type, $start, $end, $score, $strand, $phase, $attributes) = @_;
702 :    
703 :     #print "data: seqid=$seqid source=$source type=$type start=$start end=$end\n";
704 :     #print " score=$score strand=$strand phase=$phase\n";
705 :     #print " $attributes\n";
706 :    
707 :     #
708 :     # Parse this feature, creating a GFFFeature object for it.
709 :     #
710 :    
711 :     my $feature = GFFFeature->new($self->fig);
712 :    
713 :     $feature->seqid($seqid);
714 :     $feature->source($source);
715 :     $feature->type($type);
716 :     $feature->start($start);
717 :     $feature->end($end);
718 :     $feature->score($score);
719 :     $feature->strand($strand);
720 :     $feature->phase($phase);
721 :    
722 :     my $athash = {};
723 :    
724 :     for my $attr (split(/;/, $attributes))
725 :     {
726 :     my($name, $value) = split(/=/, $attr);
727 :    
728 : olson 1.4 my @values = map { uri_unescape($_) } split(/,/, $value);
729 :    
730 : olson 1.2 #
731 : olson 1.4 # This might be a little goofy for the users, but we will use it
732 :     # for now:
733 :     #
734 :     # if there is more than one value, the value is a ref to a list
735 :     # of the values.
736 :     #
737 :     # Otherwise, the value is a scalar.
738 : olson 1.2 #
739 :    
740 :     if (@values > 1)
741 :     {
742 : olson 1.4 #
743 :     # Yes, you can do this ... I had to look it up :-).
744 :     #
745 :     # It's in 'man perlfaq3'.
746 :     #
747 :    
748 :     $value = \@values;
749 : olson 1.2 }
750 :    
751 :     $athash->{$name} = $value;
752 :    
753 : olson 1.4 #
754 :     # Handle the GFF3-defined attributes.
755 :     #
756 :     # These show up as Class::Accessor's on the feature object.
757 :     #
758 :    
759 : olson 1.2 if ($GFFFeature::GFF_Tags{$name})
760 :     {
761 :     $feature->set($name, $value);
762 : olson 1.4
763 :     if ($name eq "Dbxref")
764 :     {
765 :     #
766 :     # If we have a SEED:figid DBxref, set the genome and fig_id attributes.
767 :     #
768 :    
769 :     my @seed_xref = grep /^"?SEED:/, @values;
770 :     if (@seed_xref and $seed_xref[0] =~ /^"?SEED:(fig\|(\d+\.\d+)\..*)/)
771 :     {
772 :     $feature->genome($2);
773 :     $feature->fig_id($1);
774 :     }
775 :    
776 :     }
777 : olson 1.2 }
778 :     }
779 :     $feature->attributes($athash);
780 :    
781 : olson 1.4
782 : olson 1.2 $self->current_file->add_feature($feature);
783 :     }
784 :    
785 :     #
786 :     # We come in here with the first line of the fasta already read
787 :     # in order to support the backward-compatiblity syntax that
788 :     # lets a file skip the ##FASTA directive if it wishes.
789 :     #
790 :     sub parse_fasta
791 :     {
792 :     my($self, $fh, $first_line) = @_;
793 :     my($cur, $cur_id);
794 :    
795 :     for ($_ = $first_line; $_; $_ = <$fh>, chomp)
796 :     {
797 :     if (/^>\s*(\S+)/)
798 :     {
799 :     if ($cur)
800 :     {
801 :     $self->handle_fasta_block($cur_id, $cur);
802 :     }
803 :    
804 :     $cur = '';
805 :     $cur_id = $1;
806 :     }
807 :     else
808 :     {
809 :     s/^\s*$//;
810 :     s/\s*$//;
811 :     if (/\s/)
812 :     {
813 :     die "FASTA data had embedded space: $_\n";
814 :     }
815 :     $cur .= $_;
816 :     }
817 :     }
818 :     if ($cur)
819 :     {
820 :     $self->handle_fasta_block($cur_id, $cur);
821 :     }
822 :     }
823 :    
824 :     sub handle_fasta_block
825 :     {
826 :     my($self, $id, $data) = @_;
827 :    
828 :     my $len = length($data);
829 :     $self->current_file->set_fasta_data($id, $data);
830 :     }
831 :    
832 :     package GFFFeature;
833 :    
834 :     use strict;
835 :     use base qw(Class::Accessor);
836 :    
837 :     our @GFF_Tags = qw(ID Name Alias Parent Target Gap Note Dbxref Ontology_term);
838 :     our %GFF_Tags;
839 :    
840 :     map { $GFF_Tags{$_} = 1 } @GFF_Tags;
841 :    
842 : olson 1.4 __PACKAGE__->mk_accessors(qw(fig seqid source type start end score strand phase attributes genome fig_id),
843 : olson 1.2 @GFF_Tags);
844 :    
845 :    
846 :     sub new
847 :     {
848 :     my($class, $fig) = @_;
849 :    
850 :     my $self = {
851 :     fig => $fig,
852 :     };
853 :    
854 :     return bless($self, $class);
855 :     }
856 :    
857 : olson 1.4 sub find_local_feature
858 :     {
859 :     my($self, $local_genome) = @_;
860 :     my $db = $self->fig->db_handle;
861 :    
862 :     # For debugging.
863 :     undef $local_genome;
864 :     if ($local_genome)
865 :     {
866 :     #
867 :     # It's a precise match. We need to determine if we have this
868 :     # particular feature in this SEED (it is possible for one to
869 :     # have exported an annotation for a feature that was added
870 :     # to a genome after its initial release).
871 :     #
872 :     # We do this by searching for a local feature with the same contig,
873 :     # start, and stop as this feature.
874 :     #
875 :    
876 :     my $qry = qq(SELECT id
877 :     FROM features
878 :     WHERE (genome = ? AND
879 :     contig = ? AND
880 :     minloc = ? AND
881 :     maxloc = ?));
882 :     my $res = $db->SQL($qry, undef, $local_genome, $self->seqid,
883 :     $self->start, $self->end);
884 :    
885 :     return map { $_->[0] } @$res;
886 :     }
887 :    
888 :     #
889 :     # Otherwise, we need to try a set of heuristics to match
890 :     # this id.
891 :     #
892 :    
893 :     #
894 :     # Try matching aliases first.
895 :     #
896 :    
897 :     my @aliases = grep { !/^\"?SEED/ } ref($self->Dbxref) ? @{$self->Dbxref} : ($self->Dbxref);
898 :    
899 :     my @maliases = map { FigGFF::map_dbxref_to_seed_alias($_) } @aliases;
900 :    
901 :     print "Found aliases @aliases\n";
902 :     print "Found mapped aliases @maliases\n";
903 :    
904 :     for my $malias (@maliases)
905 :     {
906 :     my $fid = $self->fig->by_alias($malias);
907 :     if ($fid)
908 :     {
909 :     print "Mapped $malias to $fid\n";
910 :     }
911 :     }
912 :    
913 :     }
914 :    
915 :    
916 : olson 1.2 package GFFFile;
917 :    
918 :     use strict;
919 :     use base qw(Class::Accessor);
920 :    
921 :     __PACKAGE__->mk_accessors(qw(fig filename features feature_index anno_start anno_end));
922 :    
923 :     #
924 :     # Package to hold the contents of a GFF file, and to hold the code
925 :     # for mapping its contents to the local SEED.
926 :     #
927 :     # Created by GFFParser->parse.
928 :     #
929 :    
930 :     sub new
931 :     {
932 :     my($class, $fig) = @_;
933 :    
934 :     my $self = {
935 :     fig => $fig,
936 :     features => [],
937 :     feature_index => {},
938 : olson 1.4 genome_checksum => {},
939 :     contig_checksum => {},
940 :     features_by_genome => {},
941 : olson 1.2 };
942 :     return bless($self, $class);
943 :     }
944 :    
945 :     sub add_feature
946 :     {
947 :     my($self, $feature) = @_;
948 :    
949 :     push(@{$self->features}, $feature);
950 :     $self->feature_index->{$feature->ID} = $feature;
951 : olson 1.4 push(@{$self->{features_by_genome}->{$feature->genome}}, $feature);
952 :     }
953 :    
954 :     sub features_for_genome
955 :     {
956 :     my($self, $genome) = @_;
957 :    
958 :     return $self->{features_by_genome}->{$genome};
959 :     }
960 :    
961 :     sub genome_checksum
962 :     {
963 :     my($self, $genome) = @_;
964 :    
965 :     return $self->{genome_checksum}->{$genome};
966 :     }
967 :    
968 :     sub genome_checksums
969 :     {
970 :     my($self) = @_;
971 :    
972 :     return [map { [$_, $self->{genome_checksum}->{$_}] } keys(%{$self->{genome_checksum}})];
973 : olson 1.2 }
974 :    
975 :     sub set_genome_checksum
976 :     {
977 :     my($self, $genome, $md5sum) = @_;
978 :     $self->{genome_checksum}->{$genome} = $md5sum;
979 :     }
980 :    
981 :     sub set_contig_checksum
982 :     {
983 :     my($self, $genome, $contig, $md5sum) = @_;
984 :     $self->{contig_checksum}->{$genome}->{$contig} = $md5sum;
985 :     }
986 :    
987 :     sub set_fasta_data
988 :     {
989 :     my($self, $id, $data) = @_;
990 :    
991 :     $self->{fasta_data}->{$id} = $data;
992 :     }
993 :    
994 :    
995 : olson 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3