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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3