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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11 - (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 : olson 1.11 my $filename = $parameters->{'filename'};
29 : olson 1.10
30 :     # set some variables
31 :     my $user = "master:master";
32 :     my $format_ending;
33 :    
34 :     # check format
35 :     if ($format eq "genbank") {
36 :     $format_ending = "gbk";
37 :     } elsif ($format eq "GTF") {
38 :     $format_ending = "gtf";
39 :     } elsif ($format eq "gff") {
40 :     $format = "GTF";
41 :     $format_ending = "gff";
42 :     } elsif ($format eq "embl") {
43 :     $format_ending = "embl";
44 :     } else {
45 :     die "unknown export format: $format\n";
46 :     }
47 :    
48 :     # initialize fig object
49 :     my $fig;
50 :     my $virt_genome;
51 :     if ($virt_genome_dir) {
52 :     $fig = new FIGV($virt_genome_dir);
53 :     $virt_genome = basename($virt_genome_dir);
54 :     } else {
55 :     $fig = new FIG;
56 :     }
57 :    
58 :     # check for genus species
59 :     my $gs = $fig->genus_species($genome);
60 :     unless ($gs) {
61 :     return (undef, "Couldn't recognize $genome.");
62 :     }
63 :    
64 :     # get taxonomy id and taxonomy
65 :     my $taxid = $genome;
66 :     $taxid =~ s/\.\d+$//;
67 :     my $taxonomy = $fig->taxonomy_of($genome);
68 :    
69 :     # get project
70 :     if ($virt_genome_dir and $genome eq $virt_genome) {
71 :     open(PROJECT, "<$virt_genome_dir/PROJECT") or die "Error opening $virt_genome_dir/PROJECT: $!\n";
72 :     } else {
73 :     open( PROJECT, "<$FIG_Config::organisms/$genome/PROJECT" ) or die "Error opening $FIG_Config::organisms/$genome/PROJECT: $!\n";
74 : paczian 1.1 }
75 : olson 1.10 my @project = <PROJECT>;
76 :     chomp(@project);
77 :     close PROJECT;
78 :     map {s/^\<a href\=\".*?\"\>(.*)\<\/a\>/$1/i} @project;
79 :    
80 :     # get md5 checksum
81 :     my $md5 = $fig->genome_md5sum($genome);
82 :    
83 :     # create the variable for the bio object
84 :     my $bio;
85 :     my $bio2;
86 :     my $gff_export;
87 :    
88 :     # get the contigs
89 :     foreach my $contig ($fig->contigs_of($genome)) {
90 :     my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
91 :     my $seq = $fig->dna_seq($genome, $cloc);
92 :     $seq =~ s/>//;
93 :     $bio->{$contig} = Bio::Seq->new( -id => $contig,
94 :     -seq => $seq,
95 :     );
96 :     $bio->{$contig}->desc("Contig $contig from $gs");
97 :    
98 :     my $feature = Bio::SeqFeature::Generic->new(
99 :     -start => 1,
100 :     -end => $fig->contig_ln($genome, $contig),
101 :     -tag => { db_xref => "taxon: $taxid",
102 :     organism => $gs,
103 :     mol_type => "genomic DNA",
104 :     genome_id => $genome,
105 :     project => join(" ", @project),
106 :     genome_md5 => $md5 },
107 :     -primary => "source" );
108 :     $bio->{$contig}->add_SeqFeature($feature);
109 : paczian 1.4 }
110 : olson 1.10
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 : paarmann 1.7 }
120 : olson 1.10 close FH;
121 :    
122 :     # get the pegs
123 :     foreach my $peg (sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome)) {
124 :     my $note;
125 :     my $func = $fig->function_of($peg, $user);
126 :    
127 :     my %ecs;
128 :     my @gos;
129 :    
130 :     # get EC / GO from role
131 :     if (defined $func) {
132 :     foreach my $role ($fig->roles_of_function($func)) {
133 :     my ($ec) = ($role =~ /\(EC (\d+\.\d+\.\d+\.\d+)\)/);
134 :     $ecs{$ec} = 1 if ($ec);
135 :     push @gos, @{$fr2go->{$role}} if ($fr2go->{$role});
136 :     }
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
157 :     foreach my $ln ($fig->fid_links($peg)) {
158 :     my ($db, $id);
159 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
160 :     $db = "HarvardClone";
161 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
162 :     ($db, $id) = ($1, $2);
163 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
164 :     ($db, $id) = ($1, $2);
165 :     }
166 :    
167 :     $db =~ s/\://;
168 :     if (!$db && !$id) {
169 :     print STDERR "Ignored link: $ln\n";
170 :     next;
171 :     }
172 :     push @{$note->{"db_xref"}}, "$db:$id";
173 :     }
174 :    
175 :     # add FIG id as a note
176 :     # push @{$note->{"db_xref"}}, "FIG_ID:$peg";
177 :    
178 :     # get the features
179 :    
180 :     my $loc_obj;
181 :     my @location = $fig->feature_location($peg);
182 :     my @loc_info;
183 :     my $contig;
184 :     foreach my $loc (@location) {
185 :     my($start, $stop);
186 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
187 :     ($contig, $start, $stop) = ($1, $2, $3);
188 :     my $original_contig = $contig;
189 :     my $strand = '1';
190 :     my $frame = $start % 3;
191 :     if ($start > $stop) {
192 :     $frame = $stop % 3;
193 :     ($start, $stop, $strand) = ($stop, $start, '-1');
194 :     } elsif ($start == $stop) {
195 :     $strand = ".";
196 :     $frame = ".";
197 :     }
198 : paarmann 1.8
199 : olson 1.10 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";
216 :     my $type = $fig->ftype($peg);
217 :     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 : paczian 1.1 }
225 : olson 1.10
226 :     if ($type eq "peg") {
227 :     $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
228 :     -primary => 'CDS',
229 :     -tag => {
230 :     product => $func_ok,
231 :     translation => $fig->get_translation($peg),
232 :     },
233 :     );
234 :    
235 :     foreach my $tagtype (keys %$note) {
236 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
237 :     }
238 : paarmann 1.7
239 : olson 1.10 # work around to get annotations into gff
240 :     # this is probably still wrong for split locations.
241 :     for my $l (@loc_info)
242 :     {
243 :     my($contig, $start, $stop, $strand, $frame) = @$l;
244 :     push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\t$func\n";
245 :     }
246 :    
247 :    
248 :     } elsif ($type eq "rna") {
249 :     my $primary;
250 :     if ( $func =~ /tRNA/ ) {
251 :     $primary = 'tRNA';
252 :     } elsif ( $func =~ /(Ribosomal RNA|5S RNA)/ ) {
253 :     $primary = 'rRNA';
254 :     } else {
255 :     $primary = 'RNA';
256 :     }
257 : paarmann 1.7
258 : olson 1.10 $feature = Bio::SeqFeature::Generic->new(-location => $loc_obj,
259 :     -primary => $primary,
260 :     -tag => {
261 :     product => $func,
262 :     },
263 :    
264 :     );
265 :     foreach my $tagtype (keys %$note) {
266 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
267 :    
268 :     # work around to get annotations into gff
269 :     for my $l (@loc_info)
270 :     {
271 :     my($contig, $start, $stop, $strand, $frame) = @$l;
272 :     push @$gff_export, "$contig\t$source\t$primary\t$start\t$stop\t.\t$strand\t.\t$func\n";
273 :     }
274 :     }
275 : paczian 1.1
276 : olson 1.10 } else {
277 :     print STDERR "unhandled feature type: $type\n";
278 : paczian 1.1 }
279 : paarmann 1.7 $bio->{$contig}->add_SeqFeature($feature);
280 :     }
281 : olson 1.10
282 :    
283 :     # generate filename
284 : olson 1.11 if (!$filename)
285 :     {
286 :     $filename = $directory . $genome . "." . $format_ending;
287 :     }
288 : paarmann 1.7
289 : olson 1.10 # check for FeatureIO or SeqIO
290 :     if ($format eq "GTF") {
291 :     #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
292 :     #foreach my $feature (@$bio2) {
293 :     #$fio->write_feature($feature);
294 :     #}
295 :     open (GTF, ">$filename") or die "Cannot open file $filename.";
296 :     print GTF "##gff-version 3\n";
297 :     foreach (@$gff_export) {
298 :     print GTF $_;
299 :     }
300 :     close(GTF);
301 :    
302 :     } else {
303 :     my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
304 :     foreach my $seq (keys %$bio) {
305 :     $sio->write_seq($bio->{$seq});
306 :     }
307 : paczian 1.1 }
308 : paczian 1.2
309 : olson 1.10 return ($filename, "Output file successfully written.");
310 : paczian 1.1 }
311 : paarmann 1.7
312 :    
313 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3