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

Annotation of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.3 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 : efrank 1.1 package gjoseqlib;
19 :    
20 : golsen 1.2 # A sequence entry is ( $id, $def, $seq )
21 :     # A list of entries is a list of references
22 :     #
23 :     # @seq_entry = read_next_fasta_seq( *FILEHANDLE )
24 :     # @seq_entries = read_fasta_seqs( *FILEHANDLE )
25 :     # $seq_ind = index_seq_list( @seq_entries ); # hash from ids to entry refs
26 :     #
27 :     # @seq_entry = seq_entry_by_id( \%seq_index, $seq_id );
28 :     # $seq_desc = seq_desc_by_id( \%seq_index, $seq_id );
29 :     # $seq = seq_data_by_id( \%seq_index, $seq_id );
30 :     #
31 :     # ( $id, $def ) = parse_fasta_title( $title )
32 :     # ( $id, $def ) = split_fasta_title( $title )
33 :     #
34 :     # print_seq_list_as_fasta( *FILEHANDLE, @seq_entry_list );
35 :     # print_seq_as_fasta( *FILEHANDLE, $id, $desc, $seq) ;
36 :     # print_seq_as_fasta( *FILEHANDLE, @seq_entry );
37 : efrank 1.1 # print_gb_locus( *FILEHANDLE, $locus, $def, $accession, $seq )
38 : golsen 1.2 #
39 :     # @entry = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );
40 :     # @entry = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );
41 :     # @entry = complement_DNA_entry( @seq_entry [, $fix_id] );
42 :     # @entry = complement_RNA_entry( @seq_entry [, $fix_id] );
43 :     # $DNAseq = complement_DNA_seq( $NA_seq );
44 :     # $RNAseq = complement_RNA_seq( $NA_seq );
45 :     # $DNAseq = to_DNA_seq( $NA_seq );
46 :     # $RNAseq = to_RNA_seq( $NA_seq );
47 :     # $seq = pack_seq( $sequence )
48 :     # $seq = clean_ae_sequence( $seq )
49 :     #
50 :     # $seq = translate_seq( $seq [, $met_start] )
51 :     # $aa = translate_codon( $triplet );
52 :     # $aa = translate_uc_DNA_codon( $upcase_DNA_triplet );
53 :     #
54 :     # User-supplied genetic code must be upper case index and match the
55 :     # DNA versus RNA type of sequence
56 :     #
57 :     # Locations (= oriented intervals) are ( id, start, end )
58 :     # Intervals are ( id, left, right )
59 :     #
60 :     # @intervals = read_intervals( *FILEHANDLE )
61 :     # @intervals = read_oriented_intervals( *FILEHANDLE )
62 :     # @intervals = standardize_intervals( @interval_refs ) # (id, left, right)
63 :     # @joined = join_intervals( @interval_refs )
64 :     # @intervals = locations_2_intervals( @locations )
65 :     # $interval = locations_2_intervals( $location )
66 :     # @reversed = reverse_intervals( @interval_refs ) # (id, end, start)
67 :    
68 :    
69 :     use strict;
70 : efrank 1.1
71 : golsen 1.2 use gjolib qw( wrap_text );
72 : efrank 1.1
73 : golsen 1.2 # Exported global variables:
74 : efrank 1.1
75 : golsen 1.2 our @aa_1_letter_order; # Alpha by 1 letter
76 :     our @aa_3_letter_order; # PAM matrix order
77 :     our @aa_n_codon_order;
78 :     our %genetic_code;
79 :     our %genetic_code_with_U;
80 :     our %amino_acid_codons_DNA;
81 :     our %amino_acid_codons_RNA;
82 :     our %n_codon_for_aa;
83 :     our %reverse_genetic_code_DNA;
84 :     our %reverse_genetic_code_RNA;
85 :     our %DNA_letter_can_be;
86 :     our %RNA_letter_can_be;
87 :     our %one_letter_to_three_letter_aa;
88 :     our %three_letter_to_one_letter_aa;
89 : efrank 1.1
90 :     require Exporter;
91 :    
92 :     our @ISA = qw(Exporter);
93 :     our @EXPORT = qw(
94 :     read_fasta_seqs
95 :     read_next_fasta_seq
96 :     parse_fasta_title
97 :     split_fasta_title
98 :     print_seq_list_as_fasta
99 :     print_seq_as_fasta
100 :     print_gb_locus
101 :    
102 :     index_seq_list
103 :     seq_entry_by_id
104 :     seq_desc_by_id
105 :     seq_data_by_id
106 :    
107 :     subseq_DNA_entry
108 :     subseq_RNA_entry
109 :     complement_DNA_entry
110 :     complement_RNA_entry
111 :     complement_DNA_seq
112 :     complement_RNA_seq
113 :     to_DNA_seq
114 :     to_RNA_seq
115 : golsen 1.2 pack_seq
116 : efrank 1.1 clean_ae_sequence
117 :    
118 :     translate_seq
119 :     translate_codon
120 :     translate_seq_with_user_code
121 :    
122 :     read_intervals
123 : golsen 1.2 standardize_intervals
124 : efrank 1.1 join_intervals
125 : golsen 1.2 locations_2_intervals
126 :     read_oriented_intervals
127 :     reverse_intervals
128 : efrank 1.1 );
129 :    
130 : golsen 1.2 our @EXPORT_OK = qw(
131 :     @aa_1_letter_order
132 :     @aa_3_letter_order
133 :     @aa_n_codon_order
134 :     %genetic_code
135 :     %genetic_code_with_U
136 :     %amino_acid_codons_DNA
137 :     %amino_acid_codons_RNA
138 :     %n_codon_for_aa
139 :     %reverse_genetic_code_DNA
140 :     %reverse_genetic_code_RNA
141 :     %DNA_letter_can_be
142 :     %RNA_letter_can_be
143 :     %one_letter_to_three_letter_aa
144 :     %three_letter_to_one_letter_aa
145 :     );
146 : efrank 1.1
147 :    
148 :     #-----------------------------------------------------------------------------
149 :     # Read fasta sequences from a file.
150 :     # Save the contents in a list of refs to arrays: (id, description, seq)
151 :     #
152 :     # @seqs = read_fasta_seqs(*FILEHANDLE)
153 :     #-----------------------------------------------------------------------------
154 :     sub read_fasta_seqs {
155 :     my $fh = shift;
156 :     wantarray || die "read_fasta_seqs requires list context\n";
157 :    
158 :     my @seqs = ();
159 :     my ($id, $desc, $seq) = ("", "", "");
160 :    
161 :     while (<$fh>) {
162 :     chomp;
163 :     if (/^>\s*(\S+)(\s+(.*))?$/) { # new id
164 :     if ($id && $seq) { push @seqs, [ $id, $desc, $seq ] }
165 :     ($id, $desc, $seq) = ($1, $3 ? $3 : "", "");
166 :     }
167 :     else {
168 : golsen 1.2 tr/ 0-9//d;
169 : efrank 1.1 $seq .= $_ ;
170 :     }
171 :     }
172 :    
173 :     if ($id && $seq) { push @seqs, [ $id, $desc, $seq ] }
174 :     return @seqs;
175 :     }
176 :    
177 :    
178 :     #-----------------------------------------------------------------------------
179 :     # Read one fasta sequence at a time from a file.
180 :     # Return the contents as an array: (id, description, seq)
181 :     #
182 :     # @seq_entry = read_next_fasta_seq(*FILEHANDLE)
183 :     #-----------------------------------------------------------------------------
184 :     # Reading always overshoots, so save next id and description
185 :    
186 : golsen 1.2 { # Use bare block to scope the header hash
187 :    
188 :     my %next_header;
189 :    
190 :     sub read_next_fasta_seq {
191 :     wantarray || die "read_next_fasta_seq requires list context\n";
192 : efrank 1.1
193 : golsen 1.2 my $fh = shift;
194 :     my ( $id, $desc );
195 : efrank 1.1
196 : golsen 1.2 if ( defined( $next_header{$fh} ) ) {
197 :     ( $id, $desc ) = parse_fasta_title( $next_header{$fh} );
198 : efrank 1.1 }
199 :     else {
200 : golsen 1.2 $next_header{$fh} = "";
201 :     ( $id, $desc ) = ( undef, "" );
202 :     }
203 :     my $seq = "";
204 :    
205 :     while ( <$fh> ) {
206 :     chomp;
207 :     if ( /^>/ ) { # new id
208 :     $next_header{$fh} = $_;
209 :     if ( defined($id) && $seq ) { return ($id, $desc, $seq) }
210 :     ( $id, $desc ) = parse_fasta_title( $next_header{$fh} );
211 :     $seq = "";
212 :     }
213 :     else {
214 :     tr/ 0-9//d;
215 :     $seq .= $_ ;
216 :     }
217 : efrank 1.1 }
218 :    
219 : golsen 1.2 # Done with file, delete "next header"
220 : efrank 1.1
221 : golsen 1.2 delete $next_header{$fh};
222 :     return (defined($id) && $seq) ? ($id, $desc, $seq) : () ;
223 :     }
224 : efrank 1.1 }
225 :    
226 :    
227 :     #-----------------------------------------------------------------------------
228 :     # Parse a fasta file header to id and definition parts
229 :     #
230 :     # ($id, $def) = parse_fasta_title( $title )
231 :     # ($id, $def) = split_fasta_title( $title )
232 :     #-----------------------------------------------------------------------------
233 :     sub parse_fasta_title {
234 :     my $title = shift;
235 :     chomp;
236 :     if ($title =~ /^>?\s*(\S+)(:?\s+(.*\S)\s*)?$/) {
237 :     return ($1, $3 ? $3 : "");
238 :     }
239 : golsen 1.2 elsif ($title =~ /^>/) {
240 :     return ("", "");
241 :     }
242 : efrank 1.1 else {
243 : golsen 1.2 return (undef, "");
244 : efrank 1.1 }
245 :     }
246 :    
247 :     sub split_fasta_title {
248 :     parse_fasta_title ( shift );
249 :     }
250 :    
251 :    
252 :     #-----------------------------------------------------------------------------
253 :     # Print list of sequence entries in fasta format
254 :     #
255 :     # print_seq_list_as_fasta(*FILEHANDLE, @seq_entry_list);
256 :     #-----------------------------------------------------------------------------
257 :     sub print_seq_list_as_fasta {
258 :     my $fh = shift;
259 :     my @seq_list = @_;
260 :    
261 :     foreach my $seq_ptr (@seq_list) {
262 :     print_seq_as_fasta($fh, @$seq_ptr);
263 :     }
264 :     }
265 :    
266 :    
267 :     #-----------------------------------------------------------------------------
268 :     # Print one sequence in fasta format to an open file
269 :     #
270 :     # print_seq_as_fasta(*FILEHANDLE, $id, $desc, $seq);
271 :     # print_seq_as_fasta(*FILEHANDLE, @seq_entry);
272 :     #-----------------------------------------------------------------------------
273 :     sub print_seq_as_fasta {
274 :     my $fh = shift;
275 :     my ($id, $desc, $seq) = @_;
276 :    
277 :     printf $fh ($desc) ? ">$id $desc\n" : ">$id\n";
278 :     my $len = length($seq);
279 :     for (my $i = 0; $i < $len; $i += 60) {
280 :     print $fh substr($seq, $i, 60) . "\n";
281 :     }
282 :     }
283 :    
284 :    
285 :     #-----------------------------------------------------------------------------
286 :     # Print one sequence in GenBank flat file format:
287 :     #
288 :     # print_gb_locus( *FILEHANDLE, $locus, $def, $accession, $seq )
289 :     #-----------------------------------------------------------------------------
290 :     sub print_gb_locus {
291 :     my ($fh, $loc, $def, $acc, $seq) = @_;
292 :     my ($len, $i, $imax);
293 :     my $istep = 10;
294 :    
295 :     $len = length($seq);
296 :     printf $fh "LOCUS %-10s%7d bp\n", substr($loc,0,10), $len;
297 :     print $fh "DEFINITION " . substr(wrap_text($def,80,12), 12) . "\n";
298 :     if ($acc) { print $fh "ACCESSION $acc\n" }
299 :     print $fh "ORIGIN\n";
300 :    
301 :     for ($i = 1; $i <= $len; ) {
302 :     printf $fh "%9d", $i;
303 :     $imax = $i + 59; if ($imax > $len) { $imax = $len }
304 :     for ( ; $i <= $imax; $i += $istep) {
305 :     print $fh " " . substr($seq, $i-1, $istep);
306 :     }
307 :     print $fh "\n";
308 :     }
309 :     print $fh "//\n";
310 :     }
311 :    
312 :    
313 :     #-----------------------------------------------------------------------------
314 :     # Build an index from seq_id to pointer to sequence entry: (id, desc, seq)
315 :     #
316 : golsen 1.2 # my \%seq_ind = index_seq_list(@seq_list);
317 : efrank 1.1 #
318 :     # Usage example:
319 :     #
320 :     # my @seq_list = read_fasta_seqs(*STDIN); # list of pointers to entries
321 : golsen 1.2 # my \%seq_ind = index_seq_list(@seq_list); # hash from names to pointers
322 : efrank 1.1 # my @chosen_seq = @{%seq_ind{"contig1"}}; # extract one entry
323 :     #
324 :     #-----------------------------------------------------------------------------
325 :     sub index_seq_list {
326 :     my %seq_index = map { @{$_}[0] => $_ } @_;
327 :     return \%seq_index;
328 :     }
329 :    
330 :    
331 :     #-----------------------------------------------------------------------------
332 :     # Three routines to access all or part of sequence entry by id:
333 :     #
334 :     # my @seq_entry = seq_entry_by_id( \%seq_index, $seq_id );
335 :     # my $seq_desc = seq_desc_by_id( \%seq_index, $seq_id );
336 :     # my $seq = seq_data_by_id( \%seq_index, $seq_id );
337 :     #
338 :     #-----------------------------------------------------------------------------
339 :     sub seq_entry_by_id {
340 :     (my $ind_ref = shift) || die "No index supplied to seq_entry_by_id\n";
341 :     (my $id = shift) || die "No id supplied to seq_entry_by_id\n";
342 :     wantarray || die "entry_by_id requires list context\n";
343 :     return @{ $ind_ref->{$id} };
344 :     }
345 :    
346 :    
347 :     sub seq_desc_by_id {
348 :     (my $ind_ref = shift) || die "No index supplied to seq_desc_by_id\n";
349 :     (my $id = shift) || die "No id supplied to seq_desc_by_id\n";
350 :     return ${ $ind_ref->{$id} }[1];
351 :     }
352 :    
353 :    
354 :     sub seq_data_by_id {
355 :     (my $ind_ref = shift) || die "No index supplied to seq_data_by_id\n";
356 :     (my $id = shift) || die "No id supplied to seq_data_by_id\n";
357 :     return ${ $ind_ref->{$id} }[2];
358 :     }
359 :    
360 :    
361 :     #-----------------------------------------------------------------------------
362 :     # Some simple sequence manipulations:
363 :     #
364 :     # @entry = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );
365 :     # @entry = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );
366 :     # @entry = complement_DNA_entry( @seq_entry [, $fix_id] );
367 :     # @entry = complement_RNA_entry( @seq_entry [, $fix_id] );
368 :     # $DNAseq = complement_DNA_seq( $NA_seq );
369 :     # $RNAseq = complement_RNA_seq( $NA_seq );
370 :     # $DNAseq = to_DNA_seq( $NA_seq );
371 :     # $RNAseq = to_RNA_seq( $NA_seq );
372 :     #
373 :     #-----------------------------------------------------------------------------
374 :    
375 :     sub subseq_DNA_entry {
376 :     my ($id, $desc, @rest) = @_;
377 :     wantarray || die "subseq_DNA_entry requires array context\n";
378 :    
379 :     my $seq;
380 :     ($id, $seq) = subseq_nt(1, $id, @rest); # 1 is for DNA, not RNA
381 :     return ($id, $desc, $seq);
382 :     }
383 :    
384 :    
385 :     sub subseq_RNA_entry {
386 :     my ($id, $desc, @rest) = @_;
387 :     wantarray || die "subseq_RNA_entry requires array context\n";
388 :    
389 :     my $seq;
390 :     ($id, $seq) = subseq_nt(0, $id, @rest); # 0 is for not DNA, i.e., RNA
391 :     return ($id, $desc, $seq);
392 :     }
393 :    
394 :    
395 :     sub subseq_nt {
396 :     my ($DNA, $id, $seq, $from, $to, $fix_id) = @_;
397 :     $fix_id ||= 0; # fix undef value
398 :    
399 :     my $len = length($seq);
400 :     if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
401 :     if (! $to || ( $to eq '$' ) || ( $to eq "" ) ) { $to = $len }
402 :    
403 :     my $left = ( $from < $to ) ? $from : $to;
404 :     my $right = ( $from < $to ) ? $to : $from;
405 :     if ( ( $right < 1 ) || ( $left > $len ) ) { return ($id, "") }
406 :     if ( $right > $len ) { $right = $len }
407 :     if ( $left < 1 ) { $left = 1 }
408 :    
409 :     $seq = substr($seq, $left-1, $right-$left+1);
410 :     if ( $from > $to ) {
411 :     $seq = reverse $seq;
412 :     if ( $DNA ) {
413 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
414 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
415 :     }
416 :     else {
417 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
418 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
419 :     }
420 :     }
421 :    
422 :     if ( $fix_id ) {
423 : golsen 1.2 if ( ( $id =~ s/_(\d+)_(\d+)$// )
424 : efrank 1.1 && ( abs($2-$1)+1 == $len ) ) {
425 :     if ( $1 <= $2 ) { $from += $1 - 1; $to += $1 - 1 }
426 :     else { $from = $1 + 1 - $from; $to = $1 + 1 - $to }
427 :     }
428 :     $id .= "_" . $from . "_" . $to;
429 :     }
430 :    
431 :     return ($id, $seq);
432 :     }
433 :    
434 :    
435 :     sub complement_DNA_entry {
436 :     my ($id, $desc, $seq, $fix_id) = @_;
437 :     $fix_id ||= 0; # fix undef values
438 :    
439 :     wantarray || die "complement_DNA_entry requires array context\n";
440 :     $seq = reverse $seq;
441 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
442 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
443 :     if ($fix_id) {
444 : golsen 1.2 if ($id =~ s/_(\d+)_(\d+)$//) {
445 : efrank 1.1 $id .= "_" . $2 . "_" . $1;
446 :     }
447 :     else {
448 :     $id .= "_" . length($seq) . "_1";
449 :     }
450 :     }
451 :    
452 :     return ($id, $desc, $seq);
453 :     }
454 :    
455 :    
456 :     sub complement_RNA_entry {
457 :     my ($id, $desc, $seq, $fix_id) = @_;
458 :     $fix_id ||= 0; # fix undef values
459 :    
460 :     wantarray || die "complement_DNA_entry requires array context\n";
461 :     $seq = reverse $seq;
462 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
463 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
464 :     if ($fix_id) {
465 : golsen 1.2 if ($id =~ s/_(\d+)_(\d+)$//) {
466 : efrank 1.1 $id .= "_" . $2 . "_" . $1;
467 :     }
468 :     else {
469 :     $id .= "_" . length($seq) . "_1";
470 :     }
471 :     }
472 :    
473 :     return ($id, $desc, $seq);
474 :     }
475 :    
476 :    
477 :     sub complement_DNA_seq {
478 :     my $seq = reverse shift;
479 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
480 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
481 :     return $seq;
482 :     }
483 :    
484 :    
485 :     sub complement_RNA_seq {
486 :     my $seq = reverse shift;
487 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
488 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
489 :     return $seq;
490 :     }
491 :    
492 :    
493 :     sub to_DNA_seq {
494 :     my $seq = shift;
495 :     $seq =~ tr/Uu/Tt/;
496 :     return $seq;
497 :     }
498 :    
499 :    
500 :     sub to_RNA_seq {
501 :     my $seq = shift;
502 :     $seq =~ tr/Tt/Uu/;
503 :     return $seq;
504 :     }
505 :    
506 :    
507 : golsen 1.2 sub pack_seq {
508 :     my $seq = shift;
509 :     $seq =~ tr/A-Za-z//cd;
510 :     return $seq;
511 :     }
512 :    
513 :    
514 : efrank 1.1 sub clean_ae_sequence {
515 :     $_ = shift;
516 :     $_ = to7bit($_);
517 :     s/[+]/1/g;
518 :     s/[^0-9A-IK-NP-Za-ik-np-z~.-]/-/g;
519 :     return $_;
520 :     }
521 :    
522 :    
523 :     sub to7bit {
524 :     $_ = shift;
525 :     my ($o, $c);
526 :     while (/\\([0-3][0-7][0-7])/) {
527 :     $o = oct($1) % 128;
528 :     $c = sprintf("%c", $o);
529 :     s/\\$1/$c/g;
530 :     }
531 :     return $_;
532 :     }
533 :    
534 :    
535 :     sub to8bit {
536 :     $_ = shift;
537 :     my ($o, $c);
538 :     while (/\\([0-3][0-7][0-7])/) {
539 :     $o = oct($1);
540 :     $c = sprintf("%c", $o);
541 :     s/\\$1/$c/g;
542 :     }
543 :     return $_;
544 :     }
545 :    
546 :    
547 :    
548 :     #-----------------------------------------------------------------------------
549 :     # Translate nucleotides to one letter protein:
550 :     #
551 :     # $seq = translate_seq( $seq [, $met_start] )
552 :     # $aa = translate_codon( $triplet );
553 :     # $aa = translate_uc_DNA_codon( $upcase_DNA_triplet );
554 :     #
555 :     # User-supplied genetic code must be upper case index and match the
556 :     # DNA versus RNA type of sequence
557 :     #
558 :     # $seq = translate_seq_with_user_code( $seq, $gen_code_hash [, $met_start] )
559 :     #
560 :     #-----------------------------------------------------------------------------
561 :    
562 : 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
563 :     @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
564 :     @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 );
565 :    
566 :     %genetic_code = (
567 :    
568 :     # DNA version
569 :    
570 : efrank 1.1 TTT => "F", TCT => "S", TAT => "Y", TGT => "C",
571 :     TTC => "F", TCC => "S", TAC => "Y", TGC => "C",
572 :     TTA => "L", TCA => "S", TAA => "*", TGA => "*",
573 :     TTG => "L", TCG => "S", TAG => "*", TGG => "W",
574 :     CTT => "L", CCT => "P", CAT => "H", CGT => "R",
575 :     CTC => "L", CCC => "P", CAC => "H", CGC => "R",
576 :     CTA => "L", CCA => "P", CAA => "Q", CGA => "R",
577 :     CTG => "L", CCG => "P", CAG => "Q", CGG => "R",
578 :     ATT => "I", ACT => "T", AAT => "N", AGT => "S",
579 :     ATC => "I", ACC => "T", AAC => "N", AGC => "S",
580 :     ATA => "I", ACA => "T", AAA => "K", AGA => "R",
581 :     ATG => "M", ACG => "T", AAG => "K", AGG => "R",
582 :     GTT => "V", GCT => "A", GAT => "D", GGT => "G",
583 :     GTC => "V", GCC => "A", GAC => "D", GGC => "G",
584 :     GTA => "V", GCA => "A", GAA => "E", GGA => "G",
585 :     GTG => "V", GCG => "A", GAG => "E", GGG => "G",
586 :    
587 : golsen 1.2 # RNA suppliment
588 :    
589 :     UUU => "F", UCU => "S", UAU => "Y", UGU => "C",
590 :     UUC => "F", UCC => "S", UAC => "Y", UGC => "C",
591 :     UUA => "L", UCA => "S", UAA => "*", UGA => "*",
592 :     UUG => "L", UCG => "S", UAG => "*", UGG => "W",
593 :     CUU => "L", CCU => "P", CAU => "H", CGU => "R",
594 :     CUC => "L",
595 :     CUA => "L",
596 :     CUG => "L",
597 :     AUU => "I", ACU => "T", AAU => "N", AGU => "S",
598 :     AUC => "I",
599 :     AUA => "I",
600 :     AUG => "M",
601 :     GUU => "V", GCU => "A", GAU => "D", GGU => "G",
602 :     GUC => "V",
603 :     GUA => "V",
604 :     GUG => "V",
605 :    
606 : efrank 1.1 # The following ambiguous encodings are not necessary, but
607 : golsen 1.2 # speed up the processing of some ambiguous triplets:
608 : efrank 1.1
609 :     TTY => "F", TCY => "S", TAY => "Y", TGY => "C",
610 :     TTR => "L", TCR => "S", TAR => "*",
611 :     TCN => "S",
612 :     CTY => "L", CCY => "P", CAY => "H", CGY => "R",
613 :     CTR => "L", CCR => "P", CAR => "Q", CGR => "R",
614 :     CTN => "L", CCN => "P", CGN => "R",
615 :     ATY => "I", ACY => "T", AAY => "N", AGY => "S",
616 :     ACR => "T", AAR => "K", AGR => "R",
617 :     ACN => "T",
618 :     GTY => "V", GCY => "A", GAY => "D", GGY => "G",
619 :     GTR => "V", GCR => "A", GAR => "E", GGR => "G",
620 : golsen 1.2 GTN => "V", GCN => "A", GGN => "G",
621 :    
622 :     UUY => "F", UCY => "S", UAY => "Y", UGY => "C",
623 :     UUR => "L", UCR => "S", UAR => "*",
624 :     UCN => "S",
625 :     CUY => "L",
626 :     CUR => "L",
627 :     CUN => "L",
628 :     AUY => "I",
629 :     GUY => "V",
630 :     GUR => "V",
631 :     GUN => "V"
632 :     );
633 :    
634 :    
635 :     # Add lower case by construction:
636 :    
637 :     foreach ( keys %genetic_code ) {
638 :     $genetic_code{ lc $_ } = lc $genetic_code{ $_ }
639 :     }
640 :    
641 :    
642 :     # Construct the genetic code with selanocysteine by difference:
643 :    
644 :     %genetic_code_with_U = map { $_ => $genetic_code{ $_ } } keys %genetic_code;
645 :     $genetic_code_with_U{ TGA } = "U";
646 :     $genetic_code_with_U{ tga } = "u";
647 :     $genetic_code_with_U{ UGA } = "U";
648 :     $genetic_code_with_U{ uga } = "u";
649 :    
650 :    
651 :     %amino_acid_codons_DNA = (
652 :     L => [ qw( TTA TTG CTA CTG CTT CTC ) ],
653 :     R => [ qw( AGA AGG CGA CGG CGT CGC ) ],
654 :     S => [ qw( AGT AGC TCA TCG TCT TCC ) ],
655 :     A => [ qw( GCA GCG GCT GCC ) ],
656 :     G => [ qw( GGA GGG GGT GGC ) ],
657 :     P => [ qw( CCA CCG CCT CCC ) ],
658 :     T => [ qw( ACA ACG ACT ACC ) ],
659 :     V => [ qw( GTA GTG GTT GTC ) ],
660 :     I => [ qw( ATA ATT ATC ) ],
661 :     C => [ qw( TGT TGC ) ],
662 :     D => [ qw( GAT GAC ) ],
663 :     E => [ qw( GAA GAG ) ],
664 :     F => [ qw( TTT TTC ) ],
665 :     H => [ qw( CAT CAC ) ],
666 :     K => [ qw( AAA AAG ) ],
667 :     N => [ qw( AAT AAC ) ],
668 :     Q => [ qw( CAA CAG ) ],
669 :     Y => [ qw( TAT TAC ) ],
670 :     M => [ qw( ATG ) ],
671 :     U => [ qw( TGA ) ],
672 :     W => [ qw( TGG ) ],
673 :    
674 :     l => [ qw( tta ttg cta ctg ctt ctc ) ],
675 :     r => [ qw( aga agg cga cgg cgt cgc ) ],
676 :     s => [ qw( agt agc tca tcg tct tcc ) ],
677 :     a => [ qw( gca gcg gct gcc ) ],
678 :     g => [ qw( gga ggg ggt ggc ) ],
679 :     p => [ qw( cca ccg cct ccc ) ],
680 :     t => [ qw( aca acg act acc ) ],
681 :     v => [ qw( gta gtg gtt gtc ) ],
682 :     i => [ qw( ata att atc ) ],
683 :     c => [ qw( tgt tgc ) ],
684 :     d => [ qw( gat gac ) ],
685 :     e => [ qw( gaa gag ) ],
686 :     f => [ qw( ttt ttc ) ],
687 :     h => [ qw( cat cac ) ],
688 :     k => [ qw( aaa aag ) ],
689 :     n => [ qw( aat aac ) ],
690 :     q => [ qw( caa cag ) ],
691 :     y => [ qw( tat tac ) ],
692 :     m => [ qw( atg ) ],
693 :     u => [ qw( tga ) ],
694 :     w => [ qw( tgg ) ],
695 :    
696 :     '*' => [ qw( TAA TAG TGA ) ]
697 : efrank 1.1 );
698 :    
699 :    
700 : golsen 1.2
701 :     %amino_acid_codons_RNA = (
702 :     L => [ qw( UUA UUG CUA CUG CUU CUC ) ],
703 :     R => [ qw( AGA AGG CGA CGG CGU CGC ) ],
704 :     S => [ qw( AGU AGC UCA UCG UCU UCC ) ],
705 :     A => [ qw( GCA GCG GCU GCC ) ],
706 :     G => [ qw( GGA GGG GGU GGC ) ],
707 :     P => [ qw( CCA CCG CCU CCC ) ],
708 :     T => [ qw( ACA ACG ACU ACC ) ],
709 :     V => [ qw( GUA GUG GUU GUC ) ],
710 :     B => [ qw( GAU GAC AAU AAC ) ],
711 :     Z => [ qw( GAA GAG CAA CAG ) ],
712 :     I => [ qw( AUA AUU AUC ) ],
713 :     C => [ qw( UGU UGC ) ],
714 :     D => [ qw( GAU GAC ) ],
715 :     E => [ qw( GAA GAG ) ],
716 :     F => [ qw( UUU UUC ) ],
717 :     H => [ qw( CAU CAC ) ],
718 :     K => [ qw( AAA AAG ) ],
719 :     N => [ qw( AAU AAC ) ],
720 :     Q => [ qw( CAA CAG ) ],
721 :     Y => [ qw( UAU UAC ) ],
722 :     M => [ qw( AUG ) ],
723 :     U => [ qw( UGA ) ],
724 :     W => [ qw( UGG ) ],
725 :    
726 :     l => [ qw( uua uug cua cug cuu cuc ) ],
727 :     r => [ qw( aga agg cga cgg cgu cgc ) ],
728 :     s => [ qw( agu agc uca ucg ucu ucc ) ],
729 :     a => [ qw( gca gcg gcu gcc ) ],
730 :     g => [ qw( gga ggg ggu ggc ) ],
731 :     p => [ qw( cca ccg ccu ccc ) ],
732 :     t => [ qw( aca acg acu acc ) ],
733 :     v => [ qw( gua gug guu guc ) ],
734 :     b => [ qw( gau gac aau aac ) ],
735 :     z => [ qw( gaa gag caa cag ) ],
736 :     i => [ qw( aua auu auc ) ],
737 :     c => [ qw( ugu ugc ) ],
738 :     d => [ qw( gau gac ) ],
739 :     e => [ qw( gaa gag ) ],
740 :     f => [ qw( uuu uuc ) ],
741 :     h => [ qw( cau cac ) ],
742 :     k => [ qw( aaa aag ) ],
743 :     n => [ qw( aau aac ) ],
744 :     q => [ qw( caa cag ) ],
745 :     y => [ qw( uau uac ) ],
746 :     m => [ qw( aug ) ],
747 :     u => [ qw( uga ) ],
748 :     w => [ qw( ugg ) ],
749 :    
750 :     '*' => [ qw( UAA UAG UGA ) ]
751 :     );
752 :    
753 :    
754 :     %n_codon_for_aa = map {
755 :     $_ => scalar @{ $amino_acid_codons_DNA{ $_ } }
756 :     } keys %amino_acid_codons_DNA;
757 :    
758 :    
759 :     %reverse_genetic_code_DNA = (
760 :     A => "GCN", a => "gcn",
761 :     C => "TGY", c => "tgy",
762 :     D => "GAY", d => "gay",
763 :     E => "GAR", e => "gar",
764 :     F => "TTY", f => "tty",
765 :     G => "GGN", g => "ggn",
766 :     H => "CAY", h => "cay",
767 :     I => "ATH", i => "ath",
768 :     K => "AAR", k => "aar",
769 :     L => "YTN", l => "ytn",
770 :     M => "ATG", m => "atg",
771 :     N => "AAY", n => "aay",
772 :     P => "CCN", p => "ccn",
773 :     Q => "CAR", q => "car",
774 :     R => "MGN", r => "mgn",
775 :     S => "WSN", s => "wsn",
776 :     T => "ACN", t => "acn",
777 :     U => "TGA", u => "tga",
778 :     V => "GTN", v => "gtn",
779 :     W => "TGG", w => "tgg",
780 :     X => "NNN", x => "nnn",
781 :     Y => "TAY", y => "tay",
782 :     '*' => "TRR"
783 :     );
784 :    
785 :     %reverse_genetic_code_RNA = (
786 :     A => "GCN", a => "gcn",
787 :     C => "UGY", c => "ugy",
788 :     D => "GAY", d => "gay",
789 :     E => "GAR", e => "gar",
790 :     F => "UUY", f => "uuy",
791 :     G => "GGN", g => "ggn",
792 :     H => "CAY", h => "cay",
793 :     I => "AUH", i => "auh",
794 :     K => "AAR", k => "aar",
795 :     L => "YUN", l => "yun",
796 :     M => "AUG", m => "aug",
797 :     N => "AAY", n => "aay",
798 :     P => "CCN", p => "ccn",
799 :     Q => "CAR", q => "car",
800 :     R => "MGN", r => "mgn",
801 :     S => "WSN", s => "wsn",
802 :     T => "ACN", t => "acn",
803 :     U => "UGA", u => "uga",
804 :     V => "GUN", v => "gun",
805 :     W => "UGG", w => "ugg",
806 :     X => "NNN", x => "nnn",
807 :     Y => "UAY", y => "uay",
808 :     '*' => "URR"
809 :     );
810 :    
811 :    
812 :     %DNA_letter_can_be = (
813 : efrank 1.1 A => ["A"], a => ["a"],
814 :     B => ["C", "G", "T"], b => ["c", "g", "t"],
815 :     C => ["C"], c => ["c"],
816 :     D => ["A", "G", "T"], d => ["a", "g", "t"],
817 :     G => ["G"], g => ["g"],
818 :     H => ["A", "C", "T"], h => ["a", "c", "t"],
819 :     K => ["G", "T"], k => ["g", "t"],
820 :     M => ["A", "C"], m => ["a", "c"],
821 :     N => ["A", "C", "G", "T"], n => ["a", "c", "g", "t"],
822 :     R => ["A", "G"], r => ["a", "g"],
823 :     S => ["C", "G"], s => ["c", "g"],
824 :     T => ["T"], t => ["t"],
825 :     U => ["T"], u => ["t"],
826 :     V => ["A", "C", "G"], v => ["a", "c", "g"],
827 :     W => ["A", "T"], w => ["a", "t"],
828 :     Y => ["C", "T"], y => ["c", "t"]
829 :     );
830 :    
831 :    
832 : golsen 1.2 %RNA_letter_can_be = (
833 : efrank 1.1 A => ["A"], a => ["a"],
834 :     B => ["C", "G", "U"], b => ["c", "g", "u"],
835 :     C => ["C"], c => ["c"],
836 :     D => ["A", "G", "U"], d => ["a", "g", "u"],
837 :     G => ["G"], g => ["g"],
838 :     H => ["A", "C", "U"], h => ["a", "c", "u"],
839 :     K => ["G", "U"], k => ["g", "u"],
840 :     M => ["A", "C"], m => ["a", "c"],
841 :     N => ["A", "C", "G", "U"], n => ["a", "c", "g", "u"],
842 :     R => ["A", "G"], r => ["a", "g"],
843 :     S => ["C", "G"], s => ["c", "g"],
844 :     T => ["U"], t => ["u"],
845 :     U => ["U"], u => ["u"],
846 :     V => ["A", "C", "G"], v => ["a", "c", "g"],
847 :     W => ["A", "U"], w => ["a", "u"],
848 :     Y => ["C", "U"], y => ["c", "u"]
849 :     );
850 :    
851 :    
852 : golsen 1.2 %one_letter_to_three_letter_aa = {
853 :     A => "Ala", a => "Ala",
854 :     B => "Asx", b => "Asx",
855 :     C => "Cys", c => "Cys",
856 :     D => "Asp", d => "Asp",
857 :     E => "Glu", e => "Glu",
858 :     F => "Phe", f => "Phe",
859 :     G => "Gly", g => "Gly",
860 :     H => "His", h => "His",
861 :     I => "Ile", i => "Ile",
862 :     K => "Lys", k => "Lys",
863 :     L => "Leu", l => "Leu",
864 :     M => "Met", m => "Met",
865 :     N => "Asn", n => "Asn",
866 :     P => "Pro", p => "Pro",
867 :     Q => "Gln", q => "Gln",
868 :     R => "Arg", r => "Arg",
869 :     S => "Ser", s => "Ser",
870 :     T => "Thr", t => "Thr",
871 :     U => "Sec", u => "Sec",
872 :     V => "Val", v => "Val",
873 :     W => "Trp", w => "Trp",
874 :     X => "Xxx", x => "Xxx",
875 :     Y => "Tyr", y => "Tyr",
876 :     Z => "Glx", z => "Glx",
877 :     '*' => "***"
878 :     };
879 :    
880 :    
881 :     %three_letter_to_one_letter_aa = (
882 :     ALA => "A", Ala => "A", ala => "a",
883 :     ARG => "R", Arg => "R", arg => "r",
884 :     ASN => "N", Asn => "N", asn => "n",
885 :     ASP => "D", Asp => "D", asp => "d",
886 :     ASX => "B", Asx => "B", asx => "b",
887 :     CYS => "C", Cys => "C", cys => "c",
888 :     GLN => "Q", Gln => "Q", gln => "q",
889 :     GLU => "E", Glu => "E", glu => "e",
890 :     GLX => "Z", Glx => "Z", glx => "z",
891 :     GLY => "G", Gly => "G", gly => "g",
892 :     HIS => "H", His => "H", his => "h",
893 :     ILE => "I", Ile => "I", ile => "i",
894 :     LEU => "L", Leu => "L", leu => "l",
895 :     LYS => "K", Lys => "K", lys => "k",
896 :     MET => "M", Met => "M", met => "m",
897 :     PHE => "F", Phe => "F", phe => "f",
898 :     PRO => "P", Pro => "P", pro => "p",
899 :     SEC => "U", Sec => "U", sec => "u",
900 :     SER => "S", Ser => "S", ser => "s",
901 :     THR => "T", Thr => "T", thr => "t",
902 :     TRP => "W", Trp => "W", trp => "w",
903 :     TYR => "Y", Tyr => "Y", tyr => "y",
904 :     VAL => "V", Val => "V", val => "v",
905 :     XAA => "X", Xaa => "X", xaa => "x",
906 :     XXX => "X", Xxx => "X", xxx => "x",
907 :     '***' => "*"
908 : efrank 1.1 );
909 :    
910 :    
911 :     #-----------------------------------------------------------------------------
912 :     # Translate nucleotides to one letter protein:
913 :     #
914 :     # $seq = translate_seq( $seq [, $met_start] )
915 :     #
916 :     #-----------------------------------------------------------------------------
917 :    
918 :     sub translate_seq {
919 :     my $seq = uc shift;
920 :     $seq =~ tr/UX/TN/; # make it DNA, and allow X
921 :     $seq =~ tr/-//d; # remove gaps
922 :    
923 :     my $met = shift || 0; # a second argument that is true
924 :     # forces first amino acid to be Met
925 :     # (note: undef is false)
926 :    
927 :     my $imax = length($seq) - 2; # will try to translate 2 nucleotides!
928 : golsen 1.2 my $pep = ( ($met && ($imax >= 0)) ? "M" : "" );
929 : efrank 1.1 for (my $i = $met ? 3 : 0; $i <= $imax; $i += 3) {
930 :     $pep .= translate_uc_DNA_codon( substr($seq,$i,3) );
931 :     }
932 :    
933 :     return $pep;
934 :     }
935 :    
936 :    
937 :     #-----------------------------------------------------------------------------
938 :     # Translate a single triplet with "universal" genetic code
939 :     # Uppercase and DNA are performed, then translate_uc_DNA_codon
940 :     # is called.
941 :     #
942 :     # $aa = translate_codon( $triplet )
943 :     #
944 :     #-----------------------------------------------------------------------------
945 :    
946 :     sub translate_codon {
947 :     my $codon = uc shift; # Make it uppercase
948 :     $codon =~ tr/UX/TN/; # Make it DNA, and allow X
949 :     return translate_uc_DNA_codon($codon);
950 :     }
951 :    
952 :    
953 :     #-----------------------------------------------------------------------------
954 :     # Translate a single triplet with "universal" genetic code
955 :     # Uppercase and DNA assumed
956 :     # Intended for private use by translate_codon and translate_seq
957 :     #
958 :     # $aa = translate_uc_DNA_codon( $triplet )
959 :     #
960 :     #-----------------------------------------------------------------------------
961 :    
962 :     sub translate_uc_DNA_codon {
963 :     my $codon = shift;
964 :     my $aa;
965 :    
966 :     # Try a simple lookup:
967 :    
968 :     if ( $aa = $genetic_code{ $codon } ) { return $aa }
969 :    
970 :     # With the expanded code defined above, this catches simple N, R
971 :     # and Y ambiguities in the third position. Other codes like
972 :     # GG[KMSWBDHV], or even GG, might be unambiguously translated by
973 :     # converting the last position to N and seeing if this is in the
974 :     # (expanded) code table:
975 :    
976 :     if ( $aa = $genetic_code{ substr($codon,0,2) . "N" } ) { return $aa }
977 :    
978 :     # Test that codon is valid and might have unambiguous aa:
979 :    
980 :     if ( $codon !~ m/^[ACGTMY][ACGT][ACGTKMRSWYBDHVN]$/ ) { return "X" }
981 :     # ^^
982 :     # |+- for leucine YTR
983 :     # +-- for arginine MGR
984 :    
985 :     # Expand all ambiguous nucleotides to see if they all yield same aa.
986 :     # Loop order tries to fail quickly with first position change.
987 :    
988 :     $aa = "";
989 :     for my $n2 ( @{ $DNA_letter_can_be{ substr($codon,1,1) } } ) {
990 :     for my $n3 ( @{ $DNA_letter_can_be{ substr($codon,2,1) } } ) {
991 :     for my $n1 ( @{ $DNA_letter_can_be{ substr($codon,0,1) } } ) {
992 :     # set the first value of $aa
993 :     if ($aa eq "") { $aa = $genetic_code{ $n1 . $n2 . $n3 } }
994 :     # or break out if any other amino acid is detected
995 :     elsif ($aa ne $genetic_code{ $n1 . $n2 . $n3 } ) { return "X" }
996 :     }
997 :     }
998 :     }
999 :    
1000 :     return $aa || "X";
1001 :     }
1002 :    
1003 :    
1004 :     #-----------------------------------------------------------------------------
1005 :     # Translate with a user-supplied genetic code to translate a sequence.
1006 :     # Diagnose the use of upper versus lower, and T versus U in the supplied
1007 :     # code, and transform the supplied nucleotide sequence to match.
1008 :     #
1009 :     # translate_seq_with_user_code($seq, \%gen_code [, $start_with_met] )
1010 :     #
1011 :     #-----------------------------------------------------------------------------
1012 :    
1013 :     sub translate_seq_with_user_code {
1014 :     my $seq = shift;
1015 :     $seq =~ tr/-//d; # remove gaps *** Why?
1016 :     $seq =~ tr/Xx/Nn/; # allow X
1017 :    
1018 :     my $gc = shift; # Reference to hash of DNA alphabet code
1019 :     if (! $gc || ref($gc) ne "HASH") {
1020 :     die "translate_seq_with_user_code needs genetic code hash as secondargument.";
1021 :     }
1022 :    
1023 :     # Test the type of code supplied: uppercase versus lowercase
1024 :    
1025 :     my ($RNA_F, $DNA_F, $M, $N, $X);
1026 :    
1027 :     if ($gc->{ "AAA" }) { # Looks like uppercase code table
1028 :     $seq = uc $seq; # Uppercase sequence
1029 :     $RNA_F = "UUU"; # Uppercase RNA Phe codon
1030 :     $DNA_F = "TTT"; # Uppercase DNA Phe codon
1031 :     $M = "M"; # Uppercase initiator
1032 :     $N = "N"; # Uppercase ambiguous nuc
1033 :     $X = "X"; # Uppercase ambiguous aa
1034 :     }
1035 :     elsif ($gc->{ "aaa" }) { # Looks like lowercase code table
1036 :     $seq = lc $seq; # Lowercase sequence
1037 :     $RNA_F = "uuu"; # Lowercase RNA Phe codon
1038 :     $DNA_F = "ttt"; # Lowercase DNA Phe codon
1039 :     $M = "m"; # Lowercase initiator
1040 :     $N = "n"; # Lowercase ambiguous nuc
1041 :     $X = "x"; # Lowercase ambiguous aa
1042 :     }
1043 :     else {
1044 :     die "User-supplied genetic code does not have aaa or AAA\n";
1045 :     }
1046 :    
1047 :     # Test the type of code supplied: UUU versus TTT
1048 :    
1049 :     my ($ambigs);
1050 :    
1051 :     if ($gc->{ $RNA_F }) { # Looks like RNA code table
1052 :     $seq =~ tr/Tt/Uu/;
1053 :     $ambigs = \%RNA_letter_can_be;
1054 :     }
1055 :     elsif ($gc->{ $DNA_F }) { # Looks like DNA code table
1056 :     $seq =~ tr/Uu/Tt/;
1057 :     $ambigs = \%DNA_letter_can_be;
1058 :     }
1059 :     else {
1060 :     die "User-supplied genetic code does not have $RNA_F or $DNA_F\n";
1061 :     }
1062 :    
1063 :     my $imax = length($seq) - 2; # will try to translate 2 nucleotides!
1064 :    
1065 :     my $met = shift; # a third argument that is true
1066 :     # forces first amino acid to be Met
1067 :     # (note: undef is false)
1068 :     my $pep = ($met && ($imax >= 0)) ? $M : "";
1069 :     my $aa;
1070 :    
1071 :     for (my $i = $met ? 3 : 0; $i <= $imax; $i += 3) {
1072 :     $pep .= translate_codon_with_user_code( substr($seq,$i,3), $gc, $N, $X, $ambigs );
1073 :     }
1074 :    
1075 :     return $pep;
1076 :     }
1077 :    
1078 :    
1079 :     #-----------------------------------------------------------------------------
1080 :     # Translate with user-supplied genetic code hash. For speed, no error
1081 :     # check on the hash. Calling programs should check for the hash at a
1082 :     # higher level.
1083 :     #
1084 :     # Should only be called through translate_seq_with_user_code
1085 :     #
1086 :     # translate_codon_with_user_code( $triplet, \%code, $N, $X, $ambig_table )
1087 :     #
1088 :     # $triplet speaks for itself
1089 :     # $code ref to the hash with the codon translations
1090 :     # $N character to use for ambiguous nucleotide
1091 :     # $X character to use for ambiguous amino acid
1092 :     # $ambig_table ref to hash with lists of nucleotides for each ambig code
1093 :     #-----------------------------------------------------------------------------
1094 :    
1095 :    
1096 :     sub translate_codon_with_user_code {
1097 :     my $codon = shift;
1098 :     my $gc = shift;
1099 :     my $aa;
1100 :    
1101 :     # Try a simple lookup:
1102 :    
1103 :     if ( $aa = $gc->{ $codon } ) { return $aa }
1104 :    
1105 :     # Test that codon is valid and might have unambiguous aa:
1106 :    
1107 :     my ($N, $X, $ambigs) = @_;
1108 :     if ( $codon =~ m/^[ACGTUMY][ACGTU]$/i ) { $codon .= $N }
1109 :     if ( $codon !~ m/^[ACGTUMY][ACGTU][ACGTUKMRSWYBDHVN]$/i ) { return $X }
1110 :     # ^^
1111 :     # |+- for leucine YTR
1112 :     # +-- for arginine MGR
1113 :    
1114 :     # Expand all ambiguous nucleotides to see if they all yield same aa.
1115 :     # Loop order tries to fail quickly with first position change.
1116 :    
1117 :     $aa = "";
1118 :     for my $n2 ( @{ $ambigs->{ substr($codon,1,1) } } ) {
1119 :     for my $n3 ( @{ $ambigs->{ substr($codon,2,1) } } ) {
1120 :     for my $n1 ( @{ $ambigs->{ substr($codon,0,1) } } ) {
1121 :     # set the first value of $aa
1122 :     if ($aa eq "") { $aa = $gc->{ $n1 . $n2 . $n3 } }
1123 :     # break out if any other amino acid is detected
1124 :     elsif ($aa ne $gc->{ $n1 . $n2 . $n3 } ) { return "X" }
1125 :     }
1126 :     }
1127 :     }
1128 :    
1129 :     return $aa || $X;
1130 :     }
1131 :    
1132 :    
1133 :     #-----------------------------------------------------------------------------
1134 :     # Read a list of intervals from a file.
1135 :     # Allow id_start_end, or id \s start \s end formats
1136 :     #
1137 : golsen 1.2 # @intervals = read_intervals( *FILEHANDLE )
1138 : efrank 1.1 #-----------------------------------------------------------------------------
1139 :     sub read_intervals {
1140 :     my $fh = shift;
1141 :     my @intervals = ();
1142 :    
1143 :     while (<$fh>) {
1144 :     chomp;
1145 :     /^(\S+)_(\d+)_(\d+)(\s.*)?$/ # id_start_end WIT2
1146 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/ # id_start-end ???
1147 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/ # id=start=end Badger
1148 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/ # id \s start \s end
1149 :     || next;
1150 :    
1151 :     # Matched a pattern. Store reference to (id, left, right):
1152 :     push @intervals, ($2 < $3) ? [ $1, $2+0, $3+0 ]
1153 :     : [ $1, $3+0, $2+0 ];
1154 :     }
1155 :     return @intervals;
1156 :     }
1157 :    
1158 :    
1159 :     #-----------------------------------------------------------------------------
1160 : golsen 1.2 # Convert a list of intervals to read [ id, left_end, right_end ].
1161 :     #
1162 :     # @intervals = standardize_intervals( @interval_refs )
1163 :     #-----------------------------------------------------------------------------
1164 :     sub standardize_intervals {
1165 :     map { ( $_->[1] < $_->[2] ) ? $_ : [ $_->[0], $_->[2], $_->[1] ] } @_;
1166 :     }
1167 :    
1168 :    
1169 :     #-----------------------------------------------------------------------------
1170 : efrank 1.1 # Take the union of a list of intervals
1171 :     #
1172 :     # @joined = join_intervals( @interval_refs )
1173 :     #-----------------------------------------------------------------------------
1174 :     sub join_intervals {
1175 :     my @ordered = sort { $a->[0] cmp $b->[0] # first by id
1176 :     || $a->[1] <=> $b->[1] # next by left end
1177 :     || $b->[2] <=> $a->[2] # finally longest first
1178 :     } @_;
1179 :    
1180 :     my @joined = ();
1181 :     my $n_int = @ordered;
1182 :    
1183 :     my ($cur_id) = "";
1184 :     my ($cur_left) = -1;
1185 :     my ($cur_right) = -1;
1186 :     my ($new_id, $new_left, $new_right);
1187 :    
1188 :     for (my $i = 0; $i < $n_int; $i++) {
1189 :     ($new_id, $new_left, $new_right) = @{$ordered[$i]}; # get the new data
1190 :    
1191 :     if ( ( $new_id ne $cur_id) # if new contig
1192 :     || ( $new_left > $cur_right + 1) # or not touching previous
1193 :     ) { # push the previous interval
1194 :     if ($cur_id) { push (@joined, [ $cur_id, $cur_left, $cur_right ]) }
1195 :     $cur_id = $new_id; # update the current interval
1196 :     $cur_left = $new_left;
1197 :     $cur_right = $new_right;
1198 :     }
1199 :    
1200 :     elsif ($new_right > $cur_right) { # extend the right end if necessary
1201 :     $cur_right = $new_right;
1202 :     }
1203 :     }
1204 :    
1205 :     if ($cur_id) { push (@joined, [$cur_id, $cur_left, $cur_right]) }
1206 :     return @joined;
1207 :     }
1208 :    
1209 : golsen 1.2
1210 :     #-----------------------------------------------------------------------------
1211 :     # Split location strings to oriented intervals.
1212 :     #
1213 :     # @intervals = locations_2_intervals( @locations )
1214 :     # $interval = locations_2_intervals( $location )
1215 :     #-----------------------------------------------------------------------------
1216 :     sub locations_2_intervals {
1217 :     my @intervals = map { /^(\S+)_(\d+)_(\d+)(\s.*)?$/
1218 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/
1219 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/
1220 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/
1221 :     ? [ $1, $2+0, $3+0 ]
1222 :     : ()
1223 :     } @_;
1224 :    
1225 :     return wantarray ? @intervals : $intervals[0];
1226 :     }
1227 :    
1228 :    
1229 :     #-----------------------------------------------------------------------------
1230 :     # Read a list of oriented intervals from a file.
1231 :     # Allow id_start_end, or id \s start \s end formats
1232 :     #
1233 :     # @intervals = read_oriented_intervals( *FILEHANDLE )
1234 :     #-----------------------------------------------------------------------------
1235 :     sub read_oriented_intervals {
1236 :     my $fh = shift;
1237 :     my @intervals = ();
1238 :    
1239 :     while (<$fh>) {
1240 :     chomp;
1241 :     /^(\S+)_(\d+)_(\d+)(\s.*)?$/ # id_start_end WIT2
1242 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/ # id_start-end ???
1243 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/ # id=start=end Badger
1244 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/ # id \s start \s end
1245 :     || next;
1246 :    
1247 :     # Matched a pattern. Store reference to (id, start, end):
1248 :     push @intervals, [ $1, $2+0, $3+0 ];
1249 :     }
1250 :     return @intervals;
1251 :     }
1252 :    
1253 :    
1254 :     #-----------------------------------------------------------------------------
1255 :     # Reverse the orientation of a list of intervals
1256 :     #
1257 :     # @reversed = reverse_intervals( @interval_refs )
1258 :     #-----------------------------------------------------------------------------
1259 :     sub reverse_intervals {
1260 :     map { [ $_->[0], $_->[2], $_->[1] ] } @_;
1261 :     }
1262 :    
1263 :    
1264 : efrank 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3