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

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.33 - (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.16
190 : overbeek 1.14 my @genome = ();
191 :     my @full_tax = ();
192 : gdpusch 1.16 for ($i=0; ($i < @lines) && ($lines[$i] !~ /;/); $i++) {
193 : overbeek 1.14 push(@genome,$lines[$i]);
194 :     }
195 : gdpusch 1.16
196 :     while ($i < @lines) {
197 : overbeek 1.14 push(@full_tax,$lines[$i]);
198 : gdpusch 1.27 ++$i;
199 : overbeek 1.14 }
200 : gdpusch 1.16
201 : gdpusch 1.18 $genome = join(qq( ),map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @genome);
202 :     $taxonomy = join(qq( ),map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @full_tax);
203 : gdpusch 1.16
204 : gdpusch 1.18 $taxonomy =~ s/\n\s+//og;
205 :     $taxonomy =~ s/ {2,}/ /og;
206 :     $taxonomy =~ s/\.$//o;
207 :     $taxonomy = $taxonomy . qq(; $genome);
208 : overbeek 1.15
209 : gdpusch 1.16 if (! $written_genome) {
210 : gdpusch 1.20 if ($force_bioname) { $genome = $force_bioname; }
211 :     print GENOME qq($genome\n);
212 : overbeek 1.4 close(GENOME);
213 : gdpusch 1.9
214 : gdpusch 1.18 if ($force_taxonomy) { $taxonomy = $force_taxonomy; }
215 :     print TAXONOMY qq($taxonomy\n);
216 :     close(TAXONOMY);
217 :    
218 : gdpusch 1.21 $written_genome = qq($taxonomy\t$genome);
219 : overbeek 1.4 }
220 : gdpusch 1.18 elsif (($written_genome ne qq($taxonomy\t$bioname)) && (!$force_bioname || !$force_taxonomy)) {
221 :     print STDERR qq(serious mismatch in GENOME/TAX for $id\n$written_genome\n$taxonomy\t$bioname\n\n);
222 : overbeek 1.4 }
223 : efrank 1.1 }
224 : gdpusch 1.16
225 :     while ($record =~ m/\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/ogs) {
226 : gdpusch 1.18 my $cds = $1;
227 : gdpusch 1.16 if (($cds !~ m/\/pseudo/o) &&
228 :     (($cds !~ m/\/exception/o) || ($cds =~ m/\/translation/o))
229 :     ) {
230 : gdpusch 1.18 &process_cds($id, \$cds, $prefixP, \$idNp, $contigs, $fh_peg_tbl, $fh_peg_fasta, $fh_assigned_funcs, $fh_annotations, $fh_ec_nums);
231 : overbeek 1.7 }
232 : efrank 1.2 }
233 :    
234 : gdpusch 1.27 while ($record =~ m/\n\s{3,6}(([tr]|misc\_)RNA)\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/ogs)
235 : efrank 1.2 {
236 : gdpusch 1.27 $type = $2;
237 :     $rna = $3;
238 : gdpusch 1.28 &process_rna($id, $type, \$rna, $prefixR, \$idNr, $contigs, $fh_rna_tbl, $fh_rna_fasta, $fh_assigned_funcs, $fh_annotations);
239 : efrank 1.1 }
240 :     }
241 :     }
242 : gdpusch 1.18 close($fh_contigs);
243 :     close($fh_peg_tbl);
244 :     close($fh_peg_fasta);
245 :     close($fh_rna_tbl);
246 :     close($fh_rna_fasta);
247 :     close($fh_assigned_funcs);
248 :     close($fh_annotations);
249 :     close($fh_ec_nums);
250 :    
251 :     if (!-s $assigned_funcs_file) {
252 : gdpusch 1.20 warn qq(WARNING: no assigned_functions in $org_dir --- deleting\n);
253 : gdpusch 1.18 unlink("$assigned_funcs_file") unless $ENV{DEBUG};
254 :     }
255 :    
256 : gdpusch 1.22 if (!-s $annotations_file) {
257 :     warn qq(WARNING: no annotations in $org_dir --- deleting\n);
258 :     unlink("$annotations_file") unless $ENV{DEBUG};
259 :     }
260 :    
261 : gdpusch 1.18 if (!-s $ec_nums_file) {
262 : gdpusch 1.20 warn qq(WARNING: no EC_numbers in $org_dir --- deleting\n);
263 : gdpusch 1.18 unlink("$ec_nums_file") unless $ENV{DEBUG};
264 :     }
265 : efrank 1.2
266 : gdpusch 1.20 if ((!-s $rna_tbl_file) || (!-s $rna_fasta_file)) {
267 :     warn qq(WARNING: no RNAs in $org_dir --- deleting\n);
268 : olson 1.29 system('rm', '-rf', "$org_dir/Features/rna") unless $ENV{DEBUG};;
269 : gdpusch 1.20 }
270 :    
271 : gdpusch 1.18 if ((!-s $peg_tbl_file) || (!-s $peg_fasta_file)) {
272 : gdpusch 1.20 warn qq(WARNING: no PEGs in $org_dir --- deleting\n);
273 : olson 1.29 system('rm', '-rf', "$org_dir/Features/peg") unless $ENV{DEBUG};
274 : gdpusch 1.18 }
275 :    
276 : gdpusch 1.20 if ((!-d qq($org_dir/Features/peg)) && (!-d qq($org_dir/Features/peg))) {
277 :     warn qq(WARNING: no Features in $org_dir --- deleting\n);
278 : olson 1.29 system("rmdir", "$org_dir/Features") unless $ENV{DEBUG};;
279 : gdpusch 1.18 }
280 : gdpusch 1.20
281 : gdpusch 1.18 if (!-s $contigs_file) {
282 :     $trouble = 1;
283 :     unlink($contigs_file);
284 : gdpusch 1.20 print STDERR qq(ERROR: No contigs in $org_dir\n);
285 : gdpusch 1.18 }
286 :    
287 :     if ($trouble) {
288 :     warn qq(Genome directory $org_dir is corrupt --- deleting\n\n);
289 : olson 1.29 system('rm', '-rf', $org_dir) unless $ENV{DEBUG};
290 : gdpusch 1.18 }
291 :    
292 :     exit($trouble);
293 : efrank 1.1
294 :     sub process_cds {
295 : gdpusch 1.20 my ($contigID, $cdsP, $prefix, $idNp, $contigs, $fh_tbl, $fh_fasta, $fh_assign, $fh_annot, $fh_ec_nums) = @_;
296 : gdpusch 1.18 my ($id, $prot);
297 : gdpusch 1.9
298 : overbeek 1.13 ++$recnum;
299 : gdpusch 1.18 my ($loc, $precise) = &get_loc($contigID,$cdsP);
300 :     my @aliases = &get_aliases($cdsP);
301 :     my @ec_nums = &get_ec_numbers($cdsP);
302 :     my ($func, $notes) = &get_func($cdsP);
303 :     my $trans = &get_trans($cdsP);
304 :    
305 : gdpusch 1.26 if (not $trans) {
306 : gdpusch 1.33 if ($no_trans) {
307 :     #...Do nothing
308 :     }
309 :     else {
310 :     warn qq(WARNING: Translation missing for CDS $recnum; attempting to generate translation\n) if $ENV{VERBOSE};
311 :    
312 :     if ($loc && $precise) {
313 :     my $dna = &SeedUtils::extract_seq($contigs,$loc);
314 :     if ($dna) {
315 :     my $genetic_code = 11;
316 :     if ($$cdsP =~ m/\/transl_table=(\d+)/o) {
317 :     $genetic_code = $1;
318 :     }
319 :    
320 :     if (not defined($trans_table{$genetic_code})) {
321 :     $trans_table{$genetic_code} = &SeedUtils::genetic_code($genetic_code);
322 :     }
323 :    
324 :     my $prot = &SeedUtils::translate($dna, $trans_table{$genetic_code}, 1);
325 :     $prot =~ s/\*$//o;
326 :    
327 :     if ($prot !~ m/\*/o) {
328 :     $trans = $prot;
329 :     }
330 :     else {
331 :     warn qq(Translation contains STOPs, changing to \'x\'s, for CDS $recnum:\n$$cdsP\n);
332 :     }
333 : gdpusch 1.9 }
334 : efrank 1.1 }
335 : gdpusch 1.33 else {
336 :     warn (qq(Could not get DNA sequence for CDS at $loc for entry:\n)
337 :     , q( CDS )
338 :     , $$cdsP
339 :     , qq(\nof record begining with:\n)
340 :     , substr($record, 0, 160)
341 :     , qq(\n\n)
342 :     );
343 :     }
344 : gdpusch 1.26 }
345 :     }
346 :    
347 :     if ($trans) {
348 : gdpusch 1.27 $id = $prefix . qq($$idNp);
349 :     ++$$idNp;
350 :    
351 : olson 1.19 print $fh_tbl qq($id\t$loc\t) . join("\t",@aliases) . "\n";
352 : gdpusch 1.16
353 :     if ($func) {
354 : gdpusch 1.20 print $fh_assign qq($id\t$func\n);
355 : efrank 1.1 }
356 : gdpusch 1.16
357 : gdpusch 1.18 if (@ec_nums) {
358 :     print $fh_ec_nums (join(qq(\t), ($id, @ec_nums)), qq(\n));
359 :     }
360 :    
361 : gdpusch 1.26 foreach my $note (@$notes) {
362 :     &make_annotation($id, $note, $fh_annot);
363 :     }
364 : gdpusch 1.18
365 : olson 1.30 &SeedUtils::display_id_and_seq($id, \$trans, $fh_fasta);
366 : efrank 1.1 }
367 : gdpusch 1.16 else {
368 : gdpusch 1.27 warn (qq(No translation for CDS --- skipping entry:\n), &Dumper($cdsP), qq(\n));
369 : efrank 1.1 }
370 :     }
371 :    
372 :     sub get_loc {
373 :     my($contigID,$cdsP) = @_;
374 :     my($beg,$end,$loc,$locS,$iter,$n,$args,@pieces);
375 :     my($precise);
376 : gdpusch 1.16
377 :     if ($$cdsP =~ m/^(([^\n]+)((\n\s{21,23}[^\/ ][^\n]+)+)?)/os) {
378 : efrank 1.1 $locS = $1;
379 :     $locS =~ s/\s//g;
380 : gdpusch 1.16 $precise = ($locS !~ m/[<>]/o);
381 : gdpusch 1.10
382 : efrank 1.1 @pieces = ();
383 :     $n = 0;
384 :     $iter = 500;
385 : gdpusch 1.16 while (($locS !~ m/^\[\d+\]$/o) && $iter)
386 : efrank 1.1 {
387 : gdpusch 1.27 --$iter;
388 : gdpusch 1.16 if ($locS =~ s/[<>]?(\d+)\.\.[<>]?(\d+)/\[$n\]/o) {
389 : efrank 1.1 push(@pieces,["loc",$1,$2]);
390 : gdpusch 1.27 ++$n;
391 : efrank 1.1 }
392 : gdpusch 1.10
393 : gdpusch 1.16 if ($locS =~ s/([,\(])(\d+)([,\)])/$1\[$n\]$3/o) {
394 : efrank 1.1 push(@pieces,["loc",$2,$2]);
395 : gdpusch 1.27 ++$n;
396 : efrank 1.1 }
397 : gdpusch 1.10
398 : gdpusch 1.16 if ($locS =~ s/join\((\[\d+\](,\[\d+\])*)\)/\[$n\]/o) {
399 : efrank 1.1 $args = $1;
400 : gdpusch 1.16 push(@pieces,["join",map { $_ =~ /^\[(\d+)\]$/; $1 } split(m/,/,$args)]);
401 : gdpusch 1.27 ++$n;
402 : efrank 1.1 }
403 : gdpusch 1.10
404 : gdpusch 1.16 if ($locS =~ s/complement\((\[\d+\](,\[\d+\])*)\)/\[$n\]/og) {
405 : efrank 1.1 $args = $1;
406 : gdpusch 1.16 push(@pieces,["complement", map { $_ =~ m/^\[(\d+)\]$/o; $1 } split(m/,/o, $args)]);
407 : gdpusch 1.27 ++$n;
408 : efrank 1.1 }
409 :     }
410 : gdpusch 1.10
411 : gdpusch 1.16 if ($locS =~ m/^\[(\d+)\]$/o) {
412 : efrank 1.1 $loc = &conv(\@pieces,$1,$contigID);
413 :     }
414 : gdpusch 1.16 else {
415 : efrank 1.1 print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
416 : gdpusch 1.18 die qq(aborted);
417 : efrank 1.1 }
418 : gdpusch 1.10
419 : gdpusch 1.16 my @locs = split(m/,/o, $loc);
420 : gdpusch 1.10 #...STOP is now included, so don't trim it off...
421 :     # &trim_stop(\@locs);
422 : efrank 1.1 $loc = join(",",@locs);
423 :     }
424 : gdpusch 1.16 else {
425 : gdpusch 1.18 print STDERR &Dumper($cdsP); die qq(could not parse location);
426 :     die qq(aborted);
427 : efrank 1.1 }
428 :     return ($loc,$precise);
429 :     }
430 :    
431 :     sub trim_stop {
432 :     my($locs) = @_;
433 :     my($to_trim,$n);
434 : gdpusch 1.16
435 : efrank 1.1 $to_trim = 3;
436 : gdpusch 1.16 while ($to_trim && (@$locs > 0)) {
437 : efrank 1.1 $n = @$locs - 1;
438 : gdpusch 1.16 if ($locs->[$n] =~ m/^(\S+)_(\d+)_(\d+)$/o) {
439 :     if ($2 <= ($3-$to_trim)) {
440 : efrank 1.1 $locs->[$n] = join("_",($1,$2,$3-$to_trim));
441 :     $to_trim = 0;
442 :     }
443 : gdpusch 1.16 elsif ($2 >= ($3 + $to_trim)) {
444 : efrank 1.1 $locs->[$n] = join("_",($1,$2,$3+$to_trim));
445 :     $to_trim = 0;
446 :     }
447 : gdpusch 1.16 else {
448 : efrank 1.1 $to_trim -= abs($3-$2) + 1;
449 :     pop @$locs;
450 :     }
451 :     }
452 : gdpusch 1.16 else {
453 : gdpusch 1.18 die qq(could not parse $locs->[$n]);
454 : efrank 1.1 }
455 :     }
456 :     }
457 :    
458 :    
459 :     sub conv {
460 :     my($pieces,$n,$contigID) = @_;
461 : gdpusch 1.16
462 : gdpusch 1.18 if ($pieces->[$n]->[0] eq qq(loc)) {
463 : efrank 1.1 return join("_",$contigID,@{$pieces->[$n]}[1..2]);
464 :     }
465 : gdpusch 1.18 elsif ($pieces->[$n]->[0] eq qq(join)) {
466 : efrank 1.1 return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
467 :     }
468 : gdpusch 1.18 elsif ($pieces->[$n]->[0] eq qq(complement)) {
469 : gdpusch 1.16 return join(",",&complement(join(",", map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
470 : efrank 1.1 }
471 :     }
472 :    
473 :     sub complement {
474 :     my($loc) = @_;
475 :     my(@locs);
476 :    
477 :     @locs = reverse split(/,/,$loc);
478 : gdpusch 1.16 foreach $loc (@locs) {
479 :     if ($loc =~ m/^(\S+)_(\d+)_(\d+)$/o) {
480 :     $loc = join("_",($1,$3,$2));
481 :     }
482 :     else {
483 : gdpusch 1.18 confess qq(Bad location: $loc);
484 : gdpusch 1.16 }
485 : efrank 1.1 }
486 :     return join(",",@locs);
487 :     }
488 :    
489 :     sub get_aliases {
490 :     my($cdsP) = @_;
491 : gdpusch 1.18 my($id, $prefix, $alias, $db_ref);
492 : gdpusch 1.16
493 : efrank 1.1 my @aliases = ();
494 : gdpusch 1.16 while ($$cdsP =~ m/\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/ogs) {
495 : gdpusch 1.20 ($type, $alias) = ($1, $2);
496 : gdpusch 1.18
497 :     # define prefixes for different types of ids
498 :     if ($type eq qq(locus_tag)){
499 :     $id = qq(locus|$alias);
500 :     }
501 :     elsif ( $type eq qq(protein_id) ) {
502 :     $id = qq(protein_id|$alias);
503 :     }
504 :     elsif ( $type eq qq(gene) ){
505 :     $id = qq(gene_name|$alias);
506 :     }
507 :     else{
508 :     $id = $alias;
509 :     }
510 :    
511 :     push(@aliases,$id);
512 : efrank 1.1 }
513 :    
514 : gdpusch 1.16 while ($$cdsP =~ m/\/db_xref=\"([^\"]+)\"/ogs) {
515 : efrank 1.1 $db_ref = $1;
516 : gdpusch 1.25 $db_ref =~ s/[\s\n]+//ogs;
517 : gdpusch 1.16 $db_ref =~ s/^GI:/gi\|/o;
518 : wilke 1.17 $db_ref =~ s/^GeneID:/geneID\|/o;
519 : gdpusch 1.16 $db_ref =~ s/SWISS-PROT:/sp\|/o;
520 : efrank 1.1 push(@aliases,$db_ref);
521 :     }
522 : gdpusch 1.16
523 : efrank 1.1 return @aliases;
524 :     }
525 :    
526 : gdpusch 1.18 sub get_ec_numbers {
527 :     my ($cdsP) = @_;
528 :     my @ec_numbers = ($$cdsP =~ m{/EC_number=\"([0-9]+\.[0-9\-]+\.[0-9\-]+\.[0-9\-]+)\"}ogs);
529 :     return @ec_numbers;
530 :     }
531 :    
532 : efrank 1.1 sub get_trans {
533 :     my($cdsP) = @_;
534 : gdpusch 1.18 my $tran = qq();
535 : gdpusch 1.26
536 :     if ($$cdsP =~ m/\/translation=\"([^\"]*)\"/os) {
537 : efrank 1.1 $tran = $1;
538 :     $tran =~ s/\s//gs;
539 :     }
540 : gdpusch 1.26
541 :     # elsif ($$cdsP =~ m/\/protein_id=\"([^\"]+)\"/o) {
542 :     # $tran = $1;
543 :     # $tran =~ s/\s//ogs;
544 :     # }
545 :    
546 : efrank 1.1 return $tran;
547 :     }
548 :    
549 :     sub get_func {
550 :     my($cdsP) = @_;
551 : overbeek 1.11
552 : gdpusch 1.18 my $functions = [];
553 :     my $products = [];
554 :     my $prot_descs = [];
555 :     my $notes = [];
556 :    
557 :     print STDERR qq(\nRecord $recnum:\n$$cdsP\n) if $ENV{VERBOSE};
558 :    
559 : gdpusch 1.20 @$functions = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/function=\"([^\"]*)\"/ogs);
560 : gdpusch 1.18 if ($ENV{VERBOSE} && @$functions) {
561 :     print STDERR (qq(Functions: ), ((@$functions > 1) ? qq(\n) : qq()));
562 :     print STDERR (join(qq(\n), @$functions), qq(\n));
563 :     }
564 :    
565 : gdpusch 1.20 @$products = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/product=\"([^\"]*)\"/ogs);
566 : gdpusch 1.18 if ($ENV{VERBOSE} && @$products) {
567 :     print STDERR (qq(Products: ), ((@$products > 1) ? qq(\n) : qq()));
568 :     print STDERR (join(qq(\n), @$products), qq(\n));
569 : overbeek 1.4 }
570 : gdpusch 1.18
571 : gdpusch 1.20 @$prot_descs = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/prot_desc=\"([^\"]*)\"/ogs);
572 : gdpusch 1.18 if ($ENV{VERBOSE} && @$prot_descs) {
573 :     print STDERR (qq(Prot_Descs: ), ((@$prot_descs > 1) ? qq(\n) : qq()));
574 :     print STDERR (join(qq(\n), @$prot_descs), qq(\n));
575 :     }
576 :    
577 :     @$notes = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/note=\"([^\"]*)\"/ogs);
578 :     if ($ENV{VERBOSE} && @$notes) {
579 :     print STDERR (qq(Notes: ), ((@$notes > 1) ? qq(\n) : qq()));
580 :     print STDERR (join(qq(\n), @$notes), qq(\n));
581 :     }
582 :    
583 :     @$products = grep { !m/^hypothetical\s+protein$/io } @$products;
584 :    
585 :     my $func = qq();
586 :     my $annotations = [];
587 :     if (@$products) {
588 :     $func = join(q( / ), @$products);
589 :     @$annotations = (@$functions, @$prot_descs, @$notes);
590 :     }
591 :     elsif (@$functions) {
592 :     $func =join(q( / ), @$functions);
593 :     @$annotations =(@$prot_descs, @$notes);
594 :     }
595 :     elsif (@$prot_descs) {
596 :     $func =join(q( / ), @$prot_descs);
597 :     @$annotations = @$notes;
598 : overbeek 1.4 }
599 : gdpusch 1.16 else {
600 : gdpusch 1.18 $func = qq(hypothetical protein);
601 :     @$annotations = @$notes;
602 : efrank 1.2 }
603 : gdpusch 1.18
604 :     return ($func, $annotations);
605 : overbeek 1.4 }
606 :    
607 :     sub fixup_func {
608 :     my($func) = @_;
609 : gdpusch 1.16
610 :     $func =~ s/^COG\d+:\s+//oi;
611 :     $func =~ s/^ORFID:\S+\s+//oi;
612 : efrank 1.2 return $func;
613 :     }
614 :    
615 : overbeek 1.4 sub ok_func {
616 :     my($func) = @_;
617 : gdpusch 1.16
618 : overbeek 1.4 return (
619 : gdpusch 1.16 ($func !~ m/^[a-zA-Z]{1,3}\d+[a-zA-Z]{0,3}(\.\d+([a-zA-Z])?)?$/o) &&
620 :     ($func !~ m/^\d+$/o) &&
621 :     ($func !~ m/ensangp/oi) &&
622 :     ($func !~ m/agr_/oi) &&
623 :     ($func !~ m/^SC.*:/o) &&
624 :     ($func !~ m/^RIKEN/o) &&
625 :     ($func !~ m/\slen: \d+/o) &&
626 :     ($func !~ m/^similar to/oi) &&
627 :     ($func !~ m/^CPX/oi) &&
628 :     ($func !~ m/^\S{3,4}( or \S+ protein)?$/oi) &&
629 :     ($func !~ m/^putative$/o) &&
630 :     ($func !~ m/\scomes from/oi) &&
631 :     ($func !~ m/^unknown( function)?$/oi) &&
632 :     ($func !~ m/^hypothetical( (function|protein))?$/oi) &&
633 :     ($func !~ m/^orf /oi) &&
634 :     ($func !~ m/^ORF\s?\d+[lr]?$/oi) &&
635 :     ($func !~ m/^[a-z]{1,3}\s?\d+$/oi)
636 :     );
637 : overbeek 1.4 }
638 :    
639 : efrank 1.2 sub process_rna {
640 : gdpusch 1.28 my ($contigID, $type, $rnaP, $prefix, $idNr, $contigs, $fh_tbl, $fh_dna, $fh_func, $fh_annot) = @_;
641 : gdpusch 1.27 my ($loc, @aliases, $func, $trans, $id, $precise, $DNA);
642 : gdpusch 1.9
643 : gdpusch 1.27 ($loc, $precise) = &get_loc($contigID,$rnaP);
644 :     $func = &get_rna_func($rnaP);
645 :    
646 :     if ($loc && $precise) {
647 : olson 1.30 if ($dna = &SeedUtils::extract_seq($contigs,$loc)) {
648 : gdpusch 1.18 $id = $prefix . qq($$idNr);
649 : gdpusch 1.27 ++$$idNr;
650 : gdpusch 1.16
651 :     if (! $func ) {
652 : gdpusch 1.27 warn qq(WARNING: could not get func $$rnaP\n) if $ENV{VERBOSE};
653 : gdpusch 1.9 }
654 : gdpusch 1.16 else {
655 : gdpusch 1.28 print $fh_func qq($id\t$func\n);
656 :     print $fh_tbl qq($id\t$loc\t\n);
657 :    
658 :     &make_annotation($id, qq(Set function to\n$func\nInitial import.), $fh_annot);
659 : olson 1.30 &SeedUtils::display_id_and_seq($id, \$dna, $fh_dna);
660 : gdpusch 1.9 }
661 : efrank 1.2 }
662 : gdpusch 1.16 else {
663 : gdpusch 1.27 warn (qq(Could not get DNA sequence for RNA at \'$loc\' for entry\n),
664 :     , q( ), $type, q( )
665 :     , $$rnaP
666 :     , qq(\nof record begining with:\n)
667 :     , substr($record, 0, 160)
668 :     , qq(\n\n)
669 :     );
670 : efrank 1.2 }
671 :     }
672 : gdpusch 1.16 else {
673 : gdpusch 1.27 warn (qq(Could not handle RNA at \'$loc\' for entry\n),
674 :     , q( ), $type, q( )
675 :     , $$rnaP
676 :     , qq(\nof record begining with:\n)
677 :     , substr($record, 0, 160)
678 :     , qq(\n\n)
679 :     );
680 : efrank 1.2 }
681 :     }
682 :    
683 :     sub get_rna_func {
684 :     my($cdsP) = @_;
685 : gdpusch 1.18 my $func = qq();
686 : gdpusch 1.16
687 : gdpusch 1.26 if ($$cdsP =~ m/\/product=\"([^\"]*)\"/os) {
688 : efrank 1.2 $func = $1;
689 :     $func =~ s/\s+/ /gs;
690 :     }
691 : gdpusch 1.26 elsif ($$cdsP =~ m/\/gene=\"([^\"]*)\"/os) {
692 : efrank 1.2 $func = $1;
693 : gdpusch 1.16 $func =~ s/\s+/ /ogs;
694 : efrank 1.2 }
695 : gdpusch 1.26 elsif ($$cdsP =~ m/\/note=\"([^\"]*)\"/os) {
696 : efrank 1.1 $func = $1;
697 : gdpusch 1.16 $func =~ s/\s+/ /ogs;
698 : efrank 1.1 }
699 :     return $func;
700 :     }
701 : gdpusch 1.18
702 :     sub make_annotation {
703 :     my ($fid, $note, $fh_ann) = @_;
704 :     # print STDERR Dumper($fid, $note, $fh_ann);
705 :    
706 :     print $fh_ann ($fid, qq(\n));
707 :     print $fh_ann (time, qq(\n));
708 :     print $fh_ann qq(parse_genbank\n);
709 :     print $fh_ann ($note, qq(\n));
710 :     print $fh_ann qq(//\n);
711 :    
712 :     return 1;
713 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3