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

Annotation of /FortyEight/SeedExport.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3