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

Annotation of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : efrank 1.1 package gjoseqlib;
2 :    
3 : golsen 1.2 # A sequence entry is ( $id, $def, $seq )
4 :     # A list of entries is a list of references
5 :     #
6 : overbeek 1.4 # @seq_entry = read_next_fasta_seq( \*FILEHANDLE )
7 :     # @seq_entries = read_fasta_seqs( \*FILEHANDLE ) # Original form
8 :     # @seq_entries = read_fasta( ) # STDIN
9 :     # @seq_entries = read_fasta( \*FILEHANDLE )
10 :     # @seq_entries = read_fasta( $filename )
11 :     # @seq_entries = read_clustal( ) # STDIN
12 :     # @seq_entries = read_clustal( \*FILEHANDLE )
13 :     # @seq_entries = read_clustal( $filename )
14 :     # @seq_entries = read_clustal_file( $filename )
15 : golsen 1.2 #
16 : overbeek 1.4 # $seq_ind = index_seq_list( @seq_entries ); # hash from ids to entries
17 : golsen 1.2 # @seq_entry = seq_entry_by_id( \%seq_index, $seq_id );
18 :     # $seq_desc = seq_desc_by_id( \%seq_index, $seq_id );
19 :     # $seq = seq_data_by_id( \%seq_index, $seq_id );
20 :     #
21 :     # ( $id, $def ) = parse_fasta_title( $title )
22 :     # ( $id, $def ) = split_fasta_title( $title )
23 :     #
24 : overbeek 1.4 # print_seq_list_as_fasta( \*FILEHANDLE, @seq_entry_list ); # Original form
25 :     # print_alignment_as_fasta( @seq_entry_list ); # STDOUT
26 :     # print_alignment_as_fasta( \@seq_entry_list ); # STDOUT
27 :     # print_alignment_as_fasta( \*FILEHANDLE, @seq_entry_list );
28 :     # print_alignment_as_fasta( \*FILEHANDLE, \@seq_entry_list );
29 :     # print_alignment_as_fasta( $filename, @seq_entry_list );
30 :     # print_alignment_as_fasta( $filename, \@seq_entry_list );
31 :     # print_alignment_as_phylip( @seq_entry_list ); # STDOUT
32 :     # print_alignment_as_phylip( \@seq_entry_list ); # STDOUT
33 :     # print_alignment_as_phylip( \*FILEHANDLE, @seq_entry_list );
34 :     # print_alignment_as_phylip( \*FILEHANDLE, \@seq_entry_list );
35 :     # print_alignment_as_phylip( $filename, @seq_entry_list );
36 :     # print_alignment_as_phylip( $filename, \@seq_entry_list );
37 :     # print_alignment_as_nexus( [ \%label_hash, ] @seq_entry_list );
38 :     # print_alignment_as_nexus( [ \%label_hash, ] \@seq_entry_list );
39 :     # print_alignment_as_nexus( \*FILEHANDLE, [ \%label_hash, ] @seq_entry_list );
40 :     # print_alignment_as_nexus( \*FILEHANDLE, [ \%label_hash, ] \@seq_entry_list );
41 :     # print_alignment_as_nexus( $filename, [ \%label_hash, ] @seq_entry_list );
42 :     # print_alignment_as_nexus( $filename, [ \%label_hash, ] \@seq_entry_list );
43 :     # print_seq_as_fasta( \*FILEHANDLE, $id, $desc, $seq) ;
44 :     # print_seq_as_fasta( \*FILEHANDLE, @seq_entry );
45 :     # print_gb_locus( \*FILEHANDLE, $locus, $def, $accession, $seq );
46 : golsen 1.2 #
47 : golsen 1.6 # @packed_seqs = pack_alignment( @seqs )
48 :     # @packed_seqs = pack_alignment( \@seqs )
49 :     # \@packed_seqs = pack_alignment( @seqs )
50 :     # \@packed_seqs = pack_alignment( \@seqs )
51 : golsen 1.5 #
52 : golsen 1.2 # @entry = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );
53 :     # @entry = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );
54 : golsen 1.9 # $DNAseq = DNA_subseq( $seq, $from, $to );
55 :     # $DNAseq = DNA_subseq( \$seq, $from, $to );
56 :     # $RNAseq = RNA_subseq( $seq, $from, $to );
57 :     # $RNAseq = RNA_subseq( \$seq, $from, $to );
58 : golsen 1.2 # @entry = complement_DNA_entry( @seq_entry [, $fix_id] );
59 :     # @entry = complement_RNA_entry( @seq_entry [, $fix_id] );
60 :     # $DNAseq = complement_DNA_seq( $NA_seq );
61 :     # $RNAseq = complement_RNA_seq( $NA_seq );
62 :     # $DNAseq = to_DNA_seq( $NA_seq );
63 :     # $RNAseq = to_RNA_seq( $NA_seq );
64 :     # $seq = pack_seq( $sequence )
65 :     # $seq = clean_ae_sequence( $seq )
66 :     #
67 : golsen 1.10 # $aa = translate_seq( $nt, $met_start )
68 :     # $aa = translate_seq( $nt )
69 :     # $aa = translate_codon( $triplet );
70 :     #
71 :     # User-supplied genetic code. The supplied code needs to be complete in
72 :     # RNA and/or DNA, and upper and/or lower case. The program guesses based
73 :     # on lysine and phenylalanine codons.
74 : golsen 1.2 #
75 : golsen 1.10 # $aa = translate_seq_with_user_code( $nt, $gen_code_hash, $met_start )
76 :     # $aa = translate_seq_with_user_code( $nt, $gen_code_hash )
77 : golsen 1.9 #
78 : golsen 1.10 # Locations (= oriented intervals) are ( id, start, end )
79 :     # Intervals are ( id, left, right )
80 : golsen 1.2 #
81 : overbeek 1.4 # @intervals = read_intervals( \*FILEHANDLE )
82 :     # @intervals = read_oriented_intervals( \*FILEHANDLE )
83 : golsen 1.2 # @intervals = standardize_intervals( @interval_refs ) # (id, left, right)
84 :     # @joined = join_intervals( @interval_refs )
85 :     # @intervals = locations_2_intervals( @locations )
86 :     # $interval = locations_2_intervals( $location )
87 :     # @reversed = reverse_intervals( @interval_refs ) # (id, end, start)
88 : overbeek 1.4 #
89 :     # Convert GenBank locations to SEED locations
90 :     #
91 :     # @seed_locs = gb_location_2_seed( $contig, @gb_locs )
92 : golsen 1.7 #
93 :     # Read quality scores from a fasta-like file:
94 :     #
95 :     # @seq_entries = read_qual( ) # STDIN
96 :     # \@seq_entries = read_qual( ) # STDIN
97 :     # @seq_entries = read_qual( \*FILEHANDLE )
98 :     # \@seq_entries = read_qual( \*FILEHANDLE )
99 :     # @seq_entries = read_qual( $filename )
100 :     # \@seq_entries = read_qual( $filename )
101 :     #
102 : golsen 1.2
103 :     use strict;
104 : golsen 1.9 use Carp;
105 : efrank 1.1
106 : golsen 1.2 # Exported global variables:
107 : efrank 1.1
108 : golsen 1.2 our @aa_1_letter_order; # Alpha by 1 letter
109 :     our @aa_3_letter_order; # PAM matrix order
110 :     our @aa_n_codon_order;
111 :     our %genetic_code;
112 :     our %genetic_code_with_U;
113 :     our %amino_acid_codons_DNA;
114 :     our %amino_acid_codons_RNA;
115 :     our %n_codon_for_aa;
116 :     our %reverse_genetic_code_DNA;
117 :     our %reverse_genetic_code_RNA;
118 :     our %DNA_letter_can_be;
119 :     our %RNA_letter_can_be;
120 :     our %one_letter_to_three_letter_aa;
121 :     our %three_letter_to_one_letter_aa;
122 : efrank 1.1
123 :     require Exporter;
124 :    
125 :     our @ISA = qw(Exporter);
126 :     our @EXPORT = qw(
127 :     read_fasta_seqs
128 : overbeek 1.4 read_fasta
129 : efrank 1.1 read_next_fasta_seq
130 : overbeek 1.4 read_clustal_file
131 :     read_clustal
132 : efrank 1.1 parse_fasta_title
133 :     split_fasta_title
134 :     print_seq_list_as_fasta
135 : overbeek 1.4 print_alignment_as_fasta
136 :     print_alignment_as_phylip
137 :     print_alignment_as_nexus
138 : efrank 1.1 print_seq_as_fasta
139 :     print_gb_locus
140 :    
141 :     index_seq_list
142 :     seq_entry_by_id
143 :     seq_desc_by_id
144 :     seq_data_by_id
145 :    
146 : golsen 1.5 pack_alignment
147 :    
148 : efrank 1.1 subseq_DNA_entry
149 :     subseq_RNA_entry
150 : golsen 1.9 DNA_subseq
151 :     RNA_subseq
152 : efrank 1.1 complement_DNA_entry
153 :     complement_RNA_entry
154 :     complement_DNA_seq
155 :     complement_RNA_seq
156 :     to_DNA_seq
157 :     to_RNA_seq
158 : golsen 1.2 pack_seq
159 : efrank 1.1 clean_ae_sequence
160 :    
161 :     translate_seq
162 :     translate_codon
163 :     translate_seq_with_user_code
164 :    
165 :     read_intervals
166 : golsen 1.2 standardize_intervals
167 : efrank 1.1 join_intervals
168 : golsen 1.2 locations_2_intervals
169 :     read_oriented_intervals
170 :     reverse_intervals
171 : overbeek 1.4
172 :     gb_location_2_seed
173 : golsen 1.7
174 :     read_qual
175 : efrank 1.1 );
176 :    
177 : golsen 1.2 our @EXPORT_OK = qw(
178 :     @aa_1_letter_order
179 :     @aa_3_letter_order
180 :     @aa_n_codon_order
181 :     %genetic_code
182 :     %genetic_code_with_U
183 :     %amino_acid_codons_DNA
184 :     %amino_acid_codons_RNA
185 :     %n_codon_for_aa
186 :     %reverse_genetic_code_DNA
187 :     %reverse_genetic_code_RNA
188 :     %DNA_letter_can_be
189 :     %RNA_letter_can_be
190 :     %one_letter_to_three_letter_aa
191 :     %three_letter_to_one_letter_aa
192 :     );
193 : efrank 1.1
194 :    
195 :     #-----------------------------------------------------------------------------
196 : overbeek 1.4 # Helper function for defining an input filehandle:
197 :     # filehandle is passed through
198 :     # string is taken as file name to be openend
199 :     # undef or "" defaults to STDOUT
200 :     #
201 :     # ( \*FH, $name, $close [, $file] ) = input_filehandle( $file );
202 :     #
203 :     #-----------------------------------------------------------------------------
204 :     sub input_filehandle
205 :     {
206 :     my $file = shift;
207 :    
208 :     # FILEHANDLE
209 :    
210 :     return ( $file, $file, 0 ) if ( ref( $file ) eq "GLOB" );
211 :    
212 :     # Null string or undef
213 :    
214 :     return ( \*STDIN, "", 0 ) if ( ! defined( $file ) || ( $file eq "" ) );
215 :    
216 :     # File name
217 :    
218 :     if ( ! ref( $file ) )
219 :     {
220 :     my $fh;
221 :     -f $file or die "Could not find input file \"$file\"\n";
222 :     open( $fh, "<$file" ) || die "Could not open \"$file\" for input\n";
223 :     return ( $fh, $file, 1 );
224 :     }
225 :    
226 :     # Some other kind of reference; return the unused value
227 :    
228 :     return ( \*STDIN, undef, 0, $file );
229 :     }
230 :    
231 :    
232 :     #-----------------------------------------------------------------------------
233 :     # Read fasta sequences from a filehandle (legacy interface; use read_fasta)
234 : efrank 1.1 # Save the contents in a list of refs to arrays: (id, description, seq)
235 :     #
236 : overbeek 1.4 # @seq_entries = read_fasta_seqs( \*FILEHANDLE )
237 :     #-----------------------------------------------------------------------------
238 :     sub read_fasta_seqs { read_fasta( @_ ) }
239 :    
240 :    
241 : efrank 1.1 #-----------------------------------------------------------------------------
242 : overbeek 1.4 # Read fasta sequences.
243 :     # Save the contents in a list of refs to arrays: (id, description, seq)
244 :     #
245 :     # @seq_entries = read_fasta( ) # STDIN
246 :     # \@seq_entries = read_fasta( ) # STDIN
247 :     # @seq_entries = read_fasta( \*FILEHANDLE )
248 :     # \@seq_entries = read_fasta( \*FILEHANDLE )
249 :     # @seq_entries = read_fasta( $filename )
250 :     # \@seq_entries = read_fasta( $filename )
251 :     #-----------------------------------------------------------------------------
252 :     sub read_fasta {
253 :     my ( $fh, $name, $close, $unused ) = input_filehandle( $_[0] );
254 :     $unused && die "Bad reference type (" . ref( $unused ) . ") passed to read_fasta\n";
255 : efrank 1.1
256 :     my @seqs = ();
257 :     my ($id, $desc, $seq) = ("", "", "");
258 :    
259 : overbeek 1.4 while ( <$fh> ) {
260 : efrank 1.1 chomp;
261 :     if (/^>\s*(\S+)(\s+(.*))?$/) { # new id
262 :     if ($id && $seq) { push @seqs, [ $id, $desc, $seq ] }
263 :     ($id, $desc, $seq) = ($1, $3 ? $3 : "", "");
264 :     }
265 :     else {
266 : golsen 1.2 tr/ 0-9//d;
267 : efrank 1.1 $seq .= $_ ;
268 :     }
269 :     }
270 : overbeek 1.4 close( $fh ) if $close;
271 : efrank 1.1
272 : overbeek 1.4 if ( $id && $seq ) { push @seqs, [ $id, $desc, $seq ] }
273 :     return wantarray ? @seqs : \@seqs;
274 : efrank 1.1 }
275 :    
276 :    
277 :     #-----------------------------------------------------------------------------
278 :     # Read one fasta sequence at a time from a file.
279 :     # Return the contents as an array: (id, description, seq)
280 :     #
281 : overbeek 1.4 # @seq_entry = read_next_fasta_seq( \*FILEHANDLE )
282 : efrank 1.1 #-----------------------------------------------------------------------------
283 :     # Reading always overshoots, so save next id and description
284 :    
285 : golsen 1.2 { # Use bare block to scope the header hash
286 :    
287 :     my %next_header;
288 :    
289 :     sub read_next_fasta_seq {
290 :     my $fh = shift;
291 :     my ( $id, $desc );
292 : efrank 1.1
293 : golsen 1.2 if ( defined( $next_header{$fh} ) ) {
294 :     ( $id, $desc ) = parse_fasta_title( $next_header{$fh} );
295 : efrank 1.1 }
296 :     else {
297 : golsen 1.2 $next_header{$fh} = "";
298 :     ( $id, $desc ) = ( undef, "" );
299 :     }
300 :     my $seq = "";
301 :    
302 :     while ( <$fh> ) {
303 :     chomp;
304 :     if ( /^>/ ) { # new id
305 :     $next_header{$fh} = $_;
306 : overbeek 1.4 if ( defined($id) && $seq )
307 :     {
308 :     return wantarray ? ($id, $desc, $seq) : [$id, $desc, $seq]
309 :     }
310 : golsen 1.2 ( $id, $desc ) = parse_fasta_title( $next_header{$fh} );
311 :     $seq = "";
312 :     }
313 :     else {
314 :     tr/ 0-9//d;
315 :     $seq .= $_ ;
316 :     }
317 : efrank 1.1 }
318 :    
319 : golsen 1.2 # Done with file, delete "next header"
320 : efrank 1.1
321 : golsen 1.2 delete $next_header{$fh};
322 : overbeek 1.4 return ( defined($id) && $seq ) ? ( wantarray ? ($id, $desc, $seq)
323 :     : [$id, $desc, $seq]
324 :     )
325 :     : () ;
326 : golsen 1.2 }
327 : efrank 1.1 }
328 :    
329 :    
330 :     #-----------------------------------------------------------------------------
331 : overbeek 1.4 # Read a clustal alignment from a file.
332 :     # Save the contents in a list of refs to arrays: (id, description, seq)
333 :     #
334 :     # @seq_entries = read_clustal_file( $filename )
335 :     #-----------------------------------------------------------------------------
336 :     sub read_clustal_file { read_clustal( @_ ) }
337 :    
338 :    
339 :     #-----------------------------------------------------------------------------
340 :     # Read a clustal alignment.
341 :     # Save the contents in a list of refs to arrays: (id, description, seq)
342 :     #
343 :     # @seq_entries = read_clustal( ) # STDIN
344 :     # @seq_entries = read_clustal( \*FILEHANDLE )
345 :     # @seq_entries = read_clustal( $filename )
346 :     #-----------------------------------------------------------------------------
347 :     sub read_clustal {
348 :     my ( $fh, undef, $close, $unused ) = input_filehandle( shift );
349 :     $unused && die "Bad reference type (" . ref( $unused ) . ") passed to read_clustal_file\n";
350 :    
351 :     my ( %seq, @ids, $line );
352 :     while ( defined( $line = <$fh> ) )
353 :     {
354 :     ( $line =~ /^[A-Za-z0-9]/ ) or next;
355 :     chomp $line;
356 :     my @flds = split /\s+/, $line;
357 :     if ( @flds == 2 )
358 :     {
359 :     $seq{ $flds[0] } or push @ids, $flds[0];
360 :     push @{ $seq{ $flds[0] } }, $flds[1];
361 :     }
362 :     }
363 :     close( $fh ) if $close;
364 :    
365 :     map { [ $_, "", join( "", @{$seq{$_}} ) ] } @ids;
366 :     }
367 :    
368 :    
369 :     #-----------------------------------------------------------------------------
370 : efrank 1.1 # Parse a fasta file header to id and definition parts
371 :     #
372 :     # ($id, $def) = parse_fasta_title( $title )
373 :     # ($id, $def) = split_fasta_title( $title )
374 :     #-----------------------------------------------------------------------------
375 :     sub parse_fasta_title {
376 :     my $title = shift;
377 :     chomp;
378 :     if ($title =~ /^>?\s*(\S+)(:?\s+(.*\S)\s*)?$/) {
379 :     return ($1, $3 ? $3 : "");
380 :     }
381 : golsen 1.2 elsif ($title =~ /^>/) {
382 :     return ("", "");
383 :     }
384 : efrank 1.1 else {
385 : golsen 1.2 return (undef, "");
386 : efrank 1.1 }
387 :     }
388 :    
389 :     sub split_fasta_title {
390 :     parse_fasta_title ( shift );
391 :     }
392 :    
393 :    
394 :     #-----------------------------------------------------------------------------
395 : overbeek 1.4 # Helper function for defining an output filehandle:
396 :     # filehandle is passed through
397 :     # string is taken as file name to be openend
398 :     # undef or "" defaults to STDOUT
399 :     #
400 :     # ( \*FH, $name, $close [, $file] ) = output_filehandle( $file );
401 : efrank 1.1 #
402 :     #-----------------------------------------------------------------------------
403 : overbeek 1.4 sub output_filehandle
404 :     {
405 :     my $file = shift;
406 :    
407 :     # FILEHANDLE
408 :    
409 :     return ( $file, $file, 0 ) if ( ref( $file ) eq "GLOB" );
410 :    
411 :     # Null string or undef
412 :    
413 :     return ( \*STDOUT, "", 0 ) if ( ! defined( $file ) || ( $file eq "" ) );
414 :    
415 :     # File name
416 :    
417 :     if ( ! ref( $file ) )
418 :     {
419 :     my $fh;
420 :     open( $fh, ">$file" ) || die "Could not open output $file\n";
421 :     return ( $fh, $file, 1 );
422 :     }
423 :    
424 :     # Some other kind of reference; return the unused value
425 :    
426 :     return ( \*STDOUT, undef, 0, $file );
427 :     }
428 :    
429 :    
430 :     #-----------------------------------------------------------------------------
431 :     # Legacy function for printing fasta sequence set:
432 :     #
433 :     # print_seq_list_as_fasta( \*FILEHANDLE, @seq_entry_list );
434 :     #-----------------------------------------------------------------------------
435 :     sub print_seq_list_as_fasta { print_alignment_as_fasta( @_ ) }
436 :    
437 :    
438 :     #-----------------------------------------------------------------------------
439 :     # Print list of sequence entries in fasta format.
440 :     # Missing, undef or "" filename defaults to STDOUT.
441 :     #
442 :     # print_alignment_as_fasta( @seq_entry_list ); # STDOUT
443 :     # print_alignment_as_fasta( \@seq_entry_list ); # STDOUT
444 :     # print_alignment_as_fasta( \*FILEHANDLE, @seq_entry_list );
445 :     # print_alignment_as_fasta( \*FILEHANDLE, \@seq_entry_list );
446 :     # print_alignment_as_fasta( $filename, @seq_entry_list );
447 :     # print_alignment_as_fasta( $filename, \@seq_entry_list );
448 :     #-----------------------------------------------------------------------------
449 :     sub print_alignment_as_fasta {
450 :     my ( $fh, undef, $close, $unused ) = output_filehandle( shift );
451 :     ( unshift @_, $unused ) if $unused;
452 :    
453 : overbeek 1.8 ( ref( $_[0] ) eq "ARRAY" ) or confess "Bad sequence entry passed to print_alignment_as_fasta\n";
454 : overbeek 1.4
455 :     # Expand the sequence entry list if necessary:
456 :    
457 :     if ( ref( $_[0]->[0] ) eq "ARRAY" ) { @_ = @{ $_[0] } }
458 :    
459 :     foreach my $seq_ptr ( @_ ) { print_seq_as_fasta( $fh, @$seq_ptr ) }
460 :    
461 :     close( $fh ) if $close;
462 :     }
463 :    
464 :    
465 :     #-----------------------------------------------------------------------------
466 :     # Print list of sequence entries in phylip format.
467 :     # Missing, undef or "" filename defaults to STDOUT.
468 :     #
469 :     # print_alignment_as_phylip( @seq_entry_list ); # STDOUT
470 :     # print_alignment_as_phylip( \@seq_entry_list ); # STDOUT
471 :     # print_alignment_as_phylip( \*FILEHANDLE, @seq_entry_list );
472 :     # print_alignment_as_phylip( \*FILEHANDLE, \@seq_entry_list );
473 :     # print_alignment_as_phylip( $filename, @seq_entry_list );
474 :     # print_alignment_as_phylip( $filename, \@seq_entry_list );
475 :     #-----------------------------------------------------------------------------
476 :     sub print_alignment_as_phylip {
477 :     my ( $fh, undef, $close, $unused ) = output_filehandle( shift );
478 :     ( unshift @_, $unused ) if $unused;
479 :    
480 :     ( ref( $_[0] ) eq "ARRAY" ) or die die "Bad sequence entry passed to print_alignment_as_phylip\n";
481 :    
482 :     my @seq_list = ( ref( $_[0]->[0] ) eq "ARRAY" ) ? @{ $_[0] } : @_;
483 :    
484 :     my ( %id2, %used );
485 :     my $maxlen = 0;
486 :     foreach ( @seq_list )
487 :     {
488 :     my ( $id, undef, $seq ) = @$_;
489 :    
490 :     # Need a name that is unique within 10 characters
491 :    
492 :     my $id2 = substr( $id, 0, 10 );
493 :     $id2 =~ s/_/ /g; # PHYLIP sequence files accept spaces
494 :     my $n = "0";
495 :     while ( $used{ $id2 } )
496 :     {
497 :     $n++;
498 :     $id2 = substr( $id, 0, 10 - length( $n ) ) . $n;
499 :     }
500 :     $used{ $id2 } = 1;
501 :     $id2{ $id } = $id2;
502 :    
503 :     # Prepare to pad sequences (should not be necessary, but ...)
504 :    
505 :     my $len = length( $seq );
506 :     $maxlen = $len if ( $len > $maxlen );
507 :     }
508 : efrank 1.1
509 : overbeek 1.4 my $nseq = @seq_list;
510 :     print $fh "$nseq $maxlen\n";
511 :     foreach ( @seq_list )
512 :     {
513 :     my ( $id, undef, $seq ) = @$_;
514 :     my $len = length( $seq );
515 :     printf $fh "%-10s %s%s\n", $id2{ $id },
516 :     $seq,
517 :     $len<$maxlen ? ("?" x ($maxlen-$len)) : "";
518 : efrank 1.1 }
519 : overbeek 1.4
520 :     close( $fh ) if $close;
521 :     }
522 :    
523 :    
524 :     #-----------------------------------------------------------------------------
525 :     # Print list of sequence entries in nexus format.
526 :     # Missing, undef or "" filename defaults to STDOUT.
527 :     #
528 :     # print_alignment_as_nexus( [ \%label_hash, ] @seq_entry_list );
529 :     # print_alignment_as_nexus( [ \%label_hash, ] \@seq_entry_list );
530 :     # print_alignment_as_nexus( \*FILEHANDLE, [ \%label_hash, ] @seq_entry_list );
531 :     # print_alignment_as_nexus( \*FILEHANDLE, [ \%label_hash, ] \@seq_entry_list );
532 :     # print_alignment_as_nexus( $filename, [ \%label_hash, ] @seq_entry_list );
533 :     # print_alignment_as_nexus( $filename, [ \%label_hash, ] \@seq_entry_list );
534 :     #-----------------------------------------------------------------------------
535 :     sub print_alignment_as_nexus {
536 :     my ( $fh, undef, $close, $unused ) = output_filehandle( shift );
537 :     ( unshift @_, $unused ) if $unused;
538 :    
539 :     my $lbls = ( ref( $_[0] ) eq "HASH" ) ? shift : undef;
540 :    
541 :     ( ref( $_[0] ) eq "ARRAY" ) or die "Bad sequence entry passed to print_alignment_as_nexus\n";
542 :    
543 :     my @seq_list = ( ref( $_[0]->[0] ) eq "ARRAY" ) ? @{ $_[0] } : @_;
544 :    
545 :     my %id2;
546 :     my ( $maxidlen, $maxseqlen ) = ( 0, 0 );
547 :     my ( $n1, $n2, $nt, $nu ) = ( 0, 0, 0, 0 );
548 :     foreach ( @seq_list )
549 :     {
550 :     my ( $id, undef, $seq ) = @$_;
551 :     my $id2 = $lbls ? ( $lbls->{ $id } || $id ) : $id;
552 :     if ( $id2 !~ /^[-+.0-9A-Za-z~_|]+$/ )
553 :     {
554 :     $id2 =~ s/'/''/g;
555 :     $id2 = qq('$id2');
556 :     }
557 :     $id2{ $id } = $id2;
558 :     my $idlen = length( $id2 );
559 :     $maxidlen = $idlen if ( $idlen > $maxidlen );
560 :    
561 :     my $seqlen = length( $seq );
562 :     $maxseqlen = $seqlen if ( $seqlen > $maxseqlen );
563 :    
564 :     $nt += $seq =~ tr/Tt//d;
565 :     $nu += $seq =~ tr/Uu//d;
566 :     $n1 += $seq =~ tr/ACGNacgn//d;
567 :     $n2 += $seq =~ tr/A-Za-z//d;
568 :     }
569 :    
570 :     my $nseq = @seq_list;
571 :     my $type = ( $n1 < 2 * $n2 ) ? 'protein' : ($nt>$nu) ? 'DNA' : 'RNA';
572 :    
573 :     print $fh <<"END_HEAD";
574 :     #NEXUS
575 :    
576 :     BEGIN Data;
577 :     Dimensions
578 :     NTax=$nseq
579 :     NChar=$maxseqlen
580 :     ;
581 :     Format
582 :     DataType=$type
583 :     Gap=-
584 :     Missing=?
585 :     ;
586 :     Matrix
587 :    
588 :     END_HEAD
589 :    
590 :     foreach ( @seq_list )
591 :     {
592 :     my ( $id, undef, $seq ) = @$_;
593 :     my $len = length( $seq );
594 :     printf $fh "%-${maxidlen}s %s%s\n",
595 :     $id2{ $id },
596 :     $seq,
597 :     $len<$maxseqlen ? ("?" x ($maxseqlen-$len)) : "";
598 :     }
599 :    
600 :     print $fh <<"END_TAIL";
601 :     ;
602 :     END;
603 :     END_TAIL
604 :    
605 :     close( $fh ) if $close;
606 : efrank 1.1 }
607 :    
608 :    
609 :     #-----------------------------------------------------------------------------
610 :     # Print one sequence in fasta format to an open file
611 :     #
612 : overbeek 1.4 # print_seq_as_fasta( \*FILEHANDLE, $id, $desc, $seq );
613 :     # print_seq_as_fasta( \*FILEHANDLE, @seq_entry );
614 : efrank 1.1 #-----------------------------------------------------------------------------
615 :     sub print_seq_as_fasta {
616 :     my $fh = shift;
617 :     my ($id, $desc, $seq) = @_;
618 :    
619 :     printf $fh ($desc) ? ">$id $desc\n" : ">$id\n";
620 :     my $len = length($seq);
621 :     for (my $i = 0; $i < $len; $i += 60) {
622 :     print $fh substr($seq, $i, 60) . "\n";
623 :     }
624 :     }
625 :    
626 :    
627 :     #-----------------------------------------------------------------------------
628 :     # Print one sequence in GenBank flat file format:
629 :     #
630 : overbeek 1.4 # print_gb_locus( \*FILEHANDLE, $locus, $def, $accession, $seq )
631 : efrank 1.1 #-----------------------------------------------------------------------------
632 :     sub print_gb_locus {
633 :     my ($fh, $loc, $def, $acc, $seq) = @_;
634 :     my ($len, $i, $imax);
635 :     my $istep = 10;
636 :    
637 :     $len = length($seq);
638 :     printf $fh "LOCUS %-10s%7d bp\n", substr($loc,0,10), $len;
639 :     print $fh "DEFINITION " . substr(wrap_text($def,80,12), 12) . "\n";
640 :     if ($acc) { print $fh "ACCESSION $acc\n" }
641 :     print $fh "ORIGIN\n";
642 :    
643 :     for ($i = 1; $i <= $len; ) {
644 :     printf $fh "%9d", $i;
645 :     $imax = $i + 59; if ($imax > $len) { $imax = $len }
646 :     for ( ; $i <= $imax; $i += $istep) {
647 :     print $fh " " . substr($seq, $i-1, $istep);
648 :     }
649 :     print $fh "\n";
650 :     }
651 :     print $fh "//\n";
652 :     }
653 :    
654 :    
655 :     #-----------------------------------------------------------------------------
656 : golsen 1.7 # Return a string with text wrapped to defined line lengths:
657 :     #
658 :     # $wrapped_text = wrap_text( $str ) # default len = 80
659 :     # $wrapped_text = wrap_text( $str, $len ) # default ind = 0
660 :     # $wrapped_text = wrap_text( $str, $len, $indent ) # default ind_n = ind
661 :     # $wrapped_text = wrap_text( $str, $len, $indent_1, $indent_n )
662 :     #-----------------------------------------------------------------------------
663 :     sub wrap_text {
664 :     my ($str, $len, $ind, $indn) = @_;
665 :    
666 :     defined($str) || die "wrap_text called without a string\n";
667 :     defined($len) || ($len = 80);
668 :     defined($ind) || ($ind = 0);
669 :     ($ind < $len) || die "wrap error: indent greater than line length\n";
670 :     defined($indn) || ($indn = $ind);
671 :     ($indn < $len) || die "wrap error: indent_n greater than line length\n";
672 :    
673 :     $str =~ s/\s+$//;
674 :     $str =~ s/^\s+//;
675 :     my ($maxchr, $maxchr1);
676 :     my (@lines) = ();
677 :    
678 :     while ($str) {
679 :     $maxchr1 = ($maxchr = $len - $ind) - 1;
680 :     if ($maxchr >= length($str)) {
681 :     push @lines, (" " x $ind) . $str;
682 :     last;
683 :     }
684 :     elsif ($str =~ /^(.{0,$maxchr1}\S)\s+(\S.*)$/) { # no expr in {}
685 :     push @lines, (" " x $ind) . $1;
686 :     $str = $2;
687 :     }
688 :     elsif ($str =~ /^(.{0,$maxchr1}-)(.*)$/) {
689 :     push @lines, (" " x $ind) . $1;
690 :     $str = $2;
691 :     }
692 :     else {
693 :     push @lines, (" " x $ind) . substr($str, 0, $maxchr);
694 :     $str = substr($str, $maxchr);
695 :     }
696 :     $ind = $indn;
697 :     }
698 :    
699 :     return join("\n", @lines);
700 :     }
701 :    
702 :    
703 :     #-----------------------------------------------------------------------------
704 : efrank 1.1 # Build an index from seq_id to pointer to sequence entry: (id, desc, seq)
705 :     #
706 : overbeek 1.4 # my \%seq_ind = index_seq_list( @seq_list );
707 :     # my \%seq_ind = index_seq_list( \@seq_list );
708 : efrank 1.1 #
709 :     # Usage example:
710 :     #
711 : overbeek 1.4 # my @seq_list = read_fasta_seqs(\*STDIN); # list of pointers to entries
712 :     # my \%seq_ind = index_seq_list(@seq_list); # hash from names to pointers
713 :     # my @chosen_seq = @{%seq_ind{"contig1"}}; # extract one entry
714 : efrank 1.1 #
715 :     #-----------------------------------------------------------------------------
716 :     sub index_seq_list {
717 : overbeek 1.4 ( ref( $_[0] ) ne 'ARRAY' ) ? {}
718 :     : ( ref( $_[0]->[0] ) ne 'ARRAY' ) ? { map { $_->[0] => $_ } @_ }
719 :     : { map { $_->[0] => $_ } @{ $_[0] } }
720 : efrank 1.1 }
721 :    
722 :    
723 :     #-----------------------------------------------------------------------------
724 :     # Three routines to access all or part of sequence entry by id:
725 :     #
726 : overbeek 1.4 # @seq_entry = seq_entry_by_id( \%seq_index, $seq_id );
727 :     # \@seq_entry = seq_entry_by_id( \%seq_index, $seq_id );
728 :     # $seq_desc = seq_desc_by_id( \%seq_index, $seq_id );
729 :     # $seq = seq_data_by_id( \%seq_index, $seq_id );
730 : efrank 1.1 #
731 :     #-----------------------------------------------------------------------------
732 :     sub seq_entry_by_id {
733 :     (my $ind_ref = shift) || die "No index supplied to seq_entry_by_id\n";
734 :     (my $id = shift) || die "No id supplied to seq_entry_by_id\n";
735 : overbeek 1.4 return wantarray ? @{ $ind_ref->{$id} } : $ind_ref->{$id};
736 : efrank 1.1 }
737 :    
738 :    
739 :     sub seq_desc_by_id {
740 :     (my $ind_ref = shift) || die "No index supplied to seq_desc_by_id\n";
741 :     (my $id = shift) || die "No id supplied to seq_desc_by_id\n";
742 :     return ${ $ind_ref->{$id} }[1];
743 :     }
744 :    
745 :    
746 :     sub seq_data_by_id {
747 :     (my $ind_ref = shift) || die "No index supplied to seq_data_by_id\n";
748 :     (my $id = shift) || die "No id supplied to seq_data_by_id\n";
749 :     return ${ $ind_ref->{$id} }[2];
750 :     }
751 :    
752 : golsen 1.5 #-----------------------------------------------------------------------------
753 :     # Remove columns of alignment gaps from sequences:
754 :     #
755 : golsen 1.6 # @packed_seqs = pack_alignment( @seqs )
756 :     # @packed_seqs = pack_alignment( \@seqs )
757 :     # \@packed_seqs = pack_alignment( @seqs )
758 :     # \@packed_seqs = pack_alignment( \@seqs )
759 : golsen 1.5 #
760 :     #-----------------------------------------------------------------------------
761 :    
762 :     sub pack_alignment
763 :     {
764 :     my @seqs = ( ref( $_[0] ) eq 'ARRAY' and ref( $_[0]->[0] ) eq 'ARRAY' ) ? @{$_[0] } : @_;
765 :     @seqs or return wantarray ? () : [];
766 :    
767 :     my $mask = pack_mask( $seqs[0]->[2] );
768 :     foreach ( @seqs[ 1 .. (@seqs-1) ] )
769 :     {
770 :     $mask |= pack_mask( $_->[2] );
771 :     }
772 :    
773 :     my $seq;
774 :     my @seqs2 = map { $seq = $_->[2] & $mask;
775 :     $seq =~ tr/\000//d;
776 :     [ $_->[0], $_->[1], $seq ]
777 :     }
778 :     @seqs;
779 :    
780 :     return wantarray ? @seqs2 : \@seqs2;
781 :     }
782 :    
783 :     sub pack_mask
784 :     {
785 :     my $mask = shift;
786 :     $mask =~ tr/-/\000/;
787 :     $mask =~ tr/\000/\377/c;
788 :     return $mask;
789 :     }
790 : efrank 1.1
791 :     #-----------------------------------------------------------------------------
792 :     # Some simple sequence manipulations:
793 :     #
794 :     # @entry = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );
795 :     # @entry = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );
796 :     # @entry = complement_DNA_entry( @seq_entry [, $fix_id] );
797 :     # @entry = complement_RNA_entry( @seq_entry [, $fix_id] );
798 :     # $DNAseq = complement_DNA_seq( $NA_seq );
799 :     # $RNAseq = complement_RNA_seq( $NA_seq );
800 :     # $DNAseq = to_DNA_seq( $NA_seq );
801 :     # $RNAseq = to_RNA_seq( $NA_seq );
802 :     #
803 :     #-----------------------------------------------------------------------------
804 :    
805 :     sub subseq_DNA_entry {
806 :     my ($id, $desc, @rest) = @_;
807 :     wantarray || die "subseq_DNA_entry requires array context\n";
808 :    
809 :     my $seq;
810 :     ($id, $seq) = subseq_nt(1, $id, @rest); # 1 is for DNA, not RNA
811 :     return ($id, $desc, $seq);
812 :     }
813 :    
814 :    
815 :     sub subseq_RNA_entry {
816 :     my ($id, $desc, @rest) = @_;
817 :     wantarray || die "subseq_RNA_entry requires array context\n";
818 :    
819 :     my $seq;
820 :     ($id, $seq) = subseq_nt(0, $id, @rest); # 0 is for not DNA, i.e., RNA
821 :     return ($id, $desc, $seq);
822 :     }
823 :    
824 :    
825 :     sub subseq_nt {
826 :     my ($DNA, $id, $seq, $from, $to, $fix_id) = @_;
827 :     $fix_id ||= 0; # fix undef value
828 :    
829 :     my $len = length($seq);
830 :     if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
831 :     if (! $to || ( $to eq '$' ) || ( $to eq "" ) ) { $to = $len }
832 :    
833 :     my $left = ( $from < $to ) ? $from : $to;
834 :     my $right = ( $from < $to ) ? $to : $from;
835 :     if ( ( $right < 1 ) || ( $left > $len ) ) { return ($id, "") }
836 :     if ( $right > $len ) { $right = $len }
837 :     if ( $left < 1 ) { $left = 1 }
838 :    
839 :     $seq = substr($seq, $left-1, $right-$left+1);
840 :     if ( $from > $to ) {
841 :     $seq = reverse $seq;
842 :     if ( $DNA ) {
843 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
844 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
845 :     }
846 :     else {
847 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
848 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
849 :     }
850 :     }
851 :    
852 :     if ( $fix_id ) {
853 : golsen 1.2 if ( ( $id =~ s/_(\d+)_(\d+)$// )
854 : efrank 1.1 && ( abs($2-$1)+1 == $len ) ) {
855 :     if ( $1 <= $2 ) { $from += $1 - 1; $to += $1 - 1 }
856 :     else { $from = $1 + 1 - $from; $to = $1 + 1 - $to }
857 :     }
858 :     $id .= "_" . $from . "_" . $to;
859 :     }
860 :    
861 :     return ($id, $seq);
862 :     }
863 :    
864 :    
865 : golsen 1.9 sub DNA_subseq
866 :     {
867 :     my ( $seq, $from, $to ) = @_;
868 :    
869 :     my $len = ref( $seq ) eq 'SCALAR' ? length( $$seq )
870 :     : length( $seq );
871 :     if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
872 :     if ( ( $to eq '$' ) || ( ! $to ) ) { $to = $len }
873 :    
874 :     my $left = ( $from < $to ) ? $from : $to;
875 :     my $right = ( $from < $to ) ? $to : $from;
876 :     if ( ( $right < 1 ) || ( $left > $len ) ) { return "" }
877 :     if ( $right > $len ) { $right = $len }
878 :     if ( $left < 1 ) { $left = 1 }
879 :    
880 :     my $subseq = ref( $seq ) eq 'SCALAR' ? substr( $$seq, $left-1, $right-$left+1 )
881 :     : substr( $seq, $left-1, $right-$left+1 );
882 :    
883 :     if ( $from > $to )
884 :     {
885 :     $subseq = reverse $subseq;
886 :     $subseq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
887 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
888 :     }
889 :    
890 :     $subseq
891 :     }
892 :    
893 :    
894 :     sub RNA_subseq
895 :     {
896 :     my ( $seq, $from, $to ) = @_;
897 :    
898 :     my $len = ref( $seq ) eq 'SCALAR' ? length( $$seq )
899 :     : length( $seq );
900 :     if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
901 :     if ( ( $to eq '$' ) || ( ! $to ) ) { $to = $len }
902 :    
903 :     my $left = ( $from < $to ) ? $from : $to;
904 :     my $right = ( $from < $to ) ? $to : $from;
905 :     if ( ( $right < 1 ) || ( $left > $len ) ) { return "" }
906 :     if ( $right > $len ) { $right = $len }
907 :     if ( $left < 1 ) { $left = 1 }
908 :    
909 :     my $subseq = ref( $seq ) eq 'SCALAR' ? substr( $$seq, $left-1, $right-$left+1 )
910 :     : substr( $seq, $left-1, $right-$left+1 );
911 :    
912 :     if ( $from > $to )
913 :     {
914 :     $subseq = reverse $subseq;
915 :     $subseq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
916 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
917 :     }
918 :    
919 :     $subseq
920 :     }
921 :    
922 :    
923 : efrank 1.1 sub complement_DNA_entry {
924 :     my ($id, $desc, $seq, $fix_id) = @_;
925 :     $fix_id ||= 0; # fix undef values
926 :    
927 :     wantarray || die "complement_DNA_entry requires array context\n";
928 :     $seq = reverse $seq;
929 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
930 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
931 :     if ($fix_id) {
932 : golsen 1.2 if ($id =~ s/_(\d+)_(\d+)$//) {
933 : efrank 1.1 $id .= "_" . $2 . "_" . $1;
934 :     }
935 :     else {
936 :     $id .= "_" . length($seq) . "_1";
937 :     }
938 :     }
939 :    
940 :     return ($id, $desc, $seq);
941 :     }
942 :    
943 :    
944 :     sub complement_RNA_entry {
945 :     my ($id, $desc, $seq, $fix_id) = @_;
946 :     $fix_id ||= 0; # fix undef values
947 :    
948 :     wantarray || die "complement_DNA_entry requires array context\n";
949 :     $seq = reverse $seq;
950 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
951 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
952 :     if ($fix_id) {
953 : golsen 1.2 if ($id =~ s/_(\d+)_(\d+)$//) {
954 : efrank 1.1 $id .= "_" . $2 . "_" . $1;
955 :     }
956 :     else {
957 :     $id .= "_" . length($seq) . "_1";
958 :     }
959 :     }
960 :    
961 :     return ($id, $desc, $seq);
962 :     }
963 :    
964 :    
965 :     sub complement_DNA_seq {
966 :     my $seq = reverse shift;
967 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
968 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
969 :     return $seq;
970 :     }
971 :    
972 :    
973 :     sub complement_RNA_seq {
974 :     my $seq = reverse shift;
975 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
976 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
977 :     return $seq;
978 :     }
979 :    
980 :    
981 :     sub to_DNA_seq {
982 :     my $seq = shift;
983 :     $seq =~ tr/Uu/Tt/;
984 :     return $seq;
985 :     }
986 :    
987 :    
988 :     sub to_RNA_seq {
989 :     my $seq = shift;
990 :     $seq =~ tr/Tt/Uu/;
991 :     return $seq;
992 :     }
993 :    
994 :    
995 : golsen 1.2 sub pack_seq {
996 :     my $seq = shift;
997 :     $seq =~ tr/A-Za-z//cd;
998 :     return $seq;
999 :     }
1000 :    
1001 :    
1002 : efrank 1.1 sub clean_ae_sequence {
1003 :     $_ = shift;
1004 :     $_ = to7bit($_);
1005 :     s/[+]/1/g;
1006 :     s/[^0-9A-IK-NP-Za-ik-np-z~.-]/-/g;
1007 :     return $_;
1008 :     }
1009 :    
1010 :    
1011 :     sub to7bit {
1012 :     $_ = shift;
1013 :     my ($o, $c);
1014 :     while (/\\([0-3][0-7][0-7])/) {
1015 :     $o = oct($1) % 128;
1016 :     $c = sprintf("%c", $o);
1017 :     s/\\$1/$c/g;
1018 :     }
1019 :     return $_;
1020 :     }
1021 :    
1022 :    
1023 :     sub to8bit {
1024 :     $_ = shift;
1025 :     my ($o, $c);
1026 :     while (/\\([0-3][0-7][0-7])/) {
1027 :     $o = oct($1);
1028 :     $c = sprintf("%c", $o);
1029 :     s/\\$1/$c/g;
1030 :     }
1031 :     return $_;
1032 :     }
1033 :    
1034 :    
1035 :    
1036 :     #-----------------------------------------------------------------------------
1037 :     # Translate nucleotides to one letter protein:
1038 :     #
1039 : golsen 1.10 # $seq = translate_seq( $seq [, $met_start] )
1040 :     # $aa = translate_codon( $triplet )
1041 :     # $aa = translate_DNA_codon( $triplet ) # Does not rely on DNA
1042 :     # $aa = translate_uc_DNA_codon( $triplet ) # Does not rely on uc or DNA
1043 : efrank 1.1 #
1044 :     # User-supplied genetic code must be upper case index and match the
1045 :     # DNA versus RNA type of sequence
1046 :     #
1047 : golsen 1.10 # $seq = translate_seq_with_user_code( $seq, $gen_code_hash [, $met_start] )
1048 : efrank 1.1 #
1049 :     #-----------------------------------------------------------------------------
1050 :    
1051 : golsen 1.2 @aa_1_letter_order = qw( A C D E F G H I K L M N P Q R S T V W Y ); # Alpha by 1 letter
1052 :     @aa_3_letter_order = qw( A R N D C Q E G H I L K M F P S T W Y V ); # PAM matrix order
1053 :     @aa_n_codon_order = qw( L R S A G P T V I C D E F H K N Q Y M W );
1054 :    
1055 :     %genetic_code = (
1056 :    
1057 :     # DNA version
1058 :    
1059 : golsen 1.10 TTT => 'F', TCT => 'S', TAT => 'Y', TGT => 'C',
1060 :     TTC => 'F', TCC => 'S', TAC => 'Y', TGC => 'C',
1061 :     TTA => 'L', TCA => 'S', TAA => '*', TGA => '*',
1062 :     TTG => 'L', TCG => 'S', TAG => '*', TGG => 'W',
1063 :     CTT => 'L', CCT => 'P', CAT => 'H', CGT => 'R',
1064 :     CTC => 'L', CCC => 'P', CAC => 'H', CGC => 'R',
1065 :     CTA => 'L', CCA => 'P', CAA => 'Q', CGA => 'R',
1066 :     CTG => 'L', CCG => 'P', CAG => 'Q', CGG => 'R',
1067 :     ATT => 'I', ACT => 'T', AAT => 'N', AGT => 'S',
1068 :     ATC => 'I', ACC => 'T', AAC => 'N', AGC => 'S',
1069 :     ATA => 'I', ACA => 'T', AAA => 'K', AGA => 'R',
1070 :     ATG => 'M', ACG => 'T', AAG => 'K', AGG => 'R',
1071 :     GTT => 'V', GCT => 'A', GAT => 'D', GGT => 'G',
1072 :     GTC => 'V', GCC => 'A', GAC => 'D', GGC => 'G',
1073 :     GTA => 'V', GCA => 'A', GAA => 'E', GGA => 'G',
1074 :     GTG => 'V', GCG => 'A', GAG => 'E', GGG => 'G',
1075 : golsen 1.2
1076 : efrank 1.1 # The following ambiguous encodings are not necessary, but
1077 : golsen 1.2 # speed up the processing of some ambiguous triplets:
1078 : efrank 1.1
1079 : golsen 1.10 TTY => 'F', TCY => 'S', TAY => 'Y', TGY => 'C',
1080 :     TTR => 'L', TCR => 'S', TAR => '*',
1081 :     TCN => 'S',
1082 :     CTY => 'L', CCY => 'P', CAY => 'H', CGY => 'R',
1083 :     CTR => 'L', CCR => 'P', CAR => 'Q', CGR => 'R',
1084 :     CTN => 'L', CCN => 'P', CGN => 'R',
1085 :     ATY => 'I', ACY => 'T', AAY => 'N', AGY => 'S',
1086 :     ACR => 'T', AAR => 'K', AGR => 'R',
1087 :     ACN => 'T',
1088 :     GTY => 'V', GCY => 'A', GAY => 'D', GGY => 'G',
1089 :     GTR => 'V', GCR => 'A', GAR => 'E', GGR => 'G',
1090 :     GTN => 'V', GCN => 'A', GGN => 'G'
1091 : golsen 1.2 );
1092 :    
1093 : golsen 1.10 # Add RNA by construction:
1094 :    
1095 :     foreach ( grep { /T/ } keys %genetic_code )
1096 :     {
1097 :     my $codon = $_;
1098 :     $codon =~ s/T/U/g;
1099 :     $genetic_code{ $codon } = lc $genetic_code{ $_ }
1100 :     }
1101 : golsen 1.2
1102 :     # Add lower case by construction:
1103 :    
1104 : golsen 1.10 foreach ( keys %genetic_code )
1105 :     {
1106 : golsen 1.2 $genetic_code{ lc $_ } = lc $genetic_code{ $_ }
1107 :     }
1108 :    
1109 :    
1110 : golsen 1.9 # Construct the genetic code with selenocysteine by difference:
1111 : golsen 1.2
1112 : golsen 1.10 %genetic_code_with_U = %genetic_code;
1113 :     $genetic_code_with_U{ TGA } = 'U';
1114 :     $genetic_code_with_U{ tga } = 'u';
1115 :     $genetic_code_with_U{ UGA } = 'U';
1116 :     $genetic_code_with_U{ uga } = 'u';
1117 : golsen 1.2
1118 :    
1119 :     %amino_acid_codons_DNA = (
1120 : overbeek 1.4 L => [ qw( TTA TTG CTA CTG CTT CTC ) ],
1121 :     R => [ qw( AGA AGG CGA CGG CGT CGC ) ],
1122 :     S => [ qw( AGT AGC TCA TCG TCT TCC ) ],
1123 :     A => [ qw( GCA GCG GCT GCC ) ],
1124 :     G => [ qw( GGA GGG GGT GGC ) ],
1125 :     P => [ qw( CCA CCG CCT CCC ) ],
1126 :     T => [ qw( ACA ACG ACT ACC ) ],
1127 :     V => [ qw( GTA GTG GTT GTC ) ],
1128 :     I => [ qw( ATA ATT ATC ) ],
1129 :     C => [ qw( TGT TGC ) ],
1130 :     D => [ qw( GAT GAC ) ],
1131 :     E => [ qw( GAA GAG ) ],
1132 :     F => [ qw( TTT TTC ) ],
1133 :     H => [ qw( CAT CAC ) ],
1134 :     K => [ qw( AAA AAG ) ],
1135 :     N => [ qw( AAT AAC ) ],
1136 :     Q => [ qw( CAA CAG ) ],
1137 :     Y => [ qw( TAT TAC ) ],
1138 :     M => [ qw( ATG ) ],
1139 :     U => [ qw( TGA ) ],
1140 :     W => [ qw( TGG ) ],
1141 :    
1142 :     l => [ qw( tta ttg cta ctg ctt ctc ) ],
1143 :     r => [ qw( aga agg cga cgg cgt cgc ) ],
1144 :     s => [ qw( agt agc tca tcg tct tcc ) ],
1145 :     a => [ qw( gca gcg gct gcc ) ],
1146 :     g => [ qw( gga ggg ggt ggc ) ],
1147 :     p => [ qw( cca ccg cct ccc ) ],
1148 :     t => [ qw( aca acg act acc ) ],
1149 :     v => [ qw( gta gtg gtt gtc ) ],
1150 :     i => [ qw( ata att atc ) ],
1151 :     c => [ qw( tgt tgc ) ],
1152 :     d => [ qw( gat gac ) ],
1153 :     e => [ qw( gaa gag ) ],
1154 :     f => [ qw( ttt ttc ) ],
1155 :     h => [ qw( cat cac ) ],
1156 :     k => [ qw( aaa aag ) ],
1157 :     n => [ qw( aat aac ) ],
1158 :     q => [ qw( caa cag ) ],
1159 :     y => [ qw( tat tac ) ],
1160 :     m => [ qw( atg ) ],
1161 :     u => [ qw( tga ) ],
1162 :     w => [ qw( tgg ) ],
1163 : golsen 1.2
1164 : overbeek 1.4 '*' => [ qw( TAA TAG TGA ) ]
1165 : efrank 1.1 );
1166 :    
1167 :    
1168 : golsen 1.2
1169 :     %amino_acid_codons_RNA = (
1170 : overbeek 1.4 L => [ qw( UUA UUG CUA CUG CUU CUC ) ],
1171 :     R => [ qw( AGA AGG CGA CGG CGU CGC ) ],
1172 :     S => [ qw( AGU AGC UCA UCG UCU UCC ) ],
1173 :     A => [ qw( GCA GCG GCU GCC ) ],
1174 :     G => [ qw( GGA GGG GGU GGC ) ],
1175 :     P => [ qw( CCA CCG CCU CCC ) ],
1176 :     T => [ qw( ACA ACG ACU ACC ) ],
1177 :     V => [ qw( GUA GUG GUU GUC ) ],
1178 :     B => [ qw( GAU GAC AAU AAC ) ],
1179 :     Z => [ qw( GAA GAG CAA CAG ) ],
1180 :     I => [ qw( AUA AUU AUC ) ],
1181 :     C => [ qw( UGU UGC ) ],
1182 :     D => [ qw( GAU GAC ) ],
1183 :     E => [ qw( GAA GAG ) ],
1184 :     F => [ qw( UUU UUC ) ],
1185 :     H => [ qw( CAU CAC ) ],
1186 :     K => [ qw( AAA AAG ) ],
1187 :     N => [ qw( AAU AAC ) ],
1188 :     Q => [ qw( CAA CAG ) ],
1189 :     Y => [ qw( UAU UAC ) ],
1190 :     M => [ qw( AUG ) ],
1191 :     U => [ qw( UGA ) ],
1192 :     W => [ qw( UGG ) ],
1193 :    
1194 :     l => [ qw( uua uug cua cug cuu cuc ) ],
1195 :     r => [ qw( aga agg cga cgg cgu cgc ) ],
1196 :     s => [ qw( agu agc uca ucg ucu ucc ) ],
1197 :     a => [ qw( gca gcg gcu gcc ) ],
1198 :     g => [ qw( gga ggg ggu ggc ) ],
1199 :     p => [ qw( cca ccg ccu ccc ) ],
1200 :     t => [ qw( aca acg acu acc ) ],
1201 :     v => [ qw( gua gug guu guc ) ],
1202 :     b => [ qw( gau gac aau aac ) ],
1203 :     z => [ qw( gaa gag caa cag ) ],
1204 :     i => [ qw( aua auu auc ) ],
1205 :     c => [ qw( ugu ugc ) ],
1206 :     d => [ qw( gau gac ) ],
1207 :     e => [ qw( gaa gag ) ],
1208 :     f => [ qw( uuu uuc ) ],
1209 :     h => [ qw( cau cac ) ],
1210 :     k => [ qw( aaa aag ) ],
1211 :     n => [ qw( aau aac ) ],
1212 :     q => [ qw( caa cag ) ],
1213 :     y => [ qw( uau uac ) ],
1214 :     m => [ qw( aug ) ],
1215 :     u => [ qw( uga ) ],
1216 :     w => [ qw( ugg ) ],
1217 : golsen 1.2
1218 : overbeek 1.4 '*' => [ qw( UAA UAG UGA ) ]
1219 : golsen 1.2 );
1220 :    
1221 :    
1222 :     %n_codon_for_aa = map {
1223 :     $_ => scalar @{ $amino_acid_codons_DNA{ $_ } }
1224 :     } keys %amino_acid_codons_DNA;
1225 :    
1226 :    
1227 :     %reverse_genetic_code_DNA = (
1228 : overbeek 1.4 A => "GCN", a => "gcn",
1229 :     C => "TGY", c => "tgy",
1230 :     D => "GAY", d => "gay",
1231 :     E => "GAR", e => "gar",
1232 :     F => "TTY", f => "tty",
1233 :     G => "GGN", g => "ggn",
1234 :     H => "CAY", h => "cay",
1235 :     I => "ATH", i => "ath",
1236 :     K => "AAR", k => "aar",
1237 :     L => "YTN", l => "ytn",
1238 :     M => "ATG", m => "atg",
1239 :     N => "AAY", n => "aay",
1240 :     P => "CCN", p => "ccn",
1241 :     Q => "CAR", q => "car",
1242 :     R => "MGN", r => "mgn",
1243 :     S => "WSN", s => "wsn",
1244 :     T => "ACN", t => "acn",
1245 :     U => "TGA", u => "tga",
1246 :     V => "GTN", v => "gtn",
1247 :     W => "TGG", w => "tgg",
1248 :     X => "NNN", x => "nnn",
1249 :     Y => "TAY", y => "tay",
1250 :     '*' => "TRR"
1251 : golsen 1.2 );
1252 :    
1253 :     %reverse_genetic_code_RNA = (
1254 : overbeek 1.4 A => "GCN", a => "gcn",
1255 :     C => "UGY", c => "ugy",
1256 :     D => "GAY", d => "gay",
1257 :     E => "GAR", e => "gar",
1258 :     F => "UUY", f => "uuy",
1259 :     G => "GGN", g => "ggn",
1260 :     H => "CAY", h => "cay",
1261 :     I => "AUH", i => "auh",
1262 :     K => "AAR", k => "aar",
1263 :     L => "YUN", l => "yun",
1264 :     M => "AUG", m => "aug",
1265 :     N => "AAY", n => "aay",
1266 :     P => "CCN", p => "ccn",
1267 :     Q => "CAR", q => "car",
1268 :     R => "MGN", r => "mgn",
1269 :     S => "WSN", s => "wsn",
1270 :     T => "ACN", t => "acn",
1271 :     U => "UGA", u => "uga",
1272 :     V => "GUN", v => "gun",
1273 :     W => "UGG", w => "ugg",
1274 :     X => "NNN", x => "nnn",
1275 :     Y => "UAY", y => "uay",
1276 :     '*' => "URR"
1277 : golsen 1.2 );
1278 :    
1279 :    
1280 :     %DNA_letter_can_be = (
1281 : efrank 1.1 A => ["A"], a => ["a"],
1282 :     B => ["C", "G", "T"], b => ["c", "g", "t"],
1283 :     C => ["C"], c => ["c"],
1284 :     D => ["A", "G", "T"], d => ["a", "g", "t"],
1285 :     G => ["G"], g => ["g"],
1286 :     H => ["A", "C", "T"], h => ["a", "c", "t"],
1287 :     K => ["G", "T"], k => ["g", "t"],
1288 :     M => ["A", "C"], m => ["a", "c"],
1289 :     N => ["A", "C", "G", "T"], n => ["a", "c", "g", "t"],
1290 :     R => ["A", "G"], r => ["a", "g"],
1291 :     S => ["C", "G"], s => ["c", "g"],
1292 :     T => ["T"], t => ["t"],
1293 :     U => ["T"], u => ["t"],
1294 :     V => ["A", "C", "G"], v => ["a", "c", "g"],
1295 :     W => ["A", "T"], w => ["a", "t"],
1296 :     Y => ["C", "T"], y => ["c", "t"]
1297 :     );
1298 :    
1299 :    
1300 : golsen 1.2 %RNA_letter_can_be = (
1301 : efrank 1.1 A => ["A"], a => ["a"],
1302 :     B => ["C", "G", "U"], b => ["c", "g", "u"],
1303 :     C => ["C"], c => ["c"],
1304 :     D => ["A", "G", "U"], d => ["a", "g", "u"],
1305 :     G => ["G"], g => ["g"],
1306 :     H => ["A", "C", "U"], h => ["a", "c", "u"],
1307 :     K => ["G", "U"], k => ["g", "u"],
1308 :     M => ["A", "C"], m => ["a", "c"],
1309 :     N => ["A", "C", "G", "U"], n => ["a", "c", "g", "u"],
1310 :     R => ["A", "G"], r => ["a", "g"],
1311 :     S => ["C", "G"], s => ["c", "g"],
1312 :     T => ["U"], t => ["u"],
1313 :     U => ["U"], u => ["u"],
1314 :     V => ["A", "C", "G"], v => ["a", "c", "g"],
1315 :     W => ["A", "U"], w => ["a", "u"],
1316 :     Y => ["C", "U"], y => ["c", "u"]
1317 :     );
1318 :    
1319 :    
1320 : overbeek 1.4 %one_letter_to_three_letter_aa = (
1321 :     A => "Ala", a => "Ala",
1322 :     B => "Asx", b => "Asx",
1323 :     C => "Cys", c => "Cys",
1324 :     D => "Asp", d => "Asp",
1325 :     E => "Glu", e => "Glu",
1326 :     F => "Phe", f => "Phe",
1327 :     G => "Gly", g => "Gly",
1328 :     H => "His", h => "His",
1329 :     I => "Ile", i => "Ile",
1330 :     K => "Lys", k => "Lys",
1331 :     L => "Leu", l => "Leu",
1332 :     M => "Met", m => "Met",
1333 :     N => "Asn", n => "Asn",
1334 :     P => "Pro", p => "Pro",
1335 :     Q => "Gln", q => "Gln",
1336 :     R => "Arg", r => "Arg",
1337 :     S => "Ser", s => "Ser",
1338 :     T => "Thr", t => "Thr",
1339 :     U => "Sec", u => "Sec",
1340 :     V => "Val", v => "Val",
1341 :     W => "Trp", w => "Trp",
1342 :     X => "Xxx", x => "Xxx",
1343 :     Y => "Tyr", y => "Tyr",
1344 :     Z => "Glx", z => "Glx",
1345 :     '*' => "***"
1346 :     );
1347 : golsen 1.2
1348 :    
1349 :     %three_letter_to_one_letter_aa = (
1350 :     ALA => "A", Ala => "A", ala => "a",
1351 :     ARG => "R", Arg => "R", arg => "r",
1352 :     ASN => "N", Asn => "N", asn => "n",
1353 :     ASP => "D", Asp => "D", asp => "d",
1354 :     ASX => "B", Asx => "B", asx => "b",
1355 :     CYS => "C", Cys => "C", cys => "c",
1356 :     GLN => "Q", Gln => "Q", gln => "q",
1357 :     GLU => "E", Glu => "E", glu => "e",
1358 :     GLX => "Z", Glx => "Z", glx => "z",
1359 :     GLY => "G", Gly => "G", gly => "g",
1360 :     HIS => "H", His => "H", his => "h",
1361 :     ILE => "I", Ile => "I", ile => "i",
1362 :     LEU => "L", Leu => "L", leu => "l",
1363 :     LYS => "K", Lys => "K", lys => "k",
1364 :     MET => "M", Met => "M", met => "m",
1365 :     PHE => "F", Phe => "F", phe => "f",
1366 :     PRO => "P", Pro => "P", pro => "p",
1367 :     SEC => "U", Sec => "U", sec => "u",
1368 :     SER => "S", Ser => "S", ser => "s",
1369 :     THR => "T", Thr => "T", thr => "t",
1370 :     TRP => "W", Trp => "W", trp => "w",
1371 :     TYR => "Y", Tyr => "Y", tyr => "y",
1372 :     VAL => "V", Val => "V", val => "v",
1373 :     XAA => "X", Xaa => "X", xaa => "x",
1374 :     XXX => "X", Xxx => "X", xxx => "x",
1375 :     '***' => "*"
1376 : efrank 1.1 );
1377 :    
1378 :    
1379 :     #-----------------------------------------------------------------------------
1380 : golsen 1.10 # Translate nucleotides to one letter protein. Respects case of the
1381 :     # nucleotide sequence.
1382 : efrank 1.1 #
1383 : golsen 1.10 # $aa = translate_seq( $nt, $met_start )
1384 :     # $aa = translate_seq( $nt )
1385 : efrank 1.1 #
1386 :     #-----------------------------------------------------------------------------
1387 :    
1388 : golsen 1.10 sub translate_seq
1389 :     {
1390 :     my $seq = shift;
1391 : efrank 1.1 $seq =~ tr/-//d; # remove gaps
1392 :    
1393 : golsen 1.10 my @codons = $seq =~ m/(...?)/g; # Will try to translate last 2 nt
1394 :    
1395 :     # A second argument that is true forces first amino acid to be Met
1396 : efrank 1.1
1397 : golsen 1.10 my @met;
1398 :     if ( ( shift @_ ) && ( my $codon1 = shift @codons ) )
1399 :     {
1400 :     push @met, ( $codon1 =~ /[a-z]/ ? 'm' : 'M' );
1401 : efrank 1.1 }
1402 :    
1403 : golsen 1.10 join( '', @met, map { translate_codon( $_ ) } @codons )
1404 : efrank 1.1 }
1405 :    
1406 :    
1407 :     #-----------------------------------------------------------------------------
1408 : golsen 1.10 # Translate a single triplet with "universal" genetic code.
1409 : efrank 1.1 #
1410 :     # $aa = translate_codon( $triplet )
1411 : golsen 1.10 # $aa = translate_DNA_codon( $triplet )
1412 :     # $aa = translate_uc_DNA_codon( $triplet )
1413 : efrank 1.1 #
1414 :     #-----------------------------------------------------------------------------
1415 :    
1416 : golsen 1.10 sub translate_DNA_codon { translate_codon( @_ ) }
1417 : efrank 1.1
1418 : golsen 1.10 sub translate_uc_DNA_codon { translate_codon( uc $_[0] ) }
1419 : efrank 1.1
1420 : golsen 1.10 sub translate_codon
1421 :     {
1422 :     my $codon = shift;
1423 :     $codon =~ tr/Uu/Tt/; # Make it DNA
1424 :    
1425 :     # Try a simple lookup:
1426 : efrank 1.1
1427 :     my $aa;
1428 : golsen 1.10 if ( $aa = $genetic_code{ $codon } ) { return $aa }
1429 : efrank 1.1
1430 : golsen 1.10 # Attempt to recover from mixed-case codons:
1431 : efrank 1.1
1432 : golsen 1.10 $codon = ( $codon =~ /[a-z]/ ) ? lc $codon : uc $codon;
1433 : efrank 1.1 if ( $aa = $genetic_code{ $codon } ) { return $aa }
1434 :    
1435 : golsen 1.10 # The code defined above catches simple N, R and Y ambiguities in the
1436 :     # third position. Other codons (e.g., GG[KMSWBDHV], or even GG) might
1437 :     # be unambiguously translated by converting the last position to N and
1438 :     # seeing if this is in the code table:
1439 :    
1440 :     my $N = ( $codon =~ /[a-z]/ ) ? 'n' : 'N';
1441 :     if ( $aa = $genetic_code{ substr($codon,0,2) . $N } ) { return $aa }
1442 :    
1443 :     # Test that codon is valid for an unambiguous aa:
1444 :    
1445 :     my $X = ( $codon =~ /[a-z]/ ) ? 'x' : 'X';
1446 :     if ( $codon !~ m/^[ACGTMY][ACGT][ACGTKMRSWYBDHVN]$/i
1447 :     && $codon !~ m/^YT[AGR]$/i # Leu YTR
1448 :     && $codon !~ m/^MG[AGR]$/i # Arg MGR
1449 :     )
1450 :     {
1451 :     return $X;
1452 :     }
1453 : efrank 1.1
1454 :     # Expand all ambiguous nucleotides to see if they all yield same aa.
1455 :    
1456 : golsen 1.10 my @n1 = @{ $DNA_letter_can_be{ substr( $codon, 0, 1 ) } };
1457 :     my $n2 = substr( $codon, 1, 1 );
1458 :     my @n3 = @{ $DNA_letter_can_be{ substr( $codon, 2, 1 ) } };
1459 :     my @triples = map { my $n12 = $_ . $n2; map { $n12 . $_ } @n3 } @n1;
1460 :    
1461 :     my $triple = shift @triples;
1462 :     $aa = $genetic_code{ $triple };
1463 :     $aa or return $X;
1464 :    
1465 :     foreach $triple ( @triples ) { return $X if $aa ne $genetic_code{$triple} }
1466 : efrank 1.1
1467 : golsen 1.10 return $aa;
1468 : efrank 1.1 }
1469 :    
1470 :    
1471 :     #-----------------------------------------------------------------------------
1472 :     # Translate with a user-supplied genetic code to translate a sequence.
1473 :     # Diagnose the use of upper versus lower, and T versus U in the supplied
1474 :     # code, and transform the supplied nucleotide sequence to match.
1475 :     #
1476 : golsen 1.10 # $aa = translate_seq_with_user_code( $nt, \%gen_code )
1477 :     # $aa = translate_seq_with_user_code( $nt, \%gen_code, $start_with_met )
1478 : efrank 1.1 #
1479 : golsen 1.10 # Modified 2007-11-22 to be less intrusive in these diagnoses by sensing
1480 :     # the presence of both versions in the user code.
1481 : efrank 1.1 #-----------------------------------------------------------------------------
1482 :    
1483 : golsen 1.10 sub translate_seq_with_user_code
1484 :     {
1485 : efrank 1.1 my $seq = shift;
1486 :     $seq =~ tr/-//d; # remove gaps *** Why?
1487 :    
1488 : golsen 1.10 my $gc = shift; # Reference to hash of code
1489 :     if (! $gc || ref($gc) ne "HASH")
1490 :     {
1491 :     print STDERR "translate_seq_with_user_code needs genetic code hash as second argument.";
1492 :     return undef;
1493 :     }
1494 :    
1495 :     # Test code support for upper vs lower case:
1496 :    
1497 :     my ( $TTT, $UUU );
1498 :     if ( $gc->{AAA} && ! $gc->{aaa} ) # Uppercase only code table
1499 :     {
1500 :     $seq = uc $seq; # Uppercase sequence
1501 :     ( $TTT, $UUU ) = ( 'TTT', 'UUU' );
1502 :     }
1503 :     elsif ( $gc->{aaa} && ! $gc->{AAA} ) # Lowercase only code table
1504 :     {
1505 :     $seq = lc $seq; # Lowercase sequence
1506 :     ( $TTT, $UUU ) = ( 'ttt', 'uuu' );
1507 :     }
1508 :     elsif ( $gc->{aaa} )
1509 :     {
1510 :     ( $TTT, $UUU ) = ( 'ttt', 'uuu' );
1511 : efrank 1.1 }
1512 : golsen 1.10 else
1513 :     {
1514 :     print STDERR "User-supplied genetic code does not have aaa or AAA\n";
1515 :     return undef;
1516 : efrank 1.1 }
1517 :    
1518 : golsen 1.10 # Test code support for U vs T:
1519 : efrank 1.1
1520 : golsen 1.10 my $ambigs;
1521 :     if ( $gc->{$UUU} && ! $gc->{$TTT} ) # RNA only code table
1522 :     {
1523 :     $seq = tr/Tt/Uu/;
1524 : efrank 1.1 $ambigs = \%RNA_letter_can_be;
1525 :     }
1526 : golsen 1.10 elsif ( $gc->{$TTT} && ! $gc->{$UUU} ) # DNA only code table
1527 :     {
1528 :     $seq = tr/Uu/Tt/;
1529 : efrank 1.1 $ambigs = \%DNA_letter_can_be;
1530 :     }
1531 : golsen 1.10 else
1532 :     {
1533 :     my $t = $seq =~ tr/Tt//;
1534 :     my $u = $seq =~ tr/Uu//;
1535 :     $ambigs = ( $t > $u ) ? \%DNA_letter_can_be : \%RNA_letter_can_be;
1536 : efrank 1.1 }
1537 :    
1538 : golsen 1.10 # We can now do the codon-by-codon translation:
1539 : efrank 1.1
1540 : golsen 1.10 my @codons = $seq =~ m/(...?)/g; # will try to translate last 2 nt
1541 :    
1542 :     # A third argument that is true forces first amino acid to be Met
1543 : efrank 1.1
1544 : golsen 1.10 my @met;
1545 :     if ( ( shift @_ ) && ( my $codon1 = shift @codons ) )
1546 :     {
1547 :     push @met, ( $codon1 =~ /[a-z]/ ? 'm' : 'M' );
1548 : efrank 1.1 }
1549 :    
1550 : golsen 1.10 join( '', @met, map { translate_codon_with_user_code( $_, $gc, $ambigs ) } @codons )
1551 : efrank 1.1 }
1552 :    
1553 :    
1554 :     #-----------------------------------------------------------------------------
1555 : golsen 1.10 # Translate with user-supplied genetic code hash. No error check on the code.
1556 :     # Should only be called through translate_seq_with_user_code.
1557 : efrank 1.1 #
1558 : golsen 1.10 # $aa = translate_codon_with_user_code( $triplet, \%code, $ambig_table )
1559 : efrank 1.1 #
1560 :     # $triplet speaks for itself
1561 :     # $code ref to the hash with the codon translations
1562 :     # $ambig_table ref to hash with lists of nucleotides for each ambig code
1563 :     #-----------------------------------------------------------------------------
1564 :    
1565 : golsen 1.10 sub translate_codon_with_user_code
1566 :     {
1567 :     my ( $codon, $gc, $N, $X, $ambigs ) = @_;
1568 :    
1569 :     # Try a simple lookup:
1570 : efrank 1.1
1571 :     my $aa;
1572 : golsen 1.10 if ( $aa = $gc->{ $codon } ) { return $aa }
1573 : efrank 1.1
1574 : golsen 1.10 # Attempt to recover from mixed-case codons:
1575 : efrank 1.1
1576 : golsen 1.10 $codon = ( $codon =~ /[a-z]/ ) ? lc $codon : uc $codon;
1577 :     if ( $aa = $genetic_code{ $codon } ) { return $aa }
1578 : efrank 1.1
1579 : golsen 1.10 # Test that codon is valid for an unambiguous aa:
1580 : efrank 1.1
1581 : golsen 1.10 my $X = ( $codon =~ /[a-z]/ ) ? 'x' : 'X';
1582 :    
1583 :     if ( $codon =~ m/^[ACGTU][ACGTU]$/i ) # Add N?
1584 :     {
1585 :     $codon .= ( $codon =~ /[a-z]/ ) ? 'n' : 'N';
1586 :     }
1587 :     elsif ( $codon !~ m/^[ACGTUMY][ACGTU][ACGTUKMRSWYBDHVN]$/i ) # Valid?
1588 :     {
1589 :     return $X;
1590 :     }
1591 : efrank 1.1
1592 :     # Expand all ambiguous nucleotides to see if they all yield same aa.
1593 :    
1594 : golsen 1.10 my @n1 = @{ $ambigs->{ substr( $codon, 0, 1 ) } };
1595 :     my $n2 = substr( $codon, 1, 1 );
1596 :     my @n3 = @{ $ambigs->{ substr( $codon, 2, 1 ) } };
1597 :     my @triples = map { my $n12 = $_ . $n2; map { $n12 . $_ } @n3 } @n1;
1598 :    
1599 :     my $triple = shift @triples;
1600 :     $aa = $gc->{ $triple } || $gc->{ lc $triple } || $gc->{ uc $triple };
1601 :     $aa or return $X;
1602 :    
1603 :     foreach $triple ( @triples )
1604 :     {
1605 :     return $X if $aa ne ( $gc->{$triple} || $gc->{lc $triple} || $gc->{uc $triple} );
1606 : efrank 1.1 }
1607 :    
1608 : golsen 1.10 return $aa;
1609 : efrank 1.1 }
1610 :    
1611 :    
1612 :     #-----------------------------------------------------------------------------
1613 :     # Read a list of intervals from a file.
1614 :     # Allow id_start_end, or id \s start \s end formats
1615 :     #
1616 : overbeek 1.4 # @intervals = read_intervals( \*FILEHANDLE )
1617 : efrank 1.1 #-----------------------------------------------------------------------------
1618 :     sub read_intervals {
1619 :     my $fh = shift;
1620 :     my @intervals = ();
1621 :    
1622 :     while (<$fh>) {
1623 :     chomp;
1624 :     /^(\S+)_(\d+)_(\d+)(\s.*)?$/ # id_start_end WIT2
1625 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/ # id_start-end ???
1626 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/ # id=start=end Badger
1627 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/ # id \s start \s end
1628 :     || next;
1629 :    
1630 :     # Matched a pattern. Store reference to (id, left, right):
1631 :     push @intervals, ($2 < $3) ? [ $1, $2+0, $3+0 ]
1632 :     : [ $1, $3+0, $2+0 ];
1633 :     }
1634 :     return @intervals;
1635 :     }
1636 :    
1637 :    
1638 :     #-----------------------------------------------------------------------------
1639 : golsen 1.2 # Convert a list of intervals to read [ id, left_end, right_end ].
1640 :     #
1641 :     # @intervals = standardize_intervals( @interval_refs )
1642 :     #-----------------------------------------------------------------------------
1643 :     sub standardize_intervals {
1644 :     map { ( $_->[1] < $_->[2] ) ? $_ : [ $_->[0], $_->[2], $_->[1] ] } @_;
1645 :     }
1646 :    
1647 :    
1648 :     #-----------------------------------------------------------------------------
1649 : efrank 1.1 # Take the union of a list of intervals
1650 :     #
1651 :     # @joined = join_intervals( @interval_refs )
1652 :     #-----------------------------------------------------------------------------
1653 :     sub join_intervals {
1654 :     my @ordered = sort { $a->[0] cmp $b->[0] # first by id
1655 :     || $a->[1] <=> $b->[1] # next by left end
1656 :     || $b->[2] <=> $a->[2] # finally longest first
1657 :     } @_;
1658 :    
1659 :     my @joined = ();
1660 :     my $n_int = @ordered;
1661 :    
1662 :     my ($cur_id) = "";
1663 :     my ($cur_left) = -1;
1664 :     my ($cur_right) = -1;
1665 :     my ($new_id, $new_left, $new_right);
1666 :    
1667 :     for (my $i = 0; $i < $n_int; $i++) {
1668 :     ($new_id, $new_left, $new_right) = @{$ordered[$i]}; # get the new data
1669 :    
1670 :     if ( ( $new_id ne $cur_id) # if new contig
1671 :     || ( $new_left > $cur_right + 1) # or not touching previous
1672 :     ) { # push the previous interval
1673 :     if ($cur_id) { push (@joined, [ $cur_id, $cur_left, $cur_right ]) }
1674 :     $cur_id = $new_id; # update the current interval
1675 :     $cur_left = $new_left;
1676 :     $cur_right = $new_right;
1677 :     }
1678 :    
1679 :     elsif ($new_right > $cur_right) { # extend the right end if necessary
1680 :     $cur_right = $new_right;
1681 :     }
1682 :     }
1683 :    
1684 :     if ($cur_id) { push (@joined, [$cur_id, $cur_left, $cur_right]) }
1685 :     return @joined;
1686 :     }
1687 :    
1688 : golsen 1.2
1689 :     #-----------------------------------------------------------------------------
1690 :     # Split location strings to oriented intervals.
1691 :     #
1692 :     # @intervals = locations_2_intervals( @locations )
1693 :     # $interval = locations_2_intervals( $location )
1694 :     #-----------------------------------------------------------------------------
1695 :     sub locations_2_intervals {
1696 :     my @intervals = map { /^(\S+)_(\d+)_(\d+)(\s.*)?$/
1697 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/
1698 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/
1699 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/
1700 :     ? [ $1, $2+0, $3+0 ]
1701 :     : ()
1702 :     } @_;
1703 :    
1704 :     return wantarray ? @intervals : $intervals[0];
1705 :     }
1706 :    
1707 :    
1708 :     #-----------------------------------------------------------------------------
1709 :     # Read a list of oriented intervals from a file.
1710 :     # Allow id_start_end, or id \s start \s end formats
1711 :     #
1712 : overbeek 1.4 # @intervals = read_oriented_intervals( \*FILEHANDLE )
1713 : golsen 1.2 #-----------------------------------------------------------------------------
1714 :     sub read_oriented_intervals {
1715 :     my $fh = shift;
1716 :     my @intervals = ();
1717 :    
1718 :     while (<$fh>) {
1719 :     chomp;
1720 :     /^(\S+)_(\d+)_(\d+)(\s.*)?$/ # id_start_end WIT2
1721 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/ # id_start-end ???
1722 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/ # id=start=end Badger
1723 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/ # id \s start \s end
1724 :     || next;
1725 :    
1726 :     # Matched a pattern. Store reference to (id, start, end):
1727 :     push @intervals, [ $1, $2+0, $3+0 ];
1728 :     }
1729 :     return @intervals;
1730 :     }
1731 :    
1732 :    
1733 :     #-----------------------------------------------------------------------------
1734 :     # Reverse the orientation of a list of intervals
1735 :     #
1736 :     # @reversed = reverse_intervals( @interval_refs )
1737 :     #-----------------------------------------------------------------------------
1738 :     sub reverse_intervals {
1739 :     map { [ $_->[0], $_->[2], $_->[1] ] } @_;
1740 :     }
1741 :    
1742 :    
1743 : overbeek 1.4 #-----------------------------------------------------------------------------
1744 :     # Convert GenBank locations to SEED locations
1745 :     #
1746 :     # @seed_locs = gb_location_2_seed( $contig, @gb_locs )
1747 :     #-----------------------------------------------------------------------------
1748 :     sub gb_location_2_seed
1749 :     {
1750 :     my $contig = shift @_;
1751 :     $contig or die "First arg of gb_location_2_seed must be contig_id\n";
1752 :    
1753 :     map { join( ',', gb_loc_2_seed_2( $contig, $_ ) ) || undef } @_
1754 :     }
1755 :    
1756 :     sub gb_loc_2_seed_2
1757 :     {
1758 :     my ( $contig, $loc ) = @_;
1759 :    
1760 :     if ( $loc =~ /^(\d+)\.\.(\d+)$/ )
1761 :     {
1762 :     join( '_', $contig, $1, $2 )
1763 :     }
1764 :    
1765 :     elsif ( $loc =~ /^join\((.*)\)$/ )
1766 :     {
1767 :     $loc = $1;
1768 :     my $lvl = 0;
1769 :     for ( my $i = length( $loc )-1; $i >= 0; $i-- )
1770 :     {
1771 :     for ( substr( $loc, $i, 1 ) )
1772 :     {
1773 :     /,/ && ! $lvl and substr( $loc, $i, 1 ) = "\t";
1774 :     /\(/ and $lvl--;
1775 :     /\)/ and $lvl++;
1776 :     }
1777 :     }
1778 :     $lvl == 0 or print STDERR "Paren matching error: $loc\n" and die;
1779 :     map { gb_loc_2_seed_2( $contig, $_ ) } split /\t/, $loc
1780 :     }
1781 :    
1782 :     elsif ( $loc =~ /^complement\((.*)\)$/ )
1783 :     {
1784 :     map { s/_(\d+)_(\d+)$/_$2_$1/; $_ }
1785 :     reverse
1786 :     gb_loc_2_seed_2( $contig, $1 )
1787 :     }
1788 :    
1789 :     else
1790 :     {
1791 :     ()
1792 :     }
1793 :     }
1794 :    
1795 :    
1796 : golsen 1.7 #-----------------------------------------------------------------------------
1797 :     # Read qual.
1798 :     #
1799 :     # Save the contents in a list of refs to arrays: [ $id, $descript, \@qual ]
1800 :     #
1801 :     # @seq_entries = read_qual( ) # STDIN
1802 :     # \@seq_entries = read_qual( ) # STDIN
1803 :     # @seq_entries = read_qual( \*FILEHANDLE )
1804 :     # \@seq_entries = read_qual( \*FILEHANDLE )
1805 :     # @seq_entries = read_qual( $filename )
1806 :     # \@seq_entries = read_qual( $filename )
1807 :     #-----------------------------------------------------------------------------
1808 :     sub read_qual {
1809 :     my ( $fh, $name, $close, $unused ) = input_filehandle( $_[0] );
1810 :     $unused && die "Bad reference type (" . ref( $unused ) . ") passed to read_qual\n";
1811 :    
1812 :     my @quals = ();
1813 :     my ($id, $desc, $qual) = ("", "", []);
1814 :    
1815 :     while ( <$fh> ) {
1816 :     chomp;
1817 :     if (/^>\s*(\S+)(\s+(.*))?$/) { # new id
1818 :     if ($id && @$qual) { push @quals, [ $id, $desc, $qual ] }
1819 :     ($id, $desc, $qual) = ($1, $3 ? $3 : "", []);
1820 :     }
1821 :     else {
1822 :     push @$qual, split;
1823 :     }
1824 :     }
1825 :     close( $fh ) if $close;
1826 :    
1827 :     if ($id && @$qual) { push @quals, [ $id, $desc, $qual ] }
1828 :     return wantarray ? @quals : \@quals;
1829 :     }
1830 :    
1831 :    
1832 : efrank 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3