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

Annotation of /FigKernelPackages/gjogenbank.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.1 #!/usr/bin/env perl
2 :    
3 :     # This is a SAS component.
4 :    
5 :     package gjogenbank;
6 :    
7 :     #===============================================================================
8 :     #
9 :     # Parse one or more GenBank files into perl structures
10 :     #
11 :     # @entries = parse_genbank( ) # \*STDIN
12 :     # \@entries = parse_genbank( ) # \*STDIN
13 :     # @entries = parse_genbank( \*FH )
14 :     # \@entries = parse_genbank( \*FH )
15 :     # @entries = parse_genbank( $file )
16 :     # \@entries = parse_genbank( $file )
17 :     #
18 :     # Access functions:
19 :     #
20 :     # Sequence of a feature, optionally including information on partial ends.
21 :     #
22 :     # $seq = ftr_dna( $dna, $ftr )
23 :     # $seq = ftr_dna( \$dna, $ftr )
24 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( $dna, $ftr )
25 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( \$dna, $ftr )
26 :     #
27 :     # $ftr_location = location( $ftr ) # Returns empty string on failure.
28 :     #
29 :     # \%ftr_qualifiers = qualifiers( $ftr ) # Returns empty hash reference on failure.
30 :     #
31 :     # $gene = gene( $ftr )
32 :     # @gene_and_synonyms = gene( $ftr )
33 :     #
34 :     # $id = CDS_id( $ftr ) # Prefer protein_id as id:
35 :     # $id = CDS_gi_or_id( $ftr ) # Prefer gi number as id:
36 :     # $gi = CDS_gi( $ftr ) # gi number or nothing:
37 :     #
38 :     # $product = product( $ftr )
39 :     #
40 :     # @EC_number = EC_number( $ftr )
41 :     # \@EC_number = EC_number( $ftr )
42 :     #
43 :     # $translation = CDS_translation( $ftr ) # Uses in situ if found
44 :     #
45 :     # Convert GenBank location to [ [ $contig, $begin, $dir, $length ], ... ]
46 :     #
47 :     # \@cbdl = genbank_loc_2_cbdl( $loc, $contig_id )
48 :     #
49 :     # Convert GenBank location to a SEED or Sapling location string.
50 :     #
51 :     # $loc = genbank_loc_2_seed( $acc, $loc )
52 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_seed( $acc, $loc )
53 :     #
54 :     # $loc = genbank_loc_2_sapling( $acc, $loc )
55 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_sapling( $acc, $loc )
56 :     #
57 :     # Convert a [ [ contig, begin, dir, length ], ... ] location to GenBank.
58 :     #
59 :     # $gb_location = cbdl_2_genbank( \@cbdl )
60 :     # ( $contig, $gb_location ) = cbdl_2_genbank( \@cbdl )
61 :     #
62 :     #===============================================================================
63 :    
64 :     use strict;
65 :     # eval { require gjoseqlib; } # required for translations
66 :    
67 :     require Exporter;
68 :     our @ISA = qw(Exporter);
69 :     our @EXPORT = qw( parse_genbank );
70 :    
71 :     #===============================================================================
72 :     #
73 :     # @entries = parse_genbank( ) # \*STDIN
74 :     # \@entries = parse_genbank( ) # \*STDIN
75 :     # @entries = parse_genbank( \*FH )
76 :     # \@entries = parse_genbank( \*FH )
77 :     # @entries = parse_genbank( $file )
78 :     # \@entries = parse_genbank( $file )
79 :     #
80 :     # Each entry is a hash with key value pairs, some of which have special
81 :     # processing.
82 :     #
83 :     # Key => Value
84 :     #
85 :     # ACCESSION => [ Accession, ... ]
86 :     # DATE => Date
87 :     # DEFINITION => Definition
88 :     # GEOMETRY => Geometry
89 :     # GI => GI_number
90 :     # KEYWORDS => { Key_phrase => 1, Key_phrase => 1, ... }
91 :     # LOCUS => Locus
92 :     # ORGANISM => Organism_name
93 :     # REFERENCES => [ { Field => Value, Field => Value, ... }, ... ]
94 :     # SEQUENCE => Sequence
95 :     # SOURCE => Source_string
96 :     # TAXONOMY => [ Taxon, Taxon, ... ]
97 :     # VERSION => [ Version, Other_information ]
98 :     #
99 :     # Feature records are merged by type. Slash is removed from qualifier name.
100 :     # Surrounding quotation marks are stripped from qualifier values.
101 :     #
102 :     # FEATURES => { Type => [ [ Location, { Qualifier => \@values } ],
103 :     # [ Location, { Qualifier => \@values } ],
104 :     # ...
105 :     # ],
106 :     # Type => ...
107 :     # }
108 :     #
109 :     #===============================================================================
110 :     sub parse_genbank
111 :     {
112 :     my $file = shift;
113 :     my ( $fh, $close ) = input_filehandle( $file );
114 :     $fh or return wantarray() ? () : [];
115 :    
116 :     my @entries;
117 :    
118 :     my $state = 0;
119 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
120 :     # 0 = Looking for LOCUS
121 :     # 1 = Header information
122 :     # 2 = Features
123 :     # 3 = Sequence
124 :     # -1 = Error
125 :    
126 :     while ( $state >= 0 )
127 :     {
128 :     my %entry = ();
129 :     my @sequence;
130 :    
131 :     # 1 2 3 4 5 6 7 8
132 :     # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
133 :     # LOCUS NC_000909 1664970 bp DNA circular BCT 03-DEC-2005
134 :     # LOCUS DGRINCAD_6 9696 BP DS-DNA SYN 22-AUG-2006
135 :     #
136 :     while ( $state == 0 )
137 :     {
138 :     if ( /^LOCUS\s+(\S+)\s.*\s(\S+)\s+\S+\s+(\S+)\s*$/ )
139 :     {
140 :     $entry{ LOCUS } = $1;
141 :     $entry{ GEOMETRY } = $2;
142 :     $entry{ DATE } = $3;
143 :     $state = 1;
144 :     }
145 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
146 :     }
147 :    
148 :     # Reading the header requires merging continuations, then dealing
149 :     # with the data:
150 :    
151 :     while ( $state == 1 )
152 :     {
153 :     if ( /^FEATURES / )
154 :     {
155 :     $state = 2;
156 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
157 :     }
158 :    
159 :     elsif ( /^ORIGIN / )
160 :     {
161 :     $state = 3;
162 :     }
163 :    
164 :     elsif ( /^REFERENCE / )
165 :     {
166 :     my $ref;
167 :     ( $ref, $state, $_ ) = read_ref( $_, $fh );
168 :     push @{ $entry{ REFERENCES } }, $ref if $ref;
169 :     defined() or $state = -1;
170 :     }
171 :    
172 :     elsif ( /^(..........) (.*\S)\s*$/ ) # Any other keyword
173 :     {
174 :     my ( $tag, $value ) = ( $1, $2 );
175 :     $tag =~ s/^ +//;
176 :     $tag =~ s/ +$//;
177 :    
178 :     # Merge continuations:
179 :    
180 :     my $sep = ( $tag eq 'ORGANISM' ) ? "\t" :
181 :     ( $tag eq 'COMMENT' ) ? "\n" : ' ';
182 :    
183 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
184 :     while ( $state >= 0 && s/^ {12}/$sep/ )
185 :     {
186 :     $value .= $_ if length() > 1;
187 :     $sep = ' ' if $sep eq "\t";
188 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
189 :     }
190 :    
191 :     # Special case formats
192 :    
193 :     if ( $tag eq 'ACCESSION' )
194 :     {
195 :     $entry{ ACCESSION } = [ split / +/, $value ];
196 :     }
197 :     elsif ( $tag eq 'VERSION' )
198 :     {
199 :     my ( $gi ) = $value =~ s/ +GI:(\d+)//;
200 :     $entry{ GI } = $gi if $gi;
201 :     $entry{ VERSION } = [ split / +/, $value ];
202 :     }
203 :     elsif ( $tag eq 'KEYWORDS' )
204 :     {
205 :     $value =~ s/\s*\.$//;
206 :     $entry{ KEYWORDS } = { map { $_ => 1 } split /; */, $value } if $value;
207 :     }
208 :     elsif ( $tag eq 'ORGANISM' )
209 :     {
210 :     $value =~ s/\.$//;
211 :     my ( $org, $tax ) = split /\t/, $value;
212 :     $entry{ ORGANISM } = $org;
213 :     $entry{ TAXONOMY } = [ split /; */, $tax ] if $tax;
214 :     }
215 :     else
216 :     {
217 :     $entry{ $tag } = $value;
218 :     }
219 :    
220 :     # To know that we are at end of continuations, we must have
221 :     # read another line.
222 :    
223 :     defined() or $state = -1;
224 :     }
225 :    
226 :     else # This is really a format error, but let's skip it.
227 :     {
228 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
229 :     }
230 :     }
231 :    
232 :     # Reading the features requires merging continuations, then dealing
233 :     # with the data:
234 :    
235 :     while ( $state == 2 )
236 :     {
237 :     if ( /^ORIGIN/ )
238 :     {
239 :     $state = 3;
240 :     }
241 :    
242 :     elsif ( /^ (\S+)\s+(\S+)/ )
243 :     {
244 :     my ( $type, $loc ) = ( $1, $2 );
245 :     my ( $qualif, $value, %qualifs );
246 :    
247 :     # Collect the rest of the location:
248 :    
249 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
250 :     while ( $state >= 0 && /^ {15}/ && ( $_ !~ /^\s*\/\w/ ) )
251 :     {
252 :     s/^ +//;
253 :     $loc .= $_;
254 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
255 :     }
256 :    
257 :     # Collect qualiiers:
258 :    
259 :     while ( $state == 2 && ( $_ =~ /^\s*\/\w+/ ) )
260 :     {
261 :     # Qualifiers without = get an undef value (intentionally)
262 :    
263 :     ( $qualif, undef, $value ) = /^\s*\/(\w+)(=(.*))?/;
264 :    
265 :     # Quoted strings can have value lines that start with /, so
266 :     # we must track quotation marks.
267 :    
268 :     my $nquote = $value ? $value =~ tr/"// : 0;
269 :    
270 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
271 :     while ( $state >= 0 && /^ {15}/ && ( $nquote % 2 || ( ! /^\s*\/\w/ ) ) )
272 :     {
273 :     s/^ +//;
274 :     $nquote += tr/"//;
275 :     $value .= " $_";
276 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
277 :     }
278 :    
279 :     if ( $nquote % 2 )
280 :     {
281 :     print STDERR "Feature quotation nesting error: $type, $loc, $qualif=$value\n";
282 :     exit;
283 :     }
284 :    
285 :     if ( $qualif )
286 :     {
287 :     if ( $value && $value =~ /^".*"$/ )
288 :     {
289 :     $value =~ s/^"//;
290 :     $value =~ s/"$//;
291 :     $value =~ s/""/"/g;
292 :     }
293 :    
294 :     if ( $qualif eq 'translation' ) { $value =~ s/ +//g }
295 :    
296 :     push @{ $qualifs{ $qualif } }, $value;
297 :     }
298 :     }
299 :    
300 :     push @{ $entry{ FEATURES }->{ $type } }, [ $loc, \%qualifs ] if ( $type && $loc );
301 :    
302 :     defined() or $state = -1;
303 :     }
304 :    
305 :     elsif ( /^\s{0,4}\S/ ) # Not feature and not origin
306 :     {
307 :     $state = -1;
308 :     }
309 :    
310 :     else
311 :     {
312 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
313 :     }
314 :     }
315 :    
316 :     # Read the sequence:
317 :    
318 :     while ( $state == 3 )
319 :     {
320 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
321 :     while ( $state >= 0 && /^[ 0-9]{10}/ )
322 :     {
323 :     s/[^A-Za-z]+//g;
324 :     push @sequence, $_;
325 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
326 :     }
327 :     $entry{ SEQUENCE } = join "", @sequence;
328 :    
329 :     push @entries, \%entry;
330 :    
331 :     $state = ( $_ eq '//' ) ? 0 : -1 if $state >= 0;
332 :     }
333 :     }
334 :    
335 :     wantarray ? @entries : \@entries;
336 :     }
337 :    
338 :    
339 :     #-------------------------------------------------------------------------------
340 :     # Parse a reference.
341 :     #-------------------------------------------------------------------------------
342 :     sub read_ref
343 :     {
344 :     my ( $line, $fh ) = @_;
345 :     my $state = 1;
346 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
347 :    
348 :     my ( $tag, $value );
349 :     my %ref = ();
350 :     while ( ( $state >= 0 ) && /^ / )
351 :     {
352 :     if ( substr( $_, 0, 10 ) =~ /\S/ )
353 :     {
354 :     ( $tag, $value ) = $_ =~ /\s*(\w+)\s+(.*)$/;
355 :     }
356 :     elsif ( /\S/ )
357 :     {
358 :     s/^ +//;
359 :     $value .= " $_" if $_;
360 :     }
361 :    
362 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
363 :    
364 :     if ( ( $state < 0 ) || ( /^ {0,9}\S/ ) )
365 :     {
366 :     if ( $tag && $value )
367 :     {
368 :     $ref{ $tag } = $value;
369 :     }
370 :     $tag = $value = undef;
371 :     }
372 :     }
373 :    
374 :     ( (keys %ref ? \%ref : undef ), $state, $_ )
375 :     }
376 :    
377 :    
378 :    
379 :     #===============================================================================
380 :     # Access methods for some features and feature data
381 :     #===============================================================================
382 :     # Sequence of a feature, optionally including information on partial ends.
383 :     #
384 :     # $seq = ftr_dna( $dna, $ftr )
385 :     # $seq = ftr_dna( \$dna, $ftr )
386 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( $dna, $ftr )
387 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( \$dna, $ftr )
388 :     #
389 :     #-------------------------------------------------------------------------------
390 :     sub ftr_dna
391 :     {
392 :     my ( $dnaR, $ftr ) = @_;
393 :     $dnaR and ( ! ref( $dnaR ) or ref( $dnaR ) eq 'SCALAR' ) or return undef;
394 :     $dnaR = \$dnaR if ! ref( $dnaR );
395 :    
396 :     my $loc = &location( $ftr );
397 :     $loc or return undef;
398 :    
399 :     my $loc0 = $loc;
400 :     my $complement = ( $loc =~ s/^complement\((.*)\)$/$1/ );
401 :     $loc =~ s/^join\((.*)\)$/$1/;
402 :     my @spans = split /,/, $loc;
403 :     if ( grep { ! /^<?\d+\.\.>?\d+$/ } @spans )
404 :     {
405 :     print STDERR "*** Feature location parse error: $loc0\n";
406 :     return undef;
407 :     }
408 :    
409 :     my $partial_5 = $spans[ 0] =~ s/^<//;
410 :     my $partial_3 = $spans[-1] =~ s/\.\.>/../;
411 :     ( $partial_5, $partial_3 ) = ( $partial_3, $partial_5 ) if $complement;
412 :    
413 :     my $seq = join( '', map { extract_span( $dnaR, $_ ) } @spans );
414 :     $seq = gjoseqlib::complement_DNA_seq( $seq ) if $complement;
415 :    
416 :     # Sequences that run off the end can start at other than the first
417 :     # nucleotide of a codon.
418 :    
419 :     my $qual = &qualifiers( $ftr );
420 :     my $codon_start = $qual->{ codon_start } ? $qual->{ codon_start }->[0] : 1;
421 :     $seq = substr( $seq, $codon_start-1 ) if $codon_start > 1;
422 :    
423 :     wantarray ? ( $seq, $partial_5, $partial_3 ) : $seq;
424 :     }
425 :    
426 :     sub extract_span
427 :     {
428 :     my ( $dnaR, $span ) = @_;
429 :     my ( $beg, $end ) = $span =~ /^<?(\d+)\.\.>?(\d+)$/;
430 :     ( $beg > 0 ) && ( $beg <= $end ) && ( $end <= length( $$dnaR ) ) or return '';
431 :    
432 :     substr( $$dnaR, $beg-1, $end-$beg+1 );
433 :     }
434 :    
435 :    
436 :     #-------------------------------------------------------------------------------
437 :     #
438 :     # $ftr_location = location( $ftr ) # Returns empty string on failure.
439 :     #
440 :     #-------------------------------------------------------------------------------
441 :     sub location
442 :     {
443 :     my ( $ftr ) = @_;
444 :    
445 :     ( defined( $ftr )
446 :     && ( ref( $ftr ) eq 'ARRAY' )
447 :     && ( @$ftr > 1 ) )
448 :     ? $ftr->[0]
449 :     : '';
450 :     }
451 :    
452 :    
453 :     #-------------------------------------------------------------------------------
454 :     #
455 :     # \%ftr_qualifiers = qualifiers( $ftr ) # Returns empty hash reference on failure.
456 :     #
457 :     #-------------------------------------------------------------------------------
458 :     sub qualifiers
459 :     {
460 :     my ( $ftr ) = @_;
461 :     my $qual;
462 :     ( defined( $ftr )
463 :     && ( ref( $ftr ) eq 'ARRAY' )
464 :     && ( @$ftr > 1 )
465 :     && defined( $qual = $ftr->[1] )
466 :     && ( ref( $qual ) eq 'HASH' ) )
467 :     ? $qual
468 :     : {};
469 :     }
470 :    
471 :    
472 :     #-------------------------------------------------------------------------------
473 :     #
474 :     # $gene = gene( $ftr )
475 :     # @gene_and_synonyms = gene( $ftr )
476 :     #
477 :     #-------------------------------------------------------------------------------
478 :     sub gene
479 :     {
480 :     my $qual = &qualifiers( @_ );
481 :     my %seen;
482 :     my @gene = grep { ! $seen{ $_ }++ }
483 :     ( ( $qual->{ gene } ? @{ $qual->{ gene } } : () ),
484 :     ( $qual->{ gene_synonym } ? @{ $qual->{ gene_synonym } } : () )
485 :     );
486 :    
487 :     wantarray ? @gene : $gene[0];
488 :     }
489 :    
490 :    
491 :     #-------------------------------------------------------------------------------
492 :     # Prefer protein_id as id:
493 :     #
494 :     # $id = CDS_id( $ftr )
495 :     #
496 :     #-------------------------------------------------------------------------------
497 :     sub CDS_id
498 :     {
499 :     my $qual = &qualifiers( @_ );
500 :     my $id;
501 :    
502 :     ( $id ) = @{ $qual->{ protein_id } } if $qual->{ protein_id };
503 :     ( $id ) = map { m/^GI:(.+)$/i ? $1 : () } @{ $qual->{ db_xref } } if ! $id && $qual->{ db_xref };
504 :     ( $id ) = @{ $qual->{ locus_tag } } if ! $id && $qual->{ locus_tag };
505 :    
506 :     $id;
507 :     }
508 :    
509 :    
510 :     #-------------------------------------------------------------------------------
511 :     # Prefer gi number as id:
512 :     #
513 :     # $id = CDS_gi_or_id( $ftr )
514 :     #
515 :     #-------------------------------------------------------------------------------
516 :     sub CDS_gi_or_id
517 :     {
518 :     my $qual = &qualifiers( @_ );
519 :     my $id;
520 :    
521 :     ( $id ) = map { m/^GI:(.+)$/i ? $1 : () } @{ $qual->{ db_xref } } if $qual->{ db_xref };
522 :     ( $id ) = @{ $qual->{ protein_id } } if ! $id && $qual->{ protein_id };
523 :     ( $id ) = @{ $qual->{ locus_tag } } if ! $id && $qual->{ locus_tag };
524 :    
525 :     $id;
526 :     }
527 :    
528 :    
529 :     #-------------------------------------------------------------------------------
530 :     # gi number or nothing:
531 :     #
532 :     # $gi = CDS_gi( $ftr )
533 :     #
534 :     #-------------------------------------------------------------------------------
535 :     sub CDS_gi
536 :     {
537 :     my $qual = &qualifiers( @_ );
538 :    
539 :     my ( $id ) = map { m/^GI:(.+)$/i ? $1 : () } @{ $qual->{ db_xref } } if $qual->{ db_xref };
540 :    
541 :     $id;
542 :     }
543 :    
544 :    
545 :     #-------------------------------------------------------------------------------
546 :     #
547 :     # $product = product( $ftr )
548 :     #
549 :     #-------------------------------------------------------------------------------
550 :     sub product
551 :     {
552 :     my $qual = &qualifiers( @_ );
553 :     my $prod;
554 :    
555 :     ( $prod ) = @{ $qual->{ product } } if $qual->{ product };
556 :     ( $prod ) = @{ $qual->{ function } } if ! $prod && $qual->{ function };
557 :     ( $prod ) = @{ $qual->{ note } } if ! $prod && $qual->{ note };
558 :    
559 :     $prod;
560 :     }
561 :    
562 :    
563 :     #-------------------------------------------------------------------------------
564 :     #
565 :     # @EC_number = EC_number( $ftr )
566 :     # \@EC_number = EC_number( $ftr )
567 :     #
568 :     #-------------------------------------------------------------------------------
569 :     sub EC_number
570 :     {
571 :     my $qual = &qualifiers( @_ );
572 :     my @EC = $qual->{ EC_number } ? @{ $qual->{ EC_number } } : ();
573 :    
574 :     wantarray ? @EC : \@EC;
575 :     }
576 :    
577 :    
578 :     #-------------------------------------------------------------------------------
579 :     # This is the in situ translation. Will extract from the DNA sequence if
580 :     # supplied.
581 :     #
582 :     # $translation = CDS_translation( $ftr )
583 :     # $translation = CDS_translation( $ftr, $dna )
584 :     # $translation = CDS_translation( $ftr, \$dna )
585 :     #
586 :     #-------------------------------------------------------------------------------
587 :     sub CDS_translation
588 :     {
589 :     my ( $ftr, $seqR ) = @_;
590 :     my $qual = &qualifiers( $ftr );
591 :    
592 :     return $qual->{ translation }->[0] if $qual->{ translation };
593 :     return undef unless $seqR;
594 :    
595 :     my $CDS_dna = ftr_dna( $seqR, $ftr ) or return undef;
596 :    
597 :     my $have_lib = 0;
598 :     eval { require gjoseqlib; $have_lib = 1; };
599 :    
600 :     return
601 :     }
602 :    
603 :    
604 :     #===============================================================================
605 :     # Utilities for locations and location strings.
606 :     #===============================================================================
607 :     # Convert GenBank location to a SEED location.
608 :     #
609 :     # $loc = genbank_loc_2_seed( $acc, $loc )
610 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_seed( $acc, $loc )
611 :     #
612 :     #-------------------------------------------------------------------------------
613 :     sub genbank_loc_2_seed
614 :     {
615 :     my ( $acc, $loc ) = @_;
616 :     $acc && $loc or return undef;
617 :     genbank_loc_2_string( $acc, $loc, 'seed' );
618 :     }
619 :    
620 :    
621 :     #-------------------------------------------------------------------------------
622 :     # Convert GenBank location to a Sapling location.
623 :     #
624 :     # $loc = genbank_loc_2_sapling( $acc, $loc )
625 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_sapling( $acc, $loc )
626 :     #
627 :     #-------------------------------------------------------------------------------
628 :     sub genbank_loc_2_sapling
629 :     {
630 :     my ( $acc, $loc ) = @_;
631 :     $acc && $loc or return undef;
632 :     genbank_loc_2_string( $acc, $loc, 'sapling' );
633 :     }
634 :    
635 :    
636 :     #-------------------------------------------------------------------------------
637 :     # Convert GenBank location to another location format.
638 :     # At present, only 'sapling' (D) and 'seed' are supported.
639 :     #
640 :     # $loc = genbank_loc_2_string( $acc, $loc, $format )
641 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_string( $acc, $loc, $format )
642 :     #
643 :     #-------------------------------------------------------------------------------
644 :     sub genbank_loc_2_string
645 :     {
646 :     my ( $acc, $loc, $format ) = @_;
647 :     $acc && $loc or return undef;
648 :    
649 :     my ( $cbdl, $partial_5, $partial_3 ) = genbank_loc_2_cbdl( $loc, $acc );
650 :     my $str = cbdl_2_string( $cbdl, $format || 'sapling' );
651 :    
652 :     wantarray ? ( $str, $partial_5, $partial_3 ) : $str;
653 :     }
654 :    
655 :    
656 :     #-------------------------------------------------------------------------------
657 :     # Convert GenBank location to a list of [contig, begin, dir, len ] locations.
658 :     # order() is treated as join(). Nesting is allowed (unlike the standard).
659 :     #
660 :     # \@cbdl = genbank_loc_2_cbdl( $loc, $accession )
661 :     # ( \@cbdl, $partial_5, $partial_3 ) = genbank_loc_2_cbdl( $loc, $accession )
662 :     #
663 :     # Elements are:
664 :     #
665 :     # (accession:)?<?\d+..>?\d+ # range of sites
666 :     # (accession:)?<?\d+^>?\d+ # site between residues
667 :     # (accession:)?\d+ # single residue
668 :     # (accession:)?complement\(element\)
669 :     # (accession:)?join\(element,element,...\)
670 :     # (accession:)?order\(element,element,...\)
671 :     #
672 :     #-------------------------------------------------------------------------------
673 :     # Paterns used in the parsing. They are in each subroutine due to a very
674 :     # strange initialization issue.
675 :     #
676 :     # Because $paranthetical is self-referential, it must be declared before it is
677 :     # defined.
678 :     #
679 :     # my $paranthetical;
680 :     # $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
681 :     # my $contigid = qr/[^\s:(),]+/;
682 :     # my $complement = qr/(?:$contigid:)?complement$paranthetical/;
683 :     # my $complement_parts = qr/(?:($contigid):)?complement($paranthetical)/;
684 :     # my $join = qr/(?:$contigid:)?join$paranthetical/;
685 :     # my $join_parts = qr/(?:($contigid):)?join($paranthetical)/;
686 :     # my $order = qr/(?:$contigid:)?order$paranthetical/;
687 :     # my $order_parts = qr/(?:($contigid):)?order($paranthetical)/;
688 :     # my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
689 :     # my $range_parts = qr/(?:($contigid):)?(<?)(\d+)\.\.(>?)(\d+)/;
690 :     # my $site = qr/(?:$contigid:)?<?\d+^>?\d+/;
691 :     # my $site_parts = qr/(?:($contigid):)?(<?)(\d+)^(>?)(\d+)/;
692 :     # my $position = qr/(?:$contigid:)?\d+/;
693 :     # my $position_parts = qr/(?:($contigid):)?(\d+)/;
694 :     # my $element = qr/$range|$position|$complement|$join|$order/;
695 :     # my $elementlist = qr/$element(?:,$element)*/;
696 :     #
697 :     #-------------------------------------------------------------------------------
698 :    
699 :     sub genbank_loc_2_cbdl
700 :     {
701 :     my ( $loc, $acc ) = @_;
702 :    
703 :     my $contigid = qr/[^\s:(),]+/;
704 :    
705 :     my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
706 :     return gb_loc_range( $loc, $acc ) if $loc =~ /^$range$/;
707 :    
708 :     # This cannot by part of any other format except complement
709 :     my $site = qr/(?:$contigid:)?<?\d+\^>?\d+/;
710 :     return gb_loc_site( $loc, $acc ) if $loc =~ /^$site$/;
711 :    
712 :     my $position = qr/(?:$contigid:)?\d+/;
713 :     return gb_loc_position( $loc, $acc ) if $loc =~ /^$position$/;
714 :    
715 :     my $paranthetical;
716 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
717 :     my $complement = qr/(?:$contigid:)?complement$paranthetical/;
718 :     return gb_loc_complement( $loc, $acc ) if $loc =~ /^$complement$/;
719 :    
720 :     my $join = qr/(?:$contigid:)?join$paranthetical/;
721 :     return gb_loc_join( $loc, $acc ) if $loc =~ /^$join$/;
722 :    
723 :     # Treated as a join
724 :     my $order = qr/(?:$contigid:)?order$paranthetical/;
725 :     return gb_loc_order( $loc, $acc ) if $loc =~ /^$order$/;
726 :    
727 :     return ();
728 :     }
729 :    
730 :    
731 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
732 :     # A range of positions, with optional accession number prefix, optional less
733 :     # than first position (5' partial), begin, .., optional greater than end
734 :     # position (3' partial), and end position.
735 :     #
736 :     # (\S+:)?<?\d+\.\.>?\d+
737 :     #
738 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_range( $loc, $acc )
739 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
740 :     sub gb_loc_range
741 :     {
742 :     my ( $loc, $acc ) = @_;
743 :    
744 :     my $contigid = qr/[^\s:(),]+/;
745 :     my $range_parts = qr/(?:($contigid):)?(<?)(\d+)\.\.(>?)(\d+)/;
746 :     my ( $acc2, $p5, $beg, $p3, $end ) = $loc =~ /^$range_parts$/;
747 :     $beg && $end or return ();
748 :     $acc2 ||= $acc;
749 :    
750 :     # GenBank standard is always $beg <= $end. We will relax that.
751 :    
752 :     ( [ [ $acc2, $beg, (($end>=$beg)?'+':'-'), abs($end-$beg)+1 ] ], $p5, $p3 );
753 :     }
754 :    
755 :    
756 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
757 :     # A range of positions, with optional accession number prefix, optional less
758 :     # than first position (5' partial), begin, ^, optional greater than end
759 :     # position (3' partial), and end position.
760 :     #
761 :     # (\S+:)?<?\d+^>?\d+
762 :     #
763 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_site( $loc, $acc )
764 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
765 :     sub gb_loc_site
766 :     {
767 :     my ( $loc, $acc ) = @_;
768 :    
769 :     my $contigid = qr/[^\s:(),]+/;
770 :     my $site_parts = qr/(?:($contigid):)?(<?)(\d+)\^(>?)(\d+)/;
771 :     my ( $acc2, $p5, $beg, $p3, $end ) = $loc =~ /^$site_parts$/;
772 :     $beg && $end or return ();
773 :     $acc2 ||= $acc;
774 :    
775 :     # GenBank standard is always $beg <= $end. We will relax that.
776 :    
777 :     ( [ [ $acc2, $beg, (($end>=$beg)?'+':'-'), abs($end-$beg)+1 ] ], $p5, $p3 );
778 :     }
779 :    
780 :    
781 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
782 :     # A singe position, with optional accession number prefix.
783 :     #
784 :     # (\S+:)?\d+
785 :     #
786 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_position( $loc, $acc )
787 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
788 :     sub gb_loc_position
789 :     {
790 :     my ( $loc, $acc ) = @_;
791 :    
792 :     my $contigid = qr/[^\s:(),]+/;
793 :     my $position_parts = qr/(?:($contigid):)?(\d+)/;
794 :     my ( $acc2, $beg ) = $loc =~ /^$position_parts$/;
795 :     $beg or return ();
796 :    
797 :     ( [ [ $acc2 || $acc, $beg, '+', 1 ] ], '', '' );
798 :     }
799 :    
800 :    
801 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
802 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_complement( $loc, $acc )
803 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
804 :     sub gb_loc_complement
805 :     {
806 :     my ( $loc, $acc ) = @_;
807 :    
808 :     my $paranthetical;
809 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
810 :     my $contigid = qr/[^\s:(),]+/;
811 :     my $complement_parts = qr/(?:($contigid):)?complement($paranthetical)/;
812 :    
813 :     my ( $acc2, $loc2 ) = $loc =~ /^$complement_parts$/;
814 :     $loc2 && $loc2 =~ s/^\(// && $loc2 =~ s/\)$// or return ();
815 :     my ( $locs, $p5, $p3 ) = genbank_loc_2_cbdl( $loc2, $acc2 || $acc );
816 :     $locs && ref( $locs ) eq 'ARRAY' && @$locs or return ();
817 :    
818 :     ( [ map { complement_cbdl( @$_ ) } reverse @$locs ], $p3, $p5 );
819 :     }
820 :    
821 :    
822 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
823 :     #
824 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_join( $loc, $acc )
825 :     #
826 :     # There is no warning about partial sequences internal to list.
827 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
828 :     sub gb_loc_join
829 :     {
830 :     my ( $loc, $acc ) = @_;
831 :    
832 :     my $paranthetical;
833 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
834 :     my $contigid = qr/[^\s:(),]+/;
835 :     my $complement = qr/(?:$contigid:)?complement$paranthetical/;
836 :     my $join = qr/(?:$contigid:)?join$paranthetical/;
837 :     my $join_parts = qr/(?:($contigid):)?join($paranthetical)/;
838 :     my $order = qr/(?:$contigid:)?order$paranthetical/;
839 :     my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
840 :     my $position = qr/(?:$contigid:)?\d+/;
841 :     my $element = qr/$range|$position|$complement|$join|$order/;
842 :     my $elementlist = qr/$element(?:,$element)*/;
843 :    
844 :     my ( $acc2, $locs ) = $loc =~ /^$join_parts$/;
845 :     $locs && $locs =~ s/^\(// && $locs =~ s/\)$//
846 :     && $locs =~ /^$elementlist$/
847 :     or return ();
848 :     $acc2 ||= $acc;
849 :    
850 :     my @elements = map { [ genbank_loc_2_cbdl( $_, $acc2 ) ] }
851 :     $locs =~ m/($element)/g;
852 :     @elements or return ();
853 :    
854 :     ( [ map { @{ $_->[0] } } @elements ], $elements[0]->[1], $elements[-1]->[2] );
855 :     }
856 :    
857 :    
858 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
859 :     # Ordered list is treated as a join:
860 :     #
861 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_order( $loc, $acc )
862 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
863 :     sub gb_loc_order
864 :     {
865 :     my ( $loc, $acc ) = @_;
866 :    
867 :     my $paranthetical;
868 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
869 :     my $contigid = qr/[^\s:(),]+/;
870 :     my $complement = qr/(?:$contigid:)?complement$paranthetical/;
871 :     my $join = qr/(?:$contigid:)?join$paranthetical/;
872 :     my $order = qr/(?:$contigid:)?order$paranthetical/;
873 :     my $order_parts = qr/(?:($contigid):)?order($paranthetical)/;
874 :     my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
875 :     my $position = qr/(?:$contigid:)?\d+/;
876 :     my $element = qr/$range|$position|$complement|$join|$order/;
877 :     my $elementlist = qr/$element(?:,$element)*/;
878 :    
879 :     my ( $acc2, $locs ) = $loc =~ /^$order_parts$/;
880 :     $locs && $locs =~ s/^\(// && $locs =~ s/\)$//
881 :     && $locs =~ /^$elementlist$/
882 :     or return ();
883 :    
884 :     gb_loc_join( "join($locs)", $acc2 || $acc );
885 :     }
886 :    
887 :    
888 :     #-------------------------------------------------------------------------------
889 :     # $cbdl = complement_cbdl( $contig, $beg, $dir, $len )
890 :     # $cbdl = complement_cbdl( [ $contig, $beg, $dir, $len ] )
891 :     #-------------------------------------------------------------------------------
892 :     sub complement_cbdl
893 :     {
894 :     defined $_[0] or return ();
895 :     my ( $contig, $beg, $dir, $len ) = ref( $_[0] ) ? @{$_[0]} : @_;
896 :    
897 :     ( $dir =~ /^-/ ) ? [ $contig, $beg -= $len - 1, '+', $len ]
898 :     : [ $contig, $beg += $len - 1, '-', $len ];
899 :     }
900 :    
901 :    
902 :     #-------------------------------------------------------------------------------
903 :     # $loc = cbdl_2_string( \@cbdl, $format )
904 :     #-------------------------------------------------------------------------------
905 :     sub cbdl_2_string
906 :     {
907 :     my ( $cbdl, $format ) = @_;
908 :     $cbdl && ( ref( $cbdl ) eq 'ARRAY' ) or return undef;
909 :     $format = 'sapling' if ! defined $format;
910 :     return cbdl_2_genbank( $cbdl ) if $foramt =~ m/genbank/i;
911 :     join( ',', map { cbdl_part_2_string( $_, $format ) } @$cbdl );
912 :     }
913 :    
914 :    
915 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
916 :     # Support function for formatting one contiguous part of location.
917 :     #
918 :     # $loc_part = cbdl_part_2_string( $cbdl_part, $format )
919 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
920 :     sub cbdl_part_2_string
921 :     {
922 :     my ( $part, $format ) = @_;
923 :     $part && ( ref( $part ) eq 'ARRAY' ) && ( @$part == 4 ) or return ();
924 :     my ( $contig, $beg, $dir, $len ) = @$part;
925 :     $dir = $dir =~ /^-/ ? '-' : '+';
926 :    
927 :     if ( $format =~ m/seed/i )
928 :     {
929 :     my $n2 = ( $dir eq '+' ) ? $beg + $len - 1 : $beg - $len + 1;
930 :     return join( '_', $contig, $beg, $n2 );
931 :     }
932 :    
933 :     # Default is sapling:
934 :    
935 :     return $contig . '_' . $beg . $dir . $len;
936 :     }
937 :    
938 :     #-------------------------------------------------------------------------------
939 :     # Convert a [ [ contig, begin, dir, length ], ... ] location to GenBank.
940 :     #
941 :     # $gb_location = cbdl_2_genbank( \@cbdl )
942 :     # ( $contig, $gb_location ) = cbdl_2_genbank( \@cbdl )
943 :     #-------------------------------------------------------------------------------
944 :     sub cbdl_2_genbank
945 :     {
946 :     my ( $cbdl, $contig ) = @_;
947 :     $cbdl && ref( $cbdl ) eq 'ARRAY' && @$cbdl or return '';
948 :     my @cbdl = ref( $cbdl->[0] ) ? @$cbdl : ( $cbdl );
949 :     @cbdl or return '';
950 :    
951 :     my $dir = $cbdl[0]->[2];
952 :     @cbdl = map { complement_cbdl( $_ ) } reverse @cbdl if $dir =~ /^-/;
953 :    
954 :     $contig = $dbdl[0]->[0];
955 :     my @gb = map { cbdl_part_2_genbank( $_, $contig ) } @cbdl;
956 :    
957 :     my $gb = ( @gb > 1 ) ? 'join(' . join( ',', @gb ) . ')' : $gb[0];
958 :    
959 :     $gb = "complement($gb)" if $dir =~ /^-/;
960 :    
961 :     return wantarray ? ( $contig, $gb ) : $gb;
962 :     }
963 :    
964 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
965 :     # Support function for formatting one contiguous part of location.
966 :     #
967 :     # $loc_part = cbdl_part_2_genbank( $cbdl_part, $contig )
968 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
969 :     sub cbdl_part_2_genbank
970 :     {
971 :     my ( $part, $contig0 ) = @_;
972 :     $part && ( ref( $part ) eq 'ARRAY' ) && ( @$part == 4 ) or return ();
973 :     my ( $contig, $beg, $dir, $len ) = @$part;
974 :    
975 :     my $gb;
976 :     if ( $dir =~ /^-/ )
977 :     {
978 :     my $end = $beg - ( $len-1);
979 :     $gb = "complement($end..$beg)";
980 :     }
981 :     else
982 :     {
983 :     my $end = $beg + ($len-1);
984 :     $gb = "$beg..$end";
985 :     }
986 :    
987 :     $gb = "$contig:$gb" if $contig0 && $contig ne $contig0;
988 :    
989 :     return $gb;
990 :     }
991 :    
992 :     #===============================================================================
993 :     # Helper function for defining an input filehandle:
994 :     #
995 :     # filehandle is passed through
996 :     # string is taken as file name to be openend
997 :     # undef or "" defaults to STDOUT
998 :     #
999 :     # \*FH = input_filehandle( $file );
1000 :     # ( \*FH, $close ) = input_filehandle( $file );
1001 :     #
1002 :     #===============================================================================
1003 :     sub input_filehandle
1004 :     {
1005 :     my $file = shift;
1006 :    
1007 :     # Null string or undef
1008 :    
1009 :     if ( ! defined( $file ) || ( $file eq '' ) )
1010 :     {
1011 :     return wantarray ? ( \*STDIN, 0 ) : \*STDIN;
1012 :     }
1013 :    
1014 :     # FILEHANDLE?
1015 :    
1016 :     if ( ref( $file ) eq "GLOB" )
1017 :     {
1018 :     return wantarray ? ( $file, 0 ) : $file;
1019 :     }
1020 :    
1021 :     # File name
1022 :    
1023 :     if ( ! ref( $file ) )
1024 :     {
1025 :     -f $file or die "Could not find input file \"$file\"\n";
1026 :     my $fh;
1027 :     open( $fh, "<$file" ) || die "Could not open \"$file\" for input\n";
1028 :     return wantarray ? ( $fh, 1 ) : $fh;
1029 :     }
1030 :    
1031 :     return wantarray ? ( \*STDIN, undef ) : \*STDIN;
1032 :     }
1033 :    
1034 :    
1035 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3