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

Annotation of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3