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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download) (as text)

1 : paczian 1.1 package SeedExport;
2 :    
3 :     1;
4 :    
5 :     use strict;
6 :     use FIG;
7 :     use FIGV;
8 : paczian 1.2 use URI::Escape;
9 : paczian 1.1 use Bio::FeatureIO;
10 :     use Bio::SeqIO;
11 :     use Bio::Seq;
12 :     use Bio::SeqFeature::Generic;
13 : paczian 1.2 use Bio::SeqFeature::Annotated;
14 : paczian 1.1 use File::Basename;
15 :    
16 :     sub export {
17 :     my ($parameters) = @_;
18 :    
19 :     # get parameters
20 :     my $virt_genome_dir = $parameters->{'virtual_genome_directory'};
21 :     my $genome = $parameters->{'genome'};
22 :     my $directory = $parameters->{'directory'};
23 :     my $format = $parameters->{'export_format'};
24 :    
25 :     # set some variables
26 :     my $user = "master:master";
27 :     my $format_ending;
28 :    
29 :     # check format
30 :     if ($format eq "genbank") {
31 : paczian 1.2 $format_ending = "gbk";
32 : paczian 1.1 } elsif ($format eq "GTF") {
33 : paczian 1.2 $format_ending = "gtf";
34 :     } elsif ($format eq "gff") {
35 :     $format = "GTF";
36 :     $format_ending = "gff";
37 :     } elsif ($format eq "embl") {
38 :     $format_ending = "embl";
39 : paczian 1.1 } else {
40 : paczian 1.2 die "unknown export format: $format\n";
41 : paczian 1.1 }
42 :    
43 :     # initialize fig object
44 :     my $fig;
45 :     my $virt_genome;
46 :     if ($virt_genome_dir) {
47 :     $fig = new FIGV($virt_genome_dir);
48 :     $virt_genome = basename($virt_genome_dir);
49 :     } else {
50 :     $fig = new FIG;
51 :     }
52 :    
53 :     # check for genus species
54 :     my $gs = $fig->genus_species($genome);
55 :     unless ($gs) {
56 :     return (undef, "Couldn't recognize $genome.");
57 :     }
58 :    
59 :     # get taxonomy id and taxonomy
60 :     my $taxid = $genome;
61 :     $taxid =~ s/\.\d+$//;
62 :     my $taxonomy = $fig->taxonomy_of($genome);
63 :    
64 :     # get project
65 :     if ($virt_genome_dir and $genome eq $virt_genome) {
66 :     open(PROJECT, "<$virt_genome_dir/PROJECT") or die "Error opening $virt_genome_dir/PROJECT: $!\n";
67 :     } else {
68 :     open( PROJECT, "<$FIG_Config::organisms/$genome/PROJECT" ) or die "Error opening $FIG_Config::organisms/$genome/PROJECT: $!\n";
69 :     }
70 :     my @project = <PROJECT>;
71 :     chomp(@project);
72 :     close PROJECT;
73 :     map {s/^\<a href\=\".*?\"\>(.*)\<\/a\>/$1/i} @project;
74 :    
75 :     # get md5 checksum
76 :     my $md5 = $fig->genome_md5sum($genome);
77 :    
78 :     # create the variable for the bio object
79 :     my $bio;
80 : paczian 1.2 my $bio2;
81 : paczian 1.1
82 :     # get the contigs
83 :     foreach my $contig ($fig->contigs_of($genome)) {
84 :     my $cloc = $contig.'_1_'.$fig->contig_ln($genome, $contig);
85 :     my $seq = $fig->dna_seq($genome, $cloc);
86 : paczian 1.2 $seq =~ s/>//;
87 : paczian 1.1 $bio->{$contig} = Bio::Seq->new(-id => "$contig", seq => $seq);
88 :     $bio->{$contig}->desc("Contig $contig from $gs");
89 :     my $feature = Bio::SeqFeature::Generic->new(
90 :     -start => 1,
91 :     -end => $fig->contig_ln($genome, $contig),
92 :     -tag => { dbxref => "taxon: $taxid",
93 :     organism => $gs,
94 :     mol_type => "genomic DNA",
95 :     genome_id => $genome,
96 :     project => join(" ", @project),
97 :     genome_md5 => $md5 },
98 :     -primary => "source" );
99 :     $bio->{$contig}->add_SeqFeature($feature);
100 :     }
101 :    
102 :     # get the pegs
103 :     foreach my $peg (sort { &FIG::by_fig_id($a,$b) } $fig->pegs_of($genome), $fig->rnas_of($genome)) {
104 :     my $note;
105 :     my $func = $fig->function_of($peg, $user);
106 :    
107 :     # parse out the ec number if there is one
108 :     my $ec;
109 :     if ($func =~ /E\.*C\.*\s+(\d+|-)\.(\d+|-)\.(\d+|-)\.(\d+|-)/) {
110 :     push @{$note->{"ec_number"}}, "$1.$2.$3.$4";
111 :     }
112 :    
113 :     # check the aliases
114 :     foreach my $alias ($fig->feature_aliases($peg))
115 :     {
116 :     if ($alias =~ /^NP_/ || $alias =~ /YP_/) {push @{$note->{"Dbxref"}}, "NCBI_genpept:$alias"}
117 :     elsif ($alias =~ /gi\|/)
118 :     {
119 :     $alias =~ s/^gi\|//;
120 :     push @{$note->{"Dbxref"}}, "NCBI_gi:$alias";
121 :     }
122 :     elsif ($alias =~ /^kegg/i)
123 :     {
124 :     $alias =~ s/kegg\|/KEGG:/i;
125 :     $alias =~ s/^(.*):/$1+/;
126 :     push @{$note->{"Dbxref"}}, $alias
127 :     }
128 :     elsif ($alias =~ /^uni/)
129 :     {
130 :     $alias =~ s/uni\|/UniProt:/;
131 :     push @{$note->{"Dbxref"}}, $alias;
132 :     }
133 :     elsif ($alias =~ /^sp\|/)
134 :     {
135 :     $alias =~ s/sp\|/Swiss-Prot:/;
136 :     push @{$note->{"Dbxref"}}, $alias
137 :     }
138 :     elsif ($alias =~ /^[a-zA-Z][a-zA-Z0-9]{2,}\_[a-zA-Z]*\d+$/)
139 :     {
140 :     # that should be the regexp for a valid locus tag. Starts with a letter,
141 :     # has at least three alphanumerics before the _
142 :     # then has a series of optional letters and ends with numbers
143 :     push @{$note->{"locus_tag"}}, $alias;
144 :     }
145 :     elsif ($alias =~ /^.{4,6}$/)
146 :     {
147 :     push @{$note->{"gene_symbol"}}, $alias;
148 :     }
149 :     elsif ($alias =~ /^[a-zA-Z]+\d+$/)
150 :     {
151 :     push @{$note->{"locus"}}, $alias;
152 :     }
153 :     else {
154 :     # don't know what it is so keep it as an alias
155 :     push @{$note->{"Alias"}}, $alias;
156 :     }
157 :     }
158 :    
159 :     # get the links
160 :     foreach my $ln ($fig->fid_links($peg))
161 :     {
162 :     my ($db, $id);
163 :     if ($ln =~ /field0=CATID/ && $ln =~ /query0=(\d+)/ && ($id=$1) && $ln =~ /pcenter=harvard/) {
164 :     $db = "HarvardClone";
165 :     } elsif ($ln =~ /(PIRSF)(\d+)/) {
166 :     ($db, $id) = ($1, $2);
167 :     } elsif ($ln =~ />(\S+)\s+(\S+.*?)</) {
168 :     ($db, $id) = ($1, $2);
169 :     }
170 :    
171 :     $db =~ s/\://;
172 :     if (!$db && !$id) {
173 :     print STDERR "Ignored link: $ln\n";
174 :     next;
175 :     }
176 :     push @{$note->{"Dbxref"}}, "$db:$id";
177 :     }
178 :    
179 :     # add FIG id as a note
180 :     push @{$note->{"Dbxref"}}, "FIG_ID:$peg";
181 :    
182 :     # get the features
183 :     my @location = $fig->feature_location($peg);
184 :     foreach my $loc (@location) {
185 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
186 : paczian 1.2 my ($contig, $start, $stop) = ($1, $2, $3);
187 :     my $original_contig = $contig;
188 :     my $strand = '1';
189 :     my $frame = $start % 3;
190 :     if ($start > $stop) {
191 :     $frame = $stop % 3;
192 :     ($start, $stop, $strand) = ($stop, $start, '-1');
193 :     }
194 :     elsif ($start == $stop) {
195 :     $strand = ".";
196 :     $frame = ".";
197 :     }
198 :     my $source = "FIG";
199 : paczian 1.1 my $type = $fig->ftype($peg);
200 :     my $feature;
201 :     if ($type eq "peg") {
202 :     $feature = Bio::SeqFeature::Generic->new(
203 :     -start => $start,
204 :     -end => $stop,
205 :     -strand => $strand,
206 :     -primary => 'CDS',
207 : paczian 1.2 -tag => {
208 :     product => $func,
209 :     translation => $fig->get_translation($peg),
210 :     },
211 : paczian 1.1 );
212 :    
213 :     foreach my $tagtype (keys %$note) {
214 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
215 :     }
216 : paczian 1.2
217 :     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 :    
226 :     push(@$bio2, $feature2);
227 :    
228 : paczian 1.1 } elsif ($type eq "rna") {
229 :     $feature = Bio::SeqFeature::Generic->new(
230 :     -start => $start,
231 :     -end => $stop,
232 :     -strand => $strand,
233 :     -primary => 'RNA',
234 :     );
235 :     foreach my $tagtype (keys %$note) {
236 :     $feature->add_tag_value($tagtype, @{$note->{$tagtype}});
237 :     }
238 :    
239 :     } else {
240 :     print STDERR "unhandled feature type: $type\n";
241 :     }
242 :    
243 :     $bio->{$contig}->add_SeqFeature($feature);
244 :     }
245 :     }
246 :    
247 :     # generate filename
248 :     my $filename = $directory . $genome . "." . $format_ending;
249 : paczian 1.2
250 : paczian 1.1 # check for FeatureIO or SeqIO
251 :     if ($format eq "GTF") {
252 : paczian 1.2 my $fio = Bio::FeatureIO->new(-file => ">$filename", -format => "GTF");
253 :     foreach my $feature (@$bio2) {
254 :     $fio->write_feature($feature);
255 :     }
256 : paczian 1.1 } else {
257 : paczian 1.2 my $sio = Bio::SeqIO->new(-file => ">$filename", -format => $format);
258 :     foreach my $seq (keys %$bio) {
259 :     $sio->write_seq($bio->{$seq});
260 :     }
261 : paczian 1.1 }
262 : paczian 1.2
263 : paczian 1.1 return ($filename, "Output file successfully written.");
264 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3