[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.21, Sun Feb 5 00:45:19 2006 UTC
# Line 1  Line 1 
1  #  #
2    # Copyright (c) 2003-2006 University of Chicago and Fellowship
3    # for Interpretations of Genomes. All Rights Reserved.
4    #
5    # This file is part of the SEED Toolkit.
6    #
7    # The SEED Toolkit is free software. You can redistribute
8    # it and/or modify it under the terms of the SEED Toolkit
9    # Public License.
10    #
11    # You should have received a copy of the SEED Toolkit Public License
12    # along with this program; if not write to the University of Chicago
13    # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14    # Genomes at veronika@thefig.info or download a copy from
15    # http://www.theseed.org/LICENSE.TXT.
16    #
17    
18    #
19  # FIG GFF utilities.  # FIG GFF utilities.
20  #  #
21    
# Line 6  Line 23 
23  # A GFFWriter handles the generation of GFF files from SEED data structures.  # A GFFWriter handles the generation of GFF files from SEED data structures.
24  #  #
25    
26    package FigGFF;
27    
28    use strict;
29    
30    use base qw(Exporter);
31    use vars qw(@EXPORT);
32    @EXPORT = qw(map_seed_alias_to_dbxref map_dbxref_to_seed_alias);
33    
34    #
35    # General GFF-related routines.
36    #
37    
38    
39    #
40    # Alias translation.
41    #
42    # These routines map between the SEED aliases and the standard
43    # dbxref names as defined here:
44    #
45    # ftp://ftp.geneontology.org/pub/go/doc/GO.xrf_abbs
46    #
47    # In particular:
48    #
49    # abbreviation: NCBI_gi
50    # database: NCBI databases.
51    # object: Identifier.
52    # example_id: NCBI_gi:10727410
53    # generic_url: http://www.ncbi.nlm.nih.gov/
54    # url_syntax:
55    # url_example: http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=nucleotide&list_uids=10727410&dopt=GenBank
56    #
57    # abbreviation: NCBI_NP
58    # database: NCBI RefSeq.
59    # object: Protein identifier.
60    # example_id: NCBI_NP:123456
61    # generic_url: http://www.ncbi.nlm.nih.gov/
62    # url_syntax:
63    #
64    # abbreviation: KEGG
65    # database: Kyoto Encyclopedia of Genes and Genomes.
66    # generic_url: http://www.genome.ad.jp/kegg/
67    
68    sub map_seed_alias_to_dbxref
69    {
70        my($alias) = @_;
71    
72        $_ = $alias;
73        if (/^NP_(\d+.*)/)
74        {
75            return "NCBI_NP:$1";
76        }
77        elsif (/^NM_(\d+.*)/)
78        {
79            return "NCBI_NM:$1";
80        }
81        elsif (/^gi\|(\d+)/)
82        {
83            return "NCBI_gi:$1";
84        }
85        elsif (/^kegg\|(\S+):(\S+)/)
86        {
87            return "KEGG:$1 $2";
88        }
89        elsif (/^uni\|(\S+)/)
90        {
91            return "UniProt:$1";
92        }
93        elsif (/^sp\|(\S+)/)
94        {
95            return "Swiss-Prot:$1";
96        }
97    
98        return undef;
99    }
100    
101    #
102    # And map back again.
103    #
104    
105    sub map_dbxref_to_seed_alias
106    {
107        my($dbxref) = @_;
108    
109        # if it is not a valid xref just return it
110        return $dbxref unless $dbxref =~ m/:/;
111    
112        my($type, $ref) = split(/:/, $dbxref, 2);
113    
114        if (lc($type) eq "ncbi_np")
115        {
116            return "$ref";
117        }
118        elsif (lc($type) eq "ncbi_nm")
119        {
120            return "$ref";
121        }
122        elsif (lc($type) eq "ncbi_pid")
123        {
124            return "$ref";
125        }
126        elsif (lc($type) eq "ncbi_gi")
127        {
128            return "gi|$ref";
129        }
130        elsif (lc($type) eq "kegg")
131        {
132            $ref =~ s/ /:/;
133            return "kegg|$ref";
134        }
135        elsif (lc($type) eq "uniprot")
136        {
137            return "uni|$ref";
138        }
139        elsif (lc($type) eq "swiss-prot")
140        {
141            return "sp|$ref";
142        }
143    
144        return $dbxref; # just return itself if we don't know what it is.
145    }
146    
147  package GFFWriter;  package GFFWriter;
148    
149  use strict;  use strict;
150  use FIG;  use FIG;
151    use FigGFF;
152    
153  use Carp;  use Carp;
154  use URI::Escape;  use URI::Escape;
# Line 28  Line 166 
166    
167      map { $default_options->{$_} = $options{$_} } keys(%options);      map { $default_options->{$_} = $options{$_} } keys(%options);
168    
169        # added contig_start_cache and contig_end_cache because we have something like
170        # sequence-region   contigname      1       10000
171        # in general we will set contig_start_cache == 1
172        # and contig_end_cache == contig_length_cache
173    
174      my $self = {      my $self = {
175          options => $default_options,          options => $default_options,
176          contig_length_cache => {},          contig_length_cache => {},
177            contig_start_cache  => {},
178            contig_end_cache    => {},
179          fig => $fig,          fig => $fig,
180      };      };
181    
# Line 52  Line 197 
197    
198  sub gff3_for_feature  sub gff3_for_feature
199  {  {
200      my($self, $fid, $user) = @_;      my($self, $fid, $user, $user_note, $in_aliases, $in_loc) = @_;
201    
202      #      #
203      # Options we need to figure out somehow.      # Options we need to figure out somehow.
# Line 61  Line 206 
206    
207      my $escapespace = $options->{escapespace};      my $escapespace = $options->{escapespace};
208      my $outputfasta = $options->{outputfasta};      my $outputfasta = $options->{outputfasta};
209      my %outputtype = (cds => 1,  
210                        protein => 1,      my %outputtype;
211                        transcript => 1,      map { $outputtype{$_} = 1 } @{$options->{outputtype}};
                       gene => 1);  
212    
213      my $fastasequences = '';      my $fastasequences = '';
214      my $contig_data;      my $contig_data;
215      my $linelength = $options->{linelength};      my $linelength = $options->{linelength};
216    
217        my $beg = $self->{options}->{beg};
218        my $end = $self->{options}->{end};
219    
220      my $fig = $self->{fig};      my $fig = $self->{fig};
221    
222      #      #
223      # Do this first to make sure that we really have a feature.      # Do this first to make sure that we really have a feature.
224      #      #
225      my @location = $fig->feature_location($fid);      my @location = ref($in_loc) ? @$in_loc : $fig->feature_location($fid);
     print Dumper(\@location);  
226      if (@location == 0 or !defined($location[0]))      if (@location == 0 or !defined($location[0]))
227      {      {
228          warn "No location found for feature $fid\n";          warn "No location found for feature $fid\n";
# Line 99  Line 245 
245      # all the aliases we are going to use      # all the aliases we are going to use
246      #      #
247      my @alias;      my @alias;
248        my @xref;
249    
250        if ($options->{with_assignments})
251        {
252      my $func = $fig->function_of($fid, $user);      my $func = $fig->function_of($fid, $user);
253      if ($func)      if ($func)
254      {      {
255          push @$note, ("Note=\"" . uri_escape($func) . '"');              push @$note, ("Note=" . uri_escape($func));
     }  
   
     # now find aliases  
     foreach my $alias ($fig->feature_aliases($fid))  
     {  
         if ($alias =~ /^NP/)  
         {  
             push @$note, "Dbxref=\"NCBI_genpept:$alias\"";  
256          }          }
         elsif ($alias =~ /gi\|/)  
         {  
             $alias =~ s/^gi\|//;  
             push @$note, "Dbxref=\"NCBI_gi:$alias\"";  
257          }          }
258          elsif ($alias =~ /^kegg/i)  
259        if ($options->{with_aliases})
260          {          {
261              $alias =~ s/kegg\|/KEGG:/i;          # now find aliases
262              $alias =~ s/^(.*):/$1+/;          my @feat_aliases = ref($in_aliases) ? @$in_aliases : $fig->feature_aliases($fid);
263              push @$note, "Dbxref=\"$alias\"";          foreach my $alias (@feat_aliases)
         }  
         elsif ($alias =~ /^uni/)  
264          {          {
265              $alias =~ s/uni\|/UniProt:/;              my $mapped = FigGFF::map_seed_alias_to_dbxref($alias);
266              # $note = check_go($note, $alias) if ($go);              if ($mapped)
             push @$note, "Dbxref=\"$alias\"";  
         }  
         elsif ($alias =~ /^sp\|/)  
267          {          {
268              $alias =~ s/sp\|/Swiss-Prot:/;                  push(@xref, $mapped);
             push @$note, "Dbxref=\"$alias\"";  
269          }          }
270          else          else
271          {          {
272              # don't know what it is so keep it as an alias                  push(@alias, $alias);
273              $alias = uri_escape($alias); # just in case!              }
             push @alias, $alias;  
274          }          }
275    }    }
276    
277      # 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
278      if (@alias)      if (@alias)
279      {      {
280          push @$note, "Alias=\"". join (",", @alias) . '"';          push @$note, "Alias=". join (",", map { uri_escape($_) } @alias);
281        }
282    
283        #
284        # If we have user note passed in, add it.
285        #
286    
287        if ($user_note)
288        {
289            push @$note, $user_note;
290      }      }
291    
292      # 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
293      push @$note, "Dbxref=\"SEED:$fid\"";      #
294        # for now, make SEED xref the first in the list so we can search for DBxref=SEEd
295        #
296    
297        unshift(@xref, "SEED:$fid");
298    
299        push(@$note, "Dbxref=" . join(",", map { uri_escape($_) } @xref));
300    
301      # 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
302      my $allnotes;      my $allnotes;
# Line 239  Line 385 
385              # gene: SO:0000704              # gene: SO:0000704
386              # cds: SO:0000316              # cds: SO:0000316
387              # protein:  NOT VALID: should be protein_coding_primary_transcript SO:0000120              # protein:  NOT VALID: should be protein_coding_primary_transcript SO:0000120
388              foreach my $type (keys %outputtype)  
389              {              #
390                  my ($id, $type)=($trunc{$type} . "." . $geneid, $type);              # For now, we will only output CDS features, and include
391                # the translation as an attribute (attribuute name is
392                # translation_id, value is a key that corresponds to a FASTA
393                # section at the end of the file).
394                #
395    
396                my $type = "cds";
397    
398                my $protein_id = "pro.$geneid";
399                my $cds_id = "cds.$geneid";
400    
401                  # we want to store some sequences to be output                  # we want to store some sequences to be output
402                  if ($outputfasta && $type eq "protein")              if ($outputfasta)
403                  {                  {
404                      my $addseq = $fig->get_translation($fid);                      my $addseq = $fig->get_translation($fid);
405    
# Line 252  Line 408 
408                      #                      #
409                      $addseq =~ s/(.{$linelength})/$1\n/g;                      $addseq =~ s/(.{$linelength})/$1\n/g;
410                      chomp($addseq);                      chomp($addseq);
411                      $fastasequences .= ">$id\n$addseq\n";                  $fastasequences .= ">$protein_id\n$addseq\n";
412                  }  
413                  if ($outputfasta && $type eq "cds")                  $addseq = uc($fig->dna_seq($genome, @location));
                 {  
                     my $addseq = uc($fig->dna_seq($genome, $fig->feature_location($fid)));  
414                      $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);                      $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
415                      $fastasequences .= ">$id\n$addseq\n";  
416                  }                  $fastasequences .= ">$cds_id\n$addseq\n";
417                  # correct the incorrect sofa term  
418                  if ($type eq "protein")                  $allnotes .= ";translation_id=$protein_id";
                 {  
                     $type = "SO:0000120";  
419                  }                  }
420    
421                  push (@{$contig_data->{$contig_key}},                  push (@{$contig_data->{$contig_key}},
422                        (join "\t",                        (join "\t",
423                         ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "Id=$id;$allnotes")));                     ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$cds_id;$allnotes")));
424              } # end the foreach my $type          }
         } # end the if type==peg  
425          elsif ($type eq "rna")          elsif ($type eq "rna")
426          {          {
427              $fid =~ /\.rna\.(\d+)/;              $fid =~ /\.rna\.(\d+)/;
# Line 281  Line 432 
432              my ($id, $type)=("rna.$geneid", "tRNA");              my ($id, $type)=("rna.$geneid", "tRNA");
433              if ($outputfasta)              if ($outputfasta)
434              {              {
435                  my $addseq = $fig->dna_seq($genome, $fig->feature_location($fid));                  my $addseq = $fig->dna_seq($genome, @location);
436                  $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);                  $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
437                  $fastasequences .= ">$id\n$addseq\n";                  $fastasequences .= ">$id\n$addseq\n";
438              }              }
439              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")));
440          } # end the if type == rna          } # end the if type == rna
441          else          else
442          {          {
# Line 367  Line 518 
518          }          }
519          for my $list (@{$contigs{$contig}})          for my $list (@{$contigs{$contig}})
520          {          {
521              print $fh join "\n", @$list;              print $fh join("\n", @$list), "\n";
522          }          }
523      }      }
524    
# Line 410  Line 561 
561      close($fh) if $close_output;      close($fh) if $close_output;
562  }  }
563    
564    package GFFParser;
565    
566    use strict;
567    use URI::Escape;
568    use Carp;
569    use Data::Dumper;
570    
571    use base qw(Class::Accessor);
572    
573    __PACKAGE__->mk_accessors(qw(fig current_file));
574    
575    my $count;
576    
577    
578    #
579    # GFF file parser. Creates GFFFiles.
580    #
581    
582    sub new
583    {
584        my($class, $fig) = @_;
585    
586        my $self = {
587            fig => $fig,
588        };
589    
590        return bless($self, $class);
591    }
592    
593    sub parse
594    {
595        my($self, $file) = @_;
596    
597        my($fh, $close_handle);
598    
599        my $fobj = GFFFile->new($self->fig);
600        $self->current_file($fobj);
601    
602        if (ref($file) ? (ref($file) eq 'GLOB'
603                           || UNIVERSAL::isa($file, 'GLOB')
604                           || UNIVERSAL::isa($file, 'IO::Handle'))
605            : (ref(\$file) eq 'GLOB'))
606        {
607            $fh = $file;
608        }
609        else
610        {
611            open($fh, "<$file") or confess "Cannot open $file: $!";
612            $fobj->filename($file);
613            $close_handle = 1;
614        }
615    
616        #
617        # Start parsing by verifying this is a gff3 file.
618        #
619    
620        $_ = <$fh>;
621    
622        if (m,^\#gff-version\t(\S+),)
623        {
624            if ($1 != 3)
625            {
626                confess "Invalid GFF File: version is not 3";
627            }
628        }
629    
630        #
631        # Now parse.
632        #
633    
634        while (<$fh>)
635        {
636            chomp;
637            next unless ($_); # ignore empty lines
638            #
639            # Check first for the fasta directive so we can run off and parse that
640            # separately.
641            #
642    
643    
644            if (/^>/)
645            {
646                $self->parse_fasta($fh, $_);
647                last;
648            }
649            elsif (/^\#\#FASTA/)
650            {
651                # print "Got fasta directive\n";
652                $_ = <$fh>;
653                chomp;
654                $self->parse_fasta($fh, $_);
655                last;
656            }
657            elsif (/^\#\s/)
658            {
659                #
660                # comment.
661                #
662                next;
663            }
664            elsif (/^\#$/)
665            {
666                #
667                # blank line starting with #
668                #
669                next;
670            }
671            elsif (/^\#\#(\S+)(?:\t(.*))?/)
672            {
673                #
674                # GFF3 directive.
675                #
676    
677                $self->parse_gff3_directive($1, $2);
678    
679            }
680            elsif (/^\#(\S+)(?:\t(.*))?/)
681            {
682                #
683                # Directive.
684                #
685    
686                if (lc($1) eq "seed")
687                {
688                    $self->parse_seed_directive($2);
689                }
690                else
691                {
692                    $self->parse_local_directive($1, $2);
693                }
694    
695            }
696            elsif (/^([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)$/)
697            {
698                $self->parse_feature($1, $2, $3, $4, $5, $6, $7, $8, $9);
699            }
700            else
701            {
702                die "bad line: '$_'\n";
703            }
704        }
705    
706        return $fobj;
707    }
708    
709    sub parse_gff3_directive
710    {
711        my($self, $directive, $rest) = @_;
712    
713        $directive = lc($directive);
714        # this should catch both #seed and ##seed :-)
715        if ($directive eq "seed")
716        {
717          return $self->parse_seed_directive($rest);
718        }
719    
720        my @rest=split /\t/, $rest;
721    
722        # removed genome, genome_md5, origin, taxnomy as they are not real gff directives. These are in seed_directives below
723        if ($directive eq "project")
724        {
725            # I am not sure if PROJECT is a seed directive or a GFF directive
726            $self->current_file->project($rest[0]);
727        }
728        elsif ($directive eq "sequence-region")
729        {
730            $self->current_file->contigs($rest[0]);
731            $self->{contig_length_cache}->{$rest[0]}=$rest[2]-$rest[1];
732            $self->{contig_start_cache}->{$rest[0]}=$rest[1];
733            $self->{contig_end_cache}->{$rest[0]}=$rest[2];
734        }
735        else
736        {
737            print STDERR "Have gff3 directive '$directive' rest='$rest'\n";
738        }
739    
740    }
741    
742    sub parse_seed_directive
743    {
744        my($self, $rest) = @_;
745    
746        my($verb, @rest) = split(/\t/, $rest);
747    
748        # are we case sensitive? I don't think so
749        $verb=lc($verb);
750    
751        if ($verb eq "genome_id")
752        {
753            $self->current_file->genome_id($rest[0]);
754        }
755        elsif ($verb eq "name")
756        {
757            $self->current_file->genome_name($rest[0]);
758        }
759        elsif ($verb eq "genome_md5")
760        {
761            $self->current_file->set_genome_checksum($rest[0]);
762        }
763        elsif ($verb eq "project")
764        {
765            # I am not sure if PROJECT is a seed directive or a GFF directive
766            $self->current_file->project($rest[0]);
767        }
768        elsif ($verb eq "taxonomy")
769        {
770            $self->current_file->taxonomy($rest[0]);
771        }
772        elsif ($verb eq "taxon_id")
773        {
774            $self->current_file->taxon_id($rest[0]);
775        }
776        elsif ($verb eq "anno_start")
777        {
778            $self->current_file->anno_start($rest[0]);
779        }
780        elsif ($verb eq "anno_end")
781        {
782            $self->current_file->anno_start($rest[0]);
783        }
784        elsif ($verb eq "contig_md5")
785        {
786            $self->current_file->set_contig_checksum(@rest[0,1,2]);
787        }
788    }
789    
790    sub parse_local_directive
791    {
792        my($self, $directive, $rest) = @_;
793    
794        print STDERR "Have local directive '$directive' rest='$rest'\n";
795    }
796    
797    sub parse_feature
798    {
799        my($self, $seqid, $source, $type, $start, $end, $score, $strand, $phase, $attributes) = @_;
800    
801        #print "data: seqid=$seqid source=$source type=$type start=$start end=$end\n";
802        #print "      score=$score strand=$strand phase=$phase\n";
803        #print "      $attributes\n";
804    
805        #
806        # Parse this feature, creating a GFFFeature object for it.
807        #
808    
809        my $feature = GFFFeature->new($self->fig);
810    
811        $feature->seqid($seqid);
812        $feature->source($source);
813        $feature->type($type);
814        $feature->start($start);
815        $feature->end($end);
816        $feature->score($score);
817        $feature->strand($strand);
818        $feature->phase($phase);
819    
820        my $athash = {};
821    
822        for my $attr (split(/;/, $attributes))
823        {
824            my($name, $value) = split(/=/, $attr);
825    
826            my @values = map { uri_unescape($_) } split(/,/, $value);
827    
828            # handle the aliases
829            if ($name eq "Alias") {
830             foreach my $val (@values)
831             {
832               $val = FigGFF::map_dbxref_to_seed_alias($val);
833             }
834            }
835    
836            #
837            # This might be a little goofy for the users, but we will use it
838            # for now:
839            #
840            # if there is more than one value, the value is a ref to a list
841            # of the values.
842            #
843            # Otherwise, the value is a scalar.
844            #
845    
846            if (@values > 1)
847            {
848                #
849                # Yes, you can do this ... I had to look it up :-).
850                #
851                # It's in 'man perlfaq3'.
852                #
853    
854                $value = \@values;
855            }
856            else
857            {
858                $value = $values[0];
859            }
860    
861    
862            $athash->{$name} = $value;
863    
864            #
865            # Handle the GFF3-defined attributes.
866            #
867            # These show up as Class::Accessor's on the feature object.
868            #
869    
870            if ($GFFFeature::GFF_Tags{$name})
871            {
872                $feature->set($name, $value);
873    
874                if ($name eq "Dbxref")
875                {
876                    #
877                    # If we have a SEED:figid DBxref, set the genome and fig_id attributes.
878                    #
879    
880                    my @seed_xref = grep /^"?SEED:/, @values;
881                    if (@seed_xref and $seed_xref[0] =~ /^"?SEED:(fig\|(\d+\.\d+)\..*)/)
882                    {
883                        $feature->genome($2);
884                        $feature->fig_id($1);
885                    }
886    
887                }
888            }
889        }
890        $feature->attributes($athash);
891    
892    
893        $self->current_file->add_feature($feature);
894    }
895    
896    #
897    # We come in here with the first line of the fasta already read
898    # in order to support the backward-compatiblity syntax that
899    # lets a file skip the ##FASTA directive if it wishes.
900    #
901    sub parse_fasta
902    {
903        my($self, $fh, $first_line) = @_;
904        my($cur, $cur_id);
905    
906        for ($_ = $first_line; defined($_);  $_ = <$fh>, chomp)
907        {
908            if (/^>\s*(\S+)/)
909            {
910                if ($cur)
911                {
912                    $self->handle_fasta_block($cur_id, $cur);
913                }
914    
915                $cur = '';
916                $cur_id = $1;
917            }
918            else
919            {
920                s/^\s*$//;
921                s/\s*$//;
922                if (/\s/)
923                {
924                    die "FASTA data had embedded space: $_\n";
925                }
926                $cur .= $_;
927            }
928        }
929        if ($cur)
930        {
931            $self->handle_fasta_block($cur_id, $cur);
932        }
933    }
934    
935    sub handle_fasta_block
936    {
937        my($self, $id, $data) = @_;
938    
939        my $len = length($data);
940        $self->current_file->fasta_data($id, $data);
941    }
942    
943    package GFFFeature;
944    
945    use strict;
946    use base qw(Class::Accessor);
947    
948    our @GFF_Tags = qw(ID Name  Alias Parent Target Gap Note Dbxref Ontology_term);
949    our %GFF_Tags;
950    
951    map { $GFF_Tags{$_} = 1 } @GFF_Tags;
952    
953    __PACKAGE__->mk_accessors(qw(fig seqid source type start end score strand phase attributes
954                                 genome fig_id),
955                              @GFF_Tags);
956    
957    
958    sub new
959    {
960        my($class, $fig) = @_;
961    
962        my $self = {
963            fig => $fig,
964        };
965    
966        return bless($self, $class);
967    }
968    
969    sub find_local_feature
970    {
971        my($self, $local_genome) = @_;
972        my $db = $self->fig->db_handle;
973    
974        # For debugging.
975        undef $local_genome;
976        if ($local_genome)
977        {
978            #
979            # It's a precise match. We need to determine if we have this
980            # particular feature in this SEED (it is possible for one to
981            # have exported an annotation for a feature that was added
982            # to a genome after its initial release).
983            #
984            # We do this by searching for a local feature with the same contig,
985            # start, and stop as this feature.
986            #
987    
988            my $qry = qq(SELECT id
989                         FROM features
990                         WHERE (genome = ? AND
991                                contig = ? AND
992                                minloc = ? AND
993                                maxloc = ?));
994            my $res = $db->SQL($qry, undef, $local_genome, $self->seqid,
995                               $self->start, $self->end);
996    
997            return map { $_->[0] } @$res;
998        }
999    
1000        #
1001        # Otherwise, we need to try a set of heuristics to match
1002        # this id.
1003        #
1004    
1005        #
1006        # Try matching aliases first.
1007        #
1008    
1009        my @aliases = grep { !/^\"?SEED/ } ref($self->Dbxref) ? @{$self->Dbxref} : ($self->Dbxref);
1010    
1011        my @maliases = map { FigGFF::map_dbxref_to_seed_alias($_) } @aliases;
1012    
1013        print "Found aliases @aliases\n";
1014        print "Found mapped aliases @maliases\n";
1015    
1016        for my $malias (@maliases)
1017        {
1018            my $fid = $self->fig->by_alias($malias);
1019            if ($fid)
1020            {
1021                print "Mapped $malias to $fid\n";
1022            }
1023        }
1024    
1025    }
1026    
1027    
1028    package GFFFile;
1029    
1030    use strict;
1031    use base qw(Class::Accessor);
1032    
1033    __PACKAGE__->mk_accessors(qw(fig filename features feature_index anno_start anno_end taxon_id genome_id));
1034    
1035    #
1036    # Package to hold the contents of a GFF file, and to hold the code
1037    # for mapping its contents to the local SEED.
1038    #
1039    # Created by GFFParser->parse.
1040    #
1041    
1042    sub new
1043    {
1044        my($class, $fig) = @_;
1045    
1046        my $self = {
1047            fig => $fig,
1048            features => [],
1049            contigs  => [],
1050            feature_index => {},
1051            genome_checksum => '',
1052            contig_checksum => {},
1053            features_by_genome => {},
1054        };
1055        return bless($self, $class);
1056    }
1057    
1058    sub add_feature
1059    {
1060        my($self, $feature) = @_;
1061    
1062        push(@{$self->features}, $feature);
1063        $self->feature_index->{$feature->ID} = $feature;
1064        push(@{$self->{features_by_genome}->{$feature->genome}}, $feature);
1065    }
1066    
1067    sub features_for_genome
1068    {
1069        my($self, $genome) = @_;
1070    
1071        return $self->{features_by_genome}->{$genome};
1072    }
1073    
1074    sub genome_checksum
1075    {
1076        my($self) = @_;
1077    
1078        return $self->{genome_checksum};
1079    }
1080    
1081    sub set_genome_checksum
1082    {
1083        my($self, $md5sum) = @_;
1084        $self->{genome_checksum} = $md5sum;
1085    }
1086    
1087    sub set_contig_checksum
1088    {
1089        my($self, $genome, $contig, $md5sum) = @_;
1090        $self->{contig_checksum}->{$genome}->{$contig} = $md5sum;
1091    }
1092    
1093    =head2 fasta_data()
1094    
1095    Get or set the fasta data. Given an id and some data will set the data for that id. Given an id will return the data for that id. Called without arguments will return a reference to a hash of sequences.
1096    
1097    This means that if you give it an id and sequence it will return that sequence. Hmmm.
1098    
1099    =cut
1100    
1101    sub fasta_data
1102    {
1103        my($self, $id, $data) = @_;
1104        $id && $data && ($self->{fasta_data}->{$id} = $data);
1105        $id && return $self->{fasta_data}->{$id};
1106        return $self->{fasta_data};
1107    }
1108    
1109    
1110    =head2 contigs()
1111    
1112    Add a contig to the list, or return a reference to an array of contigs
1113    
1114    =cut
1115    
1116    sub contigs
1117    {
1118        my($self, $contig) = @_;
1119        if ($contig && $contig =~ /\w\w\_\d+\.\d+/) {
1120          print STDERR "WARNING: $contig appears to have a version number. We should standardize on timming that somewhere\n";
1121        }
1122        $contig && (push @{$self->{contigs}}, $contig);
1123        return $self->{contigs};
1124    }
1125    
1126    =head2 contig_length()
1127    
1128    Get or set the length of a specfic contig.
1129      my $length=$fob->contig_length($contig, $length);
1130      my $length=$fob->contig_length($contig);
1131    
1132    =cut
1133    
1134    sub contig_length
1135    {
1136       my($self, $contig, $length) = @_;
1137       $length && ($self->{contig_length_cache}->{$contig}=$length);
1138       return $self->{contig_length_cache}->{$contig};
1139    }
1140    
1141    =head1 Information about the source of the sequence.
1142    
1143    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.
1144    
1145    =cut
1146    
1147    =head2 genome_id()
1148    
1149    Get or set a genome id for this file.
1150    
1151    =cut
1152    
1153    sub genome_id
1154    {
1155        my($self, $genomeid) = @_;
1156        $genomeid && ($self->{genome_id}=$genomeid);
1157        return $self->{genome_id};
1158    }
1159    
1160    =head2 genome_name()
1161    
1162    Get or set a genome id for this file.
1163    
1164    =cut
1165    
1166    sub genome_name
1167    {
1168        my($self, $genomename) = @_;
1169        $genomename && ($self->{genome_name}=$genomename);
1170        return $self->{genome_name};
1171    }
1172    
1173    =head2 project()
1174    
1175    Get or set the project.
1176    
1177    =cut
1178    
1179    sub project
1180    {
1181         my ($self, $pro) = @_;
1182         $pro && ($self->{project}=$pro);
1183         return $self->{project};
1184    }
1185    
1186    =head2 taxonomy()
1187    
1188    Get or set the taxonomy
1189    
1190    =cut
1191    
1192    sub taxonomy
1193    {
1194        my($self, $tax) = @_;
1195        $tax && ($self->{taxonomy}=$tax);
1196        return $self->{taxonomy};
1197    }
1198    
1199    
1200    
1201    
1202    
1203  1;  1;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3