[Bio] / FigKernelPackages / FigGFF.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/FigGFF.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 #
2 :     # FIG GFF utilities.
3 :     #
4 :    
5 :     #
6 :     # A GFFWriter handles the generation of GFF files from SEED data structures.
7 :     #
8 :    
9 :    
10 :     package GFFWriter;
11 :    
12 :     use strict;
13 :     use FIG;
14 :    
15 :     use Carp;
16 :     use URI::Escape;
17 :     use Data::Dumper;
18 :    
19 :     sub new
20 :     {
21 :     my($class, $fig, %options) = @_;
22 :    
23 :     my $default_options = {
24 :     escapespace => 0,
25 :     outputfasta => 1,
26 :     linelength => 60,
27 :     };
28 :    
29 :     map { $default_options->{$_} = $options{$_} } keys(%options);
30 :    
31 :     my $self = {
32 :     options => $default_options,
33 :     contig_length_cache => {},
34 :     fig => $fig,
35 :     };
36 :    
37 :     return bless $self, $class;
38 :     }
39 :    
40 :    
41 :     =head1 gff3_for_feature
42 :    
43 :     Returns the GFF3 information for a given feature.
44 :    
45 :     The return is a pair ($contig_data, $fasta_sequences) that can be passed
46 :     into write_gff3().
47 :    
48 :     $contig_data is a hashref mapping a contig name to a list of GFF3 file lines
49 :     for the sequences in that contig.
50 :    
51 :     =cut
52 :    
53 :     sub gff3_for_feature
54 :     {
55 :     my($self, $fid, $user) = @_;
56 :    
57 :     #
58 :     # Options we need to figure out somehow.
59 :     #
60 :     my $options = $self->{options};
61 :    
62 :     my $escapespace = $options->{escapespace};
63 :     my $outputfasta = $options->{outputfasta};
64 :     my %outputtype = (cds => 1,
65 :     protein => 1,
66 :     transcript => 1,
67 :     gene => 1);
68 :    
69 :     my $fastasequences = '';
70 :     my $contig_data;
71 :     my $linelength = $options->{linelength};
72 :    
73 :     my $fig = $self->{fig};
74 :    
75 :     #
76 :     # Do this first to make sure that we really have a feature.
77 :     #
78 :     my @location = $fig->feature_location($fid);
79 :     print Dumper(\@location);
80 :     if (@location == 0 or !defined($location[0]))
81 :     {
82 :     warn "No location found for feature $fid\n";
83 :     return ({}, "");
84 :     }
85 :    
86 :     ###########
87 :     #
88 :     # Begin figuring out the column 9 information about notes and aliases and GO terms
89 :     # All the information is temporarily stored in @alias or @note, and at the end is joined
90 :     # into $allnote
91 :     #
92 :     ###########
93 :    
94 :     #
95 :     # the notes for the last column
96 :     #
97 :     my $note;
98 :     #
99 :     # all the aliases we are going to use
100 :     #
101 :     my @alias;
102 :    
103 :     my $func = $fig->function_of($fid, $user);
104 :     if ($func)
105 :     {
106 :     push @$note, ("Note=\"" . uri_escape($func) . '"');
107 :     }
108 :    
109 :     # now find aliases
110 :     foreach my $alias ($fig->feature_aliases($fid))
111 :     {
112 :     if ($alias =~ /^NP/)
113 :     {
114 :     push @$note, "Dbxref=\"NCBI_genpept:$alias\"";
115 :     }
116 :     elsif ($alias =~ /gi\|/)
117 :     {
118 :     $alias =~ s/^gi\|//;
119 :     push @$note, "Dbxref=\"NCBI_gi:$alias\"";
120 :     }
121 :     elsif ($alias =~ /^kegg/i)
122 :     {
123 :     $alias =~ s/kegg\|/KEGG:/i;
124 :     $alias =~ s/^(.*):/$1+/;
125 :     push @$note, "Dbxref=\"$alias\"";
126 :     }
127 :     elsif ($alias =~ /^uni/)
128 :     {
129 :     $alias =~ s/uni\|/UniProt:/;
130 :     # $note = check_go($note, $alias) if ($go);
131 :     push @$note, "Dbxref=\"$alias\"";
132 :     }
133 :     elsif ($alias =~ /^sp\|/)
134 :     {
135 :     $alias =~ s/sp\|/Swiss-Prot:/;
136 :     push @$note, "Dbxref=\"$alias\"";
137 :     }
138 :     else
139 :     {
140 :     # don't know what it is so keep it as an alias
141 :     $alias = uri_escape($alias); # just in case!
142 :     push @alias, $alias;
143 :     }
144 :     }
145 :    
146 :     # now just join all the aliases and put them into @$note so we can add it to the array
147 :     if (@alias)
148 :     {
149 :     push @$note, "Alias=\"". join (",", @alias) . '"';
150 :     }
151 :    
152 :     # the LAST thing I am going to add as a note is the FIG id so that I can grep it out easily
153 :     push @$note, "Dbxref=\"SEED:$fid\"";
154 :    
155 :     # finally join all the notes into a long string that can be added as column 9
156 :     my $allnotes;
157 :     $allnotes = join ";", @$note;
158 :    
159 :     # do we want to convert '%20' to ' '
160 :     unless ($escapespace)
161 :     {
162 :     $allnotes =~ s/\%20/ /g;
163 :     }
164 :    
165 :     ###########
166 :     #
167 :     # End figuring out the column 9 information about notes and aliases and GO terms
168 :     #
169 :     ###########
170 :    
171 :     #
172 :     # Cache contig lengths.
173 :     #
174 :     my $len = $self->{contig_length_cache};
175 :    
176 :     my $genome = $fig->genome_of($fid);
177 :    
178 :     foreach my $loc (@location)
179 :     {
180 :     $loc =~ /^(.*)\_(\d+)\_(\d+)$/;
181 :     my ($contig, $start, $stop) = ($1, $2, $3);
182 :     my $original_contig=$contig;
183 :    
184 :     #
185 :     # the contig name must be escaped
186 :     #
187 :     $contig = uri_escape($contig);
188 :    
189 :     #my $contig_key = "$genome:$contig";
190 :     my $contig_key = $contig;
191 :    
192 :     unless (defined $len->{$contig})
193 :     {
194 :     $len->{$contig}=$fig->contig_ln($genome, $original_contig);
195 :     }
196 :     my $strand='+';
197 :    
198 :     #
199 :     # These were bounds-checking for dumping all of a genome.
200 :     #
201 :     #next if (defined $beg && ($start < $beg || $stop < $beg));
202 :     #next if (defined $end && ($start > $end || $stop > $end));
203 :    
204 :     if ($start > $stop)
205 :     {
206 :     ($start, $stop, $strand)=($stop, $start, '-');
207 :     }
208 :     elsif ($start == $stop)
209 :     {
210 :     $strand=".";
211 :     }
212 :    
213 :     my $type=$fig->ftype($fid);
214 :    
215 :     if ($type eq "peg")
216 :     {
217 :     # it is a protein coding gene
218 :     # create an artificial id that is just the fid.(\d+) information
219 :     # we will use this to create ids in the form cds.xxx; trn.xxxx; pro.xxx; gen.xxx;
220 :     $fid =~ /\.peg\.(\d+)/;
221 :     my $geneid=$1;
222 :    
223 :     ############## KLUDGE
224 :     #
225 :     # At the moment the outputs for transcript, gene, CDS, and pro are all the same.
226 :     # This is clearly a kludge and wrong, but it will work at the moment
227 :     #
228 :    
229 :     # defined some truncations
230 :     my %trunc=(
231 :     "transcript" => "trn",
232 :     "gene" => "gen",
233 :     "protein" => "pro",
234 :     "cds" => "cds",
235 :     );
236 :    
237 :     # SO terms:
238 :     # transcript: SO:0000673
239 :     # gene: SO:0000704
240 :     # cds: SO:0000316
241 :     # protein: NOT VALID: should be protein_coding_primary_transcript SO:0000120
242 :     foreach my $type (keys %outputtype)
243 :     {
244 :     my ($id, $type)=($trunc{$type} . "." . $geneid, $type);
245 :     # we want to store some sequences to be output
246 :     if ($outputfasta && $type eq "protein")
247 :     {
248 :     my $addseq = $fig->get_translation($fid);
249 :    
250 :     #
251 :     # the chomp is so that we know for sure to add the line back
252 :     #
253 :     $addseq =~ s/(.{$linelength})/$1\n/g;
254 :     chomp($addseq);
255 :     $fastasequences .= ">$id\n$addseq\n";
256 :     }
257 :     if ($outputfasta && $type eq "cds")
258 :     {
259 :     my $addseq = uc($fig->dna_seq($genome, $fig->feature_location($fid)));
260 :     $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
261 :     $fastasequences .= ">$id\n$addseq\n";
262 :     }
263 :     # correct the incorrect sofa term
264 :     if ($type eq "protein")
265 :     {
266 :     $type = "SO:0000120";
267 :     }
268 :    
269 :     push (@{$contig_data->{$contig_key}},
270 :     (join "\t",
271 :     ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "Id=$id;$allnotes")));
272 :     } # end the foreach my $type
273 :     } # end the if type==peg
274 :     elsif ($type eq "rna")
275 :     {
276 :     $fid =~ /\.rna\.(\d+)/;
277 :     my $geneid=$1;
278 :     #
279 :     # tRNA is a valid SOFA term == SO:0000253
280 :     #
281 :     my ($id, $type)=("rna.$geneid", "tRNA");
282 :     if ($outputfasta)
283 :     {
284 :     my $addseq = $fig->dna_seq($genome, $fig->feature_location($fid));
285 :     $addseq =~ s/(.{$linelength})/$1\n/g; chomp($addseq);
286 :     $fastasequences .= ">$id\n$addseq\n";
287 :     }
288 :     push (@{$contig_data->{$contig_key}}, (join "\t", ($contig, "The SEED", $type, $start, $stop, ".", $strand, ".", "Id=$id;$allnotes")));
289 :     } # end the if type == rna
290 :     else
291 :     {
292 :     die "Don't know what type: |$type| is";
293 :     }
294 :     }
295 :     return ($contig_data, $fastasequences);
296 :     }
297 :    
298 :     =head1 write_gff3
299 :    
300 :     Write a set of gff3 per-contig data and fasta sequence data to a file or filehandle.
301 :    
302 :     $genome is the genome these contigs are a part of.
303 :     $contig_list is a list of contig-data hashes as returned by gff_for_feature.
304 :     $fast_list is a list of fasta data strings.
305 :    
306 :     =cut
307 :    
308 :     sub write_gff3
309 :     {
310 :     my($self, $output, $genome, $contig_list, $fasta_list) = @_;
311 :    
312 :     my $fig = $self->{fig};
313 :    
314 :     my $len = $self->{contig_length_cache};
315 :     my $fh;
316 :    
317 :     my $beg = $self->{options}->{beg};
318 :     my $end = $self->{options}->{end};
319 :    
320 :     my $close_output;
321 :    
322 :     if (ref($output))
323 :     {
324 :     $fh = $output;
325 :     }
326 :     else
327 :     {
328 :     open($fh, ">$output") or confess "Cannot open output '$output': $!";
329 :     $close_output = 1;
330 :     }
331 :    
332 :     #
333 :     # Build a data structure from the list of contigs
334 :     # that has a list of lists of data per contig name.
335 :     # (Do this so we don't copy all of the contig data itself, as it
336 :     # could be quite large).
337 :     #
338 :     my %contigs;
339 :    
340 :     #
341 :     # iterate over the given list of contig hashes.
342 :     #
343 :     for my $chash (@$contig_list)
344 :     {
345 :     #
346 :     # Then for each contig in the individual contig hashes,
347 :     # add the data list to %contigs.
348 :     #
349 :     for my $contig (keys %$chash)
350 :     {
351 :     push(@{$contigs{$contig}}, $chash->{$contig});
352 :     }
353 :     }
354 :    
355 :     foreach my $contig (sort keys %contigs)
356 :     {
357 :     print $fh "##sequence-region\t$contig\t";
358 :     if (defined $beg) {
359 :     print $fh "$beg\t";
360 :     } else {
361 :     print $fh "1\t";
362 :     }
363 :     if (defined $end) {
364 :     print $fh "$end\n";
365 :     } else {
366 :     print $fh "$len->{$contig}\n";
367 :     }
368 :     for my $list (@{$contigs{$contig}})
369 :     {
370 :     print $fh join "\n", @$list;
371 :     }
372 :     }
373 :    
374 :     print $fh "##FASTA\n";
375 :     # print out the cds and pro if we need them
376 :    
377 :     if ($self->{options}->{outputfasta})
378 :     {
379 :     for my $fastasequences (@$fasta_list)
380 :     {
381 :     print $fh $fastasequences;
382 :     }
383 :     }
384 :    
385 :     my $ll = $self->{options}->{linelength};
386 :     foreach my $contig (sort keys %contigs)
387 :     {
388 :     my $len=$fig->contig_ln($genome, $contig);
389 :     my $dna_seq=$fig->dna_seq($genome, $contig . "_1_". $len);
390 :     if (defined $beg)
391 :     {
392 :     unless (defined $end) {
393 :     $end=$len;
394 :     }
395 :     $dna_seq = substr($dna_seq, $beg, $end);
396 :     }
397 :     elsif (defined $end)
398 :     {
399 :     $beg=1;
400 :     $dna_seq = substr($dna_seq, $beg, $end);
401 :     }
402 :    
403 :     my $contig=uri_escape($contig);
404 :    
405 :     $dna_seq =~ s/(.{$ll})/$1\n/g;
406 :     chomp($dna_seq); # just remove the last \n if there is one
407 :     print $fh ">$contig\n$dna_seq\n";
408 :     }
409 :    
410 :     close($fh) if $close_output;
411 :     }
412 :    
413 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3