[Bio] / FortyEight / SeedExport.pm Repository:
ViewVC logotype

Diff of /FortyEight/SeedExport.pm

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

revision 1.2, Tue Jan 16 20:29:55 2007 UTC revision 1.7, Thu Feb 15 22:41:39 2007 UTC
# Line 4  Line 4 
4    
5  use strict;  use strict;
6  use FIG;  use FIG;
7    use FIG_Config;
8  use FIGV;  use FIGV;
9  use URI::Escape;  use URI::Escape;
10  use Bio::FeatureIO;  use Bio::FeatureIO;
# Line 78  Line 79 
79      # create the variable for the bio object      # create the variable for the bio object
80      my $bio;      my $bio;
81      my $bio2;      my $bio2;
82      my $gff_export;
83    
84      # get the contigs      # get the contigs
85      foreach my $contig ($fig->contigs_of($genome)) {      foreach my $contig ($fig->contigs_of($genome)) {
# Line 99  Line 101 
101          $bio->{$contig}->add_SeqFeature($feature);          $bio->{$contig}->add_SeqFeature($feature);
102      }      }
103    
104      # get the functional role name -> GO file
105      open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";
106      my $fr2go = {};
107      while (<FH>) {
108        chomp;
109        my ($fr, $go) = split /\t/;
110        $fr2go->{$fr} = [] unless (exists $fr2go->{$fr});
111        push @{$fr2go->{$fr}}, $go;
112      }
113      close FH;
114    
115      # get the pegs      # get the pegs
116      foreach my $peg (sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome)) {      foreach my $peg (sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome)) {
117          my $note;          my $note;
118          my $func = $fig->function_of($peg, $user);          my $func = $fig->function_of($peg, $user);
119    
120          # parse out the ec number if there is one      my %ecs;
121          my $ec;      my @gos;
122          if ($func =~ /E\.*C\.*\s+(\d+|-)\.(\d+|-)\.(\d+|-)\.(\d+|-)/) {  
123              push @{$note->{"ec_number"}}, "$1.$2.$3.$4";      # get EC / GO from role
124          }      if (defined $func) {
125          foreach my $role ($fig->roles_of_function($func)) {
126          # check the aliases          my ($ec) = ($role =~ /\(EC (\d+\.\d+\.\d+\.\d+)\)/);
127          foreach my $alias ($fig->feature_aliases($peg))          $ecs{$ec} = 1 if ($ec);
128          {          push @gos, @{$fr2go->{$role}} if ($fr2go->{$role});
129              if ($alias =~ /^NP_/ || $alias =~ /YP_/) {push @{$note->{"Dbxref"}}, "NCBI_genpept:$alias"}        }
             elsif ($alias =~ /gi\|/)  
             {  
                 $alias =~ s/^gi\|//;  
                 push @{$note->{"Dbxref"}}, "NCBI_gi:$alias";  
             }  
             elsif ($alias =~ /^kegg/i)  
             {  
                 $alias =~ s/kegg\|/KEGG:/i;  
                 $alias =~ s/^(.*):/$1+/;  
                 push @{$note->{"Dbxref"}}, $alias  
             }  
             elsif ($alias =~ /^uni/)  
             {  
                 $alias =~ s/uni\|/UniProt:/;  
                 push @{$note->{"Dbxref"}}, $alias;  
             }  
             elsif ($alias =~ /^sp\|/)  
             {  
                 $alias =~ s/sp\|/Swiss-Prot:/;  
                 push @{$note->{"Dbxref"}}, $alias  
             }  
             elsif ($alias =~ /^[a-zA-Z][a-zA-Z0-9]{2,}\_[a-zA-Z]*\d+$/)  
             {  
                 # that should be the regexp for a valid locus tag. Starts with a letter,  
                 # has at least three alphanumerics before the _  
                 # then has a series of optional letters and ends with numbers  
                 push @{$note->{"locus_tag"}}, $alias;  
             }  
             elsif ($alias =~ /^.{4,6}$/)  
             {  
                 push @{$note->{"gene_symbol"}}, $alias;  
             }  
             elsif ($alias =~ /^[a-zA-Z]+\d+$/)  
             {  
                 push @{$note->{"locus"}}, $alias;  
             }  
             else {  
                 # don't know what it is so keep it as an alias  
                 push @{$note->{"Alias"}}, $alias;  
130              }              }
131    
132        # remove duplicate gos
133        my %gos = map { $_ => 1 } @gos if (scalar(@gos));
134    
135        # add GOs
136        push @{$note->{"Dbxref"}}, @gos;
137    
138        # add ECs
139        push @{$note->{"EC_number"}}, keys(%ecs);
140    
141        # get the aliases from principal id
142        my $pid = $fig->maps_to_id($peg);
143        my @rw_aliases = map { $fig->rewrite_db_xrefs($_->[0]) } $fig->mapped_prot_ids($pid);
144        my @aliases;
145        foreach my $a (@rw_aliases) {
146          push @{$note->{"Dbxref"}}, $a if ( $a );
147          }          }
148    
149          # get the links          # get the links
150          foreach my $ln ($fig->fid_links($peg))      foreach my $ln ($fig->fid_links($peg)) {
         {  
151              my ($db, $id);              my ($db, $id);
152              if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {              if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
153                  $db = "HarvardClone";                  $db = "HarvardClone";
# Line 190  Line 179 
179              if ($start > $stop) {              if ($start > $stop) {
180                $frame = $stop % 3;                $frame = $stop % 3;
181                ($start, $stop, $strand) = ($stop, $start, '-1');                ($start, $stop, $strand) = ($stop, $start, '-1');
182              }        } elsif ($start == $stop) {
             elsif ($start == $stop) {  
183                $strand = ".";                $strand = ".";
184                $frame = ".";                $frame = ".";
185              }              }
# Line 225  Line 213 
213    
214                  push(@$bio2, $feature2);                  push(@$bio2, $feature2);
215    
216            # work around to get annotations into gff
217            push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\t$func\n"
218    
219    
220              } elsif ($type eq "rna") {              } elsif ($type eq "rna") {
221                  $feature = Bio::SeqFeature::Generic->new(                  $feature = Bio::SeqFeature::Generic->new(
222                      -start    => $start,                      -start    => $start,
# Line 249  Line 241 
241    
242      # check for FeatureIO or SeqIO      # check for FeatureIO or SeqIO
243      if ($format eq "GTF") {      if ($format eq "GTF") {
244        my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");      #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
245        foreach my $feature (@$bio2) {      #foreach my $feature (@$bio2) {
246          $fio->write_feature($feature);      #$fio->write_feature($feature);
247        #}
248        open (GTF, ">$filename") or die "Cannot open file $filename.";
249        print GTF "##gff-version 3\n";
250        foreach (@$gff_export) {
251          print GTF $_;
252        }        }
253        close(GTF);
254    
255      } else {      } else {
256        my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);        my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
257        foreach my $seq (keys %$bio) {        foreach my $seq (keys %$bio) {
# Line 262  Line 261 
261    
262      return ($filename, "Output file successfully written.");      return ($filename, "Output file successfully written.");
263  }  }
264    
265    
266    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3