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

Annotation of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3