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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (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 :     my ($parameters) = @_;
19 :    
20 :     # 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 : paczian 1.2 $format_ending = "gbk";
33 : paczian 1.1 } elsif ($format eq "GTF") {
34 : paczian 1.2 $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 : paczian 1.1 } else {
41 : paczian 1.2 die "unknown export format: $format\n";
42 : paczian 1.1 }
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 : paczian 1.2 my $bio2;
82 : paarmann 1.5 my $gff_export;
83 : paczian 1.1
84 :     # 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 : paczian 1.2 $seq =~ s/>//;
89 : paczian 1.1 $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 : paczian 1.4 # get the functional role name -> GO file
105 :     open(FH, $FIG_Config::data . "/Ontologies/GO/fr2go") or die "could not open fr2go";
106 : paarmann 1.6 my $fr2go = {};
107 : paczian 1.4 while (<FH>) {
108 :     chomp;
109 :     my ($fr, $go) = split /\t/;
110 : paarmann 1.6 $fr2go->{$fr} = [] unless (exists $fr2go->{$fr});
111 :     push @{$fr2go->{$fr}}, $go;
112 : paczian 1.4 }
113 :     close FH;
114 :    
115 : paczian 1.1 # 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 : paczian 1.4
120 :     # check for a GO number
121 :     if (exists($fr2go->{$func})) {
122 : paarmann 1.6 push @{$note->{"Dbxref"}}, @{$fr2go->{$func}};
123 : paczian 1.4 }
124 : paczian 1.1
125 :     # parse out the ec number if there is one
126 :     my $ec;
127 :     if ($func =~ /E\.*C\.*\s+(\d+|-)\.(\d+|-)\.(\d+|-)\.(\d+|-)/) {
128 : paczian 1.3 push @{$note->{"EC_number"}}, "$1.$2.$3.$4";
129 : paczian 1.1 }
130 :    
131 :     # check the aliases
132 :     foreach my $alias ($fig->feature_aliases($peg))
133 :     {
134 : paczian 1.3 if ($alias =~ /^NP_/ || $alias =~ /YP_/) {
135 :     push @{$note->{"Dbxref"}}, "genpept:$alias"
136 :     }
137 : paczian 1.1 elsif ($alias =~ /gi\|/)
138 :     {
139 :     $alias =~ s/^gi\|//;
140 : paczian 1.3 push @{$note->{"Dbxref"}}, "GI:$alias";
141 : paczian 1.1 }
142 :     elsif ($alias =~ /^kegg/i)
143 :     {
144 :     $alias =~ s/kegg\|/KEGG:/i;
145 :     $alias =~ s/^(.*):/$1+/;
146 :     push @{$note->{"Dbxref"}}, $alias
147 :     }
148 :     elsif ($alias =~ /^uni/)
149 :     {
150 :     $alias =~ s/uni\|/UniProt:/;
151 : paczian 1.3 push @{$note->{"Dbxref"}}, "UniProtKB/TrEMBL:$alias";
152 : paczian 1.1 }
153 :     elsif ($alias =~ /^sp\|/)
154 :     {
155 :     $alias =~ s/sp\|/Swiss-Prot:/;
156 : paczian 1.3 push @{$note->{"Dbxref"}}, "UniProtKB/Swiss-Prot:$alias";
157 : paczian 1.1 }
158 :     elsif ($alias =~ /^[a-zA-Z][a-zA-Z0-9]{2,}\_[a-zA-Z]*\d+$/)
159 :     {
160 :     # that should be the regexp for a valid locus tag. Starts with a letter,
161 :     # has at least three alphanumerics before the _
162 :     # then has a series of optional letters and ends with numbers
163 :     push @{$note->{"locus_tag"}}, $alias;
164 :     }
165 :     elsif ($alias =~ /^.{4,6}$/)
166 :     {
167 :     push @{$note->{"gene_symbol"}}, $alias;
168 :     }
169 :     elsif ($alias =~ /^[a-zA-Z]+\d+$/)
170 :     {
171 :     push @{$note->{"locus"}}, $alias;
172 :     }
173 :     else {
174 :     # don't know what it is so keep it as an alias
175 :     push @{$note->{"Alias"}}, $alias;
176 :     }
177 :     }
178 :    
179 :     # get the links
180 :     foreach my $ln ($fig->fid_links($peg))
181 :     {
182 :     my ($db, $id);
183 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
184 :     $db = "HarvardClone";
185 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
186 :     ($db, $id) = ($1, $2);
187 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
188 :     ($db, $id) = ($1, $2);
189 :     }
190 :    
191 :     $db =~ s/\://;
192 :     if (!$db && !$id) {
193 :     print STDERR "Ignored link: $ln\n";
194 :     next;
195 :     }
196 :     push @{$note->{"Dbxref"}}, "$db:$id";
197 :     }
198 :    
199 :     # add FIG id as a note
200 :     push @{$note->{"Dbxref"}}, "FIG_ID:$peg";
201 :    
202 :     # get the features
203 :     my @location = $fig->feature_location($peg);
204 :     foreach my $loc (@location) {
205 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
206 : paczian 1.2 my ($contig, $start, $stop) = ($1, $2, $3);
207 :     my $original_contig = $contig;
208 :     my $strand = '1';
209 :     my $frame = $start % 3;
210 :     if ($start > $stop) {
211 :     $frame = $stop % 3;
212 :     ($start, $stop, $strand) = ($stop, $start, '-1');
213 :     }
214 :     elsif ($start == $stop) {
215 :     $strand = ".";
216 :     $frame = ".";
217 :     }
218 :     my $source = "FIG";
219 : paczian 1.1 my $type = $fig->ftype($peg);
220 :     my $feature;
221 :     if ($type eq "peg") {
222 :     $feature = Bio::SeqFeature::Generic->new(
223 :     -start => $start,
224 :     -end => $stop,
225 :     -strand => $strand,
226 :     -primary => 'CDS',
227 : paczian 1.2 -tag => {
228 :     product => $func,
229 :     translation => $fig->get_translation($peg),
230 :     },
231 : paczian 1.1 );
232 :    
233 :     foreach my $tagtype (keys %$note) {
234 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
235 :     }
236 : paczian 1.2
237 :     my $feature2 = Bio::SeqFeature::Annotated->new( -start => $start,
238 :     -end => $stop,
239 :     -strand => $strand,
240 :     -phase => 0,
241 :     -frame => $frame,
242 :     -source => $source,
243 :     -type => "CDS",
244 :     -seq_id => $contig );
245 :    
246 :     push(@$bio2, $feature2);
247 : paarmann 1.5
248 :     # work around to get annotations into gff
249 :     push @$gff_export, "$contig\t$source\tCDS\t$start\t$stop\t.\t$strand\t$frame\t$func\n"
250 :    
251 : paczian 1.2
252 : paczian 1.1 } elsif ($type eq "rna") {
253 :     $feature = Bio::SeqFeature::Generic->new(
254 :     -start => $start,
255 :     -end => $stop,
256 :     -strand => $strand,
257 :     -primary => 'RNA',
258 :     );
259 :     foreach my $tagtype (keys %$note) {
260 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
261 :     }
262 :    
263 :     } else {
264 :     print STDERR "unhandled feature type: $type\n";
265 :     }
266 :    
267 :     $bio->{$contig}->add_SeqFeature($feature);
268 :     }
269 :     }
270 :    
271 :     # generate filename
272 :     my $filename = $directory . $genome . "." . $format_ending;
273 : paczian 1.2
274 : paczian 1.1 # check for FeatureIO or SeqIO
275 :     if ($format eq "GTF") {
276 : paarmann 1.5 #my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
277 :     #foreach my $feature (@$bio2) {
278 :     #$fio->write_feature($feature);
279 :     #}
280 :     open (GTF, ">$filename") or die "Cannot open file $filename.";
281 :     print GTF "##gff-version 3\n";
282 :     foreach (@$gff_export) {
283 :     print GTF $_;
284 : paczian 1.2 }
285 : paarmann 1.5 close(GTF);
286 :    
287 : paczian 1.1 } else {
288 : paczian 1.2 my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
289 :     foreach my $seq (keys %$bio) {
290 :     $sio->write_seq($bio->{$seq});
291 :     }
292 : paczian 1.1 }
293 : paczian 1.2
294 : paczian 1.1 return ($filename, "Output file successfully written.");
295 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3