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

Annotation of /FigKernelPackages/FigGFF.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (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 : redwards 1.6 # added contig_start_cache and contig_end_cache because we have something like
138 :     # sequence-region contigname 1 10000
139 :     # in general we will set contig_start_cache == 1
140 :     # and contig_end_cache == contig_length_cache
141 :    
142 : olson 1.1 my $self = {
143 :     options => $default_options,
144 :     contig_length_cache => {},
145 : redwards 1.6 contig_start_cache => {},
146 :     contig_end_cache => {},
147 : olson 1.1 fig => $fig,
148 :     };
149 :    
150 :     return bless $self, $class;
151 :     }
152 :    
153 :    
154 :     =head1 gff3_for_feature
155 :    
156 :     Returns the GFF3 information for a given feature.
157 :    
158 :     The return is a pair ($contig_data, $fasta_sequences) that can be passed
159 :     into write_gff3().
160 :    
161 :     $contig_data is a hashref mapping a contig name to a list of GFF3 file lines
162 :     for the sequences in that contig.
163 :    
164 :     =cut
165 :    
166 :     sub gff3_for_feature
167 :     {
168 : olson 1.4 my($self, $fid, $user, $user_note, $in_aliases, $in_loc) = @_;
169 : olson 1.1
170 :     #
171 :     # Options we need to figure out somehow.
172 :     #
173 :     my $options = $self->{options};
174 :    
175 :     my $escapespace = $options->{escapespace};
176 :     my $outputfasta = $options->{outputfasta};
177 : olson 1.2
178 :     my %outputtype;
179 :     map { $outputtype{$_} = 1 } @{$options->{outputtype}};
180 : olson 1.1
181 :     my $fastasequences = '';
182 :     my $contig_data;
183 :     my $linelength = $options->{linelength};
184 :    
185 : olson 1.3 my $beg = $self->{options}->{beg};
186 :     my $end = $self->{options}->{end};
187 :    
188 : olson 1.1 my $fig = $self->{fig};
189 :    
190 :     #
191 :     # Do this first to make sure that we really have a feature.
192 :     #
193 : olson 1.4 my @location = ref($in_loc) ? @$in_loc : $fig->feature_location($fid);
194 : olson 1.1 if (@location == 0 or !defined($location[0]))
195 :     {
196 :     warn "No location found for feature $fid\n";
197 :     return ({}, "");
198 :     }
199 :    
200 :     ###########
201 :     #
202 :     # Begin figuring out the column 9 information about notes and aliases and GO terms
203 :     # All the information is temporarily stored in @alias or @note, and at the end is joined
204 :     # into $allnote
205 :     #
206 :     ###########
207 :    
208 :     #
209 :     # the notes for the last column
210 :     #
211 :     my $note;
212 :     #
213 :     # all the aliases we are going to use
214 :     #
215 :     my @alias;
216 : olson 1.4 my @xref;
217 : olson 1.1
218 : olson 1.2 if ($options->{with_assignments})
219 : olson 1.1 {
220 : olson 1.2 my $func = $fig->function_of($fid, $user);
221 :     if ($func)
222 :     {
223 : olson 1.4 push @$note, ("Note=" . uri_escape($func));
224 : olson 1.2 }
225 : olson 1.1 }
226 : olson 1.2
227 :     if ($options->{with_aliases})
228 :     {
229 : olson 1.4 # now find aliases
230 :     my @feat_aliases = ref($in_aliases) ? @$in_aliases : $fig->feature_aliases($fid);
231 :     foreach my $alias (@feat_aliases)
232 : olson 1.1 {
233 : olson 1.4 my $mapped = FigGFF::map_seed_alias_to_dbxref($alias);
234 :     if ($mapped)
235 : olson 1.2 {
236 : olson 1.4 push(@xref, $mapped);
237 : olson 1.2 }
238 :     else
239 :     {
240 : olson 1.4 push(@alias, $alias);
241 : olson 1.2 }
242 : olson 1.1 }
243 : olson 1.2 }
244 : olson 1.1
245 :     # now just join all the aliases and put them into @$note so we can add it to the array
246 :     if (@alias)
247 :     {
248 : olson 1.4 push @$note, "Alias=". join (",", map { uri_escape($_) } @alias);
249 : olson 1.1 }
250 : olson 1.2
251 :     #
252 :     # If we have user note passed in, add it.
253 :     #
254 :    
255 :     if ($user_note)
256 :     {
257 :     push @$note, $user_note;
258 :     }
259 : olson 1.1
260 :     # the LAST thing I am going to add as a note is the FIG id so that I can grep it out easily
261 : olson 1.4 #
262 :     # for now, make SEED xref the first in the list so we can search for DBxref=SEEd
263 :     #
264 :    
265 :     unshift(@xref, "SEED:$fid");
266 :    
267 :     push(@$note, "Dbxref=" . join(",", map { uri_escape($_) } @xref));
268 : olson 1.1
269 :     # finally join all the notes into a long string that can be added as column 9
270 :     my $allnotes;
271 :     $allnotes = join ";", @$note;
272 :    
273 :     # do we want to convert '%20' to ' '
274 :     unless ($escapespace)
275 :     {
276 :     $allnotes =~ s/\%20/ /g;
277 :     }
278 :    
279 :     ###########
280 :     #
281 :     # End figuring out the column 9 information about notes and aliases and GO terms
282 :     #
283 :     ###########
284 :    
285 :     #
286 :     # Cache contig lengths.
287 :     #
288 :     my $len = $self->{contig_length_cache};
289 :    
290 :     my $genome = $fig->genome_of($fid);
291 :    
292 :     foreach my $loc (@location)
293 :     {
294 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
295 :     my ($contig, $start, $stop) = ($1, $2, $3);
296 :     my $original_contig=$contig;
297 :    
298 :     #
299 :     # the contig name must be escaped
300 :     #
301 :     $contig = uri_escape($contig);
302 :    
303 :     #my $contig_key = "$genome:$contig";
304 :     my $contig_key = $contig;
305 :    
306 :     unless (defined $len->{$contig})
307 :     {
308 :     $len->{$contig}=$fig->contig_ln($genome, $original_contig);
309 :     }
310 :     my $strand='+';
311 :    
312 :     #
313 :     # These were bounds-checking for dumping all of a genome.
314 :     #
315 :     #next if (defined $beg && ($start < $beg || $stop < $beg));
316 :     #next if (defined $end && ($start > $end || $stop > $end));
317 :    
318 :     if ($start > $stop)
319 :     {
320 :     ($start, $stop, $strand)=($stop, $start, '-');
321 :     }
322 :     elsif ($start == $stop)
323 :     {
324 :     $strand=".";
325 :     }
326 :    
327 :     my $type=$fig->ftype($fid);
328 :    
329 :     if ($type eq "peg")
330 :     {
331 :     # it is a protein coding gene
332 :     # create an artificial id that is just the fid.(\d+) information
333 :     # we will use this to create ids in the form cds.xxx; trn.xxxx; pro.xxx; gen.xxx;
334 :     $fid =~ /\.peg\.(\d+)/;
335 :     my $geneid=$1;
336 :    
337 :     ############## KLUDGE
338 :     #
339 :     # At the moment the outputs for transcript, gene, CDS, and pro are all the same.
340 :     # This is clearly a kludge and wrong, but it will work at the moment
341 :     #
342 :    
343 :     # defined some truncations
344 :     my %trunc=(
345 :     "transcript" => "trn",
346 :     "gene" => "gen",
347 :     "protein" => "pro",
348 :     "cds" => "cds",
349 :     );
350 :    
351 :     # SO terms:
352 :     # transcript: SO:0000673
353 :     # gene: SO:0000704
354 :     # cds: SO:0000316
355 :     # protein: NOT VALID: should be protein_coding_primary_transcript SO:0000120
356 : olson 1.3
357 :     #
358 :     # For now, we will only output CDS features, and include
359 :     # the translation as an attribute (attribuute name is
360 :     # translation_id, value is a key that corresponds to a FASTA
361 :     # section at the end of the file).
362 :     #
363 :    
364 :     my $type = "cds";
365 :    
366 :     my $protein_id = "pro.$geneid";
367 :     my $cds_id = "cds.$geneid";
368 :    
369 :     # we want to store some sequences to be output
370 :     if ($outputfasta)
371 : olson 1.1 {
372 : olson 1.3 my $addseq = $fig->get_translation($fid);
373 :    
374 :     #
375 :     # the chomp is so that we know for sure to add the line back
376 :     #
377 :     $addseq =~ s/(.{$linelength})/$1\n/g;
378 :     chomp($addseq);
379 :     $fastasequences .= ">$protein_id\n$addseq\n";
380 :    
381 : olson 1.4 my $addseq = uc($fig->dna_seq($genome, @location));
382 : olson 1.3 $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
383 :    
384 :     $fastasequences .= ">$cds_id\n$addseq\n";
385 :    
386 :     $allnotes .= ";translation_id=$protein_id";
387 :     }
388 :    
389 :     push (@{$contig_data->{$contig_key}},
390 :     (join "\t",
391 :     ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$cds_id;$allnotes")));
392 :     }
393 : olson 1.1 elsif ($type eq "rna")
394 :     {
395 :     $fid =~ /\.rna\.(\d+)/;
396 :     my $geneid=$1;
397 :     #
398 :     # tRNA is a valid SOFA term == SO:0000253
399 :     #
400 :     my ($id, $type)=("rna.$geneid", "tRNA");
401 :     if ($outputfasta)
402 :     {
403 : olson 1.4 my $addseq = $fig->dna_seq($genome, @location);
404 : olson 1.1 $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
405 :     $fastasequences .= ">$id\n$addseq\n";
406 :     }
407 : olson 1.2 push (@{$contig_data->{$contig_key}}, (join "\t", ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$id;$allnotes")));
408 : olson 1.1 } # end the if type == rna
409 :     else
410 :     {
411 :     die "Don't know what type: |$type| is";
412 :     }
413 :     }
414 :     return ($contig_data, $fastasequences);
415 :     }
416 :    
417 :     =head1 write_gff3
418 :    
419 :     Write a set of gff3 per-contig data and fasta sequence data to a file or filehandle.
420 :    
421 :     $genome is the genome these contigs are a part of.
422 :     $contig_list is a list of contig-data hashes as returned by gff_for_feature.
423 :     $fast_list is a list of fasta data strings.
424 :    
425 :     =cut
426 :    
427 :     sub write_gff3
428 :     {
429 :     my($self, $output, $genome, $contig_list, $fasta_list) = @_;
430 :    
431 :     my $fig = $self->{fig};
432 :    
433 :     my $len = $self->{contig_length_cache};
434 :     my $fh;
435 :    
436 :     my $beg = $self->{options}->{beg};
437 :     my $end = $self->{options}->{end};
438 :    
439 :     my $close_output;
440 :    
441 :     if (ref($output))
442 :     {
443 :     $fh = $output;
444 :     }
445 :     else
446 :     {
447 :     open($fh, ">$output") or confess "Cannot open output '$output': $!";
448 :     $close_output = 1;
449 :     }
450 :    
451 :     #
452 :     # Build a data structure from the list of contigs
453 :     # that has a list of lists of data per contig name.
454 :     # (Do this so we don't copy all of the contig data itself, as it
455 :     # could be quite large).
456 :     #
457 :     my %contigs;
458 :    
459 :     #
460 :     # iterate over the given list of contig hashes.
461 :     #
462 :     for my $chash (@$contig_list)
463 :     {
464 :     #
465 :     # Then for each contig in the individual contig hashes,
466 :     # add the data list to %contigs.
467 :     #
468 :     for my $contig (keys %$chash)
469 :     {
470 :     push(@{$contigs{$contig}}, $chash->{$contig});
471 :     }
472 :     }
473 :    
474 :     foreach my $contig (sort keys %contigs)
475 :     {
476 :     print $fh "##sequence-region\t$contig\t";
477 :     if (defined $beg) {
478 :     print $fh "$beg\t";
479 :     } else {
480 :     print $fh "1\t";
481 :     }
482 :     if (defined $end) {
483 :     print $fh "$end\n";
484 :     } else {
485 :     print $fh "$len->{$contig}\n";
486 :     }
487 :     for my $list (@{$contigs{$contig}})
488 :     {
489 : olson 1.2 print $fh join("\n", @$list), "\n";
490 : olson 1.1 }
491 :     }
492 :    
493 :     print $fh "##FASTA\n";
494 :     # print out the cds and pro if we need them
495 :    
496 :     if ($self->{options}->{outputfasta})
497 :     {
498 :     for my $fastasequences (@$fasta_list)
499 :     {
500 :     print $fh $fastasequences;
501 :     }
502 :     }
503 :    
504 :     my $ll = $self->{options}->{linelength};
505 :     foreach my $contig (sort keys %contigs)
506 :     {
507 :     my $len=$fig->contig_ln($genome, $contig);
508 :     my $dna_seq=$fig->dna_seq($genome, $contig . "_1_". $len);
509 :     if (defined $beg)
510 :     {
511 :     unless (defined $end) {
512 :     $end=$len;
513 :     }
514 :     $dna_seq = substr($dna_seq, $beg, $end);
515 :     }
516 :     elsif (defined $end)
517 :     {
518 :     $beg=1;
519 :     $dna_seq = substr($dna_seq, $beg, $end);
520 :     }
521 :    
522 :     my $contig=uri_escape($contig);
523 :    
524 :     $dna_seq =~ s/(.{$ll})/$1\n/g;
525 :     chomp($dna_seq); # just remove the last \n if there is one
526 :     print $fh ">$contig\n$dna_seq\n";
527 :     }
528 :    
529 :     close($fh) if $close_output;
530 :     }
531 :    
532 : olson 1.2 package GFFParser;
533 :    
534 :     use strict;
535 :     use URI::Escape;
536 :     use Carp;
537 :     use Data::Dumper;
538 :    
539 :     use base qw(Class::Accessor);
540 :    
541 :     __PACKAGE__->mk_accessors(qw(fig current_file));
542 :    
543 :     my $count;
544 :    
545 :    
546 :     #
547 :     # GFF file parser. Creates GFFFiles.
548 :     #
549 :    
550 :     sub new
551 :     {
552 :     my($class, $fig) = @_;
553 :    
554 :     my $self = {
555 :     fig => $fig,
556 :     };
557 :    
558 :     return bless($self, $class);
559 :     }
560 :    
561 :     sub parse
562 :     {
563 :     my($self, $file) = @_;
564 :    
565 :     my($fh, $close_handle);
566 :    
567 :     my $fobj = GFFFile->new($self->fig);
568 :     $self->current_file($fobj);
569 :    
570 :     if (ref($file) ? (ref($file) eq 'GLOB'
571 :     || UNIVERSAL::isa($file, 'GLOB')
572 :     || UNIVERSAL::isa($file, 'IO::Handle'))
573 :     : (ref(\$file) eq 'GLOB'))
574 :     {
575 :     $fh = $file;
576 :     }
577 :     else
578 :     {
579 :     open($fh, "<$file") or confess "Cannot open $file: $!";
580 :     $fobj->filename($file);
581 :     $close_handle = 1;
582 :     }
583 :    
584 :     #
585 :     # Start parsing by verifying this is a gff3 file.
586 :     #
587 :    
588 :     $_ = <$fh>;
589 :    
590 :     if (m,^\#gff-version\t(\S+),)
591 :     {
592 :     if ($1 != 3)
593 :     {
594 :     confess "Invalid GFF File: version is not 3";
595 :     }
596 :     }
597 :    
598 :     #
599 :     # Now parse.
600 :     #
601 :    
602 :     while (<$fh>)
603 :     {
604 :     chomp;
605 :     #
606 :     # Check first for the fasta directive so we can run off and parse that
607 :     # separately.
608 :     #
609 :    
610 :     if (/^>/)
611 :     {
612 :     $self->parse_fasta($fh, $_);
613 :     last;
614 :     }
615 :     elsif (/^\#\#FASTA/)
616 :     {
617 :     print "Got fasta directive\n";
618 :     $_ = <$fh>;
619 :     chomp;
620 :     $self->parse_fasta($fh, $_);
621 :     last;
622 :     }
623 :     elsif (/^\#\s/)
624 :     {
625 :     #
626 :     # comment.
627 :     #
628 :     next;
629 :     }
630 :     elsif (/^\#\#(\S+)(?:\t(.*))?/)
631 :     {
632 :     #
633 :     # GFF3 directive.
634 :     #
635 :    
636 :     $self->parse_gff3_directive($1, $2);
637 :    
638 :     }
639 :     elsif (/^\#(\S+)(?:\t(.*))?/)
640 :     {
641 :     #
642 :     # Directive.
643 :     #
644 :    
645 : redwards 1.5 if (lc($1) eq "seed")
646 : olson 1.2 {
647 :     $self->parse_seed_directive($2);
648 :     }
649 :     else
650 :     {
651 :     $self->parse_local_directive($1, $2);
652 :     }
653 :    
654 :     }
655 :     elsif (/^([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)$/)
656 :     {
657 :     $self->parse_feature($1, $2, $3, $4, $5, $6, $7, $8, $9);
658 :     }
659 :     else
660 :     {
661 :     die "bad line: '$_'\n";
662 :     }
663 :     }
664 :    
665 :     return $fobj;
666 :     }
667 :    
668 :     sub parse_gff3_directive
669 :     {
670 :     my($self, $directive, $rest) = @_;
671 :    
672 : redwards 1.6 $directive = lc($directive);
673 :     my @rest=split /\t/, $rest;
674 :    
675 :     if ($directive eq "genome")
676 :     {
677 :     $self->current_file->genome_id($rest[0]);
678 :     $self->current_file->genome_name($rest[1]);
679 :     }
680 :     elsif ($directive eq "genome_md5")
681 :     {
682 :     $self->current_file->set_genome_checksum(@rest[0,1]);
683 :     }
684 :     elsif ($directive eq "origin")
685 :     {
686 :     print STDERR "We have a directive called origin but this should be changed as it will conflict with NCBI's ORIGIN indicating beginning of the sequence\n";
687 :     print STDERR "At the moment ORIGIN is returned by \$feat->project\n";
688 :     $self->current_file->project($rest[0]);
689 :     }
690 :     elsif ($directive eq "taxonomy")
691 :     {
692 :     $self->current_file->taxonomy($rest);
693 :     }
694 :     elsif ($directive eq "sequence-region")
695 :     {
696 :     $self->{contig_length_cache}->{$rest[0]}=$rest[2]-$rest[1];
697 :     $self->{contig_start_cache}->{$rest[0]}=$rest[1];
698 :     $self->{contig_end_cache}->{$rest[0]}=$rest[2];
699 :     }
700 :     else
701 :     {
702 :     print "Have gff3 directive '$directive' rest='$rest'\n";
703 :     }
704 :    
705 : olson 1.2 }
706 :    
707 :     sub parse_seed_directive
708 :     {
709 :     my($self, $rest) = @_;
710 :    
711 :     my($verb, @rest) = split(/\t/, $rest);
712 :    
713 : redwards 1.6 # are we case sensitive? I don't think so
714 :     $verb-lc($verb);
715 :    
716 : olson 1.2 if ($verb eq "anno_start")
717 :     {
718 :     $self->current_file->anno_start($rest[0]);
719 :     }
720 :     elsif ($verb eq "anno_end")
721 :     {
722 :     $self->current_file->anno_start($rest[0]);
723 :     }
724 :     elsif ($verb eq "contig_md5")
725 :     {
726 :     $self->current_file->set_contig_checksum(@rest[0,1,2]);
727 :     }
728 :     }
729 :    
730 :     sub parse_local_directive
731 :     {
732 :     my($self, $directive, $rest) = @_;
733 :    
734 :     print "Have local directive '$directive' rest='$rest'\n";
735 :     }
736 :    
737 :     sub parse_feature
738 :     {
739 :     my($self, $seqid, $source, $type, $start, $end, $score, $strand, $phase, $attributes) = @_;
740 :    
741 :     #print "data: seqid=$seqid source=$source type=$type start=$start end=$end\n";
742 :     #print " score=$score strand=$strand phase=$phase\n";
743 :     #print " $attributes\n";
744 :    
745 :     #
746 :     # Parse this feature, creating a GFFFeature object for it.
747 :     #
748 :    
749 :     my $feature = GFFFeature->new($self->fig);
750 :    
751 :     $feature->seqid($seqid);
752 :     $feature->source($source);
753 :     $feature->type($type);
754 :     $feature->start($start);
755 :     $feature->end($end);
756 :     $feature->score($score);
757 :     $feature->strand($strand);
758 :     $feature->phase($phase);
759 :    
760 :     my $athash = {};
761 :    
762 :     for my $attr (split(/;/, $attributes))
763 :     {
764 :     my($name, $value) = split(/=/, $attr);
765 :    
766 : olson 1.4 my @values = map { uri_unescape($_) } split(/,/, $value);
767 :    
768 : olson 1.2 #
769 : olson 1.4 # This might be a little goofy for the users, but we will use it
770 :     # for now:
771 :     #
772 :     # if there is more than one value, the value is a ref to a list
773 :     # of the values.
774 :     #
775 :     # Otherwise, the value is a scalar.
776 : olson 1.2 #
777 :    
778 :     if (@values > 1)
779 :     {
780 : olson 1.4 #
781 :     # Yes, you can do this ... I had to look it up :-).
782 :     #
783 :     # It's in 'man perlfaq3'.
784 :     #
785 :    
786 :     $value = \@values;
787 : olson 1.2 }
788 :    
789 :     $athash->{$name} = $value;
790 :    
791 : olson 1.4 #
792 :     # Handle the GFF3-defined attributes.
793 :     #
794 :     # These show up as Class::Accessor's on the feature object.
795 :     #
796 :    
797 : olson 1.2 if ($GFFFeature::GFF_Tags{$name})
798 :     {
799 :     $feature->set($name, $value);
800 : olson 1.4
801 :     if ($name eq "Dbxref")
802 :     {
803 :     #
804 :     # If we have a SEED:figid DBxref, set the genome and fig_id attributes.
805 :     #
806 :    
807 :     my @seed_xref = grep /^"?SEED:/, @values;
808 :     if (@seed_xref and $seed_xref[0] =~ /^"?SEED:(fig\|(\d+\.\d+)\..*)/)
809 :     {
810 :     $feature->genome($2);
811 :     $feature->fig_id($1);
812 :     }
813 :    
814 :     }
815 : olson 1.2 }
816 :     }
817 :     $feature->attributes($athash);
818 :    
819 : olson 1.4
820 : olson 1.2 $self->current_file->add_feature($feature);
821 :     }
822 :    
823 :     #
824 :     # We come in here with the first line of the fasta already read
825 :     # in order to support the backward-compatiblity syntax that
826 :     # lets a file skip the ##FASTA directive if it wishes.
827 :     #
828 :     sub parse_fasta
829 :     {
830 :     my($self, $fh, $first_line) = @_;
831 :     my($cur, $cur_id);
832 :    
833 :     for ($_ = $first_line; $_; $_ = <$fh>, chomp)
834 :     {
835 :     if (/^>\s*(\S+)/)
836 :     {
837 :     if ($cur)
838 :     {
839 :     $self->handle_fasta_block($cur_id, $cur);
840 :     }
841 :    
842 :     $cur = '';
843 :     $cur_id = $1;
844 :     }
845 :     else
846 :     {
847 :     s/^\s*$//;
848 :     s/\s*$//;
849 :     if (/\s/)
850 :     {
851 :     die "FASTA data had embedded space: $_\n";
852 :     }
853 :     $cur .= $_;
854 :     }
855 :     }
856 :     if ($cur)
857 :     {
858 :     $self->handle_fasta_block($cur_id, $cur);
859 :     }
860 :     }
861 :    
862 :     sub handle_fasta_block
863 :     {
864 :     my($self, $id, $data) = @_;
865 :    
866 :     my $len = length($data);
867 :     $self->current_file->set_fasta_data($id, $data);
868 :     }
869 :    
870 :     package GFFFeature;
871 :    
872 :     use strict;
873 :     use base qw(Class::Accessor);
874 :    
875 :     our @GFF_Tags = qw(ID Name Alias Parent Target Gap Note Dbxref Ontology_term);
876 :     our %GFF_Tags;
877 :    
878 :     map { $GFF_Tags{$_} = 1 } @GFF_Tags;
879 :    
880 : olson 1.4 __PACKAGE__->mk_accessors(qw(fig seqid source type start end score strand phase attributes genome fig_id),
881 : olson 1.2 @GFF_Tags);
882 :    
883 :    
884 :     sub new
885 :     {
886 :     my($class, $fig) = @_;
887 :    
888 :     my $self = {
889 :     fig => $fig,
890 :     };
891 :    
892 :     return bless($self, $class);
893 :     }
894 :    
895 : olson 1.4 sub find_local_feature
896 :     {
897 :     my($self, $local_genome) = @_;
898 :     my $db = $self->fig->db_handle;
899 :    
900 :     # For debugging.
901 :     undef $local_genome;
902 :     if ($local_genome)
903 :     {
904 :     #
905 :     # It's a precise match. We need to determine if we have this
906 :     # particular feature in this SEED (it is possible for one to
907 :     # have exported an annotation for a feature that was added
908 :     # to a genome after its initial release).
909 :     #
910 :     # We do this by searching for a local feature with the same contig,
911 :     # start, and stop as this feature.
912 :     #
913 :    
914 :     my $qry = qq(SELECT id
915 :     FROM features
916 :     WHERE (genome = ? AND
917 :     contig = ? AND
918 :     minloc = ? AND
919 :     maxloc = ?));
920 :     my $res = $db->SQL($qry, undef, $local_genome, $self->seqid,
921 :     $self->start, $self->end);
922 :    
923 :     return map { $_->[0] } @$res;
924 :     }
925 :    
926 :     #
927 :     # Otherwise, we need to try a set of heuristics to match
928 :     # this id.
929 :     #
930 :    
931 :     #
932 :     # Try matching aliases first.
933 :     #
934 :    
935 :     my @aliases = grep { !/^\"?SEED/ } ref($self->Dbxref) ? @{$self->Dbxref} : ($self->Dbxref);
936 :    
937 :     my @maliases = map { FigGFF::map_dbxref_to_seed_alias($_) } @aliases;
938 :    
939 :     print "Found aliases @aliases\n";
940 :     print "Found mapped aliases @maliases\n";
941 :    
942 :     for my $malias (@maliases)
943 :     {
944 :     my $fid = $self->fig->by_alias($malias);
945 :     if ($fid)
946 :     {
947 :     print "Mapped $malias to $fid\n";
948 :     }
949 :     }
950 :    
951 :     }
952 :    
953 :    
954 : olson 1.2 package GFFFile;
955 :    
956 :     use strict;
957 :     use base qw(Class::Accessor);
958 :    
959 :     __PACKAGE__->mk_accessors(qw(fig filename features feature_index anno_start anno_end));
960 :    
961 :     #
962 :     # Package to hold the contents of a GFF file, and to hold the code
963 :     # for mapping its contents to the local SEED.
964 :     #
965 :     # Created by GFFParser->parse.
966 :     #
967 :    
968 :     sub new
969 :     {
970 :     my($class, $fig) = @_;
971 :    
972 :     my $self = {
973 :     fig => $fig,
974 :     features => [],
975 : redwards 1.6 contigs => [],
976 : olson 1.2 feature_index => {},
977 : olson 1.4 genome_checksum => {},
978 :     contig_checksum => {},
979 :     features_by_genome => {},
980 : olson 1.2 };
981 :     return bless($self, $class);
982 :     }
983 :    
984 :     sub add_feature
985 :     {
986 :     my($self, $feature) = @_;
987 :    
988 :     push(@{$self->features}, $feature);
989 :     $self->feature_index->{$feature->ID} = $feature;
990 : olson 1.4 push(@{$self->{features_by_genome}->{$feature->genome}}, $feature);
991 :     }
992 :    
993 :     sub features_for_genome
994 :     {
995 :     my($self, $genome) = @_;
996 :    
997 :     return $self->{features_by_genome}->{$genome};
998 :     }
999 :    
1000 :     sub genome_checksum
1001 :     {
1002 :     my($self, $genome) = @_;
1003 :    
1004 :     return $self->{genome_checksum}->{$genome};
1005 :     }
1006 :    
1007 :     sub genome_checksums
1008 :     {
1009 :     my($self) = @_;
1010 :    
1011 :     return [map { [$_, $self->{genome_checksum}->{$_}] } keys(%{$self->{genome_checksum}})];
1012 : olson 1.2 }
1013 :    
1014 :     sub set_genome_checksum
1015 :     {
1016 :     my($self, $genome, $md5sum) = @_;
1017 :     $self->{genome_checksum}->{$genome} = $md5sum;
1018 :     }
1019 :    
1020 :     sub set_contig_checksum
1021 :     {
1022 :     my($self, $genome, $contig, $md5sum) = @_;
1023 :     $self->{contig_checksum}->{$genome}->{$contig} = $md5sum;
1024 :     }
1025 :    
1026 : redwards 1.6
1027 : olson 1.2 sub set_fasta_data
1028 :     {
1029 :     my($self, $id, $data) = @_;
1030 :    
1031 :     $self->{fasta_data}->{$id} = $data;
1032 :     }
1033 :    
1034 :    
1035 : redwards 1.6 =head2 contigs()
1036 :    
1037 :     Add a contig to the list, or return a reference to an array of contigs
1038 :    
1039 :     =cut
1040 :    
1041 :     sub contigs
1042 :     {
1043 :     my($self, $contig) = @_;
1044 :     $contig && (push @{$self->{contigs}}, $contig);
1045 :     if (!$self->{contigs} && $self->{contig_length_cache}) {$self->{contigs} = keys %{$self->{contig_length_cache}}}
1046 :     return $self->{contigs};
1047 :     }
1048 :    
1049 :     =head2 contig_length()
1050 :    
1051 :     Get or set the length of a specfic contig.
1052 :     my $length=$fob->contig_length($contig, $length);
1053 :     my $length=$fob->contig_length($contig);
1054 :    
1055 :     =cut
1056 :    
1057 :     sub contig_length
1058 :     {
1059 :     my($self, $contig, $length) = @_;
1060 :     $length && ($self->{contig_length_cache}->{$contig}=$length);
1061 :     return $self->{contig_length_cache}->{$contig};
1062 :     }
1063 :    
1064 :     =head1 Information about the source of the sequence.
1065 :    
1066 :     These are things that we have parsed out the GFF3 file, or want to add into the GFF3 file. We can use these methods to get or set them as required. In general, if a value is supplied that will be used as the new value.
1067 :    
1068 :     =cut
1069 :    
1070 :     =head2 genome_id()
1071 :    
1072 :     Get or set a genome id for this file.
1073 :    
1074 :     =cut
1075 :    
1076 :     sub genome_id
1077 :     {
1078 :     my($self, $genomeid) = @_;
1079 :     $genomeid && ($self->{genome_id}=$genomeid);
1080 :     return $self->{genome_id};
1081 :     }
1082 :    
1083 :     =head2 genome_name()
1084 :    
1085 :     Get or set a genome id for this file.
1086 :    
1087 :     =cut
1088 :    
1089 :     sub genome_name
1090 :     {
1091 :     my($self, $genomename) = @_;
1092 :     $genomename && ($self->{genome_name}=$genomename);
1093 :     return $self->{genome_name};
1094 :     }
1095 :    
1096 :     =head2 project()
1097 :    
1098 :     Get or set the project.
1099 :    
1100 :     =cut
1101 :    
1102 :     sub project
1103 :     {
1104 :     my ($self, $pro) = @_;
1105 :     $pro && ($self->{project}=$pro);
1106 :     return $self->{project};
1107 :     }
1108 :    
1109 :     =head2 taxonomy()
1110 :    
1111 :     Get or set the taxonomy
1112 :    
1113 :     =cut
1114 :    
1115 :     sub taxonomy
1116 :     {
1117 :     my($self, $tax) = @_;
1118 :     $tax && ($self->{taxonomy}=$tax);
1119 :     return $self->{taxonomy};
1120 :     }
1121 :    
1122 :    
1123 :    
1124 :    
1125 :    
1126 : olson 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3