[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.1, Thu Jan 11 17:06:22 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; # this is used to escape the sequences and comes from http://search.cpan.org/~rse/lcwa-1.0.0/  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::SeqFeature::Generic;  use Bio::SeqFeature::Generic;
14    use Bio::SeqFeature::Annotated;
15  use File::Basename;  use File::Basename;
16    
17  sub export {  sub export {
# Line 30  Line 32 
32          $format_ending = "gbk";          $format_ending = "gbk";
33      } elsif ($format eq "GTF") {      } elsif ($format eq "GTF") {
34          $format_ending = "gtf";          $format_ending = "gtf";
35      } elsif ($format eq "gff") {
36        $format = "GTF";
37        $format_ending = "gff";
38      } elsif ($format eq "embl") {
39        $format_ending = "embl";
40      } else {      } else {
41          die "unknown export format: $format\n";          die "unknown export format: $format\n";
42      }      }
# Line 71  Line 78 
78    
79      # create the variable for the bio object      # create the variable for the bio object
80      my $bio;      my $bio;
81      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)) {
86          my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);          my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
87          my $seq = $fig->dna_seq($genome, $cloc);          my $seq = $fig->dna_seq($genome, $cloc);
88        $seq =~ s/>//;
89          $bio->{$contig} = Bio::Seq->new(-id => "$contig", seq => $seq);          $bio->{$contig} = Bio::Seq->new(-id => "$contig", seq => $seq);
90          $bio->{$contig}->desc("Contig $contig from $gs");          $bio->{$contig}->desc("Contig $contig from $gs");
91          my $feature = Bio::SeqFeature::Generic->new(          my $feature = Bio::SeqFeature::Generic->new(
# Line 91  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});
             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;  
129              }              }
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 178  Line 175 
175              my ($contig, $start, $stop)=($1, $2, $3);              my ($contig, $start, $stop)=($1, $2, $3);
176              my $original_contig=$contig;              my $original_contig=$contig;
177              my $strand='1';              my $strand='1';
178          my $frame = $start % 3;
179              next if (defined $beg && ($start < $beg || $stop < $beg));        if ($start > $stop) {
180              next if (defined $end && ($start > $end || $stop > $end));          $frame = $stop % 3;
181            ($start, $stop, $strand) = ($stop, $start, '-1');
182              if ($start > $stop) {($start, $stop, $strand)=($stop, $start, '-1')}        } elsif ($start == $stop) {
183              elsif ($start == $stop) {$strand="."}          $strand = ".";
184            $frame = ".";
185          }
186          my $source = "FIG";
187              my $type = $fig->ftype($peg);              my $type = $fig->ftype($peg);
188              my $feature;              my $feature;
189              if ($type eq "peg") {              if ($type eq "peg") {
# Line 202  Line 201 
201                  foreach my $tagtype (keys %$note) {                  foreach my $tagtype (keys %$note) {
202                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});                      $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
203                  }                  }
204    
205            my $feature2 = Bio::SeqFeature::Annotated->new( -start    => $start,
206                                                            -end      => $stop,
207                                                            -strand   => $strand,
208                                                            -phase    => 0,
209                                                            -frame    => $frame,
210                                                            -source   => $source,
211                                                            -type     => "CDS",
212                                                            -seq_id   => $contig );
213    
214            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 226  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 => $format);      #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
245          foreach my $seq (keys %$bio) {      #foreach my $feature (@$bio2) {
246              my $feat = Bio::FeatureIO->new(-seq => $seq , -format => 'features');      #$fio->write_feature($feature);
247              while ( my $feature = $feat->next_feature() ) {      #}
248                  $fio->write_feature($feature);      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);
# Line 243  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.1  
changed lines
  Added in v.1.7

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3