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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : paczian 1.1 package SeedExport;
2 :    
3 :     1;
4 :    
5 :     use strict;
6 :     use FIG;
7 : paczian 1.4 use FIG_Config;
8 : paczian 1.1 use FIGV;
9 : paczian 1.2 use URI::Escape;
10 : paczian 1.1 use Bio::FeatureIO;
11 :     use Bio::SeqIO;
12 :     use Bio::Seq;
13 : olson 1.10 use Bio::Location::Split;
14 :     use Bio::Location::Simple;
15 : paczian 1.1 use Bio::SeqFeature::Generic;
16 : paczian 1.2 use Bio::SeqFeature::Annotated;
17 : paczian 1.1 use File::Basename;
18 :    
19 :     sub export {
20 : olson 1.10 my ($parameters) = @_;
21 :    
22 :     # get parameters
23 :     my $virt_genome_dir = $parameters->{'virtual_genome_directory'};
24 :     my $genome = $parameters->{'genome'};
25 :     my $directory = $parameters->{'directory'};
26 :     my $format = $parameters->{'export_format'};
27 :     my $strip_ec = $parameters->{'strip_ec'};
28 :    
29 :     # set some variables
30 :     my $user = "master:master";
31 :     my $format_ending;
32 :    
33 :     # check format
34 :     if ($format eq "genbank") {
35 :     $format_ending = "gbk";
36 :     } elsif ($format eq "GTF") {
37 :     $format_ending = "gtf";
38 :     } elsif ($format eq "gff") {
39 :     $format = "GTF";
40 :     $format_ending = "gff";
41 :     } elsif ($format eq "embl") {
42 :     $format_ending = "embl";
43 :     } else {
44 :     die "unknown export format: $format\n";
45 :     }
46 :    
47 :     # initialize fig object
48 :     my $fig;
49 :     my $virt_genome;
50 :     if ($virt_genome_dir) {
51 :     $fig = new FIGV($virt_genome_dir);
52 :     $virt_genome = basename($virt_genome_dir);
53 :     } else {
54 :     $fig = new FIG;
55 :     }
56 :    
57 :     # check for genus species
58 :     my $gs = $fig->genus_species($genome);
59 :     unless ($gs) {
60 :     return (undef, "Couldn't recognize $genome.");
61 :     }
62 :    
63 :     # get taxonomy id and taxonomy
64 :     my $taxid = $genome;
65 :     $taxid =~ s/\.\d+$//;
66 :     my $taxonomy = $fig->taxonomy_of($genome);
67 :    
68 :     # get project
69 :     if ($virt_genome_dir and $genome eq $virt_genome) {
70 :     open(PROJECT, "<$virt_genome_dir/PROJECT") or die "Error opening $virt_genome_dir/PROJECT: $!\n";
71 :     } else {
72 :     open( PROJECT, "<$FIG_Config::organisms/$genome/PROJECT" ) or die "Error opening $FIG_Config::organisms/$genome/PROJECT: $!\n";
73 : paczian 1.1 }
74 : olson 1.10 my @project = <PROJECT>;
75 :     chomp(@project);
76 :     close PROJECT;
77 :     map {s/^\<a href\=\".*?\"\>(.*)\<\/a\>/$1/i} @project;
78 :    
79 :     # get md5 checksum
80 :     my $md5 = $fig->genome_md5sum($genome);
81 :    
82 :     # create the variable for the bio object
83 :     my $bio;
84 :     my $bio2;
85 :     my $gff_export;
86 :    
87 :     # get the contigs
88 :     foreach my $contig ($fig->contigs_of($genome)) {
89 :     my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
90 :     my $seq = $fig->dna_seq($genome, $cloc);
91 :     $seq =~ s/>//;
92 :     $bio->{$contig} = Bio::Seq->new( -id => $contig,
93 :     -seq => $seq,
94 :     );
95 :     $bio->{$contig}->desc("Contig $contig from $gs");
96 :    
97 :     my $feature = Bio::SeqFeature::Generic->new(
98 :     -start => 1,
99 :     -end => $fig->contig_ln($genome, $contig),
100 :     -tag => { db_xref => "taxon: $taxid",
101 :     organism => $gs,
102 :     mol_type => "genomic DNA",
103 :     genome_id => $genome,
104 :     project => join(" ", @project),
105 :     genome_md5 => $md5 },
106 :     -primary => "source" );
107 :     $bio->{$contig}->add_SeqFeature($feature);
108 : paczian 1.4 }
109 : olson 1.10
110 :     # get the functional role name -> GO file
111 :     open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";
112 :     my $fr2go = {};
113 :     while (<FH>) {
114 :     chomp;
115 :     my ($fr, $go) = split /\t/;
116 :     $fr2go->{$fr} = [] unless (exists $fr2go->{$fr});
117 :     push @{$fr2go->{$fr}}, $go;
118 : paarmann 1.7 }
119 : olson 1.10 close FH;
120 :    
121 :     # get the pegs
122 :     foreach my $peg (sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome)) {
123 :     my $note;
124 :     my $func = $fig->function_of($peg, $user);
125 :    
126 :     my %ecs;
127 :     my @gos;
128 :    
129 :     # get EC / GO from role
130 :     if (defined $func) {
131 :     foreach my $role ($fig->roles_of_function($func)) {
132 :     my ($ec) = ($role =~ /\(EC (\d+\.\d+\.\d+\.\d+)\)/);
133 :     $ecs{$ec} = 1 if ($ec);
134 :     push @gos, @{$fr2go->{$role}} if ($fr2go->{$role});
135 :     }
136 :     }
137 :    
138 :     # remove duplicate gos
139 :     my %gos = map { $_ => 1 } @gos if (scalar(@gos));
140 :    
141 :     # add GOs
142 :     push @{$note->{"db_xref"}}, @gos;
143 :    
144 :     # add ECs
145 :     push @{$note->{"EC_number"}}, keys(%ecs);
146 :    
147 :     # get the aliases from principal id
148 :     my $pid = $fig->maps_to_id($peg);
149 :     my @rw_aliases = map { $fig->rewrite_db_xrefs($_->[0]) } $fig->mapped_prot_ids($pid);
150 :     my @aliases;
151 :     foreach my $a (@rw_aliases) {
152 :     push @{$note->{"db_xref"}}, $a if ( $a );
153 :     }
154 :    
155 :     # get the links
156 :     foreach my $ln ($fig->fid_links($peg)) {
157 :     my ($db, $id);
158 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
159 :     $db = "HarvardClone";
160 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
161 :     ($db, $id) = ($1, $2);
162 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
163 :     ($db, $id) = ($1, $2);
164 :     }
165 :    
166 :     $db =~ s/\://;
167 :     if (!$db && !$id) {
168 :     print STDERR "Ignored link: $ln\n";
169 :     next;
170 :     }
171 :     push @{$note->{"db_xref"}}, "$db:$id";
172 :     }
173 :    
174 :     # add FIG id as a note
175 :     # push @{$note->{"db_xref"}}, "FIG_ID:$peg";
176 :    
177 :     # get the features
178 :    
179 :     my $loc_obj;
180 :     my @location = $fig->feature_location($peg);
181 :     my @loc_info;
182 :     my $contig;
183 :     foreach my $loc (@location) {
184 :     my($start, $stop);
185 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
186 :     ($contig, $start, $stop) = ($1, $2, $3);
187 :     my $original_contig = $contig;
188 :     my $strand = '1';
189 :     my $frame = $start % 3;
190 :     if ($start > $stop) {
191 :     $frame = $stop % 3;
192 :     ($start, $stop, $strand) = ($stop, $start, '-1');
193 :     } elsif ($start == $stop) {
194 :     $strand = ".";
195 :     $frame = ".";
196 :     }
197 : paarmann 1.8
198 : olson 1.10 push(@loc_info, [$contig, $start, $stop, $strand, $frame]);
199 :    
200 :     my $sloc = new Bio::Location::Simple(-start => $start,
201 :     -end => $stop,
202 :     -strand => $strand);
203 :     if (@location == 1)
204 :     {
205 :     $loc_obj = $sloc;
206 :     }
207 :     else
208 :     {
209 :     $loc_obj = new Bio::Location::Split() if !$loc_obj;
210 :     $loc_obj->add_sub_Location($sloc);
211 :     }
212 :     }
213 :    
214 :     my $source = "FIG";
215 :     my $type = $fig->ftype($peg);
216 :     my $feature;
217 :    
218 :     # strip EC from function
219 :     my $func_ok = $func;
220 :     if ($strip_ec) {
221 :     $func_ok =~ s/\s+\(EC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
222 :     $func_ok =~ s/\s+\(TC \d+\.(\d+|-)\.(\d+|-)\.(\d+|-)\)//g;
223 : paczian 1.1 }
224 : olson 1.10
225 :     if ($type eq "peg") {
226 :     $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
227 :     -primary => 'CDS',
228 :     -tag => {
229 :     product => $func_ok,
230 :     translation => $fig->get_translation($peg),
231 :     },
232 :     );
233 :    
234 :     foreach my $tagtype (keys %$note) {
235 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
236 :     }
237 : paarmann 1.7
238 : olson 1.10 # work around to get annotations into gff
239 :     # this is probably still wrong for split locations.
240 :     for my $l (@loc_info)
241 :     {
242 :     my($contig, $start, $stop, $strand, $frame) = @$l;
243 :     push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\t$func\n";
244 :     }
245 :    
246 :    
247 :     } elsif ($type eq "rna") {
248 :     my $primary;
249 :     if ( $func =~ /tRNA/ ) {
250 :     $primary = 'tRNA';
251 :     } elsif ( $func =~ /(Ribosomal RNA|5S RNA)/ ) {
252 :     $primary = 'rRNA';
253 :     } else {
254 :     $primary = 'RNA';
255 :     }
256 : paarmann 1.7
257 : olson 1.10 $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
258 :     -primary => $primary,
259 :     -tag => {
260 :     product => $func,
261 :     },
262 :    
263 :     );
264 :     foreach my $tagtype (keys %$note) {
265 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
266 :    
267 :     # work around to get annotations into gff
268 :     for my $l (@loc_info)
269 :     {
270 :     my($contig, $start, $stop, $strand, $frame) = @$l;
271 :     push @$gff_export, "$contig\t$source\t$primary\t$start\t$stop\t.\t$strand\t.\t$func\n";
272 :     }
273 :     }
274 : paczian 1.1
275 : olson 1.10 } else {
276 :     print STDERR "unhandled feature type: $type\n";
277 : paczian 1.1 }
278 : paarmann 1.7 $bio->{$contig}->add_SeqFeature($feature);
279 :     }
280 : olson 1.10
281 :    
282 :     # generate filename
283 :     my $filename = $directory . $genome . "." . $format_ending;
284 : paarmann 1.7
285 : olson 1.10 # check for FeatureIO or SeqIO
286 :     if ($format eq "GTF") {
287 :     #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
288 :     #foreach my $feature (@$bio2) {
289 :     #$fio->write_feature($feature);
290 :     #}
291 :     open (GTF, ">$filename") or die "Cannot open file $filename.";
292 :     print GTF "##gff-version 3\n";
293 :     foreach (@$gff_export) {
294 :     print GTF $_;
295 :     }
296 :     close(GTF);
297 :    
298 :     } else {
299 :     my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
300 :     foreach my $seq (keys %$bio) {
301 :     $sio->write_seq($bio->{$seq});
302 :     }
303 : paczian 1.1 }
304 : paczian 1.2
305 : olson 1.10 return ($filename, "Output file successfully written.");
306 : paczian 1.1 }
307 : paarmann 1.7
308 :    
309 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3