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

Annotation of /FigKernelPackages/find_special_proteins.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.11
2 :     #
3 : olson 1.12 # This is a SAS component
4 :     #
5 : olson 1.11
6 : golsen 1.3 #
7 :     # Copyright (c) 2003-2007 University of Chicago and Fellowship
8 :     # for Interpretations of Genomes. All Rights Reserved.
9 :     #
10 :     # This file is part of the SEED Toolkit.
11 :     #
12 :     # The SEED Toolkit is free software. You can redistribute
13 :     # it and/or modify it under the terms of the SEED Toolkit
14 :     # Public License.
15 :     #
16 :     # You should have received a copy of the SEED Toolkit Public License
17 :     # along with this program; if not write to the University of Chicago
18 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
19 :     # Genomes at veronika@thefig.info or download a copy from
20 :     # http://www.theseed.org/LICENSE.TXT.
21 :     #
22 :    
23 : golsen 1.1 package find_special_proteins;
24 :    
25 :     use strict;
26 :     use gjoseqlib;
27 :     use gjoparseblast;
28 : golsen 1.5 use NCBI_genetic_code;
29 : golsen 1.1
30 : golsen 1.5 # use Data::Dumper;
31 : golsen 1.1
32 :     #===============================================================================
33 :     # Use a database of reference sequences to locate selenoproteins in a genome.
34 :     # This program will not (by design) find new instances of selenocysteine.
35 : golsen 1.4 # With the pyrrolysine option, find pyrrolysyl proteins.
36 : golsen 1.1 #
37 :     # @locs = find_selenoproteins( \%params )
38 :     #
39 :     # Required parameters:
40 :     #
41 : golsen 1.6 # contigs => \@contigs # genome is a synonym
42 : golsen 1.1 #
43 :     # Optional paramaters:
44 :     #
45 : golsen 1.5 # allow_C => 1 # allow C in reference to align with SeC (D = 0)
46 :     # allow_K => 1 # allow K in reference to align with PyK (D = 0)
47 : golsen 1.6 # comment => $text # comment for output hash
48 : golsen 1.4 # is_init => \@initiator_codons # D = [ATG,GTG]
49 :     # is_alt => \@second_choice_init_codons # D = [TTG]
50 :     # is_term => \@terminator_codons # D = [TAA,TAG,TGA]
51 : golsen 1.6 # note => $text # same as 'comment'
52 : golsen 1.4 # pyrrolysine => 1 # find pyrrolysylproteins
53 :     # references => \@selenoproteins # selenocysteine must be U or X
54 :     # references => \@pyrrolysylproteins # pyrrolysine should be X
55 :     # tmp => $directory # directory for tmp_dir
56 :     # tmp_dir => $directory # directory for temporary files
57 : golsen 1.1 #
58 : golsen 1.4 # Some keys can be shortened.
59 : golsen 1.1 #===============================================================================
60 :     sub find_selenoproteins
61 :     {
62 :     my ( $params ) = @_;
63 :    
64 :     my ( $tmp_dir, $save_tmp ) = temporary_directory( $params );
65 :    
66 : golsen 1.4 my ( $contig_key ) = grep { /^contig/ || /^gen/ } keys %$params;
67 :     my $contigs = $contig_key ? $params->{ $contig_key } : undef;
68 :     ( ref( $contigs ) eq 'ARRAY' ) && @$contigs
69 : golsen 1.1 or return ();
70 :    
71 : golsen 1.4 my %contigR = map { $_->[0] => \$_->[2] } @$contigs; # Index sequence ref by id
72 :    
73 :     my ( $pyrro_key ) = grep { /^pyrrolys/ } keys %$params;
74 :     my $pyrrolys = $pyrro_key ? $params->{ $pyrro_key } : 0;
75 :    
76 : golsen 1.6 my ( $comment_key ) = grep { /^comment/ || /^note/ } keys %$params;
77 :     my $comment = $comment_key ? $params->{ $comment_key } : '';
78 :    
79 : golsen 1.4 # Here are most of the differences between selenocysteine and
80 :     # pyrrolysine processing:
81 :    
82 :     my $my_stop = $pyrrolys ? 'TAG' : 'TGA';
83 :     my $my_symb = $pyrrolys ? 'X' : 'U';
84 : golsen 1.6 $comment ||= $pyrrolys ? 'pyrrolysylprotein' : 'selenoprotein';
85 : golsen 1.1
86 : golsen 1.4 my %gencode = %gjoseqlib::genetic_code;
87 :     $gencode{ $my_stop } = $my_symb;
88 : golsen 1.1
89 : golsen 1.4 # Get reference sequences as parameter, or get them from a module:
90 : golsen 1.1
91 : golsen 1.4 my ( $ref_key ) = grep { /^ref/ } keys %$params;
92 :     my $refs = $ref_key ? $params->{ $ref_key } : undef;
93 : golsen 1.3 if ( $refs )
94 :     {
95 :     ( ref( $refs ) eq 'ARRAY' ) && @$refs
96 : golsen 1.4 or print STDERR "No reference sequences supplied to find_selenoproteins().\n"
97 : golsen 1.3 and return ();
98 :     }
99 : golsen 1.4 elsif ( $pyrrolys )
100 :     {
101 :     eval
102 :     {
103 :     require pyrrolysylprotein_ref_seq;
104 :     $refs = $pyrrolysylprotein_ref_seq::ref_seqs;
105 :     };
106 :     $refs && ( ref( $refs ) eq 'ARRAY' ) && @$refs
107 :     or print STDERR "Unable to get reference sequences from pyrrolysylprotein_ref_seq.pm.\n"
108 :     and return ();
109 :     }
110 : golsen 1.3 else
111 :     {
112 :     eval
113 :     {
114 :     require selenoprotein_ref_seq;
115 :     $refs = $selenoprotein_ref_seq::ref_seqs;
116 :     };
117 :     $refs && ( ref( $refs ) eq 'ARRAY' ) && @$refs
118 :     or print STDERR "Unable to get reference sequences from selenoprotein_ref_seq.pm.\n"
119 :     and return ();
120 :     }
121 : golsen 1.1
122 :     # Normally, the query (reference) amino acid must be U or X. To avoid
123 :     # annoying formatdb messages, U is transliterated to X before writing
124 :     # the reference database. To find a new selenoprotein, one option is
125 :     # to use a genome from a nonselenocysteine containing organism as the
126 : golsen 1.3 # reference, and allow TGA codons to align with C in the reference
127 : golsen 1.1 # sequence. This option enables that function. Otherwise C in the
128 :     # reference is not allowed to align with a proposed selenocysteine.
129 : golsen 1.4 # Ditto for pyrrolysine.
130 : golsen 1.1
131 : golsen 1.4 my ( $allow_key ) = grep { /^allow/ } keys %$params;
132 :     my $ref_aa = $params->{ $allow_key } && $pyrrolys ? qr/[KX]/
133 :     : $params->{ $allow_key } ? qr/[CX]/
134 :     : qr/X/;
135 : golsen 1.1
136 :     my $suffix = $$ . "_" . sprintf( "%09d", int( 1e9 * rand() ) );
137 :    
138 :     # Print and format the contig sequences:
139 :    
140 :     my $contig_file = "$tmp_dir/contigs_$suffix";
141 :     gjoseqlib::print_alignment_as_fasta( $contig_file, $contigs );
142 :     -f $contig_file or return ();
143 :    
144 :     my $ext_bin = $FIG_Config::ext_bin;
145 :     my $formatdb = $params->{ formatdb };
146 :     $formatdb ||= $ext_bin ? "$ext_bin/formatdb" : 'formatdb';
147 :     system "cd '$tmp_dir'; $formatdb -p f -i 'contigs_$suffix'";
148 :    
149 : golsen 1.4 # Make a clean copy of ref sequences to avoid BLAST warnings.
150 : golsen 1.1
151 : golsen 1.4 $refs = [ map { my $s = uc $_->[2]; # uppercase copy
152 :     $s =~ s/[BJOUZ]/X/g; # U etc. -> X
153 :     $s =~ s/[^A-Z]//g; # delete nonletters
154 :     [ @$_[0,1], $s ]
155 :     } @$refs
156 :     ];
157 : golsen 1.1 my $ref_file = "$tmp_dir/refs_$suffix";
158 :     gjoseqlib::print_alignment_as_fasta( $ref_file, $refs );
159 : golsen 1.4 -f $ref_file or print STDERR "Failed to write reference sequences to $ref_file.\n"
160 :     and return ();
161 : golsen 1.1
162 :     # Search the contigs for the reference sequences:
163 :    
164 :     my $blastall = $params->{ blastall };
165 :     $blastall ||= $ext_bin ? "$ext_bin/blastall" : 'blastall';
166 : golsen 1.4 my $blastcmd = "$blastall -p tblastn -d '$contig_file' -i '$ref_file' "
167 :     . "-e 1e-10 -F f -b 200 -v 200 -a 2 |";
168 :     open( BLAST, $blastcmd )
169 :     or print STDERR "Could not open pipe '$blastcmd'.\n"
170 :     and return ();
171 : golsen 1.1 my @hsps = gjoparseblast::blast_hsp_list( \*BLAST );
172 :     close BLAST;
173 :    
174 :     # Delete the temporary directory unless it already existed:
175 :    
176 :     system "/bin/rm -r $tmp_dir" if ! $save_tmp;
177 :    
178 :     #
179 :     # Stucture of the output records:
180 :     #
181 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
182 :     # qid qdef qlen sid sdef slen scr e_val p_n p_val n_mat n_id n_pos n_gap dir q1 q2 qseq s1 s2 sseq
183 :     #
184 :    
185 :     @hsps = sort { $b->[6] <=> $a->[6] } @hsps;
186 :    
187 :     my $covered = {};
188 :     my %hit = ();
189 : golsen 1.6 my %def = ();
190 : golsen 1.1 my $hsp;
191 :     foreach $hsp ( @hsps )
192 :     {
193 :     # This is a chance to filter out things that are already done.
194 :     # Currently, is_covered() returns 0
195 :    
196 :     next if is_covered( $hsp, $covered );
197 :    
198 :     # List of locations with stops:
199 :    
200 :     my $i = 0;
201 : golsen 1.6 my ( $qid, $qdef, $sid, $scr, $qseq, $s1, $s2, $sseq ) = @$hsp[ 0, 1, 3, 6, 17, 18 .. 20 ];
202 : golsen 1.1 my @stops = map { $i++; $_ eq '*' ? $i : () }
203 :     split //, $sseq;
204 :     next if ! @stops;
205 :    
206 : golsen 1.6 # Changing the next statement might allow finding new instances.
207 :     # Unless a whole genome is used as a reference set, this is more
208 :     # likely to yield false positives than authentic selenoproteins,
209 :     # but it might be a useful exploratory tool.
210 : golsen 1.1
211 :     next if grep { substr( $qseq, $_-1, 1 ) !~ /$ref_aa/o } @stops;
212 :    
213 :     my $dir = $s2 <=> $s1;
214 : golsen 1.4 my $contigR = $contigR{ $sid };
215 : golsen 1.1 my $contig_len = length( $$contigR );
216 :    
217 :     my $stop;
218 :     my $is_okay = 1;
219 :     foreach $stop ( @stops )
220 :     {
221 :     my $prefix = substr( $sseq, 0, $stop-1 );
222 :     $prefix =~ s/-//g;
223 :     my $n1 = $s1 + 3 * $dir * length( $prefix );
224 :     my $n2 = $n1 + 2 * $dir;
225 :     my $codon = uc gjoseqlib::DNA_subseq( $contigR, $n1, $n2 );
226 : golsen 1.4 $is_okay = ( ( $codon eq $my_stop ) ? 1 : 0 ) or last;
227 : golsen 1.1 }
228 :     $is_okay or next;
229 :    
230 :     # Follow the orf to a start:
231 :    
232 :     my ( $from, $partial_5 ) = find_orf_start( $contigR, $s1, $s2, $params );
233 :     $from or next;
234 :    
235 :     # Follow the orf to its end:
236 :    
237 :     my ( $to, $partial_3 ) = find_orf_end( $contigR, $s1, $s2, $params );
238 :     $to or next;
239 :    
240 :     my $ntseq = uc DNA_subseq( $contigR, $from, $partial_3 ? $to : $to - 3*$dir );
241 : golsen 1.4 my $aaseq = translate_seq_with_user_code( $ntseq, \%gencode, ! $partial_5 );
242 : golsen 1.1
243 :     # Save the protein unless we already have a longer hit to the same orf.
244 :    
245 :     my $key = "$sid\t$to\t$dir";
246 :     my $len = length( $aaseq );
247 : golsen 1.4 $hit{$key} = [ $sid, $from, $to, $aaseq, $len ] unless ( $hit{$key} && $hit{$key}->[4] >= $len );
248 : golsen 1.6 $def{$key} = [ $qid, $qdef, $scr ] unless ( $def{$key} && $def{$key}->[2] >= $scr );
249 : golsen 1.1 }
250 :    
251 : golsen 1.8 # Sort by contig and midpoint location, and return a hash of location,
252 : golsen 1.1 # sequence and comment for each:
253 :    
254 :     my @prots = map { scalar { location => join( '_', @$_[0..2] ),
255 :     sequence => $_->[3],
256 : golsen 1.6 comment => $comment,
257 :     ( $_->[5] ? ( reference_id => $_->[5] ) : () ),
258 :     ( $_->[6] ? ( reference_def => $_->[6] ) : () )
259 : golsen 1.1 }
260 :     }
261 :     sort { $a->[0] cmp $b->[0] || ( ( $a->[1]+$a->[2] ) <=> ( $b->[1]+$b->[2] ) ) }
262 : golsen 1.8 map { [ @{$hit{$_}}, @{$def{$_}} ] } # [ $sid, $from, $to, $aaseq, $len, $qid, $qdef, $scr ]
263 : golsen 1.1 keys %hit;
264 :    
265 :     wantarray ? @prots : \@prots;
266 :     }
267 :    
268 :    
269 :    
270 : golsen 1.5 #===============================================================================
271 :     # Use a database of reference sequences to locate proteins in a genome.
272 :     # Options allow alteration of intiator and terminator.
273 :     #
274 :     # @locs = find_protein_homologs( \%params )
275 :     #
276 :     # Required parameters:
277 :     #
278 :     # contigs => \@contigs # genome is a synonym
279 :     # references => \@proteins # selenocysteine must be U or X
280 :     #
281 :     # Optional paramaters:
282 :     #
283 :     # code => \%genetic_code # D = standard
284 :     # code => $NCBI_code_number # D = standard
285 :     # comment => $comment_text # comment is attached to each protein
286 :     # expect => $blast_e_value
287 :     # is_init => \@initiator_codons # D = [ATG,GTG]
288 :     # is_alt => \@second_choice_init_codons # D = [TTG]
289 :     # is_term => \@terminator_codons # D = [TAA,TAG,TGA]
290 :     # tmp => $directory # directory for tmp_dir
291 :     # tmp_dir => $directory # directory for temporary files
292 :     #
293 :     # Some keys can be shortened.
294 :     #===============================================================================
295 :     sub find_protein_homologs
296 :     {
297 :     my ( $params ) = @_;
298 :    
299 :     my ( $tmp_dir, $save_tmp ) = temporary_directory( $params );
300 :    
301 :     my ( $contig_key ) = grep { $_ !~ /code/ } # Don't get "genetic_code"
302 :     grep { /^contig/ || /^gen/ }
303 :     keys %$params;
304 :     my $contigs = $contig_key ? $params->{ $contig_key } : undef;
305 :     ( ref( $contigs ) eq 'ARRAY' ) && @$contigs
306 :     or return ();
307 :    
308 :     my %contigR = map { $_->[0] => \$_->[2] } @$contigs; # Index sequence ref by id
309 :    
310 :     # Get reference sequences as parameter:
311 :    
312 :     my ( $ref_key ) = grep { /^ref/ } keys %$params;
313 :     my $refs;
314 :     if ( ! ( $ref_key && ( $refs = $params->{ $ref_key } )
315 :     && ( ref( $refs ) eq 'ARRAY' )
316 :     && @$refs
317 :     )
318 :     )
319 :     {
320 :     print STDERR "No reference sequences supplied to find_protein_homologs().\n";
321 :     return ();
322 :     }
323 :    
324 :     # Optional parameters
325 :    
326 :     my $comment = $params->{ comment } || 'Extracted by find_protein_homologs';
327 :    
328 :     my ( $exp_key ) = grep { /^exp/ } keys %$params;
329 :     my $expect = $exp_key ? ( $params->{ $exp_key } + 0 ) : 1e-10;
330 :    
331 :     my ( $code_key ) = grep { /code/ } keys %$params;
332 :     my %gencode;
333 :     if ( $code_key )
334 :     {
335 :     my $code = $params->{ $code_key };
336 :     if ( ref( $code ) eq 'HASH' )
337 :     {
338 :     %gencode = %$code;
339 :     }
340 :     elsif ( ref NCBI_genetic_code::genetic_code( $code ) eq 'HASH' )
341 :     {
342 :     %gencode = %{ NCBI_genetic_code::genetic_code( $code ) }
343 :     }
344 :     else
345 :     {
346 :     print STDERR "find_protein_homologs genetic code not a HASH reference.\n";
347 :     return ();
348 :     }
349 :     }
350 :     else
351 :     {
352 :     %gencode = %gjoseqlib::genetic_code;
353 :     }
354 :    
355 :     # Do the analysis
356 :    
357 :     # Print and format the contig sequences:
358 :    
359 :     my $suffix = $$ . "_" . sprintf( "%09d", int( 1e9 * rand() ) );
360 :     my $contig_file = "$tmp_dir/contigs_$suffix";
361 :     gjoseqlib::print_alignment_as_fasta( $contig_file, $contigs );
362 :     -f $contig_file or return ();
363 :    
364 :     my $ext_bin = $FIG_Config::ext_bin;
365 :     my $formatdb = $params->{ formatdb };
366 :     $formatdb ||= $ext_bin ? "$ext_bin/formatdb" : 'formatdb';
367 :     system "cd '$tmp_dir'; $formatdb -p f -i 'contigs_$suffix'";
368 :    
369 :     # Make a clean copy of ref sequences to avoid BLAST warnings.
370 :    
371 :     $refs = [ map { my $s = uc $_->[2]; # uppercase copy
372 :     $s =~ s/[BJOUZ]/X/g; # U etc. -> X
373 :     $s =~ s/[^A-Z]//g; # delete nonletters
374 :     [ @$_[0,1], $s ]
375 :     } @$refs
376 :     ];
377 :     my $ref_file = "$tmp_dir/refs_$suffix";
378 :     gjoseqlib::print_alignment_as_fasta( $ref_file, $refs );
379 :     -f $ref_file or print STDERR "Failed to write reference sequences to $ref_file.\n"
380 :     and return ();
381 :    
382 :     # Search the contigs for the reference sequences:
383 :    
384 :     my $blastall = $params->{ blastall };
385 :     $blastall ||= $ext_bin ? "$ext_bin/blastall" : 'blastall';
386 :     my $blastcmd = "$blastall -p tblastn -d '$contig_file' -i '$ref_file' "
387 :     . "-e $expect -F f -b 200 -v 200 -a 2 |";
388 :     open( BLAST, $blastcmd )
389 :     or print STDERR "Could not open pipe '$blastcmd'.\n"
390 :     and return ();
391 :     my @hsps = gjoparseblast::blast_hsp_list( \*BLAST );
392 :     close BLAST;
393 :    
394 :     # Delete the temporary directory unless it already existed:
395 :    
396 :     system "/bin/rm -r $tmp_dir" if ! $save_tmp;
397 :    
398 :     #
399 :     # Stucture of the output records:
400 :     #
401 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
402 :     # qid qdef qlen sid sdef slen scr e_val p_n p_val n_mat n_id n_pos n_gap dir q1 q2 qseq s1 s2 sseq
403 :     #
404 :    
405 :     @hsps = sort { $b->[6] <=> $a->[6] } @hsps;
406 :    
407 :     my $covered = {};
408 :     my %hit = ();
409 :     my $hsp;
410 :     foreach $hsp ( @hsps )
411 :     {
412 :     # This is a chance to filter out things that are already done.
413 :     # Currently, is_covered() returns 0
414 :    
415 :     next if is_covered( $hsp, $covered );
416 :    
417 :     # For the moment, no stops allowed:
418 :    
419 :     my ( $sid, $qseq, $s1, $s2, $sseq ) = @$hsp[ 3, 17, 18 .. 20 ];
420 :     next if $sseq =~ /\*/;
421 :    
422 :     my $dir = $s2 <=> $s1;
423 :     my $contigR = $contigR{ $sid };
424 :     my $contig_len = length( $$contigR );
425 :    
426 :     # Follow the orf to a start:
427 :    
428 :     my ( $from, $partial_5 ) = find_orf_start( $contigR, $s1, $s2, $params );
429 :     $from or next;
430 :    
431 :     # Follow the orf to its end:
432 :    
433 :     my ( $to, $partial_3 ) = find_orf_end( $contigR, $s1, $s2, $params );
434 :     $to or next;
435 :    
436 :     my $ntseq = uc DNA_subseq( $contigR, $from, $partial_3 ? $to : $to - 3*$dir );
437 :     my $aaseq = translate_seq_with_user_code( $ntseq, \%gencode, ! $partial_5 );
438 :    
439 :     # Save the protein unless we already have a longer hit to the same orf.
440 :    
441 :     my $key = "$sid\t$to\t$dir";
442 :     my $len = length( $aaseq );
443 :     $hit{$key} = [ $sid, $from, $to, $aaseq, $len ] unless ( $hit{$key} && $hit{$key}->[4] >= $len );
444 :     }
445 :    
446 :     # Sort by contig and midpoint location, and return a hash or location,
447 :     # sequence and comment for each:
448 :    
449 :     my @prots = map { scalar { location => join( '_', @$_[0..2] ),
450 :     sequence => $_->[3],
451 :     comment => $comment
452 :     }
453 :     }
454 :     sort { $a->[0] cmp $b->[0] || ( ( $a->[1]+$a->[2] ) <=> ( $b->[1]+$b->[2] ) ) }
455 :     map { $hit{$_} }
456 :     keys %hit;
457 :    
458 :     wantarray ? @prots : \@prots;
459 :     }
460 :    
461 :    
462 : golsen 1.1 #-------------------------------------------------------------------------------
463 :     # This is a place holder for minimizing duplicate analysis of orfs that
464 :     # have already been found. Currently, duplicates are removed by storing
465 :     # proteins by end location and direction, saving the longest.
466 :     #-------------------------------------------------------------------------------
467 :     sub is_covered
468 :     {
469 :     my ( $hsp, $covered ) = @_;
470 :     0;
471 :     }
472 :    
473 :    
474 :     #-------------------------------------------------------------------------------
475 :     # find_orf_start
476 :     #
477 :     # ( $start, $is_partial ) = find_orf_start( \$contig, $s1, $s2, $options )
478 :     #
479 :     # The start is defined by the first AUG or GUG upstream of $s1, or the
480 : gdpusch 1.7 # first UUG upstream of $s1. The triplet starting with $s1 is
481 :     # the initial search position. If the search falls off an end of the
482 : golsen 1.1 # contig, then the $start marks the start of the last complete codon,
483 :     # and $is_partial is true. If extending the orf hits a stop before a
484 :     # start, then we we truncate the orf to the start closest to $s1. If
485 :     # this fails, we return the empty list.
486 :     #
487 :     # Options:
488 :     #
489 :     # is_init => \@initiator_triplets ( D = [ATG,GTG] )
490 :     # is_alt => \@second_choice_initiator_triplets ( D = [TTG] )
491 :     # is_term => \@terminator_triplets ( D = [TAA,TAG,TGA] )
492 :     #
493 :     #-------------------------------------------------------------------------------
494 :     sub find_orf_start
495 :     {
496 :     my ( $contigR, $s1, $s2, $options ) = @_;
497 :     my $contig_len = length( $$contigR );
498 :     my $dir = $s2 <=> $s1;
499 :    
500 :     # Alternative start codons are only returned if essential to
501 :     # avoid a terminator:
502 :    
503 :     my %is_init = opt_hash( $options, 'is_init', [ qw( ATG GTG ) ] );
504 :     my %is_alt = opt_hash( $options, 'is_alt', [ qw( TTG ) ] );
505 :     my $alt_init;
506 :    
507 :     my %is_term = opt_hash( $options, 'is_term', [ qw( TAA TAG TGA ) ] );
508 :     my $n1 = $s1;
509 :     my $n2 = $n1 + 2 * $dir;
510 :    
511 :     my $extend = 1;
512 :     while ( $extend )
513 :     {
514 :     my $codon = uc gjoseqlib::DNA_subseq( $contigR, $n1, $n2 );
515 :     if ( $is_init{ $codon } )
516 :     {
517 :     return ( $n1, 0 );
518 :     }
519 :     elsif ( $is_term{ $codon } && ( $n1 != $s1 ) )
520 :     {
521 :     return ( $alt_init, 0 ) if $alt_init;
522 :     $extend = 0;
523 :     }
524 :     # Running off end of contig?
525 :     elsif ( $dir > 0 ? $n1 - 3 < 1 : $n1 + 3 > $contig_len )
526 :     {
527 :     return ( $n1, 1 );
528 :     }
529 :     else
530 :     {
531 : gdpusch 1.7 $alt_init = $n1 if ! $alt_init && $is_alt{ $codon };
532 : golsen 1.1 $n2 = $n1 - $dir;
533 :     $n1 = $n2 - 2 * $dir;
534 :     }
535 :     }
536 :    
537 :     # We failed to extend the match region to an initiator,
538 :     # now we cut into it:
539 :    
540 : gdpusch 1.9 $n1 = $s1;
541 :     $n2 = $n1 + 2 * $dir;
542 : golsen 1.1 while ( 1 )
543 :     {
544 :     # Running off end of match region?
545 :     if ( $dir > 0 ? $n2 + 3 > $s2 : $n2 - 3 < $s2 )
546 :     {
547 :     return ();
548 :     }
549 :     $n1 = $n2 + $dir;
550 :     $n2 = $n1 + 2 * $dir;
551 :     my $codon = uc gjoseqlib::DNA_subseq( $contigR, $n1, $n2 );
552 :     if ( $is_init{ $codon } || $is_alt{ $codon } )
553 :     {
554 :     return ( $n1, 0 );
555 :     }
556 :     }
557 :     }
558 :    
559 :    
560 :     #-------------------------------------------------------------------------------
561 :     # ( $end, $is_partial ) = find_orf_end( \$contig, $s1, $s2, $options )
562 :     #
563 :     # The end is defined as the first UAA, UAG or UGA downstream of $s2. If
564 :     # not end is found before the end of the contig, then $is_partial is true.
565 :     #
566 :     # Options:
567 :     #
568 :     # is_term => \@terminator_triplets ( D = [TAA,TAG,TGA] )
569 :     #
570 :     #-------------------------------------------------------------------------------
571 :     sub find_orf_end
572 :     {
573 :     my ( $contigR, $s1, $s2, $options ) = @_;
574 :     my $contig_len = length( $$contigR );
575 :     my $dir = $s2 <=> $s1;
576 :    
577 :     my %is_term = opt_hash( $options, 'is_term', [ qw( TAA TAG TGA ) ] );
578 :    
579 :     my $n2 = $s2;
580 :     my $n1 = $n2 - 2 * $dir;
581 :     while ( 1 )
582 :     {
583 :     # Running off end of contig?
584 :     if ( $dir > 0 ? $n2 + 3 > $contig_len : $n2 - 3 < 1 )
585 :     {
586 :     return ( $n2, 1 );
587 :     }
588 :     $n1 = $n2 + $dir;
589 :     $n2 = $n1 + 2 * $dir;
590 :     my $codon = uc gjoseqlib::DNA_subseq( $contigR, $n1, $n2 );
591 :     if ( $is_term{ $codon } )
592 :     {
593 :     return ( $n2, 0 );
594 :     }
595 :     }
596 :     }
597 :    
598 :    
599 :     #-------------------------------------------------------------------------------
600 :     # Fill a hash with key => value pairs based on a has or an array in \%options,
601 :     # or a supplied default list or hash:
602 :     #
603 :     # %hash = opt_hash( \%options, $opt_name, \@defaults )
604 :     # %hash = opt_hash( \%options, $opt_name, \%defaults )
605 :     #
606 :     #-------------------------------------------------------------------------------
607 :     sub opt_hash
608 :     {
609 :     my ( $options, $opt_name, $defaults ) = @_;
610 :     ( ref( $options ) eq 'HASH' && $opt_name ) or ref( $defaults ) or return ();
611 :    
612 :     $defaults = $options->{ $opt_name } if ref( $options ) && ref( $options->{ $opt_name } );
613 :    
614 :     ref( $defaults ) eq 'ARRAY' ? map { $_ => 1 } @$defaults :
615 :     ref( $defaults ) eq 'HASH' ? %$defaults :
616 :     ()
617 :     }
618 :    
619 :    
620 :     #-------------------------------------------------------------------------------
621 :     # ( $tmp_dir, $save_tmp ) = temporary_directory( \%options )
622 :     #-------------------------------------------------------------------------------
623 :     sub temporary_directory
624 :     {
625 :     my $options = ( shift ) || {};
626 :    
627 :     # Accept these option names with or without an underscore
628 :     my $tmp_dir = $options->{ tmpdir } || $options->{ tmp_dir };
629 :     my $save_tmp = $options->{ savetmp } || $options->{ save_tmp } || '';
630 :    
631 :     if ( $tmp_dir )
632 :     {
633 :     # User-supplied directory? Don't blow it away when we're done
634 :     if ( -d $tmp_dir ) { $options->{ savetmp } = $save_tmp = 1 }
635 :     }
636 :     else
637 :     {
638 :     my $tmp = $options->{ tmp } && -d $options->{ tmp } ? $options->{ tmp }
639 :     : $FIG_Config::temp && -d $FIG_Config::temp ? $FIG_Config::temp
640 :     : -d '/tmp' ? '/tmp'
641 :     : '.';
642 :     $tmp_dir = sprintf( "$tmp/special_proteins_tmp_dir.%05d.%09d", $$, int(1000000000*rand) );
643 :     # We named the directoy, let's add it in the options hash
644 :     $options->{ tmpdir } = $tmp_dir;
645 :     }
646 :    
647 :     if ( $tmp_dir && ! -d $tmp_dir )
648 :     {
649 :     mkdir $tmp_dir;
650 :     if ( ! -d $tmp_dir )
651 :     {
652 :     print STDERR "find_special_proteins::temporary_directory could not create '$tmp_dir'\n";
653 :     $options->{ tmpdir } = $tmp_dir = undef;
654 :     }
655 :     }
656 :    
657 :     return ( $tmp_dir, $save_tmp );
658 :     }
659 :    
660 : overbeek 1.10 #===============================================================================
661 :     # Use a database of reference sequences to locate proteins with unusual starts in a genome.
662 :     #
663 :     # @locs = find_odd_starts( \%params )
664 :     #
665 :     # Required parameters:
666 :     #
667 :     # contigs => \@contigs # genome is a synonym
668 :     #
669 :     # Optional paramaters:
670 :     #
671 :     # comment => $text # comment for output hash
672 :     # is_term => \@terminator_codons # D = [TAA,TAG,TGA]
673 :     # note => $text # same as 'comment'
674 :     # references => \@odd_start_proteins # user-supplied references
675 :     # tmp => $directory # directory for tmp_dir
676 :     # tmp_dir => $directory # directory for temporary files
677 :     #
678 :     # Some keys can be shortened.
679 :     #===============================================================================
680 :     sub find_odd_starts
681 :     {
682 :     my ( $params ) = @_;
683 :     my %locs;
684 :     my ( $tmp_dir, $save_tmp ) = temporary_directory( $params );
685 :    
686 :     my ( $contig_key ) = grep { /^contig/ || /^gen/ } keys %$params;
687 :     my $contigs = $contig_key ? $params->{ $contig_key } : undef;
688 :     ( ref( $contigs ) eq 'ARRAY' ) && @$contigs
689 :     or return ();
690 :    
691 :     my %contigR = map { $_->[0] => \$_->[2] } @$contigs; # Index sequence ref by id
692 :    
693 :     my ( $comment_key ) = grep { /^comment/ || /^note/ } keys %$params;
694 :     my $comment = $comment_key ? $params->{ $comment_key } : '';
695 :    
696 :     # Get reference sequences as parameter, or get them from a module:
697 :    
698 :     my ( $ref_key ) = grep { /^ref/ } keys %$params;
699 :     my $refs = $ref_key ? $params->{ $ref_key } : undef;
700 :     if ( $refs )
701 :     {
702 :     ( ref( $refs ) eq 'ARRAY' ) && @$refs
703 :     or print STDERR "No reference sequences supplied to find_odd_starts().\n"
704 :     and return ();
705 :     }
706 :     else
707 :     {
708 :     eval
709 :     {
710 :     require OddStarts_ref;
711 :     $refs = $OddStarts_ref::ref_seqs;
712 :     };
713 :     $refs && ( ref( $refs ) eq 'ARRAY' ) && @$refs
714 :     or print STDERR "Unable to get reference sequences from OddStarts_ref.pm.\n"
715 :     and return ();
716 :     }
717 :     my $suffix = $$ . "_" . sprintf( "%09d", int( 1e9 * rand() ) );
718 :    
719 :     # Print and format the contig sequences:
720 :    
721 :     my $contig_file = "$tmp_dir/contigs_$suffix";
722 :     gjoseqlib::print_alignment_as_fasta( $contig_file, $contigs );
723 :     -f $contig_file or return ();
724 :    
725 :     my $ext_bin = $FIG_Config::ext_bin;
726 :     my $formatdb = $params->{ formatdb };
727 :     $formatdb ||= $ext_bin ? "$ext_bin/formatdb" : 'formatdb';
728 :     system "cd '$tmp_dir'; $formatdb -p f -i 'contigs_$suffix'";
729 :    
730 :     # Make a clean copy of ref sequences to avoid BLAST warnings.
731 :    
732 :     $refs = [ map { my $s = uc $_->[2]; # uppercase copy
733 :     $s =~ s/[BJOUZ]/X/g; # U etc. -> X
734 :     $s =~ s/[^A-Z]//g; # delete nonletters
735 :     [ @$_[0,1], $s ]
736 :     } @$refs
737 :     ];
738 :     my $ref_file = "$tmp_dir/refs_$suffix";
739 :     gjoseqlib::print_alignment_as_fasta( $ref_file, $refs );
740 :     -f $ref_file or print STDERR "Failed to write reference sequences to $ref_file.\n"
741 :     and return ();
742 :    
743 :     # Search the contigs for the reference sequences:
744 :    
745 :     my $blastall = $params->{ blastall };
746 :     $blastall ||= $ext_bin ? "$ext_bin/blastall" : 'blastall';
747 :     my $blastcmd = "$blastall -p tblastn -d '$contig_file' -i '$ref_file' "
748 :     . "-e 1e-10 -F f -b 200 -v 200 -a 2 |";
749 :     open( BLAST, $blastcmd )
750 :     or print STDERR "Could not open pipe '$blastcmd'.\n"
751 :     and return ();
752 :     my @hsps = gjoparseblast::blast_hsp_list( \*BLAST );
753 :     close BLAST;
754 :    
755 :     # Delete the temporary directory unless it already existed:
756 :    
757 :     system "/bin/rm -r $tmp_dir" if ! $save_tmp;
758 :    
759 :     #
760 :     # Stucture of the output records:
761 :     #
762 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
763 :     # qid qdef qlen sid sdef slen scr e_val p_n p_val n_mat n_id n_pos n_gap dir q1 q2 qseq s1 s2 sseq
764 :     #
765 :    
766 :     @hsps = sort { $b->[6] <=> $a->[6] } @hsps;
767 :    
768 :     my $covered = {};
769 :     my %hit = ();
770 :     my %def = ();
771 :     my $hsp;
772 :    
773 :     foreach $hsp ( @hsps )
774 :     {
775 :    
776 :     if ((($hsp->[11] / $hsp->[10]) > 0.6) && (($hsp->[10] / $hsp->[2]) > 0.9))
777 :     {
778 :     my $contig = $hsp->[3];
779 :     my $contigP = $contigR{$contig};
780 :     my $beg = $hsp->[18];
781 :     my $end = $hsp->[19];
782 :     my ($predE,$partial) = &find_orf_end($contigP,$beg,$end,$params);
783 :     if ($predE && (! $partial))
784 :     {
785 :     # my($predS,$partial) = &find_orf_start($contigP,$begA,$end,{is_init => ['ATT']});
786 :     my($predS,$partial) = (($beg < $end) ? ($beg - (($hsp->[15] - 1) * 3)) :
787 :     ($beg + (($hsp->[15] - 1) * 3)),0);
788 :     if ($predS && (! $partial))
789 :     {
790 :     my $loc = join("_",($contig,$predS,$predE));
791 :     $locs{$loc}++;
792 :     }
793 :     }
794 :     }
795 :     else
796 :     {
797 :     # print STDERR &Dumper($hsp);
798 :     }
799 :     }
800 :     my(@locs) = sort keys(%locs);
801 :     wantarray ? @locs : \@locs;
802 :     }
803 : golsen 1.1
804 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3