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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (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 :    
26 :     # set some variables
27 :     my $user = "master:master";
28 :     my $format_ending;
29 :    
30 :     # check format
31 :     if ($format eq "genbank") {
32 :     $format_ending = "gbk";
33 :     } elsif ($format eq "GTF") {
34 :     $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 {
41 :     die "unknown export format: $format\n";
42 :     }
43 :    
44 :     # initialize fig object
45 :     my $fig;
46 :     my $virt_genome;
47 :     if ($virt_genome_dir) {
48 :     $fig = new FIGV($virt_genome_dir);
49 :     $virt_genome = basename($virt_genome_dir);
50 :     } else {
51 :     $fig = new FIG;
52 :     }
53 :    
54 :     # check for genus species
55 :     my $gs = $fig->genus_species($genome);
56 :     unless ($gs) {
57 :     return (undef, "Couldn't recognize $genome.");
58 :     }
59 :    
60 :     # get taxonomy id and taxonomy
61 :     my $taxid = $genome;
62 :     $taxid =~ s/\.\d+$//;
63 :     my $taxonomy = $fig->taxonomy_of($genome);
64 :    
65 :     # get project
66 :     if ($virt_genome_dir and $genome eq $virt_genome) {
67 :     open(PROJECT, "<$virt_genome_dir/PROJECT") or die "Error opening $virt_genome_dir/PROJECT: $!\n";
68 :     } else {
69 :     open( PROJECT, "<$FIG_Config::organisms/$genome/PROJECT" ) or die "Error opening $FIG_Config::organisms/$genome/PROJECT: $!\n";
70 :     }
71 :     my @project = <PROJECT>;
72 :     chomp(@project);
73 :     close PROJECT;
74 :     map {s/^\<a href\=\".*?\"\>(.*)\<\/a\>/$1/i} @project;
75 :    
76 :     # get md5 checksum
77 :     my $md5 = $fig->genome_md5sum($genome);
78 :    
79 :     # create the variable for the bio object
80 :     my $bio;
81 :     my $bio2;
82 :     my $gff_export;
83 : paczian 1.1
84 : paarmann 1.7 # get the contigs
85 :     foreach my $contig ($fig->contigs_of($genome)) {
86 :     my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
87 :     my $seq = $fig->dna_seq($genome, $cloc);
88 :     $seq =~ s/>//;
89 :     $bio->{$contig} = Bio::Seq->new(-id => "$contig", seq => $seq);
90 :     $bio->{$contig}->desc("Contig $contig from $gs");
91 :     my $feature = Bio::SeqFeature::Generic->new(
92 :     -start => 1,
93 :     -end => $fig->contig_ln($genome, $contig),
94 :     -tag => { dbxref => "taxon: $taxid",
95 :     organism => $gs,
96 :     mol_type => "genomic DNA",
97 :     genome_id => $genome,
98 :     project => join(" ", @project),
99 :     genome_md5 => $md5 },
100 :     -primary => "source" );
101 :     $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
116 :     foreach my $peg (sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome)) {
117 :     my $note;
118 :     my $func = $fig->function_of($peg, $user);
119 :    
120 :     my %ecs;
121 :     my @gos;
122 :    
123 :     # get EC / GO from role
124 :     if (defined $func) {
125 :     foreach my $role ($fig->roles_of_function($func)) {
126 :     my ($ec) = ($role =~ /\(EC (\d+\.\d+\.\d+\.\d+)\)/);
127 :     $ecs{$ec} = 1 if ($ec);
128 :     push @gos, @{$fr2go->{$role}} if ($fr2go->{$role});
129 :     }
130 : paczian 1.1 }
131 : paarmann 1.7
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 : paczian 1.4 }
148 : paarmann 1.7
149 :     # get the links
150 :     foreach my $ln ($fig->fid_links($peg)) {
151 :     my ($db, $id);
152 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
153 :     $db = "HarvardClone";
154 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
155 :     ($db, $id) = ($1, $2);
156 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
157 :     ($db, $id) = ($1, $2);
158 :     }
159 : paczian 1.1
160 : paarmann 1.7 $db =~ s/\://;
161 :     if (!$db && !$id) {
162 :     print STDERR "Ignored link: $ln\n";
163 :     next;
164 :     }
165 :     push @{$note->{"Dbxref"}}, "$db:$id";
166 :     }
167 :    
168 :     # add FIG id as a note
169 :     push @{$note->{"Dbxref"}}, "FIG_ID:$peg";
170 :    
171 :     # get the features
172 :     my @location = $fig->feature_location($peg);
173 :     foreach my $loc (@location) {
174 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
175 :     my ($contig, $start, $stop) = ($1, $2, $3);
176 :     my $original_contig = $contig;
177 :     my $strand = '1';
178 :     my $frame = $start % 3;
179 :     if ($start > $stop) {
180 :     $frame = $stop % 3;
181 :     ($start, $stop, $strand) = ($stop, $start, '-1');
182 :     } elsif ($start == $stop) {
183 :     $strand = ".";
184 :     $frame = ".";
185 :     }
186 :     my $source = "FIG";
187 :     my $type = $fig->ftype($peg);
188 :     my $feature;
189 :     if ($type eq "peg") {
190 :     $feature = Bio::SeqFeature::Generic->new(
191 :     -start => $start,
192 :     -end => $stop,
193 :     -strand => $strand,
194 :     -primary => 'CDS',
195 :     -tag => {
196 :     product => $func,
197 :     translation => $fig->get_translation($peg),
198 :     },
199 :     );
200 :    
201 :     foreach my $tagtype (keys %$note) {
202 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
203 : paczian 1.1 }
204 : paarmann 1.7
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 : paczian 1.1
220 : paarmann 1.7 } elsif ($type eq "rna") {
221 :     $feature = Bio::SeqFeature::Generic->new(
222 :     -start => $start,
223 :     -end => $stop,
224 :     -strand => $strand,
225 :     -primary => 'RNA',
226 :     );
227 :     foreach my $tagtype (keys %$note) {
228 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
229 : paczian 1.1 }
230 : paarmann 1.7
231 :     } else {
232 :     print STDERR "unhandled feature type: $type\n";
233 :     }
234 : paczian 1.1
235 : paarmann 1.7 $bio->{$contig}->add_SeqFeature($feature);
236 :     }
237 :     }
238 :    
239 :     # generate filename
240 :     my $filename = $directory . $genome . "." . $format_ending;
241 :    
242 :     # check for FeatureIO or SeqIO
243 :     if ($format eq "GTF") {
244 :     #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
245 :     #foreach my $feature (@$bio2) {
246 :     #$fio->write_feature($feature);
247 :     #}
248 :     open (GTF, ">$filename") or die "Cannot open file $filename.";
249 :     print GTF "##gff-version 3\n";
250 :     foreach (@$gff_export) {
251 :     print GTF $_;
252 : paczian 1.1 }
253 : paarmann 1.7 close(GTF);
254 : paarmann 1.5
255 : paarmann 1.7 } else {
256 :     my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
257 :     foreach my $seq (keys %$bio) {
258 :     $sio->write_seq($bio->{$seq});
259 : paczian 1.1 }
260 : paarmann 1.7 }
261 : paczian 1.2
262 : paarmann 1.7 return ($filename, "Output file successfully written.");
263 : paczian 1.1 }
264 : paarmann 1.7
265 :    
266 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3