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

Annotation of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3