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

Annotation of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : efrank 1.1 package gjoseqlib;
2 :    
3 :     # read_fasta_seqs( *FILEHANDLE )
4 :     # read_next_fasta_seq( *FILEHANDLE )
5 :     # parse_fasta_title( $title )
6 :     # split_fasta_title( $title )
7 :     # print_seq_list_as_fasta( *FILEHANDLE, @seq_entry_list )
8 :     # print_seq_as_fasta( *FILEHANDLE, $id, $desc, $seq )
9 :     # print_gb_locus( *FILEHANDLE, $locus, $def, $accession, $seq )
10 :    
11 :     # index_seq_list( @seq_entry_list )
12 :     # seq_entry_by_id( \%seq_index, $seq_id )
13 :     # seq_desc_by_id( \%seq_index, $seq_id )
14 :     # seq_data_by_id( \%seq_index, $seq_id )
15 :    
16 :     # subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] )
17 :     # subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] )
18 :     # complement_DNA_entry( @seq_entry [, $fix_id] )
19 :     # complement_RNA_entry( @seq_entry [, $fix_id] )
20 :     # complement_DNA_seq( $NA_seq )
21 :     # complement_RNA_seq( $NA_seq )
22 :     # to_DNA_seq( $NA_seq )
23 :     # to_RNA_seq( $NA_seq )
24 :     # clean_ae_sequence( $seq )
25 :    
26 :     # translate_seq( $seq [, $met_start] )
27 :     # translate_codon( $triplet )
28 :     # translate_seq_with_user_code( $seq, \%gen_code [, $met_start] )
29 :    
30 :     # read_intervals( *FILEHANDLE )
31 :     # join_intervals( @interval_refs )
32 :    
33 :     use gjolib qw(wrap_text);
34 :    
35 :     require Exporter;
36 :    
37 :     our @ISA = qw(Exporter);
38 :     our @EXPORT = qw(
39 :     read_fasta_seqs
40 :     read_next_fasta_seq
41 :     parse_fasta_title
42 :     split_fasta_title
43 :     print_seq_list_as_fasta
44 :     print_seq_as_fasta
45 :     print_gb_locus
46 :    
47 :     index_seq_list
48 :     seq_entry_by_id
49 :     seq_desc_by_id
50 :     seq_data_by_id
51 :    
52 :     subseq_DNA_entry
53 :     subseq_RNA_entry
54 :     complement_DNA_entry
55 :     complement_RNA_entry
56 :     complement_DNA_seq
57 :     complement_RNA_seq
58 :     to_DNA_seq
59 :     to_RNA_seq
60 :     clean_ae_sequence
61 :    
62 :     translate_seq
63 :     translate_codon
64 :     translate_seq_with_user_code
65 :    
66 :     read_intervals
67 :     join_intervals
68 :     );
69 :    
70 :     use strict;
71 :    
72 :    
73 :     #-----------------------------------------------------------------------------
74 :     # Read fasta sequences from a file.
75 :     # Save the contents in a list of refs to arrays: (id, description, seq)
76 :     #
77 :     # @seqs = read_fasta_seqs(*FILEHANDLE)
78 :     #-----------------------------------------------------------------------------
79 :     sub read_fasta_seqs {
80 :     my $fh = shift;
81 :     wantarray || die "read_fasta_seqs requires list context\n";
82 :    
83 :     my @seqs = ();
84 :     my ($id, $desc, $seq) = ("", "", "");
85 :    
86 :     while (<$fh>) {
87 :     chomp;
88 :     if (/^>\s*(\S+)(\s+(.*))?$/) { # new id
89 :     if ($id && $seq) { push @seqs, [ $id, $desc, $seq ] }
90 :     ($id, $desc, $seq) = ($1, $3 ? $3 : "", "");
91 :     }
92 :     else {
93 :     tr/\t 0-9//d;
94 :     $seq .= $_ ;
95 :     }
96 :     }
97 :    
98 :     if ($id && $seq) { push @seqs, [ $id, $desc, $seq ] }
99 :     return @seqs;
100 :     }
101 :    
102 :    
103 :     #-----------------------------------------------------------------------------
104 :     # Read one fasta sequence at a time from a file.
105 :     # Return the contents as an array: (id, description, seq)
106 :     #
107 :     # @seq_entry = read_next_fasta_seq(*FILEHANDLE)
108 :     #-----------------------------------------------------------------------------
109 :     # Reading always overshoots, so save next id and description
110 :     # (these should really be hashes indexed by file handle):
111 :     #
112 :     my $next_fasta_seq_header = "";
113 :    
114 :     sub read_next_fasta_seq {
115 :     my $fh = shift;
116 :     wantarray || die "read_next_fasta_seq requires list context\n";
117 :    
118 :     my ( $id, $desc ) = parse_fasta_title( $next_fasta_seq_header );
119 :     my $seq = "";
120 :    
121 :     while ( <$fh> ) {
122 :     chomp;
123 :     if ( /^>/ ) { # new id
124 :     $next_fasta_seq_header = $_;
125 :     if ( $id && $seq ) { return ($id, $desc, $seq) }
126 :     ( $id, $desc ) = parse_fasta_title( $next_fasta_seq_header );
127 :     $seq = "";
128 :     }
129 :     else {
130 :     tr/\t 0-9//d;
131 :     $seq .= $_ ;
132 :     }
133 :     }
134 :    
135 :     # Done with file, clear "next header"
136 :    
137 :     $next_fasta_seq_header = "";
138 :     return ($id && $seq) ? ($id, $desc, $seq) : () ;
139 :     }
140 :    
141 :    
142 :     #-----------------------------------------------------------------------------
143 :     # Parse a fasta file header to id and definition parts
144 :     #
145 :     # ($id, $def) = parse_fasta_title( $title )
146 :     # ($id, $def) = split_fasta_title( $title )
147 :     #-----------------------------------------------------------------------------
148 :     sub parse_fasta_title {
149 :     my $title = shift;
150 :     chomp;
151 :     if ($title =~ /^>?\s*(\S+)(:?\s+(.*\S)\s*)?$/) {
152 :     return ($1, $3 ? $3 : "");
153 :     }
154 :     else {
155 :     return ("", "");
156 :     }
157 :     }
158 :    
159 :     sub split_fasta_title {
160 :     parse_fasta_title ( shift );
161 :     }
162 :    
163 :    
164 :     #-----------------------------------------------------------------------------
165 :     # Print list of sequence entries in fasta format
166 :     #
167 :     # print_seq_list_as_fasta(*FILEHANDLE, @seq_entry_list);
168 :     #-----------------------------------------------------------------------------
169 :     sub print_seq_list_as_fasta {
170 :     my $fh = shift;
171 :     my @seq_list = @_;
172 :    
173 :     foreach my $seq_ptr (@seq_list) {
174 :     print_seq_as_fasta($fh, @$seq_ptr);
175 :     }
176 :     }
177 :    
178 :    
179 :     #-----------------------------------------------------------------------------
180 :     # Print one sequence in fasta format to an open file
181 :     #
182 :     # print_seq_as_fasta(*FILEHANDLE, $id, $desc, $seq);
183 :     # print_seq_as_fasta(*FILEHANDLE, @seq_entry);
184 :     #-----------------------------------------------------------------------------
185 :     sub print_seq_as_fasta {
186 :     my $fh = shift;
187 :     my ($id, $desc, $seq) = @_;
188 :    
189 :     printf $fh ($desc) ? ">$id $desc\n" : ">$id\n";
190 :     my $len = length($seq);
191 :     for (my $i = 0; $i < $len; $i += 60) {
192 :     print $fh substr($seq, $i, 60) . "\n";
193 :     }
194 :     }
195 :    
196 :    
197 :     #-----------------------------------------------------------------------------
198 :     # Print one sequence in GenBank flat file format:
199 :     #
200 :     # print_gb_locus( *FILEHANDLE, $locus, $def, $accession, $seq )
201 :     #-----------------------------------------------------------------------------
202 :     sub print_gb_locus {
203 :     my ($fh, $loc, $def, $acc, $seq) = @_;
204 :     my ($len, $i, $imax);
205 :     my $istep = 10;
206 :    
207 :     $len = length($seq);
208 :     printf $fh "LOCUS %-10s%7d bp\n", substr($loc,0,10), $len;
209 :     print $fh "DEFINITION " . substr(wrap_text($def,80,12), 12) . "\n";
210 :     if ($acc) { print $fh "ACCESSION $acc\n" }
211 :     print $fh "ORIGIN\n";
212 :    
213 :     for ($i = 1; $i <= $len; ) {
214 :     printf $fh "%9d", $i;
215 :     $imax = $i + 59; if ($imax > $len) { $imax = $len }
216 :     for ( ; $i <= $imax; $i += $istep) {
217 :     print $fh " " . substr($seq, $i-1, $istep);
218 :     }
219 :     print $fh "\n";
220 :     }
221 :     print $fh "//\n";
222 :     }
223 :    
224 :    
225 :     #-----------------------------------------------------------------------------
226 :     # Build an index from seq_id to pointer to sequence entry: (id, desc, seq)
227 :     #
228 :     # my %seq_ind = index_seq_list(@seq_list);
229 :     #
230 :     # Usage example:
231 :     #
232 :     # my @seq_list = read_fasta_seqs(*STDIN); # list of pointers to entries
233 :     # my %seq_ind = index_seq_list(@seq_list); # hash from names to pointers
234 :     # my @chosen_seq = @{%seq_ind{"contig1"}}; # extract one entry
235 :     #
236 :     #-----------------------------------------------------------------------------
237 :     sub index_seq_list {
238 :     my %seq_index = map { @{$_}[0] => $_ } @_;
239 :     return \%seq_index;
240 :     }
241 :    
242 :    
243 :     #-----------------------------------------------------------------------------
244 :     # Three routines to access all or part of sequence entry by id:
245 :     #
246 :     # my @seq_entry = seq_entry_by_id( \%seq_index, $seq_id );
247 :     # my $seq_desc = seq_desc_by_id( \%seq_index, $seq_id );
248 :     # my $seq = seq_data_by_id( \%seq_index, $seq_id );
249 :     #
250 :     #-----------------------------------------------------------------------------
251 :     sub seq_entry_by_id {
252 :     (my $ind_ref = shift) || die "No index supplied to seq_entry_by_id\n";
253 :     (my $id = shift) || die "No id supplied to seq_entry_by_id\n";
254 :     wantarray || die "entry_by_id requires list context\n";
255 :     return @{ $ind_ref->{$id} };
256 :     }
257 :    
258 :    
259 :     sub seq_desc_by_id {
260 :     (my $ind_ref = shift) || die "No index supplied to seq_desc_by_id\n";
261 :     (my $id = shift) || die "No id supplied to seq_desc_by_id\n";
262 :     return ${ $ind_ref->{$id} }[1];
263 :     }
264 :    
265 :    
266 :     sub seq_data_by_id {
267 :     (my $ind_ref = shift) || die "No index supplied to seq_data_by_id\n";
268 :     (my $id = shift) || die "No id supplied to seq_data_by_id\n";
269 :     return ${ $ind_ref->{$id} }[2];
270 :     }
271 :    
272 :    
273 :     #-----------------------------------------------------------------------------
274 :     # Some simple sequence manipulations:
275 :     #
276 :     # @entry = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );
277 :     # @entry = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );
278 :     # @entry = complement_DNA_entry( @seq_entry [, $fix_id] );
279 :     # @entry = complement_RNA_entry( @seq_entry [, $fix_id] );
280 :     # $DNAseq = complement_DNA_seq( $NA_seq );
281 :     # $RNAseq = complement_RNA_seq( $NA_seq );
282 :     # $DNAseq = to_DNA_seq( $NA_seq );
283 :     # $RNAseq = to_RNA_seq( $NA_seq );
284 :     #
285 :     #-----------------------------------------------------------------------------
286 :    
287 :     sub subseq_DNA_entry {
288 :     my ($id, $desc, @rest) = @_;
289 :     wantarray || die "subseq_DNA_entry requires array context\n";
290 :    
291 :     my $seq;
292 :     ($id, $seq) = subseq_nt(1, $id, @rest); # 1 is for DNA, not RNA
293 :     return ($id, $desc, $seq);
294 :     }
295 :    
296 :    
297 :     sub subseq_RNA_entry {
298 :     my ($id, $desc, @rest) = @_;
299 :     wantarray || die "subseq_RNA_entry requires array context\n";
300 :    
301 :     my $seq;
302 :     ($id, $seq) = subseq_nt(0, $id, @rest); # 0 is for not DNA, i.e., RNA
303 :     return ($id, $desc, $seq);
304 :     }
305 :    
306 :    
307 :     sub subseq_nt {
308 :     my ($DNA, $id, $seq, $from, $to, $fix_id) = @_;
309 :     $fix_id ||= 0; # fix undef value
310 :    
311 :     my $len = length($seq);
312 :     if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
313 :     if (! $to || ( $to eq '$' ) || ( $to eq "" ) ) { $to = $len }
314 :    
315 :     my $left = ( $from < $to ) ? $from : $to;
316 :     my $right = ( $from < $to ) ? $to : $from;
317 :     if ( ( $right < 1 ) || ( $left > $len ) ) { return ($id, "") }
318 :     if ( $right > $len ) { $right = $len }
319 :     if ( $left < 1 ) { $left = 1 }
320 :    
321 :     $seq = substr($seq, $left-1, $right-$left+1);
322 :     if ( $from > $to ) {
323 :     $seq = reverse $seq;
324 :     if ( $DNA ) {
325 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
326 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
327 :     }
328 :     else {
329 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
330 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
331 :     }
332 :     }
333 :    
334 :     if ( $fix_id ) {
335 :     if ( ( $id =~ s/_([0-9]+)_([0-9]+)$// )
336 :     && ( abs($2-$1)+1 == $len ) ) {
337 :     if ( $1 <= $2 ) { $from += $1 - 1; $to += $1 - 1 }
338 :     else { $from = $1 + 1 - $from; $to = $1 + 1 - $to }
339 :     }
340 :     $id .= "_" . $from . "_" . $to;
341 :     }
342 :    
343 :     return ($id, $seq);
344 :     }
345 :    
346 :    
347 :     sub complement_DNA_entry {
348 :     my ($id, $desc, $seq, $fix_id) = @_;
349 :     $fix_id ||= 0; # fix undef values
350 :    
351 :     wantarray || die "complement_DNA_entry requires array context\n";
352 :     $seq = reverse $seq;
353 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
354 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
355 :     if ($fix_id) {
356 :     if ($id =~ s/_([0-9]+)_([0-9]+)$//) {
357 :     $id .= "_" . $2 . "_" . $1;
358 :     }
359 :     else {
360 :     $id .= "_" . length($seq) . "_1";
361 :     }
362 :     }
363 :    
364 :     return ($id, $desc, $seq);
365 :     }
366 :    
367 :    
368 :     sub complement_RNA_entry {
369 :     my ($id, $desc, $seq, $fix_id) = @_;
370 :     $fix_id ||= 0; # fix undef values
371 :    
372 :     wantarray || die "complement_DNA_entry requires array context\n";
373 :     $seq = reverse $seq;
374 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
375 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
376 :     if ($fix_id) {
377 :     if ($id =~ s/_([0-9]+)_([0-9]+)$//) {
378 :     $id .= "_" . $2 . "_" . $1;
379 :     }
380 :     else {
381 :     $id .= "_" . length($seq) . "_1";
382 :     }
383 :     }
384 :    
385 :     return ($id, $desc, $seq);
386 :     }
387 :    
388 :    
389 :     sub complement_DNA_seq {
390 :     my $seq = reverse shift;
391 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
392 :     [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
393 :     return $seq;
394 :     }
395 :    
396 :    
397 :     sub complement_RNA_seq {
398 :     my $seq = reverse shift;
399 :     $seq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
400 :     [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
401 :     return $seq;
402 :     }
403 :    
404 :    
405 :     sub to_DNA_seq {
406 :     my $seq = shift;
407 :     $seq =~ tr/Uu/Tt/;
408 :     return $seq;
409 :     }
410 :    
411 :    
412 :     sub to_RNA_seq {
413 :     my $seq = shift;
414 :     $seq =~ tr/Tt/Uu/;
415 :     return $seq;
416 :     }
417 :    
418 :    
419 :     sub clean_ae_sequence {
420 :     $_ = shift;
421 :     $_ = to7bit($_);
422 :     s/[+]/1/g;
423 :     s/[^0-9A-IK-NP-Za-ik-np-z~.-]/-/g;
424 :     return $_;
425 :     }
426 :    
427 :    
428 :     sub to7bit {
429 :     $_ = shift;
430 :     my ($o, $c);
431 :     while (/\\([0-3][0-7][0-7])/) {
432 :     $o = oct($1) % 128;
433 :     $c = sprintf("%c", $o);
434 :     s/\\$1/$c/g;
435 :     }
436 :     return $_;
437 :     }
438 :    
439 :    
440 :     sub to8bit {
441 :     $_ = shift;
442 :     my ($o, $c);
443 :     while (/\\([0-3][0-7][0-7])/) {
444 :     $o = oct($1);
445 :     $c = sprintf("%c", $o);
446 :     s/\\$1/$c/g;
447 :     }
448 :     return $_;
449 :     }
450 :    
451 :    
452 :    
453 :     #-----------------------------------------------------------------------------
454 :     # Translate nucleotides to one letter protein:
455 :     #
456 :     # $seq = translate_seq( $seq [, $met_start] )
457 :     # $aa = translate_codon( $triplet );
458 :     # $aa = translate_uc_DNA_codon( $upcase_DNA_triplet );
459 :     #
460 :     # User-supplied genetic code must be upper case index and match the
461 :     # DNA versus RNA type of sequence
462 :     #
463 :     # $seq = translate_seq_with_user_code( $seq, $gen_code_hash [, $met_start] )
464 :     #
465 :     #-----------------------------------------------------------------------------
466 :    
467 :     our %genetic_code = ( # work as DNA
468 :     TTT => "F", TCT => "S", TAT => "Y", TGT => "C",
469 :     TTC => "F", TCC => "S", TAC => "Y", TGC => "C",
470 :     TTA => "L", TCA => "S", TAA => "*", TGA => "*",
471 :     TTG => "L", TCG => "S", TAG => "*", TGG => "W",
472 :     CTT => "L", CCT => "P", CAT => "H", CGT => "R",
473 :     CTC => "L", CCC => "P", CAC => "H", CGC => "R",
474 :     CTA => "L", CCA => "P", CAA => "Q", CGA => "R",
475 :     CTG => "L", CCG => "P", CAG => "Q", CGG => "R",
476 :     ATT => "I", ACT => "T", AAT => "N", AGT => "S",
477 :     ATC => "I", ACC => "T", AAC => "N", AGC => "S",
478 :     ATA => "I", ACA => "T", AAA => "K", AGA => "R",
479 :     ATG => "M", ACG => "T", AAG => "K", AGG => "R",
480 :     GTT => "V", GCT => "A", GAT => "D", GGT => "G",
481 :     GTC => "V", GCC => "A", GAC => "D", GGC => "G",
482 :     GTA => "V", GCA => "A", GAA => "E", GGA => "G",
483 :     GTG => "V", GCG => "A", GAG => "E", GGG => "G",
484 :    
485 :     # The following ambiguous encodings are not necessary, but
486 :     # speed up the processing of some common ambiguous triplets:
487 :    
488 :     TTY => "F", TCY => "S", TAY => "Y", TGY => "C",
489 :     TTR => "L", TCR => "S", TAR => "*",
490 :     TCN => "S",
491 :     CTY => "L", CCY => "P", CAY => "H", CGY => "R",
492 :     CTR => "L", CCR => "P", CAR => "Q", CGR => "R",
493 :     CTN => "L", CCN => "P", CGN => "R",
494 :     ATY => "I", ACY => "T", AAY => "N", AGY => "S",
495 :     ACR => "T", AAR => "K", AGR => "R",
496 :     ACN => "T",
497 :     GTY => "V", GCY => "A", GAY => "D", GGY => "G",
498 :     GTR => "V", GCR => "A", GAR => "E", GGR => "G",
499 :     GTN => "V", GCN => "A", GGN => "G"
500 :     );
501 :    
502 :    
503 :     my %DNA_letter_can_be = (
504 :     A => ["A"], a => ["a"],
505 :     B => ["C", "G", "T"], b => ["c", "g", "t"],
506 :     C => ["C"], c => ["c"],
507 :     D => ["A", "G", "T"], d => ["a", "g", "t"],
508 :     G => ["G"], g => ["g"],
509 :     H => ["A", "C", "T"], h => ["a", "c", "t"],
510 :     K => ["G", "T"], k => ["g", "t"],
511 :     M => ["A", "C"], m => ["a", "c"],
512 :     N => ["A", "C", "G", "T"], n => ["a", "c", "g", "t"],
513 :     R => ["A", "G"], r => ["a", "g"],
514 :     S => ["C", "G"], s => ["c", "g"],
515 :     T => ["T"], t => ["t"],
516 :     U => ["T"], u => ["t"],
517 :     V => ["A", "C", "G"], v => ["a", "c", "g"],
518 :     W => ["A", "T"], w => ["a", "t"],
519 :     Y => ["C", "T"], y => ["c", "t"]
520 :     );
521 :    
522 :    
523 :     my %RNA_letter_can_be = (
524 :     A => ["A"], a => ["a"],
525 :     B => ["C", "G", "U"], b => ["c", "g", "u"],
526 :     C => ["C"], c => ["c"],
527 :     D => ["A", "G", "U"], d => ["a", "g", "u"],
528 :     G => ["G"], g => ["g"],
529 :     H => ["A", "C", "U"], h => ["a", "c", "u"],
530 :     K => ["G", "U"], k => ["g", "u"],
531 :     M => ["A", "C"], m => ["a", "c"],
532 :     N => ["A", "C", "G", "U"], n => ["a", "c", "g", "u"],
533 :     R => ["A", "G"], r => ["a", "g"],
534 :     S => ["C", "G"], s => ["c", "g"],
535 :     T => ["U"], t => ["u"],
536 :     U => ["U"], u => ["u"],
537 :     V => ["A", "C", "G"], v => ["a", "c", "g"],
538 :     W => ["A", "U"], w => ["a", "u"],
539 :     Y => ["C", "U"], y => ["c", "u"]
540 :     );
541 :    
542 :    
543 :     my %one_letter_to_three_letter_aa = (
544 :     A => "Ala",
545 :     B => "Asx",
546 :     C => "Cys",
547 :     D => "Asp",
548 :     E => "Glu",
549 :     F => "Phe",
550 :     G => "Gly",
551 :     H => "His",
552 :     I => "Ile",
553 :     K => "Lys",
554 :     L => "Leu",
555 :     M => "Met",
556 :     N => "Asn",
557 :     P => "Pro",
558 :     Q => "Gln",
559 :     R => "Arg",
560 :     S => "Ser",
561 :     T => "Thr",
562 :     U => "Sec",
563 :     V => "Val",
564 :     W => "Trp",
565 :     X => "Xxx",
566 :     Y => "Tyr",
567 :     Z => "Glx",
568 :     "*" => "***"
569 :     );
570 :    
571 :    
572 :     #-----------------------------------------------------------------------------
573 :     # Translate nucleotides to one letter protein:
574 :     #
575 :     # $seq = translate_seq( $seq [, $met_start] )
576 :     #
577 :     #-----------------------------------------------------------------------------
578 :    
579 :     sub translate_seq {
580 :     my $seq = uc shift;
581 :     $seq =~ tr/UX/TN/; # make it DNA, and allow X
582 :     $seq =~ tr/-//d; # remove gaps
583 :    
584 :     my $met = shift || 0; # a second argument that is true
585 :     # forces first amino acid to be Met
586 :     # (note: undef is false)
587 :    
588 :     my $imax = length($seq) - 2; # will try to translate 2 nucleotides!
589 :     my $pep = ($met && ($imax >= 0)) ? "M" : "";
590 :     my $aa;
591 :     for (my $i = $met ? 3 : 0; $i <= $imax; $i += 3) {
592 :     $pep .= translate_uc_DNA_codon( substr($seq,$i,3) );
593 :     }
594 :    
595 :     return $pep;
596 :     }
597 :    
598 :    
599 :     #-----------------------------------------------------------------------------
600 :     # Translate a single triplet with "universal" genetic code
601 :     # Uppercase and DNA are performed, then translate_uc_DNA_codon
602 :     # is called.
603 :     #
604 :     # $aa = translate_codon( $triplet )
605 :     #
606 :     #-----------------------------------------------------------------------------
607 :    
608 :     sub translate_codon {
609 :     my $codon = uc shift; # Make it uppercase
610 :     $codon =~ tr/UX/TN/; # Make it DNA, and allow X
611 :     return translate_uc_DNA_codon($codon);
612 :     }
613 :    
614 :    
615 :     #-----------------------------------------------------------------------------
616 :     # Translate a single triplet with "universal" genetic code
617 :     # Uppercase and DNA assumed
618 :     # Intended for private use by translate_codon and translate_seq
619 :     #
620 :     # $aa = translate_uc_DNA_codon( $triplet )
621 :     #
622 :     #-----------------------------------------------------------------------------
623 :    
624 :     sub translate_uc_DNA_codon {
625 :     my $codon = shift;
626 :     my $aa;
627 :    
628 :     # Try a simple lookup:
629 :    
630 :     if ( $aa = $genetic_code{ $codon } ) { return $aa }
631 :    
632 :     # With the expanded code defined above, this catches simple N, R
633 :     # and Y ambiguities in the third position. Other codes like
634 :     # GG[KMSWBDHV], or even GG, might be unambiguously translated by
635 :     # converting the last position to N and seeing if this is in the
636 :     # (expanded) code table:
637 :    
638 :     if ( $aa = $genetic_code{ substr($codon,0,2) . "N" } ) { return $aa }
639 :    
640 :     # Test that codon is valid and might have unambiguous aa:
641 :    
642 :     if ( $codon !~ m/^[ACGTMY][ACGT][ACGTKMRSWYBDHVN]$/ ) { return "X" }
643 :     # ^^
644 :     # |+- for leucine YTR
645 :     # +-- for arginine MGR
646 :    
647 :     # Expand all ambiguous nucleotides to see if they all yield same aa.
648 :     # Loop order tries to fail quickly with first position change.
649 :    
650 :     $aa = "";
651 :     for my $n2 ( @{ $DNA_letter_can_be{ substr($codon,1,1) } } ) {
652 :     for my $n3 ( @{ $DNA_letter_can_be{ substr($codon,2,1) } } ) {
653 :     for my $n1 ( @{ $DNA_letter_can_be{ substr($codon,0,1) } } ) {
654 :     # set the first value of $aa
655 :     if ($aa eq "") { $aa = $genetic_code{ $n1 . $n2 . $n3 } }
656 :     # or break out if any other amino acid is detected
657 :     elsif ($aa ne $genetic_code{ $n1 . $n2 . $n3 } ) { return "X" }
658 :     }
659 :     }
660 :     }
661 :    
662 :     return $aa || "X";
663 :     }
664 :    
665 :    
666 :     #-----------------------------------------------------------------------------
667 :     # Translate with a user-supplied genetic code to translate a sequence.
668 :     # Diagnose the use of upper versus lower, and T versus U in the supplied
669 :     # code, and transform the supplied nucleotide sequence to match.
670 :     #
671 :     # translate_seq_with_user_code($seq, \%gen_code [, $start_with_met] )
672 :     #
673 :     #-----------------------------------------------------------------------------
674 :    
675 :     sub translate_seq_with_user_code {
676 :     my $seq = shift;
677 :     $seq =~ tr/-//d; # remove gaps *** Why?
678 :     $seq =~ tr/Xx/Nn/; # allow X
679 :    
680 :     my $gc = shift; # Reference to hash of DNA alphabet code
681 :     if (! $gc || ref($gc) ne "HASH") {
682 :     die "translate_seq_with_user_code needs genetic code hash as secondargument.";
683 :     }
684 :    
685 :     # Test the type of code supplied: uppercase versus lowercase
686 :    
687 :     my ($RNA_F, $DNA_F, $M, $N, $X);
688 :    
689 :     if ($gc->{ "AAA" }) { # Looks like uppercase code table
690 :     $seq = uc $seq; # Uppercase sequence
691 :     $RNA_F = "UUU"; # Uppercase RNA Phe codon
692 :     $DNA_F = "TTT"; # Uppercase DNA Phe codon
693 :     $M = "M"; # Uppercase initiator
694 :     $N = "N"; # Uppercase ambiguous nuc
695 :     $X = "X"; # Uppercase ambiguous aa
696 :     }
697 :     elsif ($gc->{ "aaa" }) { # Looks like lowercase code table
698 :     $seq = lc $seq; # Lowercase sequence
699 :     $RNA_F = "uuu"; # Lowercase RNA Phe codon
700 :     $DNA_F = "ttt"; # Lowercase DNA Phe codon
701 :     $M = "m"; # Lowercase initiator
702 :     $N = "n"; # Lowercase ambiguous nuc
703 :     $X = "x"; # Lowercase ambiguous aa
704 :     }
705 :     else {
706 :     die "User-supplied genetic code does not have aaa or AAA\n";
707 :     }
708 :    
709 :     # Test the type of code supplied: UUU versus TTT
710 :    
711 :     my ($ambigs);
712 :    
713 :     if ($gc->{ $RNA_F }) { # Looks like RNA code table
714 :     $seq =~ tr/Tt/Uu/;
715 :     $ambigs = \%RNA_letter_can_be;
716 :     }
717 :     elsif ($gc->{ $DNA_F }) { # Looks like DNA code table
718 :     $seq =~ tr/Uu/Tt/;
719 :     $ambigs = \%DNA_letter_can_be;
720 :     }
721 :     else {
722 :     die "User-supplied genetic code does not have $RNA_F or $DNA_F\n";
723 :     }
724 :    
725 :     my $imax = length($seq) - 2; # will try to translate 2 nucleotides!
726 :    
727 :     my $met = shift; # a third argument that is true
728 :     # forces first amino acid to be Met
729 :     # (note: undef is false)
730 :     my $pep = ($met && ($imax >= 0)) ? $M : "";
731 :     my $aa;
732 :    
733 :     for (my $i = $met ? 3 : 0; $i <= $imax; $i += 3) {
734 :     $pep .= translate_codon_with_user_code( substr($seq,$i,3), $gc, $N, $X, $ambigs );
735 :     }
736 :    
737 :     return $pep;
738 :     }
739 :    
740 :    
741 :     #-----------------------------------------------------------------------------
742 :     # Translate with user-supplied genetic code hash. For speed, no error
743 :     # check on the hash. Calling programs should check for the hash at a
744 :     # higher level.
745 :     #
746 :     # Should only be called through translate_seq_with_user_code
747 :     #
748 :     # translate_codon_with_user_code( $triplet, \%code, $N, $X, $ambig_table )
749 :     #
750 :     # $triplet speaks for itself
751 :     # $code ref to the hash with the codon translations
752 :     # $N character to use for ambiguous nucleotide
753 :     # $X character to use for ambiguous amino acid
754 :     # $ambig_table ref to hash with lists of nucleotides for each ambig code
755 :     #-----------------------------------------------------------------------------
756 :    
757 :    
758 :     sub translate_codon_with_user_code {
759 :     my $codon = shift;
760 :     my $gc = shift;
761 :     my $aa;
762 :    
763 :     # Try a simple lookup:
764 :    
765 :     if ( $aa = $gc->{ $codon } ) { return $aa }
766 :    
767 :     # Test that codon is valid and might have unambiguous aa:
768 :    
769 :     my ($N, $X, $ambigs) = @_;
770 :     if ( $codon =~ m/^[ACGTUMY][ACGTU]$/i ) { $codon .= $N }
771 :     if ( $codon !~ m/^[ACGTUMY][ACGTU][ACGTUKMRSWYBDHVN]$/i ) { return $X }
772 :     # ^^
773 :     # |+- for leucine YTR
774 :     # +-- for arginine MGR
775 :    
776 :     # Expand all ambiguous nucleotides to see if they all yield same aa.
777 :     # Loop order tries to fail quickly with first position change.
778 :    
779 :     $aa = "";
780 :     for my $n2 ( @{ $ambigs->{ substr($codon,1,1) } } ) {
781 :     for my $n3 ( @{ $ambigs->{ substr($codon,2,1) } } ) {
782 :     for my $n1 ( @{ $ambigs->{ substr($codon,0,1) } } ) {
783 :     # set the first value of $aa
784 :     if ($aa eq "") { $aa = $gc->{ $n1 . $n2 . $n3 } }
785 :     # break out if any other amino acid is detected
786 :     elsif ($aa ne $gc->{ $n1 . $n2 . $n3 } ) { return "X" }
787 :     }
788 :     }
789 :     }
790 :    
791 :     return $aa || $X;
792 :     }
793 :    
794 :    
795 :     #-----------------------------------------------------------------------------
796 :     # Read a list of intervals from a file.
797 :     # Allow id_start_end, or id \s start \s end formats
798 :     #
799 :     # @intervals = read_intervals(*FILEHANDLE)
800 :     #-----------------------------------------------------------------------------
801 :     sub read_intervals {
802 :     my $fh = shift;
803 :     my @intervals = ();
804 :    
805 :     while (<$fh>) {
806 :     chomp;
807 :     /^(\S+)_(\d+)_(\d+)(\s.*)?$/ # id_start_end WIT2
808 :     || /^(\S+)_(\d+)-(\d+)(\s.*)?$/ # id_start-end ???
809 :     || /^(\S+)=(\d+)=(\d+)(\s.*)?$/ # id=start=end Badger
810 :     || /^(\S+)\s+(\d+)\s+(\d+)(\s.*)?$/ # id \s start \s end
811 :     || next;
812 :    
813 :     # Matched a pattern. Store reference to (id, left, right):
814 :     push @intervals, ($2 < $3) ? [ $1, $2+0, $3+0 ]
815 :     : [ $1, $3+0, $2+0 ];
816 :     }
817 :     return @intervals;
818 :     }
819 :    
820 :    
821 :     #-----------------------------------------------------------------------------
822 :     # Take the union of a list of intervals
823 :     #
824 :     # @joined = join_intervals( @interval_refs )
825 :     #-----------------------------------------------------------------------------
826 :     sub join_intervals {
827 :     my @ordered = sort { $a->[0] cmp $b->[0] # first by id
828 :     || $a->[1] <=> $b->[1] # next by left end
829 :     || $b->[2] <=> $a->[2] # finally longest first
830 :     } @_;
831 :    
832 :     my @joined = ();
833 :     my $n_int = @ordered;
834 :    
835 :     my ($cur_id) = "";
836 :     my ($cur_left) = -1;
837 :     my ($cur_right) = -1;
838 :     my ($new_id, $new_left, $new_right);
839 :    
840 :     for (my $i = 0; $i < $n_int; $i++) {
841 :     ($new_id, $new_left, $new_right) = @{$ordered[$i]}; # get the new data
842 :    
843 :     if ( ( $new_id ne $cur_id) # if new contig
844 :     || ( $new_left > $cur_right + 1) # or not touching previous
845 :     ) { # push the previous interval
846 :     if ($cur_id) { push (@joined, [ $cur_id, $cur_left, $cur_right ]) }
847 :     $cur_id = $new_id; # update the current interval
848 :     $cur_left = $new_left;
849 :     $cur_right = $new_right;
850 :     }
851 :    
852 :     elsif ($new_right > $cur_right) { # extend the right end if necessary
853 :     $cur_right = $new_right;
854 :     }
855 :     }
856 :    
857 :     if ($cur_id) { push (@joined, [$cur_id, $cur_left, $cur_right]) }
858 :     return @joined;
859 :     }
860 :    
861 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3