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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3