[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.2, Tue Mar 15 19:47:15 2005 UTC revision 1.8, Sat Apr 9 21:09:30 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        # if it is not a valid xref just return it
89        return $dbxref unless (m/:/);
90    
91        my($type, $ref) = split(/:/, $dbxref, 2);
92    
93        if ($type eq "NCBI_NP")
94        {
95            return "NP_$ref";
96        }
97        elsif ($type eq "NCBI_gi")
98        {
99            return "gi|$ref";
100        }
101        elsif ($type eq "KEGG")
102        {
103            $ref =~ s/ /:/;
104            return "kegg|$ref";
105        }
106        elsif ($type eq "UniProt")
107        {
108            return "uni|$ref";
109        }
110        elsif ($type eq "Swiss-Prot")
111        {
112            return "sp|$ref";
113        }
114    
115        return $dbxref; # just return itself if we don't know what it is.
116    }
117    
118  package GFFWriter;  package GFFWriter;
119    
120  use strict;  use strict;
121  use FIG;  use FIG;
122    use FigGFF;
123    
124  use Carp;  use Carp;
125  use URI::Escape;  use URI::Escape;
# Line 28  Line 137 
137    
138      map { $default_options->{$_} = $options{$_} } keys(%options);      map { $default_options->{$_} = $options{$_} } keys(%options);
139    
140        # added contig_start_cache and contig_end_cache because we have something like
141        # sequence-region   contigname      1       10000
142        # in general we will set contig_start_cache == 1
143        # and contig_end_cache == contig_length_cache
144    
145      my $self = {      my $self = {
146          options => $default_options,          options => $default_options,
147          contig_length_cache => {},          contig_length_cache => {},
148            contig_start_cache  => {},
149            contig_end_cache    => {},
150          fig => $fig,          fig => $fig,
151      };      };
152    
# Line 52  Line 168 
168    
169  sub gff3_for_feature  sub gff3_for_feature
170  {  {
171      my($self, $fid, $user, $user_note) = @_;      my($self, $fid, $user, $user_note, $in_aliases, $in_loc) = @_;
172    
173      #      #
174      # Options we need to figure out somehow.      # Options we need to figure out somehow.
# Line 69  Line 185 
185      my $contig_data;      my $contig_data;
186      my $linelength = $options->{linelength};      my $linelength = $options->{linelength};
187    
188        my $beg = $self->{options}->{beg};
189        my $end = $self->{options}->{end};
190    
191      my $fig = $self->{fig};      my $fig = $self->{fig};
192    
193      #      #
194      # Do this first to make sure that we really have a feature.      # Do this first to make sure that we really have a feature.
195      #      #
196      my @location = $fig->feature_location($fid);      my @location = ref($in_loc) ? @$in_loc : $fig->feature_location($fid);
197      if (@location == 0 or !defined($location[0]))      if (@location == 0 or !defined($location[0]))
198      {      {
199          warn "No location found for feature $fid\n";          warn "No location found for feature $fid\n";
# Line 97  Line 216 
216      # all the aliases we are going to use      # all the aliases we are going to use
217      #      #
218      my @alias;      my @alias;
219        my @xref;
220    
221      if ($options->{with_assignments})      if ($options->{with_assignments})
222      {      {
223          my $func = $fig->function_of($fid, $user);          my $func = $fig->function_of($fid, $user);
224          if ($func)          if ($func)
225          {          {
226              push @$note, ("Note=\"" . uri_escape($func) . '"');              push @$note, ("Note=" . uri_escape($func));
227          }          }
228      }      }
229    
230      if ($options->{with_aliases})      if ($options->{with_aliases})
231      {      {
232      # now find aliases      # now find aliases
233          foreach my $alias ($fig->feature_aliases($fid))          my @feat_aliases = ref($in_aliases) ? @$in_aliases : $fig->feature_aliases($fid);
234          {          foreach my $alias (@feat_aliases)
             if ($alias =~ /^NP/)  
             {  
                 push @$note, "Dbxref=\"NCBI_genpept:$alias\"";  
             }  
             elsif ($alias =~ /gi\|/)  
             {  
                 $alias =~ s/^gi\|//;  
                 push @$note, "Dbxref=\"NCBI_gi:$alias\"";  
             }  
             elsif ($alias =~ /^kegg/i)  
235              {              {
236                  $alias =~ s/kegg\|/KEGG:/i;              my $mapped = FigGFF::map_seed_alias_to_dbxref($alias);
237                  $alias =~ s/^(.*):/$1+/;              if ($mapped)
                 push @$note, "Dbxref=\"$alias\"";  
             }  
             elsif ($alias =~ /^uni/)  
238              {              {
239                  $alias =~ s/uni\|/UniProt:/;                  push(@xref, $mapped);
                 # $note = check_go($note, $alias) if ($go);  
                 push @$note, "Dbxref=\"$alias\"";  
             }  
             elsif ($alias =~ /^sp\|/)  
             {  
                 $alias =~ s/sp\|/Swiss-Prot:/;  
                 push @$note, "Dbxref=\"$alias\"";  
240              }              }
241              else              else
242              {              {
243                  # don't know what it is so keep it as an alias                  push(@alias, $alias);
                 $alias = uri_escape($alias); # just in case!  
                 push @alias, $alias;  
244              }              }
245          }          }
246      }      }
# Line 150  Line 248 
248      # 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
249      if (@alias)      if (@alias)
250      {      {
251          push @$note, "Alias=\"". join (",", @alias) . '"';          push @$note, "Alias=". join (",", map { uri_escape($_) } @alias);
252      }      }
253    
254      #      #
# Line 163  Line 261 
261      }      }
262    
263      # 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
264      push @$note, "Dbxref=\"SEED:$fid\"";      #
265        # for now, make SEED xref the first in the list so we can search for DBxref=SEEd
266        #
267    
268        unshift(@xref, "SEED:$fid");
269    
270        push(@$note, "Dbxref=" . join(",", map { uri_escape($_) } @xref));
271    
272      # 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
273      my $allnotes;      my $allnotes;
# Line 252  Line 356 
356              # gene: SO:0000704              # gene: SO:0000704
357              # cds: SO:0000316              # cds: SO:0000316
358              # protein:  NOT VALID: should be protein_coding_primary_transcript SO:0000120              # protein:  NOT VALID: should be protein_coding_primary_transcript SO:0000120
359              foreach my $type (keys %outputtype)  
360              {              #
361                  my ($id, $type)=($trunc{$type} . "." . $geneid, $type);              # For now, we will only output CDS features, and include
362                # the translation as an attribute (attribuute name is
363                # translation_id, value is a key that corresponds to a FASTA
364                # section at the end of the file).
365                #
366    
367                my $type = "cds";
368    
369                my $protein_id = "pro.$geneid";
370                my $cds_id = "cds.$geneid";
371    
372                  # we want to store some sequences to be output                  # we want to store some sequences to be output
373                  if ($outputfasta && $type eq "protein")              if ($outputfasta)
374                  {                  {
375                      my $addseq = $fig->get_translation($fid);                      my $addseq = $fig->get_translation($fid);
376    
# Line 265  Line 379 
379                      #                      #
380                      $addseq =~ s/(.{$linelength})/$1\n/g;                      $addseq =~ s/(.{$linelength})/$1\n/g;
381                      chomp($addseq);                      chomp($addseq);
382                      $fastasequences .= ">$id\n$addseq\n";                  $fastasequences .= ">$protein_id\n$addseq\n";
383                  }  
384                  if ($outputfasta && $type eq "cds")                  my $addseq = uc($fig->dna_seq($genome, @location));
                 {  
                     my $addseq = uc($fig->dna_seq($genome, $fig->feature_location($fid)));  
385                      $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);                      $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
386                      $fastasequences .= ">$id\n$addseq\n";  
387                  }                  $fastasequences .= ">$cds_id\n$addseq\n";
388                  # correct the incorrect sofa term  
389                  if ($type eq "protein")                  $allnotes .= ";translation_id=$protein_id";
                 {  
                     $type = "SO:0000120";  
390                  }                  }
391    
392                  push (@{$contig_data->{$contig_key}},                  push (@{$contig_data->{$contig_key}},
393                        (join "\t",                        (join "\t",
394                         ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$id;$allnotes")));                     ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "ID=$cds_id;$allnotes")));
395              } # end the foreach my $type          }
         } # end the if type==peg  
396          elsif ($type eq "rna")          elsif ($type eq "rna")
397          {          {
398              $fid =~ /\.rna\.(\d+)/;              $fid =~ /\.rna\.(\d+)/;
# Line 294  Line 403 
403              my ($id, $type)=("rna.$geneid", "tRNA");              my ($id, $type)=("rna.$geneid", "tRNA");
404              if ($outputfasta)              if ($outputfasta)
405              {              {
406                  my $addseq = $fig->dna_seq($genome, $fig->feature_location($fid));                  my $addseq = $fig->dna_seq($genome, @location);
407                  $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);                  $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
408                  $fastasequences .= ">$id\n$addseq\n";                  $fastasequences .= ">$id\n$addseq\n";
409              }              }
# Line 536  Line 645 
645              # Directive.              # Directive.
646              #              #
647    
648              if ($1 eq "seed")              if (lc($1) eq "seed")
649              {              {
650                  $self->parse_seed_directive($2);                  $self->parse_seed_directive($2);
651              }              }
# Line 563  Line 672 
672  {  {
673      my($self, $directive, $rest) = @_;      my($self, $directive, $rest) = @_;
674    
675        $directive = lc($directive);
676        my @rest=split /\t/, $rest;
677    
678        if ($directive eq "genome")
679        {
680            $self->current_file->genome_id($rest[0]);
681            $self->current_file->genome_name($rest[1]);
682        }
683        elsif ($directive eq "genome_md5")
684        {
685            $self->current_file->set_genome_checksum(@rest[0,1]);
686        }
687        elsif ($directive eq "origin")
688        {
689            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";
690            print STDERR "At the moment ORIGIN is returned by \$feat->project\n";
691            $self->current_file->project($rest[0]);
692        }
693        elsif ($directive eq "taxonomy")
694        {
695            $self->current_file->taxonomy($rest);
696        }
697        elsif ($directive eq "sequence-region")
698        {
699            $self->current_file->contigs($rest[0]);
700            $self->{contig_length_cache}->{$rest[0]}=$rest[2]-$rest[1];
701            $self->{contig_start_cache}->{$rest[0]}=$rest[1];
702            $self->{contig_end_cache}->{$rest[0]}=$rest[2];
703        }
704        else
705        {
706      print "Have gff3 directive '$directive' rest='$rest'\n";      print "Have gff3 directive '$directive' rest='$rest'\n";
707  }  }
708    
709    }
710    
711  sub parse_seed_directive  sub parse_seed_directive
712  {  {
713      my($self, $rest) = @_;      my($self, $rest) = @_;
714    
715      my($verb, @rest) = split(/\t/, $rest);      my($verb, @rest) = split(/\t/, $rest);
716    
717        # are we case sensitive? I don't think so
718        $verb-lc($verb);
719    
720      if ($verb eq "anno_start")      if ($verb eq "anno_start")
721      {      {
722          $self->current_file->anno_start($rest[0]);          $self->current_file->anno_start($rest[0]);
# Line 580  Line 725 
725      {      {
726          $self->current_file->anno_start($rest[0]);          $self->current_file->anno_start($rest[0]);
727      }      }
     elsif ($verb eq "genome_md5")  
     {  
         $self->current_file->set_genome_checksum(@rest[0,1]);  
     }  
728      elsif ($verb eq "contig_md5")      elsif ($verb eq "contig_md5")
729      {      {
730          $self->current_file->set_contig_checksum(@rest[0,1,2]);          $self->current_file->set_contig_checksum(@rest[0,1,2]);
# Line 626  Line 767 
767      {      {
768          my($name, $value) = split(/=/, $attr);          my($name, $value) = split(/=/, $attr);
769    
770            my @values = map { uri_unescape($_) } split(/,/, $value);
771    
772            # handle the aliases
773            if ($name eq "Alias") {
774             foreach my $val (@values)
775             {
776               $val = FigGFF::map_dbxref_to_seed_alias($val);
777             }
778            }
779    
780          #          #
781          # Handle the GFF3-defined attributes          # This might be a little goofy for the users, but we will use it
782            # for now:
783            #
784            # if there is more than one value, the value is a ref to a list
785            # of the values.
786            #
787            # Otherwise, the value is a scalar.
788          #          #
   
         my @values = split(/,/, $value);  
789    
790          if (@values > 1)          if (@values > 1)
791          {          {
792              my $vlist = [];              #
793              for my $value (@values)              # Yes, you can do this ... I had to look it up :-).
794              {              #
795                  $value = uri_unescape($value);              # It's in 'man perlfaq3'.
796                  push(@$vlist, $value);              #
797              }  
798              $value = $vlist;              $value = \@values;
799          }          }
800    
801          $athash->{$name} = $value;          $athash->{$name} = $value;
802    
803            #
804            # Handle the GFF3-defined attributes.
805            #
806            # These show up as Class::Accessor's on the feature object.
807            #
808    
809          if ($GFFFeature::GFF_Tags{$name})          if ($GFFFeature::GFF_Tags{$name})
810          {          {
811              $feature->set($name, $value);              $feature->set($name, $value);
812    
813                if ($name eq "Dbxref")
814                {
815                    #
816                    # If we have a SEED:figid DBxref, set the genome and fig_id attributes.
817                    #
818    
819                    my @seed_xref = grep /^"?SEED:/, @values;
820                    if (@seed_xref and $seed_xref[0] =~ /^"?SEED:(fig\|(\d+\.\d+)\..*)/)
821                    {
822                        $feature->genome($2);
823                        $feature->fig_id($1);
824          }          }
825    
826      }      }
827            }
828        }
829      $feature->attributes($athash);      $feature->attributes($athash);
830    
831    
832      $self->current_file->add_feature($feature);      $self->current_file->add_feature($feature);
833  }  }
834    
# Line 700  Line 876 
876      my($self, $id, $data) = @_;      my($self, $id, $data) = @_;
877    
878      my $len = length($data);      my $len = length($data);
879      $self->current_file->set_fasta_data($id, $data);      $self->current_file->fasta_data($id, $data);
880  }  }
881    
882  package GFFFeature;  package GFFFeature;
# Line 713  Line 889 
889    
890  map { $GFF_Tags{$_} = 1 } @GFF_Tags;  map { $GFF_Tags{$_} = 1 } @GFF_Tags;
891    
892  __PACKAGE__->mk_accessors(qw(fig seqid source type start end score strand phase attributes),  __PACKAGE__->mk_accessors(qw(fig seqid source type start end score strand phase attributes genome fig_id),
893                            @GFF_Tags);                            @GFF_Tags);
894    
895    
# Line 728  Line 904 
904      return bless($self, $class);      return bless($self, $class);
905  }  }
906    
907    sub find_local_feature
908    {
909        my($self, $local_genome) = @_;
910        my $db = $self->fig->db_handle;
911    
912        # For debugging.
913        undef $local_genome;
914        if ($local_genome)
915        {
916            #
917            # It's a precise match. We need to determine if we have this
918            # particular feature in this SEED (it is possible for one to
919            # have exported an annotation for a feature that was added
920            # to a genome after its initial release).
921            #
922            # We do this by searching for a local feature with the same contig,
923            # start, and stop as this feature.
924            #
925    
926            my $qry = qq(SELECT id
927                         FROM features
928                         WHERE (genome = ? AND
929                                contig = ? AND
930                                minloc = ? AND
931                                maxloc = ?));
932            my $res = $db->SQL($qry, undef, $local_genome, $self->seqid,
933                               $self->start, $self->end);
934    
935            return map { $_->[0] } @$res;
936        }
937    
938        #
939        # Otherwise, we need to try a set of heuristics to match
940        # this id.
941        #
942    
943        #
944        # Try matching aliases first.
945        #
946    
947        my @aliases = grep { !/^\"?SEED/ } ref($self->Dbxref) ? @{$self->Dbxref} : ($self->Dbxref);
948    
949        my @maliases = map { FigGFF::map_dbxref_to_seed_alias($_) } @aliases;
950    
951        print "Found aliases @aliases\n";
952        print "Found mapped aliases @maliases\n";
953    
954        for my $malias (@maliases)
955        {
956            my $fid = $self->fig->by_alias($malias);
957            if ($fid)
958            {
959                print "Mapped $malias to $fid\n";
960            }
961        }
962    
963    }
964    
965    
966  package GFFFile;  package GFFFile;
967    
968  use strict;  use strict;
# Line 749  Line 984 
984      my $self = {      my $self = {
985          fig => $fig,          fig => $fig,
986          features => [],          features => [],
987            contigs  => [],
988          feature_index => {},          feature_index => {},
989            genome_checksum => {},
990            contig_checksum => {},
991            features_by_genome => {},
992      };      };
993      return bless($self, $class);      return bless($self, $class);
994  }  }
# Line 760  Line 999 
999    
1000      push(@{$self->features}, $feature);      push(@{$self->features}, $feature);
1001      $self->feature_index->{$feature->ID} = $feature;      $self->feature_index->{$feature->ID} = $feature;
1002        push(@{$self->{features_by_genome}->{$feature->genome}}, $feature);
1003    }
1004    
1005    sub features_for_genome
1006    {
1007        my($self, $genome) = @_;
1008    
1009        return $self->{features_by_genome}->{$genome};
1010    }
1011    
1012    sub genome_checksum
1013    {
1014        my($self, $genome) = @_;
1015    
1016        return $self->{genome_checksum}->{$genome};
1017    }
1018    
1019    sub genome_checksums
1020    {
1021        my($self) = @_;
1022    
1023        return [map { [$_, $self->{genome_checksum}->{$_}] } keys(%{$self->{genome_checksum}})];
1024  }  }
1025    
1026  sub set_genome_checksum  sub set_genome_checksum
# Line 774  Line 1035 
1035      $self->{contig_checksum}->{$genome}->{$contig} = $md5sum;      $self->{contig_checksum}->{$genome}->{$contig} = $md5sum;
1036  }  }
1037    
1038  sub set_fasta_data  =head2 fasta_data()
1039    
1040    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.
1041    
1042    This means that if you give it an id and sequence it will return that sequence. Hmmm.
1043    
1044    =cut
1045    
1046    sub fasta_data
1047  {  {
1048      my($self, $id, $data) = @_;      my($self, $id, $data) = @_;
1049        $id && $data && ($self->{fasta_data}->{$id} = $data);
1050        $id && return $self->{fasta_data}->{$id};
1051        return $self->{fasta_data};
1052    }
1053    
1054    
1055    =head2 contigs()
1056    
1057    Add a contig to the list, or return a reference to an array of contigs
1058    
1059    =cut
1060    
1061    sub contigs
1062    {
1063        my($self, $contig) = @_;
1064        $contig && (push @{$self->{contigs}}, $contig);
1065        return $self->{contigs};
1066    }
1067    
1068    =head2 contig_length()
1069    
1070    Get or set the length of a specfic contig.
1071      my $length=$fob->contig_length($contig, $length);
1072      my $length=$fob->contig_length($contig);
1073    
1074    =cut
1075    
1076    sub contig_length
1077    {
1078       my($self, $contig, $length) = @_;
1079       $length && ($self->{contig_length_cache}->{$contig}=$length);
1080       return $self->{contig_length_cache}->{$contig};
1081    }
1082    
1083    =head1 Information about the source of the sequence.
1084    
1085    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.
1086    
1087    =cut
1088    
1089    =head2 genome_id()
1090    
1091      $self->{fasta_data}->{$id} = $data;  Get or set a genome id for this file.
1092    
1093    =cut
1094    
1095    sub genome_id
1096    {
1097        my($self, $genomeid) = @_;
1098        $genomeid && ($self->{genome_id}=$genomeid);
1099        return $self->{genome_id};
1100    }
1101    
1102    =head2 genome_name()
1103    
1104    Get or set a genome id for this file.
1105    
1106    =cut
1107    
1108    sub genome_name
1109    {
1110        my($self, $genomename) = @_;
1111        $genomename && ($self->{genome_name}=$genomename);
1112        return $self->{genome_name};
1113    }
1114    
1115    =head2 project()
1116    
1117    Get or set the project.
1118    
1119    =cut
1120    
1121    sub project
1122    {
1123         my ($self, $pro) = @_;
1124         $pro && ($self->{project}=$pro);
1125         return $self->{project};
1126    }
1127    
1128    =head2 taxonomy()
1129    
1130    Get or set the taxonomy
1131    
1132    =cut
1133    
1134    sub taxonomy
1135    {
1136        my($self, $tax) = @_;
1137        $tax && ($self->{taxonomy}=$tax);
1138        return $self->{taxonomy};
1139  }  }
1140    
1141    
1142    
1143    
1144    
1145  1;  1;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3