[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.3, Wed Jan 31 22:10:06 2007 UTC revision 1.12, Fri Sep 25 13:40:30 2009 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;
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 21  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 78  Line 83 
83      # create the variable for the bio object      # create the variable for the bio object
84      my $bio;      my $bio;
85      my $bio2;      my $bio2;
86        my $gff_export;
87    
88      # get the contigs      # get the contigs
89      foreach my $contig ($fig->contigs_of($genome)) {      foreach my $contig ($fig->contigs_of($genome)) {
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 99  Line 108 
108          $bio->{$contig}->add_SeqFeature($feature);          $bio->{$contig}->add_SeqFeature($feature);
109      }      }
110    
111        # get the functional role name -> GO file
112        open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";
113        my $fr2go = {};
114        while (<FH>) {
115            chomp;
116            my ($fr, $go) = split /\t/;
117            $fr2go->{$fr} = [] unless (exists $fr2go->{$fr});
118            push @{$fr2go->{$fr}}, $go;
119        }
120        close FH;
121    
122      # get the pegs      # get the pegs
123      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)) {
124          my $note;          my $note;
125          my $func = $fig->function_of($peg, $user);          my $func = $fig->function_of($peg, $user);
126    
127          # parse out the ec number if there is one          my %ecs;
128          my $ec;          my @gos;
         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";  
136              }              }
             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;  
             }  
             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 175  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 216  Line 236 
236                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
237                  }                  }
238    
239                  my $feature2 = Bio::SeqFeature::Annotated->new( -start    => $start,              # work around to get annotations into gff
240                                                                  -end      => $stop,              # this is probably still wrong for split locations.
241                                                                  -strand   => $strand,              $func_ok =~ s/ #.+//;
242                                                                  -phase    => 0,              $func_ok =~ s/;/%3B/g;
243                                                                  -frame    => $frame,              $func_ok =~ s/,/2C/g;
244                                                                  -source   => $source,              $func_ok =~ s/=//g;
245                                                                  -type     => "CDS",              for my $l (@loc_info)
246                                                                  -seq_id   => $contig );              {
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    
                 push(@$bio2, $feature2);  
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") {
304        my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");          #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
305        foreach my $feature (@$bio2) {          #foreach my $feature (@$bio2) {
306          $fio->write_feature($feature);          #$fio->write_feature($feature);
307            #}
308            open (GTF, ">$filename") or die "Cannot open file $filename.";
309            print GTF "##gff-version 3\n";
310            foreach (@$gff_export) {
311                print GTF $_;
312        }        }
313            close(GTF);
314    
315      } else {      } else {
316        my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);        my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
317        foreach my $seq (keys %$bio) {        foreach my $seq (keys %$bio) {
# Line 264  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.3  
changed lines
  Added in v.1.12

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3