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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.20 - (view) (download) (as text)

1 : paczian 1.1 package SeedExport;
2 :    
3 :     1;
4 :    
5 : olson 1.15 use Data::Dumper;
6 : paczian 1.1 use strict;
7 :     use FIG;
8 : paczian 1.4 use FIG_Config;
9 : paczian 1.1 use FIGV;
10 : paczian 1.2 use URI::Escape;
11 : paczian 1.1 use Bio::FeatureIO;
12 :     use Bio::SeqIO;
13 :     use Bio::Seq;
14 : olson 1.19 use Bio::Seq::RichSeq;
15 : olson 1.10 use Bio::Location::Split;
16 :     use Bio::Location::Simple;
17 : paczian 1.1 use Bio::SeqFeature::Generic;
18 : paczian 1.2 use Bio::SeqFeature::Annotated;
19 : olson 1.19 use Bio::Species;
20 : paczian 1.1 use File::Basename;
21 :    
22 :     sub export {
23 : olson 1.10 my ($parameters) = @_;
24 :    
25 :     # get parameters
26 :     my $virt_genome_dir = $parameters->{'virtual_genome_directory'};
27 :     my $genome = $parameters->{'genome'};
28 :     my $directory = $parameters->{'directory'};
29 :     my $format = $parameters->{'export_format'};
30 :     my $strip_ec = $parameters->{'strip_ec'};
31 : olson 1.11 my $filename = $parameters->{'filename'};
32 : olson 1.10
33 :     # set some variables
34 :     my $user = "master:master";
35 :     my $format_ending;
36 :    
37 :     # check format
38 :     if ($format eq "genbank") {
39 :     $format_ending = "gbk";
40 :     } elsif ($format eq "GTF") {
41 :     $format_ending = "gtf";
42 :     } elsif ($format eq "gff") {
43 :     $format = "GTF";
44 :     $format_ending = "gff";
45 :     } elsif ($format eq "embl") {
46 :     $format_ending = "embl";
47 :     } else {
48 :     die "unknown export format: $format\n";
49 :     }
50 :    
51 :     # initialize fig object
52 :     my $fig;
53 :     my $virt_genome;
54 :     if ($virt_genome_dir) {
55 :     $fig = new FIGV($virt_genome_dir);
56 :     $virt_genome = basename($virt_genome_dir);
57 :     } else {
58 :     $fig = new FIG;
59 :     }
60 :    
61 :     # check for genus species
62 :     my $gs = $fig->genus_species($genome);
63 :     unless ($gs) {
64 : olson 1.16 warn "No genome name set for $genome\n";
65 :     $gs = "Unknown sp.";
66 : olson 1.10 }
67 :    
68 :     # get taxonomy id and taxonomy
69 :     my $taxid = $genome;
70 :     $taxid =~ s/\.\d+$//;
71 :     my $taxonomy = $fig->taxonomy_of($genome);
72 :    
73 :     # get project
74 :     if ($virt_genome_dir and $genome eq $virt_genome) {
75 :     open(PROJECT, "<$virt_genome_dir/PROJECT") or die "Error opening $virt_genome_dir/PROJECT: $!\n";
76 :     } else {
77 :     open( PROJECT, "<$FIG_Config::organisms/$genome/PROJECT" ) or die "Error opening $FIG_Config::organisms/$genome/PROJECT: $!\n";
78 : paczian 1.1 }
79 : olson 1.10 my @project = <PROJECT>;
80 :     chomp(@project);
81 :     close PROJECT;
82 :     map {s/^\<a href\=\".*?\"\>(.*)\<\/a\>/$1/i} @project;
83 :    
84 :     # get md5 checksum
85 :     my $md5 = $fig->genome_md5sum($genome);
86 :    
87 :     # create the variable for the bio object
88 :     my $bio;
89 :     my $bio2;
90 :     my $gff_export;
91 : olson 1.19
92 :     my @tax = split(/;\s+/, $taxonomy);
93 : olson 1.20
94 :     #
95 :     # To avoid issues with EMBL exports on long taxonomies, truncate
96 :     # any field in the taxonomy to 74 chars.
97 :     #
98 :     if (lc($format) eq 'embl')
99 :     {
100 :     for my $t (@tax)
101 :     {
102 :     if ((my $l = length($t)) > 73)
103 :     {
104 :    
105 :     $t = substr($t, 0, 73);
106 :     print STDERR "EMBL export: truncating long taxonomy string (length=$l) to $t\n";
107 :     }
108 :     }
109 :     }
110 :    
111 : olson 1.19 my $species = Bio::Species->new(-classification => [reverse @tax]);
112 : olson 1.10
113 : olson 1.19 my $gc = 11;
114 :     if (open(my $gcfh, "<", $fig->organism_directory($genome) . "/GENETIC_CODE"))
115 :     {
116 :     $gc = <$gcfh>;
117 :     chomp $gc;
118 :     }
119 :    
120 :    
121 : olson 1.10 # get the contigs
122 :     foreach my $contig ($fig->contigs_of($genome)) {
123 :     my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
124 :     my $seq = $fig->dna_seq($genome, $cloc);
125 :     $seq =~ s/>//;
126 : olson 1.19 $bio->{$contig} = Bio::Seq::RichSeq->new( -id => $contig,
127 :     -accession_number => $contig,
128 :     -display_name => $gs,
129 :     -seq => $seq,
130 :     -species => $species,
131 : olson 1.10 );
132 : olson 1.19 $bio->{$contig}->desc($gs);
133 :    
134 : olson 1.10 my $feature = Bio::SeqFeature::Generic->new(
135 :     -start => 1,
136 :     -end => $fig->contig_ln($genome, $contig),
137 : olson 1.19 -tag => { db_xref => "taxon:$taxid",
138 : olson 1.10 organism => $gs,
139 :     mol_type => "genomic DNA",
140 :     genome_id => $genome,
141 :     project => join(" ", @project),
142 :     genome_md5 => $md5 },
143 :     -primary => "source" );
144 :     $bio->{$contig}->add_SeqFeature($feature);
145 : paczian 1.4 }
146 : olson 1.10
147 :     # get the functional role name -> GO file
148 :     open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";
149 :     my $fr2go = {};
150 :     while (<FH>) {
151 :     chomp;
152 :     my ($fr, $go) = split /\t/;
153 :     $fr2go->{$fr} = [] unless (exists $fr2go->{$fr});
154 :     push @{$fr2go->{$fr}}, $go;
155 : paarmann 1.7 }
156 : olson 1.10 close FH;
157 :    
158 : olson 1.14 my @all_features = sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome);
159 :     my $feature_fns = $fig->function_of_bulk(\@all_features);
160 :     my @all_locations = $fig->feature_location_bulk(\@all_features);
161 :     my %all_locations;
162 :     $all_locations{$_->[0]} = $_->[1] foreach @all_locations;
163 :     my $all_aliases = $fig->feature_aliases_bulk(\@all_features);
164 :    
165 : olson 1.10 # get the pegs
166 : olson 1.14 foreach my $peg (@all_features) {
167 : olson 1.10 my $note;
168 : olson 1.14 # my $func = $fig->function_of($peg, $user);
169 :     my $func = $feature_fns->{$peg};
170 : olson 1.10
171 :     my %ecs;
172 :     my @gos;
173 :    
174 :     # get EC / GO from role
175 :     if (defined $func) {
176 :     foreach my $role ($fig->roles_of_function($func)) {
177 :     my ($ec) = ($role =~ /\(EC (\d+\.\d+\.\d+\.\d+)\)/);
178 :     $ecs{$ec} = 1 if ($ec);
179 :     push @gos, @{$fr2go->{$role}} if ($fr2go->{$role});
180 :     }
181 :     }
182 : olson 1.14
183 :     push @{$note->{db_xref}}, "SEED:$peg";
184 : olson 1.19 push @{$note->{transl_table}}, $gc;
185 : olson 1.10
186 :     # remove duplicate gos
187 :     my %gos = map { $_ => 1 } @gos if (scalar(@gos));
188 :    
189 :     # add GOs
190 :     push @{$note->{"db_xref"}}, @gos;
191 :    
192 :     # add ECs
193 :     push @{$note->{"EC_number"}}, keys(%ecs);
194 :    
195 :     # get the aliases from principal id
196 :     my $pid = $fig->maps_to_id($peg);
197 :     my @rw_aliases = map { $fig->rewrite_db_xrefs($_->[0]) } $fig->mapped_prot_ids($pid);
198 :     my @aliases;
199 :     foreach my $a (@rw_aliases) {
200 :     push @{$note->{"db_xref"}}, $a if ( $a );
201 :     }
202 :    
203 :     # get the links
204 : olson 1.14 foreach my $ln ($fig->fid_links($peg, $all_aliases->{$peg})) {
205 : olson 1.10 my ($db, $id);
206 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
207 :     $db = "HarvardClone";
208 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
209 :     ($db, $id) = ($1, $2);
210 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
211 :     ($db, $id) = ($1, $2);
212 :     }
213 :    
214 :     $db =~ s/\://;
215 :     if (!$db && !$id) {
216 :     print STDERR "Ignored link: $ln\n";
217 :     next;
218 :     }
219 :     push @{$note->{"db_xref"}}, "$db:$id";
220 :     }
221 :    
222 :     # add FIG id as a note
223 :     # push @{$note->{"db_xref"}}, "FIG_ID:$peg";
224 :    
225 :     # get the features
226 :    
227 :     my $loc_obj;
228 : olson 1.14 # my @location = $fig->feature_location($peg);
229 :     my @location = split(/,/, $all_locations{$peg});
230 : olson 1.10 my @loc_info;
231 :     my $contig;
232 :     foreach my $loc (@location) {
233 :     my($start, $stop);
234 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
235 :     ($contig, $start, $stop) = ($1, $2, $3);
236 :     my $original_contig = $contig;
237 : paczian 1.12 my $strand = '+';
238 : olson 1.10 my $frame = $start % 3;
239 :     if ($start > $stop) {
240 :     $frame = $stop % 3;
241 : paczian 1.12 ($start, $stop, $strand) = ($stop, $start, '-');
242 : olson 1.10 } elsif ($start == $stop) {
243 :     $strand = ".";
244 :     $frame = ".";
245 :     }
246 : paarmann 1.8
247 : olson 1.10 push(@loc_info, [$contig, $start, $stop, $strand, $frame]);
248 :    
249 :     my $sloc = new Bio::Location::Simple(-start => $start,
250 :     -end => $stop,
251 :     -strand => $strand);
252 :     if (@location == 1)
253 :     {
254 :     $loc_obj = $sloc;
255 :     }
256 :     else
257 :     {
258 :     $loc_obj = new Bio::Location::Split() if !$loc_obj;
259 :     $loc_obj->add_sub_Location($sloc);
260 :     }
261 :     }
262 :    
263 :     my $source = "FIG";
264 :     my $type = $fig->ftype($peg);
265 :     my $feature;
266 :    
267 :     # strip EC from function
268 :     my $func_ok = $func;
269 :     if ($strip_ec) {
270 :     $func_ok =~ s/\s+\(EC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
271 :     $func_ok =~ s/\s+\(TC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
272 : paczian 1.1 }
273 : olson 1.10
274 :     if ($type eq "peg") {
275 :     $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
276 :     -primary => 'CDS',
277 :     -tag => {
278 :     product => $func_ok,
279 :     translation => $fig->get_translation($peg),
280 :     },
281 :     );
282 :    
283 :     foreach my $tagtype (keys %$note) {
284 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
285 :     }
286 : paarmann 1.7
287 : olson 1.10 # work around to get annotations into gff
288 :     # this is probably still wrong for split locations.
289 : paczian 1.12 $func_ok =~ s/ #.+//;
290 :     $func_ok =~ s/;/%3B/g;
291 : olson 1.18 $func_ok =~ s/,/%2C/g;
292 : paczian 1.12 $func_ok =~ s/=//g;
293 : olson 1.10 for my $l (@loc_info)
294 :     {
295 : paczian 1.12 my $ec = "";
296 :     my @ecs = ($func =~ /[\(\[]*EC[\s:]?(\d+\.[\d-]+\.[\d-]+\.[\d-]+)[\)\]]*/ig);
297 :     if (scalar(@ecs)) {
298 :     $ec = ";Ontology_term=".join(',', map { "KEGG_ENZYME:" . $_ } @ecs);
299 :     }
300 :     my($contig, $start, $stop, $strand, $frame) = @$l;
301 :     push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\tID=".$peg.";Name=".$func_ok.$ec."\n";
302 : olson 1.10 }
303 :    
304 :    
305 :     } elsif ($type eq "rna") {
306 :     my $primary;
307 :     if ( $func =~ /tRNA/ ) {
308 :     $primary = 'tRNA';
309 :     } elsif ( $func =~ /(Ribosomal RNA|5S RNA)/ ) {
310 :     $primary = 'rRNA';
311 :     } else {
312 :     $primary = 'RNA';
313 :     }
314 : paarmann 1.7
315 : olson 1.10 $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
316 :     -primary => $primary,
317 :     -tag => {
318 :     product => $func,
319 :     },
320 :    
321 :     );
322 : paczian 1.12 $func_ok =~ s/ #.+//;
323 :     $func_ok =~ s/;/%3B/g;
324 :     $func_ok =~ s/,/2C/g;
325 :     $func_ok =~ s/=//g;
326 : olson 1.10 foreach my $tagtype (keys %$note) {
327 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
328 :    
329 :     # work around to get annotations into gff
330 :     for my $l (@loc_info)
331 :     {
332 :     my($contig, $start, $stop, $strand, $frame) = @$l;
333 : paczian 1.12 push @$gff_export, "$contig\t$source\t$primary\t$start\t$stop\t.\t$strand\t.\tID=$peg;Name=$func_ok\n";
334 : olson 1.10 }
335 :     }
336 : paczian 1.1
337 : olson 1.10 } else {
338 :     print STDERR "unhandled feature type: $type\n";
339 : paczian 1.1 }
340 : olson 1.14
341 :     my $bc = $bio->{$contig};
342 :     if (ref($bc))
343 :     {
344 :     $bc->add_SeqFeature($feature);
345 :     }
346 :     else
347 :     {
348 :     print STDERR "No contig found for $contig on $feature\n";
349 :     }
350 : paarmann 1.7 }
351 : olson 1.10
352 :    
353 :     # generate filename
354 : olson 1.11 if (!$filename)
355 :     {
356 :     $filename = $directory . $genome . "." . $format_ending;
357 :     }
358 : paarmann 1.7
359 : olson 1.10 # check for FeatureIO or SeqIO
360 :     if ($format eq "GTF") {
361 :     #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
362 :     #foreach my $feature (@$bio2) {
363 :     #$fio->write_feature($feature);
364 :     #}
365 :     open (GTF, ">$filename") or die "Cannot open file $filename.";
366 :     print GTF "##gff-version 3\n";
367 :     foreach (@$gff_export) {
368 :     print GTF $_;
369 :     }
370 :     close(GTF);
371 :    
372 :     } else {
373 : olson 1.13 # my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
374 :     #
375 :     # bioperl writes lowercase dna. We want uppercase for biophython happiness.
376 :     #
377 :     my $sio = Bio::SeqIO->new(-file => "| sed '/^LOCUS/s/dna/DNA/' >$filename", -format => $format);
378 : olson 1.10 foreach my $seq (keys %$bio) {
379 :     $sio->write_seq($bio->{$seq});
380 :     }
381 : paczian 1.1 }
382 : paczian 1.2
383 : olson 1.10 return ($filename, "Output file successfully written.");
384 : paczian 1.1 }
385 : paarmann 1.7
386 : olson 1.15 sub export_fids_as_GB
387 :     {
388 :     my($fig, $fids, $file, $strip_ec) = @_;
389 : paarmann 1.7
390 : olson 1.15 my @locs = $fig->feature_location_bulk($fids);
391 :     my $all_aliases = $fig->feature_aliases_bulk($fids);
392 :    
393 :     my %contigs;
394 :     my %all_locations;
395 :     for my $loct (@locs)
396 :     {
397 :     my($fid, $loc) = @$loct;
398 :     my $genome = FIG::genome_of($fid);
399 :     my ($contig,$beg,$end) = $fig->boundaries_of($loc);
400 :     $contigs{$genome, $contig}++;
401 :     $all_locations{$fid} = $loc;
402 :     }
403 :    
404 :     my $bio;
405 :    
406 :     for my $ckey (keys %contigs)
407 :     {
408 :     my($genome, $contig) = split(/$;/, $ckey);
409 :     my $gs = $fig->genus_species($genome);
410 :    
411 :     my $taxid = $genome;
412 :     $taxid =~ s/\.\d+$//;
413 :    
414 :     my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
415 :    
416 :     # my $seq = $fig->dna_seq($genome, $cloc);
417 :     my $seq = '';
418 :     $seq =~ s/>//;
419 :     $bio->{$contig} = Bio::Seq->new( -id => $contig,
420 :     -seq => $seq,
421 :     );
422 :     $bio->{$contig}->desc("Contig $contig from $gs");
423 :    
424 :     my $feature = Bio::SeqFeature::Generic->new(
425 :     -start => 1,
426 :     -end => $fig->contig_ln($genome, $contig),
427 : olson 1.19 -tag => { db_xref => "taxon:$taxid",
428 : olson 1.15 organism => $gs,
429 :     mol_type => "genomic DNA",
430 :     genome_id => $genome,
431 :     },
432 :     -primary => "source" );
433 :     $bio->{$contig}->add_SeqFeature($feature);
434 :     }
435 :    
436 :     my $feature_fns = $fig->function_of_bulk($fids);
437 :    
438 :     # get the pegs
439 :     foreach my $peg (@$fids)
440 :     {
441 :     my $note;
442 :     my $func = $feature_fns->{$peg};
443 :    
444 :     my %ecs;
445 :     my @gos;
446 :    
447 :     push @{$note->{db_xref}}, "SEED:$peg";
448 :    
449 :     # remove duplicate gos
450 :     my %gos = map { $_ => 1 } @gos if (scalar(@gos));
451 :    
452 :     # add GOs
453 :     push @{$note->{"db_xref"}}, @gos;
454 :    
455 :     # add ECs
456 :     push @{$note->{"EC_number"}}, keys(%ecs);
457 :    
458 :     # get the aliases from principal id
459 :     my $pid = $fig->maps_to_id($peg);
460 :     my @rw_aliases = map { $fig->rewrite_db_xrefs($_->[0]) } $fig->mapped_prot_ids($pid);
461 :     my @aliases;
462 :     foreach my $a (@rw_aliases) {
463 :     push @{$note->{"db_xref"}}, $a if ( $a );
464 :     }
465 :    
466 :     # get the links
467 :     foreach my $ln ($fig->fid_links($peg, $all_aliases->{$peg})) {
468 :     my ($db, $id);
469 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
470 :     $db = "HarvardClone";
471 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
472 :     ($db, $id) = ($1, $2);
473 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
474 :     ($db, $id) = ($1, $2);
475 :     }
476 :    
477 :     $db =~ s/\://;
478 :     if (!$db && !$id) {
479 :     print STDERR "Ignored link: $ln\n";
480 :     next;
481 :     }
482 :     push @{$note->{"db_xref"}}, "$db:$id";
483 :     }
484 :    
485 :     # add FIG id as a note
486 :     # push @{$note->{"db_xref"}}, "FIG_ID:$peg";
487 :    
488 :     # get the features
489 :    
490 :     my $loc_obj;
491 :     my @location = split(/,/, $all_locations{$peg});
492 :     my @loc_info;
493 :     my $contig;
494 :     foreach my $loc (@location) {
495 :     my($start, $stop);
496 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
497 :     ($contig, $start, $stop) = ($1, $2, $3);
498 :     my $original_contig = $contig;
499 :     my $strand = '+';
500 :     my $frame = $start % 3;
501 :     if ($start > $stop) {
502 :     $frame = $stop % 3;
503 :     ($start, $stop, $strand) = ($stop, $start, '-');
504 :     } elsif ($start == $stop) {
505 :     $strand = ".";
506 :     $frame = ".";
507 :     }
508 :    
509 :     push(@loc_info, [$contig, $start, $stop, $strand, $frame]);
510 :    
511 :     my $sloc = new Bio::Location::Simple(-start => $start,
512 :     -end => $stop,
513 :     -strand => $strand);
514 :     if (@location == 1)
515 :     {
516 :     $loc_obj = $sloc;
517 :     }
518 :     else
519 :     {
520 :     $loc_obj = new Bio::Location::Split() if !$loc_obj;
521 :     $loc_obj->add_sub_Location($sloc);
522 :     }
523 :     }
524 :    
525 :     my $source = "FIG";
526 :     my $type = $fig->ftype($peg);
527 :     my $feature;
528 :    
529 :     # strip EC from function
530 :     my $func_ok = $func;
531 :     if ($strip_ec) {
532 :     $func_ok =~ s/\s+\(EC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
533 :     $func_ok =~ s/\s+\(TC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
534 :     }
535 :    
536 :     if ($type eq "peg") {
537 :     $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
538 :     -primary => 'CDS',
539 :     -tag => {
540 :     product => $func_ok,
541 :     translation => $fig->get_translation($peg),
542 :     },
543 :     );
544 :    
545 :     foreach my $tagtype (keys %$note) {
546 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
547 :     }
548 :    
549 :     # work around to get annotations into gff
550 :     # this is probably still wrong for split locations.
551 :     $func_ok =~ s/ #.+//;
552 :     $func_ok =~ s/;/%3B/g;
553 :     $func_ok =~ s/,/2C/g;
554 :     $func_ok =~ s/=//g;
555 :     for my $l (@loc_info)
556 :     {
557 :     my $ec = "";
558 :     my @ecs = ($func =~ /[\(\[]*EC[\s:]?(\d+\.[\d-]+\.[\d-]+\.[\d-]+)[\)\]]*/ig);
559 :     if (scalar(@ecs)) {
560 :     $ec = ";Ontology_term=".join(',', map { "KEGG_ENZYME:" . $_ } @ecs);
561 :     }
562 :     my($contig, $start, $stop, $strand, $frame) = @$l;
563 :     }
564 :    
565 :    
566 :     } elsif ($type eq "rna") {
567 :     my $primary;
568 :     if ( $func =~ /tRNA/ ) {
569 :     $primary = 'tRNA';
570 :     } elsif ( $func =~ /(Ribosomal RNA|5S RNA)/ ) {
571 :     $primary = 'rRNA';
572 :     } else {
573 :     $primary = 'RNA';
574 :     }
575 :    
576 :     $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
577 :     -primary => $primary,
578 :     -tag => {
579 :     product => $func,
580 :     },
581 :    
582 :     );
583 :     $func_ok =~ s/ #.+//;
584 :     $func_ok =~ s/;/%3B/g;
585 :     $func_ok =~ s/,/2C/g;
586 :     $func_ok =~ s/=//g;
587 :     foreach my $tagtype (keys %$note) {
588 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
589 :    
590 :     # work around to get annotations into gff
591 :     for my $l (@loc_info)
592 :     {
593 :     my($contig, $start, $stop, $strand, $frame) = @$l;
594 :     }
595 :     }
596 :    
597 :     } else {
598 :     print STDERR "unhandled feature type: $type\n";
599 :     }
600 :    
601 :     my $bc = $bio->{$contig};
602 :     if (ref($bc))
603 :     {
604 :     $bc->add_SeqFeature($feature);
605 :     }
606 :     else
607 :     {
608 :     print STDERR "No contig found for $contig on $feature\n";
609 :     }
610 :     }
611 :    
612 :     my @file;
613 :     if (ref($file))
614 :     {
615 :     @file = (-file => $file);
616 :     }
617 :     elsif ($file)
618 :     {
619 :     @file = (-file => ">$file");
620 :     }
621 :     else
622 :     {
623 :     @file = (-file => \*STDOUT);
624 :     }
625 :     print STDERR Dumper(\@file);
626 :     my $sio = Bio::SeqIO->new(-file => ">$file", -format => 'genbank');
627 :     foreach my $seq (keys %$bio) {
628 :     $sio->write_seq($bio->{$seq});
629 :     }
630 :     }
631 : paarmann 1.7
632 : olson 1.17

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3