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

Diff of /FigKernelPackages/FigGFF.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Mon Mar 14 19:17:37 2005 UTC revision 1.6, Sat Apr 9 17:21:15 2005 UTC
# Line 6  Line 6 
6  # A GFFWriter handles the generation of GFF files from SEED data structures.  # A GFFWriter handles the generation of GFF files from SEED data structures.
7  #  #
8    
9    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    
115  package GFFWriter;  package GFFWriter;
116    
117  use strict;  use strict;
118  use FIG;  use FIG;
119    use FigGFF;
120    
121  use Carp;  use Carp;
122  use URI::Escape;  use URI::Escape;
# Line 28  Line 134 
134    
135      map { $default_options->{$_} = $options{$_} } keys(%options);      map { $default_options->{$_} = $options{$_} } keys(%options);
136    
137        # 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      my $self = {      my $self = {
143          options => $default_options,          options => $default_options,
144          contig_length_cache => {},          contig_length_cache => {},
145            contig_start_cache  => {},
146            contig_end_cache    => {},
147          fig => $fig,          fig => $fig,
148      };      };
149    
# Line 52  Line 165 
165    
166  sub gff3_for_feature  sub gff3_for_feature
167  {  {
168      my($self, $fid, $user) = @_;      my($self, $fid, $user, $user_note, $in_aliases, $in_loc) = @_;
169    
170      #      #
171      # Options we need to figure out somehow.      # Options we need to figure out somehow.
# Line 61  Line 174 
174    
175      my $escapespace = $options->{escapespace};      my $escapespace = $options->{escapespace};
176      my $outputfasta = $options->{outputfasta};      my $outputfasta = $options->{outputfasta};
177      my %outputtype = (cds => 1,  
178                        protein => 1,      my %outputtype;
179                        transcript => 1,      map { $outputtype{$_} = 1 } @{$options->{outputtype}};
                       gene => 1);  
180    
181      my $fastasequences = '';      my $fastasequences = '';
182      my $contig_data;      my $contig_data;
183      my $linelength = $options->{linelength};      my $linelength = $options->{linelength};
184    
185        my $beg = $self->{options}->{beg};
186        my $end = $self->{options}->{end};
187    
188      my $fig = $self->{fig};      my $fig = $self->{fig};
189    
190      #      #
191      # Do this first to make sure that we really have a feature.      # Do this first to make sure that we really have a feature.
192      #      #
193      my @location = $fig->feature_location($fid);      my @location = ref($in_loc) ? @$in_loc : $fig->feature_location($fid);
     print Dumper(\@location);  
194      if (@location == 0 or !defined($location[0]))      if (@location == 0 or !defined($location[0]))
195      {      {
196          warn "No location found for feature $fid\n";          warn "No location found for feature $fid\n";
# Line 99  Line 213 
213      # all the aliases we are going to use      # all the aliases we are going to use
214      #      #
215      my @alias;      my @alias;
216        my @xref;
217    
218        if ($options->{with_assignments})
219        {
220      my $func = $fig->function_of($fid, $user);      my $func = $fig->function_of($fid, $user);
221      if ($func)      if ($func)
222      {      {
223          push @$note, ("Note=\"" . uri_escape($func) . '"');              push @$note, ("Note=" . uri_escape($func));
224      }      }
   
     # now find aliases  
     foreach my $alias ($fig->feature_aliases($fid))  
     {  
         if ($alias =~ /^NP/)  
         {  
             push @$note, "Dbxref=\"NCBI_genpept:$alias\"";  
225          }          }
226          elsif ($alias =~ /gi\|/)  
227          {      if ($options->{with_aliases})
             $alias =~ s/^gi\|//;  
             push @$note, "Dbxref=\"NCBI_gi:$alias\"";  
         }  
         elsif ($alias =~ /^kegg/i)  
228          {          {
229              $alias =~ s/kegg\|/KEGG:/i;          # now find aliases
230              $alias =~ s/^(.*):/$1+/;          my @feat_aliases = ref($in_aliases) ? @$in_aliases : $fig->feature_aliases($fid);
231              push @$note, "Dbxref=\"$alias\"";          foreach my $alias (@feat_aliases)
         }  
         elsif ($alias =~ /^uni/)  
232          {          {
233              $alias =~ s/uni\|/UniProt:/;              my $mapped = FigGFF::map_seed_alias_to_dbxref($alias);
234              # $note = check_go($note, $alias) if ($go);              if ($mapped)
             push @$note, "Dbxref=\"$alias\"";  
         }  
         elsif ($alias =~ /^sp\|/)  
235          {          {
236              $alias =~ s/sp\|/Swiss-Prot:/;                  push(@xref, $mapped);
             push @$note, "Dbxref=\"$alias\"";  
237          }          }
238          else          else
239          {          {
240              # don't know what it is so keep it as an alias                  push(@alias, $alias);
241              $alias = uri_escape($alias); # just in case!              }
             push @alias, $alias;  
242          }          }
243    }    }
244    
245      # now just join all the aliases and put them into @$note so we can add it to the array      # now just join all the aliases and put them into @$note so we can add it to the array
246      if (@alias)      if (@alias)
247      {      {
248          push @$note, "Alias=\"". join (",", @alias) . '"';          push @$note, "Alias=". join (",", map { uri_escape($_) } @alias);
249        }
250    
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    
260      # the LAST thing I am going to add as a note is the FIG id so that I can grep it out easily      # the LAST thing I am going to add as a note is the FIG id so that I can grep it out easily
261      push @$note, "Dbxref=\"SEED:$fid\"";      #
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    
269      # finally join all the notes into a long string that can be added as column 9      # finally join all the notes into a long string that can be added as column 9
270      my $allnotes;      my $allnotes;
# Line 239  Line 353 
353              # gene: SO:0000704              # gene: SO:0000704
354              # cds: SO:0000316              # cds: SO:0000316
355              # protein:  NOT VALID: should be protein_coding_primary_transcript SO:0000120              # protein:  NOT VALID: should be protein_coding_primary_transcript SO:0000120
356              foreach my $type (keys %outputtype)  
357              {              #
358                  my ($id, $type)=($trunc{$type} . "." . $geneid, $type);              # 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                  # we want to store some sequences to be output
370                  if ($outputfasta && $type eq "protein")              if ($outputfasta)
371                  {                  {
372                      my $addseq = $fig->get_translation($fid);                      my $addseq = $fig->get_translation($fid);
373    
# Line 252  Line 376 
376                      #                      #
377                      $addseq =~ s/(.{$linelength})/$1\n/g;                      $addseq =~ s/(.{$linelength})/$1\n/g;
378                      chomp($addseq);                      chomp($addseq);
379                      $fastasequences .= ">$id\n$addseq\n";                  $fastasequences .= ">$protein_id\n$addseq\n";
380                  }  
381                  if ($outputfasta && $type eq "cds")                  my $addseq = uc($fig->dna_seq($genome, @location));
                 {  
                     my $addseq = uc($fig->dna_seq($genome, $fig->feature_location($fid)));  
382                      $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);                      $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
383                      $fastasequences .= ">$id\n$addseq\n";  
384                  }                  $fastasequences .= ">$cds_id\n$addseq\n";
385                  # correct the incorrect sofa term  
386                  if ($type eq "protein")                  $allnotes .= ";translation_id=$protein_id";
                 {  
                     $type = "SO:0000120";  
387                  }                  }
388    
389                  push (@{$contig_data->{$contig_key}},                  push (@{$contig_data->{$contig_key}},
390                        (join "\t",                        (join "\t",
391                         ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "Id=$id;$allnotes")));                     ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$cds_id;$allnotes")));
392              } # end the foreach my $type          }
         } # end the if type==peg  
393          elsif ($type eq "rna")          elsif ($type eq "rna")
394          {          {
395              $fid =~ /\.rna\.(\d+)/;              $fid =~ /\.rna\.(\d+)/;
# Line 281  Line 400 
400              my ($id, $type)=("rna.$geneid", "tRNA");              my ($id, $type)=("rna.$geneid", "tRNA");
401              if ($outputfasta)              if ($outputfasta)
402              {              {
403                  my $addseq = $fig->dna_seq($genome, $fig->feature_location($fid));                  my $addseq = $fig->dna_seq($genome, @location);
404                  $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);                  $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
405                  $fastasequences .= ">$id\n$addseq\n";                  $fastasequences .= ">$id\n$addseq\n";
406              }              }
407              push (@{$contig_data->{$contig_key}}, (join "\t", ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "Id=$id;$allnotes")));              push (@{$contig_data->{$contig_key}}, (join "\t", ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$id;$allnotes")));
408          } # end the if type == rna          } # end the if type == rna
409          else          else
410          {          {
# Line 367  Line 486 
486          }          }
487          for my $list (@{$contigs{$contig}})          for my $list (@{$contigs{$contig}})
488          {          {
489              print $fh join "\n", @$list;              print $fh join("\n", @$list), "\n";
490          }          }
491      }      }
492    
# Line 410  Line 529 
529      close($fh) if $close_output;      close($fh) if $close_output;
530  }  }
531    
532    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                if (lc($1) eq "seed")
646                {
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        $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    }
706    
707    sub parse_seed_directive
708    {
709        my($self, $rest) = @_;
710    
711        my($verb, @rest) = split(/\t/, $rest);
712    
713        # are we case sensitive? I don't think so
714        $verb-lc($verb);
715    
716        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            my @values = map { uri_unescape($_) } split(/,/, $value);
767    
768            #
769            # 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            #
777    
778            if (@values > 1)
779            {
780                #
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            }
788    
789            $athash->{$name} = $value;
790    
791            #
792            # Handle the GFF3-defined attributes.
793            #
794            # These show up as Class::Accessor's on the feature object.
795            #
796    
797            if ($GFFFeature::GFF_Tags{$name})
798            {
799                $feature->set($name, $value);
800    
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            }
816        }
817        $feature->attributes($athash);
818    
819    
820        $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    __PACKAGE__->mk_accessors(qw(fig seqid source type start end score strand phase attributes genome fig_id),
881                              @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    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    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            contigs  => [],
976            feature_index => {},
977            genome_checksum => {},
978            contig_checksum => {},
979            features_by_genome => {},
980        };
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        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    }
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    
1027    sub set_fasta_data
1028    {
1029        my($self, $id, $data) = @_;
1030    
1031        $self->{fasta_data}->{$id} = $data;
1032    }
1033    
1034    
1035    =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  1;  1;

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3