[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.5, Thu Feb 8 21:11:01 2007 UTC revision 1.12, Fri Sep 25 13:40:30 2009 UTC
# Line 10  Line 10 
10  use Bio::FeatureIO;  use Bio::FeatureIO;
11  use Bio::SeqIO;  use Bio::SeqIO;
12  use Bio::Seq;  use Bio::Seq;
13    use Bio::Location::Split;
14    use Bio::Location::Simple;
15  use Bio::SeqFeature::Generic;  use Bio::SeqFeature::Generic;
16  use Bio::SeqFeature::Annotated;  use Bio::SeqFeature::Annotated;
17  use File::Basename;  use File::Basename;
# Line 22  Line 24 
24      my $genome          = $parameters->{'genome'};      my $genome          = $parameters->{'genome'};
25      my $directory       = $parameters->{'directory'};      my $directory       = $parameters->{'directory'};
26      my $format          = $parameters->{'export_format'};      my $format          = $parameters->{'export_format'};
27        my $strip_ec        = $parameters->{'strip_ec'};
28        my $filename        = $parameters->{'filename'};
29    
30      # set some variables      # set some variables
31      my $user = "master:master";      my $user = "master:master";
# Line 86  Line 90 
90          my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);          my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
91          my $seq = $fig->dna_seq($genome, $cloc);          my $seq = $fig->dna_seq($genome, $cloc);
92          $seq =~ s/>//;          $seq =~ s/>//;
93          $bio->{$contig} = Bio::Seq->new(-id => "$contig", seq => $seq);          $bio->{$contig} = Bio::Seq->new( -id => $contig,
94                                            -seq => $seq,
95                                           );
96          $bio->{$contig}->desc("Contig $contig from $gs");          $bio->{$contig}->desc("Contig $contig from $gs");
97    
98          my $feature = Bio::SeqFeature::Generic->new(          my $feature = Bio::SeqFeature::Generic->new(
99              -start      => 1,              -start      => 1,
100              -end        => $fig->contig_ln($genome, $contig),              -end        => $fig->contig_ln($genome, $contig),
101              -tag        => { dbxref     => "taxon: $taxid",                                                      -tag        => { db_xref     => "taxon: $taxid",
102                               organism   => $gs,                               organism   => $gs,
103                               mol_type   => "genomic DNA",                               mol_type   => "genomic DNA",
104                               genome_id  => $genome,                               genome_id  => $genome,
# Line 103  Line 110 
110    
111      # get the functional role name -> GO file      # get the functional role name -> GO file
112      open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";      open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";
113      my $fr2go;      my $fr2go = {};
114      while (<FH>) {      while (<FH>) {
115        chomp;        chomp;
116        my ($fr, $go) = split /\t/;        my ($fr, $go) = split /\t/;
117        $fr2go->{$fr} = $go;          $fr2go->{$fr} = [] unless (exists $fr2go->{$fr});
118            push @{$fr2go->{$fr}}, $go;
119      }      }
120      close FH;      close FH;
121    
# Line 116  Line 124 
124          my $note;          my $note;
125          my $func = $fig->function_of($peg, $user);          my $func = $fig->function_of($peg, $user);
126    
127          # check for a GO number          my %ecs;
128          if (exists($fr2go->{$func})) {          my @gos;
           push @{$note->{"Dbxref"}}, $fr2go->{$func};  
         }  
   
         # parse out the ec number if there is one  
         my $ec;  
         if ($func =~ /E\.*C\.*\s+(\d+|-)\.(\d+|-)\.(\d+|-)\.(\d+|-)/) {  
             push @{$note->{"EC_number"}}, "$1.$2.$3.$4";  
         }  
129    
130          # check the aliases          # get EC / GO from role
131          foreach my $alias ($fig->feature_aliases($peg))          if (defined $func) {
132          {              foreach my $role ($fig->roles_of_function($func)) {
133              if ($alias =~ /^NP_/ || $alias =~ /YP_/) {                  my ($ec) = ($role =~ /\(EC (\d+\.\d+\.\d+\.\d+)\)/);
134                push @{$note->{"Dbxref"}}, "genpept:$alias"                  $ecs{$ec} = 1 if ($ec);
135              }                  push @gos, @{$fr2go->{$role}} if ($fr2go->{$role});
             elsif ($alias =~ /gi\|/)  
             {  
                 $alias =~ s/^gi\|//;  
                 push @{$note->{"Dbxref"}}, "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"}}, "UniProtKB/TrEMBL:$alias";  
             }  
             elsif ($alias =~ /^sp\|/)  
             {  
                 $alias =~ s/sp\|/Swiss-Prot:/;  
                 push @{$note->{"Dbxref"}}, "UniProtKB/Swiss-Prot:$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;  
136              }              }
             else {  
                 # don't know what it is so keep it as an alias  
                 push @{$note->{"Alias"}}, $alias;  
137              }              }
138    
139            # remove duplicate gos
140            my %gos = map { $_ => 1 } @gos if (scalar(@gos));
141    
142            # add GOs
143            push @{$note->{"db_xref"}}, @gos;
144    
145            # add ECs
146            push @{$note->{"EC_number"}}, keys(%ecs);
147    
148            # get the aliases from principal id
149            my $pid = $fig->maps_to_id($peg);
150            my @rw_aliases = map { $fig->rewrite_db_xrefs($_->[0]) } $fig->mapped_prot_ids($pid);
151            my @aliases;
152            foreach my $a (@rw_aliases) {
153                push @{$note->{"db_xref"}}, $a if ( $a );
154          }          }
155    
156          # get the links          # get the links
157          foreach my $ln ($fig->fid_links($peg))          foreach my $ln ($fig->fid_links($peg)) {
         {  
158              my ($db, $id);              my ($db, $id);
159              if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {              if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
160                  $db = "HarvardClone";                  $db = "HarvardClone";
# Line 192  Line 169 
169                  print STDERR "Ignored link: $ln\n";                  print STDERR "Ignored link: $ln\n";
170                  next;                  next;
171              }              }
172              push @{$note->{"Dbxref"}}, "$db:$id";              push @{$note->{"db_xref"}}, "$db:$id";
173          }          }
174    
175          # add FIG id as a note          # add FIG id as a note
176          push @{$note->{"Dbxref"}}, "FIG_ID:$peg";          # push @{$note->{"db_xref"}}, "FIG_ID:$peg";
177    
178          # get the features          # get the features
179    
180            my $loc_obj;
181          my @location = $fig->feature_location($peg);          my @location = $fig->feature_location($peg);
182            my @loc_info;
183            my $contig;
184          foreach my $loc (@location) {          foreach my $loc (@location) {
185                my($start, $stop);
186              $loc =~ /^(.*)\_(\d+)\_(\d+)$/;              $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
187              my ($contig, $start, $stop) = ($1, $2, $3);              ($contig, $start, $stop) = ($1, $2, $3);
188              my $original_contig = $contig;              my $original_contig = $contig;
189              my $strand = '1';              my $strand = '+';
190              my $frame = $start % 3;              my $frame = $start % 3;
191              if ($start > $stop) {              if ($start > $stop) {
192                $frame = $stop % 3;                $frame = $stop % 3;
193                ($start, $stop, $strand) = ($stop, $start, '-1');                  ($start, $stop, $strand) = ($stop, $start, '-');
194              }              } elsif ($start == $stop) {
             elsif ($start == $stop) {  
195                $strand = ".";                $strand = ".";
196                $frame = ".";                $frame = ".";
197              }              }
198    
199                push(@loc_info, [$contig, $start, $stop, $strand, $frame]);
200    
201                my $sloc = new Bio::Location::Simple(-start => $start,
202                                                     -end => $stop,
203                                                     -strand => $strand);
204                if (@location == 1)
205                {
206                    $loc_obj = $sloc;
207                }
208                else
209                {
210                    $loc_obj = new Bio::Location::Split() if !$loc_obj;
211                    $loc_obj->add_sub_Location($sloc);
212                }
213            }
214    
215              my $source = "FIG";              my $source = "FIG";
216              my $type = $fig->ftype($peg);              my $type = $fig->ftype($peg);
217              my $feature;              my $feature;
218    
219            # strip EC from function
220            my $func_ok = $func;
221            if ($strip_ec) {
222                $func_ok =~ s/\s+\(EC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
223                $func_ok =~ s/\s+\(TC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
224            }
225    
226              if ($type eq "peg") {              if ($type eq "peg") {
227                  $feature = Bio::SeqFeature::Generic->new(              $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
                     -start    => $start,  
                     -end      => $stop,  
                     -strand   => $strand,  
228                      -primary  => 'CDS',                      -primary  => 'CDS',
229                      -tag      => {                      -tag      => {
230                                    product     => $func,                                                           product     => $func_ok,
231                                    translation => $fig->get_translation($peg),                                    translation => $fig->get_translation($peg),
232                                   },                                   },
233                      );                      );
# Line 233  Line 236 
236                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
237                  }                  }
238    
                 my $feature2 = Bio::SeqFeature::Annotated->new( -start    => $start,  
                                                                 -end      => $stop,  
                                                                 -strand   => $strand,  
                                                                 -phase    => 0,  
                                                                 -frame    => $frame,  
                                                                 -source   => $source,  
                                                                 -type     => "CDS",  
                                                                 -seq_id   => $contig );  
   
                 push(@$bio2, $feature2);  
   
239                  # work around to get annotations into gff                  # work around to get annotations into gff
240                  push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\t$func\n"              # this is probably still wrong for split locations.
241                $func_ok =~ s/ #.+//;
242                $func_ok =~ s/;/%3B/g;
243                $func_ok =~ s/,/2C/g;
244                $func_ok =~ s/=//g;
245                for my $l (@loc_info)
246                {
247                  my $ec = "";
248                  my @ecs = ($func =~ /[\(\[]*EC[\s:]?(\d+\.[\d-]+\.[\d-]+\.[\d-]+)[\)\]]*/ig);
249                  if (scalar(@ecs)) {
250                    $ec = ";Ontology_term=".join(',', map { "KEGG_ENZYME:" . $_ } @ecs);
251                  }
252                  my($contig, $start, $stop, $strand, $frame) = @$l;
253                  push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\tID=".$peg.";Name=".$func_ok.$ec."\n";
254                }
255    
256    
257              } elsif ($type eq "rna") {              } elsif ($type eq "rna") {
258                  $feature = Bio::SeqFeature::Generic->new(              my $primary;
259                      -start    => $start,              if ( $func =~ /tRNA/ ) {
260                      -end      => $stop,                  $primary = 'tRNA';
261                      -strand   => $strand,              } elsif ( $func =~ /(Ribosomal RNA|5S RNA)/ ) {
262                      -primary  => 'RNA',                  $primary = 'rRNA';
263                } else {
264                    $primary = 'RNA';
265                }
266    
267                $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
268                                                         -primary  => $primary,
269                                                         -tag      => {
270                                                             product => $func,
271                                                         },
272    
273                      );                      );
274                $func_ok =~ s/ #.+//;
275                $func_ok =~ s/;/%3B/g;
276                $func_ok =~ s/,/2C/g;
277                $func_ok =~ s/=//g;
278                  foreach my $tagtype (keys %$note) {                  foreach my $tagtype (keys %$note) {
279                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
280    
281                    # work around to get annotations into gff
282                    for my $l (@loc_info)
283                    {
284                        my($contig, $start, $stop, $strand, $frame) = @$l;
285                        push @$gff_export, "$contig\t$source\t$primary\t$start\t$stop\t.\t$strand\t.\tID=$peg;Name=$func_ok\n";
286                    }
287                  }                  }
288    
289              } else {              } else {
290                  print STDERR "unhandled feature type: $type\n";                  print STDERR "unhandled feature type: $type\n";
291              }              }
   
292              $bio->{$contig}->add_SeqFeature($feature);              $bio->{$contig}->add_SeqFeature($feature);
293          }          }
294      }  
295    
296      # generate filename      # generate filename
297      my $filename = $directory . $genome . "." . $format_ending;      if (!$filename)
298        {
299            $filename = $directory . $genome . "." . $format_ending;
300        }
301    
302      # check for FeatureIO or SeqIO      # check for FeatureIO or SeqIO
303      if ($format eq "GTF") {      if ($format eq "GTF") {
# Line 292  Line 321 
321    
322      return ($filename, "Output file successfully written.");      return ($filename, "Output file successfully written.");
323  }  }
324    
325    
326    

Legend:
Removed from v.1.5  
changed lines
  Added in v.1.12

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3