[Bio] / FigKernelScripts / parse_genbank.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.27 # -*- perl -*-
2 : olson 1.30
3 :     #
4 :     # This is a SAS component.
5 :     #
6 :    
7 : overbeek 1.14 ########################################################################
8 : olson 1.12 #
9 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
10 :     # for Interpretations of Genomes. All Rights Reserved.
11 :     #
12 :     # This file is part of the SEED Toolkit.
13 :     #
14 :     # The SEED Toolkit is free software. You can redistribute
15 :     # it and/or modify it under the terms of the SEED Toolkit
16 :     # Public License.
17 :     #
18 :     # You should have received a copy of the SEED Toolkit Public License
19 :     # along with this program; if not write to the University of Chicago
20 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
21 :     # Genomes at veronika@thefig.info or download a copy from
22 :     # http://www.theseed.org/LICENSE.TXT.
23 :     #
24 : gdpusch 1.18 ########################################################################
25 : olson 1.12
26 : olson 1.30 use SeedUtils;
27 : gdpusch 1.16 use Carp;
28 : gdpusch 1.18 use Data::Dumper;
29 :    
30 :     my %trans_table;
31 :    
32 : gdpusch 1.33 my $usage = qq(usage: parse_genbank [-user_trans_only] [-i=genbank.entry] [-bioname=GENOME_BIONAME] [-taxonomy=TAXONOMY] [-source=DATA_SOURCE] Taxon_ID Org_Dir < genbank.entry);
33 : efrank 1.1
34 : gdpusch 1.18 my $trouble = 0;
35 : gdpusch 1.33 my $no_trans = 0;
36 : gdpusch 1.18 my $force_bioname;
37 :     my $force_taxonomy;
38 : gdpusch 1.22 my $data_source = qq(Parsed GenBank File);
39 : olson 1.31 my $genbank_entry;
40 : gdpusch 1.9 while ($ARGV[0] =~ m/^-/)
41 :     {
42 : gdpusch 1.33 if ($ARGV[0] =~ m/^-{1,2}user_trans_only/o) {
43 :     $no_trans = 1;
44 :     }
45 :     elsif ($ARGV[0] =~ m/^-{1,2}bioname=(\S+)/) {
46 : gdpusch 1.23 $force_bioname = $1;
47 :     $force_bioname =~ s/^[\'\"]//o;
48 :     $force_bioname =~ s/[\'\"]$//o;
49 : gdpusch 1.18 $force_bioname =~ s/\s+/ /g;
50 :    
51 :     print STDERR qq(\nGenome bioname will be taken as: \"$force_bioname\"\n);
52 :     }
53 : gdpusch 1.23 elsif ($ARGV[0] =~ m/^-{1,2}taxonomy=(\S+)/) {
54 :     $force_taxonomy = $1;
55 :     $force_taxonomy =~ s/^[\'\"]//o;
56 :     $force_taxonomy =~ s/[\'\"]$//o;
57 : gdpusch 1.18 $force_taxonomy =~ s/\s+/ /g;
58 :     if ($force_taxonomy !~ m/\.$/) { $force_taxonomy .= qq(.); }
59 :    
60 :     print STDERR qq(\nTaxonomy will be taken as: \"$force_taxonomy\"\n);
61 :     }
62 : gdpusch 1.23 elsif ($ARGV[0] =~ m/^-{1,2}source=(\S+)/) {
63 :     $data_source = $1;
64 :     $data_source =~ s/^[\'\"]//o;
65 :     $data_source =~ s/[\'\"]$//o;
66 : gdpusch 1.18
67 :     print STDERR qq(Data source string is: \"$data_source\"\n);
68 :     }
69 : olson 1.32 elsif ($ARGV[0] =~ m/^-{1,2}i=(.*)/) {
70 : olson 1.31 $genbank_entry = $1;
71 :     }
72 : gdpusch 1.18 else {
73 :     $trouble = 1;
74 :     die qq(Could not handle $ARGV[0]);
75 : gdpusch 1.9 }
76 : gdpusch 1.23 shift @ARGV;
77 : gdpusch 1.9 }
78 :    
79 : olson 1.31 my $input_fh;
80 :     if (defined($genbank_entry))
81 :     {
82 :     if (!open($input_fh, "<", $genbank_entry))
83 :     {
84 :     die "Could not open genbank entry $genbank_entry: $!";
85 :     }
86 :     }
87 :     else
88 :     {
89 :     $input_fh = \*STDIN;
90 :     }
91 :    
92 : gdpusch 1.18 my $taxon_ID;
93 :     my $org_dir;
94 : efrank 1.2 (
95 : gdpusch 1.18 ($taxon_ID = shift(@ARGV)) &&
96 :     ($org_dir = shift(@ARGV))
97 : efrank 1.1 )
98 :     || die $usage;
99 :    
100 : gdpusch 1.18 $prefixP = qq(fig|$taxon_ID.peg.);
101 :     $prefixR = qq(fig|$taxon_ID.rna.);
102 :    
103 :     #...Create paths down to PEG and RNA dirs
104 : olson 1.30 &SeedUtils::verify_dir("$org_dir/Features/peg");
105 :     &SeedUtils::verify_dir("$org_dir/Features/rna");
106 : gdpusch 1.18
107 :     my $fh_contigs;
108 :     my $contigs_file = qq($org_dir/contigs);
109 :     open($fh_contigs, qq(>$contigs_file))
110 :     || die qq(could not open $contigs_file);
111 :    
112 :     my $fh_assigned_funcs;
113 :     my $assigned_funcs_file = qq($org_dir/assigned_functions);
114 :     open($fh_assigned_funcs, qq(>$assigned_funcs_file))
115 :     || die qq(could not write-open $assigned_funcs_file);
116 :    
117 :     my $fh_annotations;
118 :     my $annotations_file = qq($org_dir/annotations);
119 :     open($fh_annotations, qq(>$annotations_file))
120 :     || die qq(could not write-open $annotations_file);
121 :    
122 :     my $fh_ec_nums;
123 :     my $ec_nums_file = qq($org_dir/EC_numbers);
124 :     open($fh_ec_nums, qq(>$ec_nums_file))
125 :     || die" Could not write-open $ec_nums_file";
126 :    
127 :     my $fh_peg_tbl;
128 :     my $peg_tbl_file = qq($org_dir/Features/peg/tbl);
129 :     open($fh_peg_tbl, qq(>$peg_tbl_file))
130 :     || die qq(could not open $peg_tbl_file);
131 :    
132 :     my $fh_peg_fasta;
133 :     my $peg_fasta_file = qq($org_dir/Features/peg/fasta);
134 :     open($fh_peg_fasta, qq(>$peg_fasta_file))
135 :     || die qq(could not open $peg_fasta_file);
136 :    
137 :     my $fh_rna_tbl;
138 :     my $rna_tbl_file = qq($org_dir/Features/rna/tbl);
139 :     open($fh_rna_tbl, qq(>$rna_tbl_file))
140 :     || die qq(could not write-open $rna_tbl_file);
141 :    
142 :     my $fh_rna_fasta;
143 :     my $rna_fasta_file = qq($org_dir/Features/rna/fasta);
144 :     open($fh_rna_fasta, qq(>$rna_fasta_file))
145 :     || die qq(could not write-open $rna_fasta_file);
146 : efrank 1.2
147 : gdpusch 1.18
148 :     open( TAXONOMY, qq(>$org_dir/TAXONOMY)) || die qq(could not open $org_dir/TAXONOMY);
149 :     open( GENOME, qq(>$org_dir/GENOME)) || die qq(could not open $org_dir/GENOME);
150 :    
151 :     open( PROJECT, qq(>$org_dir/PROJECT)) || die qq(could not open $org_dir/PROJECT);
152 :     print PROJECT qq($data_source\n);
153 : overbeek 1.4 close(PROJECT);
154 :    
155 : gdpusch 1.18 my $idNp = 1;
156 :     my $idNr = 1;
157 : overbeek 1.4
158 : gdpusch 1.18 $/ = qq(\n//\n);
159 :     my $record;
160 :     my $contigs = {};
161 : olson 1.31
162 :     while (defined($record = <$input_fh>)) {
163 : gdpusch 1.16 $record =~ s/^\s+//os;
164 :     last unless $record;
165 :    
166 :     if ($record !~ m/LOCUS\s+(\S+)/os) {
167 : gdpusch 1.18 print STDERR qq(No LOCUS line for record begining with:\n)
168 :     , substr($record, 0, 160), qq(\n\n);
169 : gdpusch 1.9 }
170 : gdpusch 1.16 else {
171 : gdpusch 1.18 my $id = $1;
172 : redwards 1.24 if ($record =~ m/\nORIGIN[^\n]*\n(.*?)(\/\/|LOCUS)/os) {
173 : gdpusch 1.18 my $seq = $1;
174 : gdpusch 1.16 $seq =~ s/\s//ogs;
175 :     $seq =~ s/\d//og;
176 : efrank 1.1 undef $contigs;
177 :     $contigs->{$id} = $seq;
178 : olson 1.30 &SeedUtils::display_id_and_seq($id, \$seq, $fh_contigs);
179 : efrank 1.1 }
180 : gdpusch 1.16 else {
181 : gdpusch 1.18 warn qq(could not find contig sequence for $id\n);
182 : overbeek 1.4 next;
183 :     }
184 : gdpusch 1.16
185 : gdpusch 1.18 my ($taxon_ID, $taxonomy, $written_genome);
186 : gdpusch 1.16 if ($record =~ /\n {0,4}ORGANISM\s+(\S[^\n]+(\n\s{10,14}\S[^\n]+)*)/os) {
187 : overbeek 1.14 my $block = $1;
188 :     my @lines = split(/\n/,$block);
189 : gdpusch 1.34 # print STDERR Dumper($block, \@lines);
190 : gdpusch 1.16
191 : overbeek 1.14 my @genome = ();
192 :     my @full_tax = ();
193 : gdpusch 1.34 for ($i=0; ($i < $#lines) && ($lines[$i] !~ /;/); $i++) {
194 : overbeek 1.14 push(@genome,$lines[$i]);
195 :     }
196 : gdpusch 1.34 # print STDERR "i=$i\n";
197 : gdpusch 1.16
198 :     while ($i < @lines) {
199 : overbeek 1.14 push(@full_tax,$lines[$i]);
200 : gdpusch 1.27 ++$i;
201 : overbeek 1.14 }
202 : gdpusch 1.34 # print STDERR Dumper(\@genome, @full_tax);
203 : gdpusch 1.16
204 : gdpusch 1.34 $genome = join(qq( ), map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @genome);
205 :     $taxonomy = join(qq( ), map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @full_tax);
206 : gdpusch 1.16
207 : gdpusch 1.18 $taxonomy =~ s/\n\s+//og;
208 :     $taxonomy =~ s/ {2,}/ /og;
209 :     $taxonomy =~ s/\.$//o;
210 : gdpusch 1.34 $taxonomy = $taxonomy . qq(; $genome\.);
211 :     $taxonomy =~ s/\.\.$/./o;
212 :     # print STDERR "$genome\n$taxonomy\n";
213 :    
214 : gdpusch 1.16 if (! $written_genome) {
215 : gdpusch 1.20 if ($force_bioname) { $genome = $force_bioname; }
216 :     print GENOME qq($genome\n);
217 : overbeek 1.4 close(GENOME);
218 : gdpusch 1.9
219 : gdpusch 1.18 if ($force_taxonomy) { $taxonomy = $force_taxonomy; }
220 :     print TAXONOMY qq($taxonomy\n);
221 :     close(TAXONOMY);
222 :    
223 : gdpusch 1.21 $written_genome = qq($taxonomy\t$genome);
224 : overbeek 1.4 }
225 : gdpusch 1.18 elsif (($written_genome ne qq($taxonomy\t$bioname)) && (!$force_bioname || !$force_taxonomy)) {
226 :     print STDERR qq(serious mismatch in GENOME/TAX for $id\n$written_genome\n$taxonomy\t$bioname\n\n);
227 : overbeek 1.4 }
228 : efrank 1.1 }
229 : gdpusch 1.16
230 :     while ($record =~ m/\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/ogs) {
231 : gdpusch 1.18 my $cds = $1;
232 : gdpusch 1.16 if (($cds !~ m/\/pseudo/o) &&
233 :     (($cds !~ m/\/exception/o) || ($cds =~ m/\/translation/o))
234 :     ) {
235 : gdpusch 1.18 &process_cds($id, \$cds, $prefixP, \$idNp, $contigs, $fh_peg_tbl, $fh_peg_fasta, $fh_assigned_funcs, $fh_annotations, $fh_ec_nums);
236 : overbeek 1.7 }
237 : efrank 1.2 }
238 :    
239 : gdpusch 1.27 while ($record =~ m/\n\s{3,6}(([tr]|misc\_)RNA)\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/ogs)
240 : efrank 1.2 {
241 : gdpusch 1.27 $type = $2;
242 :     $rna = $3;
243 : gdpusch 1.28 &process_rna($id, $type, \$rna, $prefixR, \$idNr, $contigs, $fh_rna_tbl, $fh_rna_fasta, $fh_assigned_funcs, $fh_annotations);
244 : efrank 1.1 }
245 :     }
246 :     }
247 : gdpusch 1.18 close($fh_contigs);
248 :     close($fh_peg_tbl);
249 :     close($fh_peg_fasta);
250 :     close($fh_rna_tbl);
251 :     close($fh_rna_fasta);
252 :     close($fh_assigned_funcs);
253 :     close($fh_annotations);
254 :     close($fh_ec_nums);
255 :    
256 :     if (!-s $assigned_funcs_file) {
257 : gdpusch 1.20 warn qq(WARNING: no assigned_functions in $org_dir --- deleting\n);
258 : gdpusch 1.18 unlink("$assigned_funcs_file") unless $ENV{DEBUG};
259 :     }
260 :    
261 : gdpusch 1.22 if (!-s $annotations_file) {
262 :     warn qq(WARNING: no annotations in $org_dir --- deleting\n);
263 :     unlink("$annotations_file") unless $ENV{DEBUG};
264 :     }
265 :    
266 : gdpusch 1.18 if (!-s $ec_nums_file) {
267 : gdpusch 1.20 warn qq(WARNING: no EC_numbers in $org_dir --- deleting\n);
268 : gdpusch 1.18 unlink("$ec_nums_file") unless $ENV{DEBUG};
269 :     }
270 : efrank 1.2
271 : gdpusch 1.20 if ((!-s $rna_tbl_file) || (!-s $rna_fasta_file)) {
272 :     warn qq(WARNING: no RNAs in $org_dir --- deleting\n);
273 : olson 1.29 system('rm', '-rf', "$org_dir/Features/rna") unless $ENV{DEBUG};;
274 : gdpusch 1.20 }
275 :    
276 : gdpusch 1.18 if ((!-s $peg_tbl_file) || (!-s $peg_fasta_file)) {
277 : gdpusch 1.20 warn qq(WARNING: no PEGs in $org_dir --- deleting\n);
278 : olson 1.29 system('rm', '-rf', "$org_dir/Features/peg") unless $ENV{DEBUG};
279 : gdpusch 1.18 }
280 :    
281 : gdpusch 1.20 if ((!-d qq($org_dir/Features/peg)) && (!-d qq($org_dir/Features/peg))) {
282 :     warn qq(WARNING: no Features in $org_dir --- deleting\n);
283 : olson 1.29 system("rmdir", "$org_dir/Features") unless $ENV{DEBUG};;
284 : gdpusch 1.18 }
285 : gdpusch 1.20
286 : gdpusch 1.18 if (!-s $contigs_file) {
287 :     $trouble = 1;
288 :     unlink($contigs_file);
289 : gdpusch 1.20 print STDERR qq(ERROR: No contigs in $org_dir\n);
290 : gdpusch 1.18 }
291 :    
292 :     if ($trouble) {
293 :     warn qq(Genome directory $org_dir is corrupt --- deleting\n\n);
294 : olson 1.29 system('rm', '-rf', $org_dir) unless $ENV{DEBUG};
295 : gdpusch 1.18 }
296 :    
297 :     exit($trouble);
298 : efrank 1.1
299 :     sub process_cds {
300 : gdpusch 1.20 my ($contigID, $cdsP, $prefix, $idNp, $contigs, $fh_tbl, $fh_fasta, $fh_assign, $fh_annot, $fh_ec_nums) = @_;
301 : gdpusch 1.18 my ($id, $prot);
302 : gdpusch 1.9
303 : overbeek 1.13 ++$recnum;
304 : gdpusch 1.18 my ($loc, $precise) = &get_loc($contigID,$cdsP);
305 :     my @aliases = &get_aliases($cdsP);
306 :     my @ec_nums = &get_ec_numbers($cdsP);
307 :     my ($func, $notes) = &get_func($cdsP);
308 :     my $trans = &get_trans($cdsP);
309 :    
310 : gdpusch 1.26 if (not $trans) {
311 : gdpusch 1.33 if ($no_trans) {
312 :     #...Do nothing
313 :     }
314 :     else {
315 :     warn qq(WARNING: Translation missing for CDS $recnum; attempting to generate translation\n) if $ENV{VERBOSE};
316 :    
317 :     if ($loc && $precise) {
318 :     my $dna = &SeedUtils::extract_seq($contigs,$loc);
319 :     if ($dna) {
320 :     my $genetic_code = 11;
321 :     if ($$cdsP =~ m/\/transl_table=(\d+)/o) {
322 :     $genetic_code = $1;
323 :     }
324 :    
325 :     if (not defined($trans_table{$genetic_code})) {
326 :     $trans_table{$genetic_code} = &SeedUtils::genetic_code($genetic_code);
327 :     }
328 :    
329 :     my $prot = &SeedUtils::translate($dna, $trans_table{$genetic_code}, 1);
330 :     $prot =~ s/\*$//o;
331 :    
332 :     if ($prot !~ m/\*/o) {
333 :     $trans = $prot;
334 :     }
335 :     else {
336 :     warn qq(Translation contains STOPs, changing to \'x\'s, for CDS $recnum:\n$$cdsP\n);
337 :     }
338 : gdpusch 1.9 }
339 : efrank 1.1 }
340 : gdpusch 1.33 else {
341 :     warn (qq(Could not get DNA sequence for CDS at $loc for entry:\n)
342 :     , q( CDS )
343 :     , $$cdsP
344 :     , qq(\nof record begining with:\n)
345 :     , substr($record, 0, 160)
346 :     , qq(\n\n)
347 :     );
348 :     }
349 : gdpusch 1.26 }
350 :     }
351 :    
352 :     if ($trans) {
353 : gdpusch 1.27 $id = $prefix . qq($$idNp);
354 :     ++$$idNp;
355 :    
356 : olson 1.19 print $fh_tbl qq($id\t$loc\t) . join("\t",@aliases) . "\n";
357 : gdpusch 1.16
358 :     if ($func) {
359 : gdpusch 1.20 print $fh_assign qq($id\t$func\n);
360 : efrank 1.1 }
361 : gdpusch 1.16
362 : gdpusch 1.18 if (@ec_nums) {
363 :     print $fh_ec_nums (join(qq(\t), ($id, @ec_nums)), qq(\n));
364 :     }
365 :    
366 : gdpusch 1.26 foreach my $note (@$notes) {
367 :     &make_annotation($id, $note, $fh_annot);
368 :     }
369 : gdpusch 1.18
370 : olson 1.30 &SeedUtils::display_id_and_seq($id, \$trans, $fh_fasta);
371 : efrank 1.1 }
372 : gdpusch 1.16 else {
373 : gdpusch 1.27 warn (qq(No translation for CDS --- skipping entry:\n), &Dumper($cdsP), qq(\n));
374 : efrank 1.1 }
375 :     }
376 :    
377 :     sub get_loc {
378 :     my($contigID,$cdsP) = @_;
379 :     my($beg,$end,$loc,$locS,$iter,$n,$args,@pieces);
380 :     my($precise);
381 : gdpusch 1.16
382 :     if ($$cdsP =~ m/^(([^\n]+)((\n\s{21,23}[^\/ ][^\n]+)+)?)/os) {
383 : efrank 1.1 $locS = $1;
384 :     $locS =~ s/\s//g;
385 : gdpusch 1.16 $precise = ($locS !~ m/[<>]/o);
386 : gdpusch 1.10
387 : efrank 1.1 @pieces = ();
388 :     $n = 0;
389 :     $iter = 500;
390 : gdpusch 1.16 while (($locS !~ m/^\[\d+\]$/o) && $iter)
391 : efrank 1.1 {
392 : gdpusch 1.27 --$iter;
393 : gdpusch 1.16 if ($locS =~ s/[<>]?(\d+)\.\.[<>]?(\d+)/\[$n\]/o) {
394 : efrank 1.1 push(@pieces,["loc",$1,$2]);
395 : gdpusch 1.27 ++$n;
396 : efrank 1.1 }
397 : gdpusch 1.10
398 : gdpusch 1.16 if ($locS =~ s/([,\(])(\d+)([,\)])/$1\[$n\]$3/o) {
399 : efrank 1.1 push(@pieces,["loc",$2,$2]);
400 : gdpusch 1.27 ++$n;
401 : efrank 1.1 }
402 : gdpusch 1.10
403 : gdpusch 1.16 if ($locS =~ s/join\((\[\d+\](,\[\d+\])*)\)/\[$n\]/o) {
404 : efrank 1.1 $args = $1;
405 : gdpusch 1.16 push(@pieces,["join",map { $_ =~ /^\[(\d+)\]$/; $1 } split(m/,/,$args)]);
406 : gdpusch 1.27 ++$n;
407 : efrank 1.1 }
408 : gdpusch 1.10
409 : gdpusch 1.16 if ($locS =~ s/complement\((\[\d+\](,\[\d+\])*)\)/\[$n\]/og) {
410 : efrank 1.1 $args = $1;
411 : gdpusch 1.16 push(@pieces,["complement", map { $_ =~ m/^\[(\d+)\]$/o; $1 } split(m/,/o, $args)]);
412 : gdpusch 1.27 ++$n;
413 : efrank 1.1 }
414 :     }
415 : gdpusch 1.10
416 : gdpusch 1.16 if ($locS =~ m/^\[(\d+)\]$/o) {
417 : efrank 1.1 $loc = &conv(\@pieces,$1,$contigID);
418 :     }
419 : gdpusch 1.16 else {
420 : efrank 1.1 print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
421 : gdpusch 1.18 die qq(aborted);
422 : efrank 1.1 }
423 : gdpusch 1.10
424 : gdpusch 1.16 my @locs = split(m/,/o, $loc);
425 : gdpusch 1.10 #...STOP is now included, so don't trim it off...
426 :     # &trim_stop(\@locs);
427 : efrank 1.1 $loc = join(",",@locs);
428 :     }
429 : gdpusch 1.16 else {
430 : gdpusch 1.18 print STDERR &Dumper($cdsP); die qq(could not parse location);
431 :     die qq(aborted);
432 : efrank 1.1 }
433 :     return ($loc,$precise);
434 :     }
435 :    
436 :     sub trim_stop {
437 :     my($locs) = @_;
438 :     my($to_trim,$n);
439 : gdpusch 1.16
440 : efrank 1.1 $to_trim = 3;
441 : gdpusch 1.16 while ($to_trim && (@$locs > 0)) {
442 : efrank 1.1 $n = @$locs - 1;
443 : gdpusch 1.16 if ($locs->[$n] =~ m/^(\S+)_(\d+)_(\d+)$/o) {
444 :     if ($2 <= ($3-$to_trim)) {
445 : efrank 1.1 $locs->[$n] = join("_",($1,$2,$3-$to_trim));
446 :     $to_trim = 0;
447 :     }
448 : gdpusch 1.16 elsif ($2 >= ($3 + $to_trim)) {
449 : efrank 1.1 $locs->[$n] = join("_",($1,$2,$3+$to_trim));
450 :     $to_trim = 0;
451 :     }
452 : gdpusch 1.16 else {
453 : efrank 1.1 $to_trim -= abs($3-$2) + 1;
454 :     pop @$locs;
455 :     }
456 :     }
457 : gdpusch 1.16 else {
458 : gdpusch 1.18 die qq(could not parse $locs->[$n]);
459 : efrank 1.1 }
460 :     }
461 :     }
462 :    
463 :    
464 :     sub conv {
465 :     my($pieces,$n,$contigID) = @_;
466 : gdpusch 1.16
467 : gdpusch 1.18 if ($pieces->[$n]->[0] eq qq(loc)) {
468 : efrank 1.1 return join("_",$contigID,@{$pieces->[$n]}[1..2]);
469 :     }
470 : gdpusch 1.18 elsif ($pieces->[$n]->[0] eq qq(join)) {
471 : efrank 1.1 return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
472 :     }
473 : gdpusch 1.18 elsif ($pieces->[$n]->[0] eq qq(complement)) {
474 : gdpusch 1.16 return join(",",&complement(join(",", map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
475 : efrank 1.1 }
476 :     }
477 :    
478 :     sub complement {
479 :     my($loc) = @_;
480 :     my(@locs);
481 :    
482 :     @locs = reverse split(/,/,$loc);
483 : gdpusch 1.16 foreach $loc (@locs) {
484 :     if ($loc =~ m/^(\S+)_(\d+)_(\d+)$/o) {
485 :     $loc = join("_",($1,$3,$2));
486 :     }
487 :     else {
488 : gdpusch 1.18 confess qq(Bad location: $loc);
489 : gdpusch 1.16 }
490 : efrank 1.1 }
491 :     return join(",",@locs);
492 :     }
493 :    
494 :     sub get_aliases {
495 :     my($cdsP) = @_;
496 : gdpusch 1.18 my($id, $prefix, $alias, $db_ref);
497 : gdpusch 1.16
498 : efrank 1.1 my @aliases = ();
499 : gdpusch 1.16 while ($$cdsP =~ m/\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/ogs) {
500 : gdpusch 1.20 ($type, $alias) = ($1, $2);
501 : gdpusch 1.18
502 :     # define prefixes for different types of ids
503 :     if ($type eq qq(locus_tag)){
504 :     $id = qq(locus|$alias);
505 :     }
506 :     elsif ( $type eq qq(protein_id) ) {
507 :     $id = qq(protein_id|$alias);
508 :     }
509 :     elsif ( $type eq qq(gene) ){
510 :     $id = qq(gene_name|$alias);
511 :     }
512 :     else{
513 :     $id = $alias;
514 :     }
515 :    
516 :     push(@aliases,$id);
517 : efrank 1.1 }
518 :    
519 : gdpusch 1.16 while ($$cdsP =~ m/\/db_xref=\"([^\"]+)\"/ogs) {
520 : efrank 1.1 $db_ref = $1;
521 : gdpusch 1.25 $db_ref =~ s/[\s\n]+//ogs;
522 : gdpusch 1.16 $db_ref =~ s/^GI:/gi\|/o;
523 : wilke 1.17 $db_ref =~ s/^GeneID:/geneID\|/o;
524 : gdpusch 1.16 $db_ref =~ s/SWISS-PROT:/sp\|/o;
525 : efrank 1.1 push(@aliases,$db_ref);
526 :     }
527 : gdpusch 1.16
528 : efrank 1.1 return @aliases;
529 :     }
530 :    
531 : gdpusch 1.18 sub get_ec_numbers {
532 :     my ($cdsP) = @_;
533 :     my @ec_numbers = ($$cdsP =~ m{/EC_number=\"([0-9]+\.[0-9\-]+\.[0-9\-]+\.[0-9\-]+)\"}ogs);
534 :     return @ec_numbers;
535 :     }
536 :    
537 : efrank 1.1 sub get_trans {
538 :     my($cdsP) = @_;
539 : gdpusch 1.18 my $tran = qq();
540 : gdpusch 1.26
541 :     if ($$cdsP =~ m/\/translation=\"([^\"]*)\"/os) {
542 : efrank 1.1 $tran = $1;
543 :     $tran =~ s/\s//gs;
544 :     }
545 : gdpusch 1.26
546 :     # elsif ($$cdsP =~ m/\/protein_id=\"([^\"]+)\"/o) {
547 :     # $tran = $1;
548 :     # $tran =~ s/\s//ogs;
549 :     # }
550 :    
551 : efrank 1.1 return $tran;
552 :     }
553 :    
554 :     sub get_func {
555 :     my($cdsP) = @_;
556 : overbeek 1.11
557 : gdpusch 1.18 my $functions = [];
558 :     my $products = [];
559 :     my $prot_descs = [];
560 :     my $notes = [];
561 :    
562 :     print STDERR qq(\nRecord $recnum:\n$$cdsP\n) if $ENV{VERBOSE};
563 :    
564 : gdpusch 1.20 @$functions = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/function=\"([^\"]*)\"/ogs);
565 : gdpusch 1.18 if ($ENV{VERBOSE} && @$functions) {
566 :     print STDERR (qq(Functions: ), ((@$functions > 1) ? qq(\n) : qq()));
567 :     print STDERR (join(qq(\n), @$functions), qq(\n));
568 :     }
569 :    
570 : gdpusch 1.20 @$products = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/product=\"([^\"]*)\"/ogs);
571 : gdpusch 1.18 if ($ENV{VERBOSE} && @$products) {
572 :     print STDERR (qq(Products: ), ((@$products > 1) ? qq(\n) : qq()));
573 :     print STDERR (join(qq(\n), @$products), qq(\n));
574 : overbeek 1.4 }
575 : gdpusch 1.18
576 : gdpusch 1.20 @$prot_descs = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/prot_desc=\"([^\"]*)\"/ogs);
577 : gdpusch 1.18 if ($ENV{VERBOSE} && @$prot_descs) {
578 :     print STDERR (qq(Prot_Descs: ), ((@$prot_descs > 1) ? qq(\n) : qq()));
579 :     print STDERR (join(qq(\n), @$prot_descs), qq(\n));
580 :     }
581 :    
582 :     @$notes = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/note=\"([^\"]*)\"/ogs);
583 :     if ($ENV{VERBOSE} && @$notes) {
584 :     print STDERR (qq(Notes: ), ((@$notes > 1) ? qq(\n) : qq()));
585 :     print STDERR (join(qq(\n), @$notes), qq(\n));
586 :     }
587 :    
588 :     @$products = grep { !m/^hypothetical\s+protein$/io } @$products;
589 :    
590 :     my $func = qq();
591 :     my $annotations = [];
592 :     if (@$products) {
593 :     $func = join(q( / ), @$products);
594 :     @$annotations = (@$functions, @$prot_descs, @$notes);
595 :     }
596 :     elsif (@$functions) {
597 :     $func =join(q( / ), @$functions);
598 :     @$annotations =(@$prot_descs, @$notes);
599 :     }
600 :     elsif (@$prot_descs) {
601 :     $func =join(q( / ), @$prot_descs);
602 :     @$annotations = @$notes;
603 : overbeek 1.4 }
604 : gdpusch 1.16 else {
605 : gdpusch 1.18 $func = qq(hypothetical protein);
606 :     @$annotations = @$notes;
607 : efrank 1.2 }
608 : gdpusch 1.18
609 :     return ($func, $annotations);
610 : overbeek 1.4 }
611 :    
612 :     sub fixup_func {
613 :     my($func) = @_;
614 : gdpusch 1.16
615 :     $func =~ s/^COG\d+:\s+//oi;
616 :     $func =~ s/^ORFID:\S+\s+//oi;
617 : efrank 1.2 return $func;
618 :     }
619 :    
620 : overbeek 1.4 sub ok_func {
621 :     my($func) = @_;
622 : gdpusch 1.16
623 : overbeek 1.4 return (
624 : gdpusch 1.16 ($func !~ m/^[a-zA-Z]{1,3}\d+[a-zA-Z]{0,3}(\.\d+([a-zA-Z])?)?$/o) &&
625 :     ($func !~ m/^\d+$/o) &&
626 :     ($func !~ m/ensangp/oi) &&
627 :     ($func !~ m/agr_/oi) &&
628 :     ($func !~ m/^SC.*:/o) &&
629 :     ($func !~ m/^RIKEN/o) &&
630 :     ($func !~ m/\slen: \d+/o) &&
631 :     ($func !~ m/^similar to/oi) &&
632 :     ($func !~ m/^CPX/oi) &&
633 :     ($func !~ m/^\S{3,4}( or \S+ protein)?$/oi) &&
634 :     ($func !~ m/^putative$/o) &&
635 :     ($func !~ m/\scomes from/oi) &&
636 :     ($func !~ m/^unknown( function)?$/oi) &&
637 :     ($func !~ m/^hypothetical( (function|protein))?$/oi) &&
638 :     ($func !~ m/^orf /oi) &&
639 :     ($func !~ m/^ORF\s?\d+[lr]?$/oi) &&
640 :     ($func !~ m/^[a-z]{1,3}\s?\d+$/oi)
641 :     );
642 : overbeek 1.4 }
643 :    
644 : efrank 1.2 sub process_rna {
645 : gdpusch 1.28 my ($contigID, $type, $rnaP, $prefix, $idNr, $contigs, $fh_tbl, $fh_dna, $fh_func, $fh_annot) = @_;
646 : gdpusch 1.27 my ($loc, @aliases, $func, $trans, $id, $precise, $DNA);
647 : gdpusch 1.9
648 : gdpusch 1.27 ($loc, $precise) = &get_loc($contigID,$rnaP);
649 :     $func = &get_rna_func($rnaP);
650 :    
651 :     if ($loc && $precise) {
652 : olson 1.30 if ($dna = &SeedUtils::extract_seq($contigs,$loc)) {
653 : gdpusch 1.18 $id = $prefix . qq($$idNr);
654 : gdpusch 1.27 ++$$idNr;
655 : gdpusch 1.16
656 :     if (! $func ) {
657 : gdpusch 1.27 warn qq(WARNING: could not get func $$rnaP\n) if $ENV{VERBOSE};
658 : gdpusch 1.9 }
659 : gdpusch 1.16 else {
660 : gdpusch 1.28 print $fh_func qq($id\t$func\n);
661 :     print $fh_tbl qq($id\t$loc\t\n);
662 :    
663 :     &make_annotation($id, qq(Set function to\n$func\nInitial import.), $fh_annot);
664 : olson 1.30 &SeedUtils::display_id_and_seq($id, \$dna, $fh_dna);
665 : gdpusch 1.9 }
666 : efrank 1.2 }
667 : gdpusch 1.16 else {
668 : gdpusch 1.27 warn (qq(Could not get DNA sequence for RNA at \'$loc\' for entry\n),
669 :     , q( ), $type, q( )
670 :     , $$rnaP
671 :     , qq(\nof record begining with:\n)
672 :     , substr($record, 0, 160)
673 :     , qq(\n\n)
674 :     );
675 : efrank 1.2 }
676 :     }
677 : gdpusch 1.16 else {
678 : gdpusch 1.27 warn (qq(Could not handle RNA at \'$loc\' for entry\n),
679 :     , q( ), $type, q( )
680 :     , $$rnaP
681 :     , qq(\nof record begining with:\n)
682 :     , substr($record, 0, 160)
683 :     , qq(\n\n)
684 :     );
685 : efrank 1.2 }
686 :     }
687 :    
688 :     sub get_rna_func {
689 :     my($cdsP) = @_;
690 : gdpusch 1.18 my $func = qq();
691 : gdpusch 1.16
692 : gdpusch 1.26 if ($$cdsP =~ m/\/product=\"([^\"]*)\"/os) {
693 : efrank 1.2 $func = $1;
694 :     $func =~ s/\s+/ /gs;
695 :     }
696 : gdpusch 1.26 elsif ($$cdsP =~ m/\/gene=\"([^\"]*)\"/os) {
697 : efrank 1.2 $func = $1;
698 : gdpusch 1.16 $func =~ s/\s+/ /ogs;
699 : efrank 1.2 }
700 : gdpusch 1.26 elsif ($$cdsP =~ m/\/note=\"([^\"]*)\"/os) {
701 : efrank 1.1 $func = $1;
702 : gdpusch 1.16 $func =~ s/\s+/ /ogs;
703 : efrank 1.1 }
704 :     return $func;
705 :     }
706 : gdpusch 1.18
707 :     sub make_annotation {
708 :     my ($fid, $note, $fh_ann) = @_;
709 :     # print STDERR Dumper($fid, $note, $fh_ann);
710 :    
711 :     print $fh_ann ($fid, qq(\n));
712 :     print $fh_ann (time, qq(\n));
713 :     print $fh_ann qq(parse_genbank\n);
714 :     print $fh_ann ($note, qq(\n));
715 :     print $fh_ann qq(//\n);
716 :    
717 :     return 1;
718 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3