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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package gjoalignandtree;
2 :    
3 : golsen 1.3 # This is a SAS component.
4 :    
5 :     #-------------------------------------------------------------------------------
6 :     #
7 :     # Order a set of sequences:
8 :     #
9 :     # @seqs = order_sequences_by_length( \@seqs, \%opts )
10 :     # \@seqs = order_sequences_by_length( \@seqs, \%opts )
11 :     #
12 :     # Write representative seqeunce groups to a file
13 :     #
14 :     # write_representative_groups( \%rep, $file )
15 : overbeek 1.1 #
16 : golsen 1.3 # @align = trim_align_to_median_ends( \@align, \%opts )
17 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
18 :     #
19 :     # print_alignment_as_pseudoclustal( $align, \%opts )
20 :     #
21 :     # $prof_rep = representative_for_profile( $align )
22 :     #
23 :     # \@trimmed_seq = extract_with_psiblast( $db, $profile, \%opts )
24 :     # ( \@trimmed_seq, \%report ) = extract_with_psiblast( $db, $profile, \%opts )
25 :     #
26 : golsen 1.8 # Profile blast against protein database:
27 :     #
28 : golsen 1.3 # $structured_blast = blastpgp( $dbfile, $profilefile, \%options )
29 :     # $structured_blast = blastpgp( \@dbseq, $profilefile, \%options )
30 :     # $structured_blast = blastpgp( $dbfile, \@profileseq, \%options )
31 :     # $structured_blast = blastpgp( \@dbseq, \@profileseq, \%options )
32 :     #
33 : golsen 1.8 # Profile blast against DNA database (PSI-tblastn):
34 :     #
35 :     # $structured_blast = blastpgpn( $dbfile, $profilefile, \%options )
36 :     # $structured_blast = blastpgpn( \@dbseq, $profilefile, \%options )
37 :     # $structured_blast = blastpgpn( $dbfile, \@profileseq, \%options )
38 :     # $structured_blast = blastpgpn( \@dbseq, \@profileseq, \%options )
39 :     #
40 : golsen 1.3 # print_blast_as_records( $file, \@queries )
41 :     # print_blast_as_records( \*FH, \@queries )
42 :     # print_blast_as_records( \@queries ) # STDOUT
43 :     #
44 :     # \@queries = read_blast_from_records( $file )
45 :     # \@queries = read_blast_from_records( \*FH )
46 :     # \@queries = read_blast_from_records( ) # STDIN
47 :     #
48 :     # $exists = verify_db( $db, $type );
49 :     #
50 :     #-------------------------------------------------------------------------------
51 : overbeek 1.1
52 :     use strict;
53 : golsen 1.3 use SeedAware;
54 : overbeek 1.1 use gjoseqlib;
55 : golsen 1.3 use gjoparseblast;
56 : overbeek 1.1
57 :     #-------------------------------------------------------------------------------
58 :     # Order a set of sequences:
59 :     #
60 :     # @seqs = order_sequences_by_length( \@seqs, \%opts )
61 :     # \@seqs = order_sequences_by_length( \@seqs, \%opts )
62 :     #
63 :     # Options:
64 :     #
65 : golsen 1.3 # order => $key # increasing (D), decreasing, median_outward,
66 :     # # closest_to_median, closest_to_length
67 :     # length => $prefered_length
68 : overbeek 1.1 #
69 :     #-------------------------------------------------------------------------------
70 :    
71 :     sub order_sequences_by_length
72 :     {
73 :     my ( $seqs, $opts ) = @_;
74 :     $seqs && ref $seqs eq 'ARRAY' && @$seqs
75 :     or print STDERR "order_sequences_by_length called with invalid sequence list.\n"
76 :     and return wantarray ? () : [];
77 :    
78 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
79 :     my $order = $opts->{ order } || 'increasing';
80 :    
81 :     my @seqs = ();
82 :     if ( $order =~ /^dec/i )
83 :     {
84 :     $opts->{ order } = 'decreaseing';
85 :     @seqs = sort { length( $b->[2] ) <=> length( $a->[2] ) } @$seqs;
86 :     }
87 :     elsif ( $order =~ /^med/i )
88 :     {
89 :     $opts->{ order } = 'median_outward';
90 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
91 :     while ( $_ = shift @seqs0 )
92 :     {
93 :     push @seqs, $_;
94 :     push @seqs, pop @seqs0 if @seqs0;
95 :     }
96 :     @seqs = reverse @seqs;
97 :     }
98 :     elsif ( $order =~ /^close/i )
99 :     {
100 :     my $pref_len;
101 :     if ( defined $opts->{ length } ) # caller supplies preferred length?
102 :     {
103 :     $opts->{ order } = 'closest_to_length';
104 :     $pref_len = $opts->{ length } + 0;
105 :     }
106 :     else # preferred length is median
107 :     {
108 :     $opts->{ order } = 'closest_to_median';
109 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
110 :     my $i = int( @seqs0 / 2 );
111 :     $pref_len = ( @seqs0 % 2 ) ? length( $seqs0[$i]->[2] )
112 :     : ( length( $seqs0[$i-1]->[2] ) + length( $seqs0[$i]->[2] ) ) / 2;
113 :     }
114 :    
115 :     @seqs = map { $_->[0] }
116 :     sort { $a->[1] <=> $b->[1] } # closest to preferred first
117 :     map { [ $_, abs( length( $_->[2] ) - $pref_len ) ] } # dist from pref?
118 :     @$seqs;
119 :     }
120 :     else
121 :     {
122 :     $opts->{ order } = 'increasing';
123 :     @seqs = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
124 :     }
125 :    
126 :     wantarray ? @seqs : \@seqs;
127 :     }
128 :    
129 :    
130 :     #-------------------------------------------------------------------------------
131 :     # Write representative seqeunce groups to a file
132 :     #
133 :     # write_representative_groups( \%rep, $file )
134 :     #
135 :     #-------------------------------------------------------------------------------
136 :    
137 :     sub write_representative_groups
138 :     {
139 :     my ( $repH, $file ) = @_;
140 :     $repH && (ref $repH eq 'HASH') && keys %$repH
141 :     or print STDERR "write_representative_groups called with invalid rep hash.\n"
142 :     and return;
143 :    
144 :     my ( $fh, $close ) = &output_file_handle( $file );
145 :    
146 :     # Order keys from largest to smallest set, then alphabetical
147 :    
148 :     my @keys = sort { @{ $repH->{$b} } <=> @{ $repH->{$a} } || lc $a cmp lc $b } keys %$repH;
149 :    
150 :     foreach ( @keys ) { print $fh join( "\t", ( $_, @{ $repH->{ $_ } } ) ), "\n" }
151 :    
152 :     close $fh if $close;
153 :     }
154 :    
155 :    
156 :     #-------------------------------------------------------------------------------
157 :     #
158 :     # @align = trim_align_to_median_ends( \@align, \%opts )
159 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
160 :     #
161 :     # Options:
162 :     #
163 : fangfang 1.5 # begin => bool # Trim start (specifically)
164 :     # end => bool # Trim end (specifically)
165 :     # fract_cov => fract # Fraction of sequences to be covered (D: 0.75)
166 : overbeek 1.1 #
167 :     #-------------------------------------------------------------------------------
168 :     sub trim_align_to_median_ends
169 :     {
170 :     my ( $align, $opts ) = @_;
171 : golsen 1.8
172 : overbeek 1.1 $align && ref $align eq 'ARRAY' && @$align
173 :     or print STDERR "trim_align_to_median_ends called with invalid sequence list.\n"
174 :     and return wantarray ? () : [];
175 :    
176 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
177 :     my $tr_beg = $opts->{ begin } || $opts->{ beg } || $opts->{ start } ? 1 : 0;
178 :     my $tr_end = $opts->{ end } ? 1 : 0;
179 :     $tr_beg = $tr_end = 1 if ! ( $tr_beg || $tr_end );
180 : fangfang 1.5 my $frac = $opts->{ fract_cov } || 0.75;
181 : overbeek 1.1
182 : fangfang 1.5 my( @ngap1, @ngap2);
183 : overbeek 1.1 foreach my $seq ( @$align )
184 :     {
185 :     my( $b, $e ) = $seq->[2] =~ /^(-*).*[^-](-*)$/;
186 : fangfang 1.5 push @ngap1, length( $b || '' );
187 :     push @ngap2, length( $e || '' );
188 : overbeek 1.1 }
189 :    
190 : fangfang 1.5 @ngap1 = sort { $a <=> $b } @ngap1;
191 :     @ngap2 = sort { $a <=> $b } @ngap2;
192 : overbeek 1.1
193 : fangfang 1.5 my $ngap1 = $tr_beg ? $ngap1[ int( @ngap1 * $frac ) ] : 0;
194 :     my $ngap2 = $tr_end ? $ngap2[ int( @ngap2 * $frac ) ] : 0;
195 : overbeek 1.1
196 :     my $ori_len = length( $align->[0]->[2] );
197 : fangfang 1.5 my $new_len = $ori_len - ( $ngap1 + $ngap2 );
198 :     my @align2 = map { [ @$_[0,1], substr( $_->[2], $ngap1, $new_len ) ] }
199 : overbeek 1.1 @$align;
200 :    
201 :     wantarray ? @align2 : \@align2;
202 :     }
203 :    
204 :    
205 :     #-------------------------------------------------------------------------------
206 :     #
207 : golsen 1.3 # print_alignment_as_pseudoclustal( $align, \%opts )
208 :     #
209 : overbeek 1.1 # Options:
210 :     #
211 :     # file => $filename # supply a file name to open and write
212 :     # file => \*FH # supply an open file handle (D = STDOUT)
213 :     # line => $linelen # residues per line (D = 60)
214 :     # lower => $bool # all lower case sequence
215 :     # upper => $bool # all upper case sequence
216 :     #
217 :     #-------------------------------------------------------------------------------
218 :    
219 :     sub print_alignment_as_pseudoclustal
220 :     {
221 :     my ( $align, $opts ) = @_;
222 :     $align && ref $align eq 'ARRAY' && @$align
223 :     or print STDERR "print_alignment_as_pseudoclustal called with invalid sequence list.\n"
224 :     and return wantarray ? () : [];
225 :    
226 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
227 :     my $line_len = $opts->{ line } || 60;
228 :     my $case = $opts->{ upper } ? 1 : $opts->{ lower } ? -1 : 0;
229 :    
230 :     my ( $fh, $close ) = &output_file_handle( $opts->{ file } );
231 :    
232 :     my $namelen = 0;
233 :     foreach ( @$align ) { $namelen = length $_->[0] if $namelen < length $_->[0] }
234 :     my $fmt = "%-${namelen}s %s\n";
235 :    
236 :     my $id;
237 :     my @lines = map { $id = $_->[0]; [ map { sprintf $fmt, $id, $_ }
238 : golsen 1.3 map { $case < 0 ? lc $_ : $case > 0 ? uc $_ : $_ } # map sequence only
239 : overbeek 1.1 $_->[2] =~ m/.{1,$line_len}/g
240 :     ] }
241 :    
242 :     @$align;
243 :    
244 :     my $ngroup = @{ $lines[0] };
245 :     for ( my $i = 0; $i < $ngroup; $i++ )
246 :     {
247 :     foreach ( @lines ) { print $fh $_->[$i] if $_->[$i] }
248 :     print $fh "\n";
249 :     }
250 :    
251 :     close $fh if $close;
252 :     }
253 :    
254 :    
255 :     #-------------------------------------------------------------------------------
256 : golsen 1.3 #
257 :     # @seqs = read_pseudoclustal( ) # D = STDIN
258 :     # \@seqs = read_pseudoclustal( ) # D = STDIN
259 :     # @seqs = read_pseudoclustal( $file_name )
260 :     # \@seqs = read_pseudoclustal( $file_name )
261 :     # @seqs = read_pseudoclustal( \*FH )
262 :     # \@seqs = read_pseudoclustal( \*FH )
263 :     #
264 :     #-------------------------------------------------------------------------------
265 :    
266 :     sub read_pseudoclustal
267 :     {
268 :     my ( $file ) = @_;
269 :     my ( $fh, $close ) = input_file_handle( $file );
270 :     my %seq;
271 :     my @ids;
272 :     while ( <$fh> )
273 :     {
274 :     chomp;
275 :     my ( $id, $data ) = /^(\S+)\s+(\S.*)$/;
276 :     if ( defined $id && defined $data )
277 :     {
278 :     push @ids, $id if ! $seq{ $id };
279 :     $data =~ s/\s+//g;
280 :     push @{ $seq{ $id } }, $data;
281 :     }
282 :     }
283 :     close $fh if $close;
284 :    
285 :     my @seq = map { [ $_, '', join( '', @{ $seq{ $_ } } ) ] } @ids;
286 :     wantarray ? @seq : \@seq;
287 :     }
288 :    
289 :    
290 :     #-------------------------------------------------------------------------------
291 : overbeek 1.1 # The profile 'query' sequence:
292 :     # 1. No terminal gaps, if possible
293 :     # 1.1 Otherwise, no gaps at start
294 :     # 2. Longest sequence passing above
295 :     #
296 :     # $prof_rep = representative_for_profile( $align )
297 :     #-------------------------------------------------------------------------------
298 :     sub representative_for_profile
299 :     {
300 :     my ( $align ) = @_;
301 :     $align && ref $align eq 'ARRAY' && @$align
302 :     or die "representative_for_profile called with invalid sequence list.\n";
303 :    
304 :     my @cand;
305 :     @cand = grep { $_->[2] =~ /^[^-].*[^-]$/ } @$align; # No end gaps
306 :     @cand = grep { $_->[2] =~ /^[^-]/ } @$align if ! @cand; # No start gaps
307 : fangfang 1.4 @cand = @$align if ! @cand;
308 : overbeek 1.1
309 : golsen 1.3 my ( $rep ) = map { $_->[0] } # sequence entry
310 :     sort { $b->[1] <=> $a->[1] } # max nongaps
311 :     map { [ $_, $_->[2] =~ tr/-//c ] } # count nongaps
312 : overbeek 1.1 @cand;
313 :    
314 :     $rep = [ @$rep ]; # Make a copy
315 :     $rep->[2] =~ s/-+//g; # Compress sequence
316 :    
317 :     return $rep;
318 :     }
319 :    
320 :    
321 :     #-------------------------------------------------------------------------------
322 :     #
323 : golsen 1.3 # \@seq = extract_with_psiblast( $db, $profile, \%opts )
324 :     # ( \@seq, \%report ) = extract_with_psiblast( $db, $profile, \%opts )
325 :     #
326 :     # $db can be $db_file_name or \@db_seqs
327 :     # $profile can be $profile_file_name or \@profile_seqs
328 :     #
329 :     # If supplied as file, profile is pseudoclustal
330 :     #
331 :     # Report records:
332 :     #
333 :     # [ $sid, $slen, $status, $frac_id, $frac_pos,
334 :     # $q_uncov_n_term, $q_uncov_c_term,
335 :     # $s_uncov_n_term, $s_uncov_c_term ]
336 :     #
337 :     # Options:
338 :     #
339 :     # e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
340 : golsen 1.8 # max_e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
341 : golsen 1.3 # max_q_uncov => $aa_per_end # maximum unmatched query (D = 20)
342 :     # max_q_uncov_c => $aa_per_end # maximum unmatched query, c-term (D = 20)
343 :     # max_q_uncov_n => $aa_per_end # maximum unmatched query, n-term (D = 20)
344 :     # min_ident => $frac_ident # minimum fraction identity (D = 0.15)
345 :     # min_positive => $frac_positive # minimum fraction positive scoring (D = 0.20)
346 : golsen 1.8 # n_result => $max_results # maxiumn restuls returned (D = 1000)
347 :     # n_thread => $n_thread # number of blastpgp threads (D = 2)
348 : golsen 1.3 # query => $q_file_name # query sequence file (D = most complete)
349 :     # query => \@q_seq_entry # query sequence (D = most complete)
350 :     # stderr => $file # blastpgp stderr (D = /dev/stderr)
351 :     #
352 :     # If supplied, query must be identical to a profile sequence, but no gaps.
353 :     #
354 :     #-------------------------------------------------------------------------------
355 :    
356 :     sub extract_with_psiblast
357 :     {
358 : fangfang 1.6 my ( $seqs, $profile, $opts ) = @_;
359 : golsen 1.3
360 :     $opts && ref $opts eq 'HASH' or $opts = {};
361 :    
362 : golsen 1.8 $opts->{ e_value } ||= $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
363 :     $opts->{ n_result } ||= 1000;
364 :     $opts->{ n_thread } ||= 2;
365 : golsen 1.3 my $max_q_uncov_c = $opts->{ max_q_uncov_c } || $opts->{ max_q_uncov } || 20;
366 :     my $max_q_uncov_n = $opts->{ max_q_uncov_n } || $opts->{ max_q_uncov } || 20;
367 :     my $min_ident = $opts->{ min_ident } || 0.15;
368 :     my $min_pos = $opts->{ min_positive } || 0.20;
369 :    
370 : fangfang 1.6 my $blast = blastpgp( $seqs, $profile, $opts );
371 : golsen 1.3 $blast && @$blast or return ();
372 :    
373 :     my ( $qid, $qdef, $qlen, $qhits ) = @{ $blast->[0] };
374 :    
375 :     my @trimmed;
376 :     my %report;
377 :    
378 :     foreach my $sdata ( @$qhits )
379 :     {
380 :     my ( $sid, $sdef, $slen, $hsps ) = @$sdata;
381 :    
382 :     if ( one_real_hsp( $hsps ) )
383 :     {
384 :     # [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
385 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
386 :    
387 :     my $hsp0 = $hsps->[0];
388 :     my ( $scr, $nmat, $nid, $npos, $q1, $q2, $qseq, $s1, $s2, $sseq ) = ( @$hsp0 )[ 0, 4, 5, 6, 9, 10, 11, 12, 13, 14 ];
389 :    
390 :     my $status;
391 :     if ( $q1-1 > $max_q_uncov_n ) { $status = 'missing start' }
392 :     elsif ( $qlen-$q2 > $max_q_uncov_c ) { $status = 'missing end' }
393 :     elsif ( $nid / $nmat < $min_ident ) { $status = 'low identity' }
394 :     elsif ( $npos / $nmat < $min_pos ) { $status = 'low positives' }
395 :     else
396 :     {
397 :     $sseq =~ s/-+//g;
398 :     $sdef .= " ($s1-$s2/$slen)" if ( ( $s1 > 1 ) || ( $s2 < $slen ) );
399 :     push @trimmed, [ $sid, $sdef, $sseq ];
400 :     $status = 'included';
401 :     }
402 :    
403 :     $report{ $sid } = [ $sid, $slen, $status, $nid/$nmat, $npos/$nmat, $q1-1, $qlen-$q2, $s1-1, $slen-$s2 ];
404 :     }
405 :     }
406 :    
407 :     wantarray ? ( \@trimmed, \%report ) : \@trimmed;
408 :     }
409 :    
410 :    
411 :     #-------------------------------------------------------------------------------
412 :     #
413 :     # Allow fragmentary matches inside of the highest-scoring hsp:
414 :     #
415 :     #-------------------------------------------------------------------------------
416 :    
417 :     sub one_real_hsp
418 :     {
419 :     my ( $hsps ) = @_;
420 :     return 0 if ! ( $hsps && ( ref( $hsps ) eq 'ARRAY' ) && @$hsps );
421 :     return 1 if @$hsps == 1;
422 :    
423 :     my ( $q1_0, $q2_0 ) = ( @{ $hsps->[0] } )[9, 10];
424 :     for ( my $i = 1; $i < @$hsps; $i++ )
425 :     {
426 :     my ($q1, $q2) = ( @{ $hsps->[$i] } )[9, 10];
427 :     return 0 if $q1 < $q1_0 || $q2 > $q2_0;
428 :     }
429 :    
430 :     return 1;
431 :     }
432 :    
433 :    
434 :     #-------------------------------------------------------------------------------
435 :     #
436 : golsen 1.8 # $structured_blast = blastpgp( $dbfile, $profilefile, \%options )
437 :     # $structured_blast = blastpgp( \@dbseq, $profilefile, \%options )
438 :     # $structured_blast = blastpgp( $dbfile, \@profileseq, \%options )
439 :     # $structured_blast = blastpgp( \@dbseq, \@profileseq, \%options )
440 : overbeek 1.1 #
441 :     # Required:
442 :     #
443 : golsen 1.3 # $db_file or \@db_seq
444 :     # $profile_file or \@profile_seq
445 :     #
446 :     # Options:
447 :     #
448 : golsen 1.8 # e_value => $max_e_value # maximum E-value of matches (D = 0.01)
449 :     # max_e_value => $max_e_value # maximum E-value of matches (D = 0.01)
450 :     # n_result => $max_results # maximum matches returned (D = 1000)
451 :     # n_thread => $n_thread # number of blastpgp threads (D = 2)
452 :     # stderr => $file # place to send blast stderr (D = /dev/null)
453 :     # query => $q_file_name # most complete
454 :     # query => \@q_seq_entry # most complete
455 : overbeek 1.1 #
456 :     #-------------------------------------------------------------------------------
457 :    
458 :     sub blastpgp
459 :     {
460 :     my ( $db, $profile, $opts ) = @_;
461 :     $opts ||= {};
462 : golsen 1.8 my $tmp = SeedAware::location_of_tmp( $opts );
463 : overbeek 1.1
464 :     my ( $dbfile, $rm_db );
465 :     if ( defined $db && ref $db )
466 :     {
467 :     ref $db eq 'ARRAY'
468 :     && @$db
469 :     || print STDERR "blastpgp requires one or more database sequences.\n"
470 :     && return undef;
471 : golsen 1.8 $dbfile = SeedAware::new_file_name( "$tmp/blastpgp_db" );
472 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $dbfile, $db );
473 :     $rm_db = 1;
474 :     }
475 :     elsif ( defined $db && -f $db )
476 :     {
477 :     $dbfile = $db;
478 :     }
479 :     else
480 :     {
481 :     die "blastpgp requires database.";
482 :     }
483 : golsen 1.3 verify_db( $dbfile, 'P' ); # protein
484 : overbeek 1.1
485 :     my ( $proffile, $rm_profile );
486 :     if ( defined $profile && ref $profile )
487 :     {
488 :     ref $profile eq 'ARRAY'
489 :     && @$profile
490 :     || print STDERR "blastpgp requires one or more profile sequences.\n"
491 :     && return undef;
492 : golsen 1.8 $proffile = SeedAware::new_file_name( "$tmp/blastpgp_profile" );
493 : overbeek 1.1 &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
494 :     $rm_profile = 1;
495 :     }
496 :     elsif ( defined $profile && -f $profile )
497 :     {
498 :     $proffile = $profile;
499 :     }
500 :     else
501 :     {
502 :     die "blastpgp requires profile.";
503 :     }
504 :    
505 :     my ( $qfile, $rm_query );
506 :     my $query = $opts->{ query };
507 :     if ( defined $query && ref $query )
508 :     {
509 : golsen 1.3 ref $query eq 'ARRAY' && @$query == 3
510 :     or print STDERR "blastpgp invalid query sequence.\n"
511 :     and return undef;
512 :    
513 : golsen 1.8 $qfile = SeedAware::new_file_name( "$tmp/blastpgp_query" );
514 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
515 :     $rm_query = 1;
516 :     }
517 :     elsif ( defined $query && -f $query )
518 :     {
519 :     $qfile = $query;
520 :     }
521 :     elsif ( $profile ) # Build it from profile
522 :     {
523 :     $query = &representative_for_profile( $profile );
524 : golsen 1.8 $qfile = SeedAware::new_file_name( "$tmp/blastpgp_query" );
525 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
526 :     $rm_query = 1;
527 :     }
528 :     else
529 : golsen 1.3 {
530 : overbeek 1.1 die "blastpgp requires database.";
531 :     }
532 :    
533 : golsen 1.8 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
534 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 2;
535 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
536 :    
537 : golsen 1.3 my @cmd = ( SeedAware::executable_for( 'blastpgp' ),
538 : golsen 1.8 '-j', 1, # function is profile alignment
539 :     '-B', $proffile, # location of profile
540 : golsen 1.3 '-d', $dbfile,
541 : golsen 1.8 '-i', $qfile,
542 :     '-e', $e_val,
543 : golsen 1.3 '-b', $nkeep,
544 :     '-v', $nkeep,
545 : golsen 1.8 '-t', 1 # issues warning if not set
546 : golsen 1.3 );
547 : golsen 1.8 push @cmd, ( '-a', $n_cpu ) if $n_cpu;
548 : golsen 1.3
549 : golsen 1.8 my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
550 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
551 : golsen 1.3 or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
552 :     and return undef;
553 :     my $out = &gjoparseblast::structured_blast_output( $blastfh, 1 ); # selfmatches okay
554 :     close $blastfh;
555 : overbeek 1.1
556 :     if ( $rm_db )
557 :     {
558 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
559 :     unlink @files if @files;
560 :     }
561 :     unlink $proffile if $rm_profile;
562 :     unlink $qfile if $rm_query;
563 :    
564 :     return $out;
565 :     }
566 :    
567 :    
568 :     #-------------------------------------------------------------------------------
569 : golsen 1.8 # Do psiblast against tranlated genomic DNA
570 :     #
571 :     # $structured_blast = blastpgpn( $dbfile, $profilefile, \%options )
572 :     # $structured_blast = blastpgpn( \@dbseq, $profilefile, \%options )
573 :     # $structured_blast = blastpgpn( $dbfile, \@profileseq, \%options )
574 :     # $structured_blast = blastpgpn( \@dbseq, \@profileseq, \%options )
575 :     #
576 :     # Required:
577 :     #
578 :     # $db_file or \@db_seq
579 :     # $profile_file or \@profile_seq
580 :     #
581 :     # Options:
582 :     #
583 :     # aa_db_file => $trans_file # put translation db here
584 :     # e_value => $max_e_value # D = 0.01
585 :     # max_e_val => $max_e_value # D = 0.01
586 :     # n_cpu => $n_cpu # synonym of n_thread
587 :     # n_result => $max_seq # D = 1000
588 :     # n_thread => $n_cpu # D = 1
589 :     # query => $q_file_name # most complete
590 :     # query => \@q_seq_entry # most complete
591 :     # stderr => $file # place to send blast stderr (D = /dev/null)
592 :     # tmp => $temp_dir # location for temp files
593 :     #
594 :     # depricated alternatives:
595 :     #
596 :     # prot_file => $trans_file # put translation db here
597 :     #-------------------------------------------------------------------------------
598 :    
599 :     sub blastpgpn
600 :     {
601 :     my ( $ndb, $profile, $opts ) = @_;
602 :     $opts ||= {};
603 :     my $tmp = SeedAware::location_of_tmp( $opts );
604 :    
605 :     $opts->{ aa_db_file } ||= $opts->{ prot_file } if defined $opts->{ prot_file };
606 :     my $dbfile = $opts->{ aa_db_file } || SeedAware::new_file_name( "$tmp/blastpgpn_db" );
607 :     my $rm_db = ! $opts->{ aa_db_file };
608 :     if ( defined $dbfile && -f $dbfile && -s $dbfile )
609 :     {
610 :     # The tranaslated sequence database exists
611 :     }
612 :     elsif ( defined $ndb )
613 :     {
614 :     if ( ref $ndb eq 'ARRAY' && @$ndb )
615 :     {
616 :     ref $ndb eq 'ARRAY'
617 :     && @$ndb
618 :     || print STDERR "Bad sequence reference passed to blastpgpn.\n"
619 :     && return undef;
620 :     my @pdb = map { six_translations( $_ ) } @$ndb;
621 :     gjoseqlib::print_alignment_as_fasta( $dbfile, \@pdb );
622 :     }
623 :     elsif ( -f $ndb && -s $ndb )
624 :     {
625 :     open( PDB, ">$dbfile" )
626 :     or print STDERR "Could not open protein database file '$dbfile'.\n"
627 :     and return undef;
628 :     my $entry;
629 :     while ( $entry = gjoseqlib::next_fasta_entry( $ndb ) )
630 :     {
631 :     gjoseqlib::write_alignment_as_fasta( \*PDB, six_translations( $entry ) );
632 :     }
633 :     close PDB;
634 :     }
635 :     else
636 :     {
637 :     die "blastpgpn requires a sequence database.";
638 :     }
639 :     }
640 :     verify_db( $dbfile, 'P' ); # protein
641 :    
642 :     my ( $proffile, $rm_profile );
643 :     if ( defined $profile && ref $profile )
644 :     {
645 :     ref $profile eq 'ARRAY'
646 :     && @$profile
647 :     || print STDERR "blastpgpn requires one or more profile sequences.\n"
648 :     && return undef;
649 :     $proffile = SeedAware::new_file_name( "$tmp/blastpgpn_profile" );
650 :     &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
651 :     $rm_profile = 1;
652 :     }
653 :     elsif ( defined $profile && -f $profile )
654 :     {
655 :     $proffile = $profile;
656 :     }
657 :     else
658 :     {
659 :     die "blastpgpn requires profile.";
660 :     }
661 :    
662 :     my ( $qfile, $rm_query );
663 :     my $query = $opts->{ query };
664 :     if ( defined $query && ref $query )
665 :     {
666 :     ref $query eq 'ARRAY' && @$query == 3
667 :     or print STDERR "blastpgpn invalid query sequence.\n"
668 :     and return undef;
669 :    
670 :     $qfile = SeedAware::new_file_name( "$tmp/blastpgpn_query" );
671 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
672 :     $rm_query = 1;
673 :     }
674 :     elsif ( defined $query && -f $query )
675 :     {
676 :     $qfile = $query;
677 :     }
678 :     elsif ( $profile ) # Build it from profile
679 :     {
680 :     $query = &representative_for_profile( $profile );
681 :     $qfile = SeedAware::new_file_name( "$tmp/blastpgpn_query" );
682 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
683 :     $rm_query = 1;
684 :     }
685 :     else
686 :     {
687 :     die "blastpgpn requires database.";
688 :     }
689 :    
690 :     my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
691 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 2;
692 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
693 :    
694 :     my @cmd = ( SeedAware::executable_for( 'blastpgp' ),
695 :     '-j', 1, # function is profile alignment
696 :     '-B', $proffile, # location of protein profile
697 :     '-d', $dbfile, # protein database
698 :     '-i', $qfile, # one protein from the profile
699 :     '-e', $e_val,
700 :     '-b', $nkeep,
701 :     '-v', $nkeep,
702 :     '-t', 1, # issues warning if not set
703 :     );
704 :     push @cmd, ( '-a', $n_cpu ) if $n_cpu;
705 :    
706 :     my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
707 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
708 :     or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
709 :     and return undef;
710 :     my $out = &gjoparseblast::structured_blast_output( $blastfh );
711 :     close $blastfh;
712 :    
713 :     if ( $rm_db )
714 :     {
715 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
716 :     unlink @files if @files;
717 :     }
718 :     unlink $proffile if $rm_profile;
719 :     unlink $qfile if $rm_query;
720 :    
721 :     # Fix the blastp output to look like tblastn output
722 :    
723 :     foreach my $qdata ( @$out )
724 :     {
725 :     my %sdata; # There are now multiple sequences per subject DNA
726 :     foreach my $sdata ( @{ $qdata->[3] } )
727 :     {
728 :     my ( $sid, $sdef, undef, $hsps ) = @$sdata;
729 :     my $fr;
730 :     ( $sid, $fr ) = $sid =~ m/^(.*)\.([-+]\d)$/;
731 :     my ( $b, $e, $slen ) = $sdef =~ m/(\d+)-(\d+)\/(\d+)$/;
732 :     $sdef =~ s/ ?\S+$//;
733 :     my $sd2 = $sdata{ $sid } ||= [ $sid, $sdef, $slen, [] ];
734 :     foreach ( @$hsps ) { adjust_hsp( $_, $fr, $b ); push @{$sd2->[3]}, $_ }
735 :     }
736 :    
737 :     # Order the hsps for each subject
738 :    
739 :     foreach my $sid ( keys %sdata )
740 :     {
741 :     my $hsps = $sdata{ $sid }->[3];
742 :     @$hsps = sort { $b->[0] <=> $a->[0] } @$hsps;
743 :     }
744 :    
745 :     # Order the subjects for the query
746 :    
747 :     @{$qdata->[3]} = map { $_->[0] } # remove score
748 :     sort { $b->[1] <=> $a->[1] } # sort by score
749 :     map { [ $_, $_->[3]->[0]->[0] ] } # tag with top score
750 :     map { $sdata{ $_ } } # subject data
751 :     keys %sdata; # subject ids
752 :     }
753 :    
754 :     return $out;
755 :     }
756 :    
757 :    
758 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
759 :     # When search is versus six frame translation, there is a need to adjust
760 :     # the frame and location information in the hsp back to the DNA coordinates.
761 :     #
762 :     # adjust_hsp( $hsp, $frame, $begin )
763 :     #
764 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
765 :     # scr Eval nseg Eval naln nid npos ngap fr q1 q2 qseq s1 s2 sseq
766 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
767 :     sub adjust_hsp
768 :     {
769 :     my ( $hsp, $fr, $b ) = @_;
770 :     $hsp->[8] = $fr;
771 :     if ( $fr > 0 )
772 :     {
773 :     $hsp->[12] = $b + 3 * ( $hsp->[12] - 1 );
774 :     $hsp->[13] = $b + 3 * ( $hsp->[13] - 1 ) + 2;
775 :     }
776 :     else
777 :     {
778 :     $hsp->[12] = $b - 3 * ( $hsp->[12] - 1 );
779 :     $hsp->[13] = $b - 3 * ( $hsp->[13] - 1 ) - 2;
780 :     }
781 :     }
782 :    
783 :    
784 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
785 :     # Do a six frame translation for use by blastpgpn. This modifications of
786 :     # the identifiers and definitions are essential to the interpretation of
787 :     # the blast results. The program 'translate_fasta_6' produces the same
788 :     # output format, and is much faster.
789 :     #
790 :     # @translations = six_translations( $nucleotide_entry )
791 :     #
792 :     # The ids are modified by adding ".frame" (+1, +2, +3, -1, -2, -3).
793 :     # The definition is mofidified by adding " begin-end/of_length".
794 :     # NCBI frame numbers reverse strand translation frames from the end of the
795 :     # sequence (i.e., the beginning of the complement of the strand).
796 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
797 :    
798 :     sub six_translations
799 :     {
800 :     my ( $id, $def, $seq ) = map { defined($_) ? $_ : '' } @{ $_[0] };
801 :     my $l = length( $seq );
802 :    
803 :     return () if $l < 15;
804 :    
805 :     # fr beg end
806 :     my @intervals = ( [ '+1', 1, $l - ( $l % 3 ) ],
807 :     [ '+2', 2, $l - ( ($l-1) % 3 ) ],
808 :     [ '+3', 3, $l - ( ($l-2) % 3 ) ],
809 :     [ '-1', $l, 1 + ( $l % 3 ) ],
810 :     [ '-2', $l-1, 1 + ( ($l-1) % 3 ) ],
811 :     [ '-3', $l-2, 1 + ( ($l-2) % 3 ) ]
812 :     );
813 :     my ( $fr, $b, $e );
814 :    
815 :     map { ( $fr, $b, $e ) = @$_;
816 :     [ "$id.$fr", "$def $b-$e/$l", gjoseqlib::translate_seq( gjoseqlib::DNA_subseq( \$seq, $b, $e ) ) ]
817 :     } @intervals;
818 :     }
819 :    
820 :    
821 :     #-------------------------------------------------------------------------------
822 : golsen 1.3 # Get an input file handle, and boolean on whether to close or not:
823 :     #
824 :     # ( \*FH, $close ) = input_file_handle( $filename );
825 :     # ( \*FH, $close ) = input_file_handle( \*FH );
826 :     # ( \*FH, $close ) = input_file_handle( ); # D = STDIN
827 :     #
828 :     #-------------------------------------------------------------------------------
829 :    
830 :     sub input_file_handle
831 :     {
832 :     my ( $file, $umask ) = @_;
833 :    
834 :     my ( $fh, $close );
835 :     if ( defined $file )
836 :     {
837 :     if ( ref $file eq 'GLOB' )
838 :     {
839 :     $fh = $file;
840 :     $close = 0;
841 :     }
842 :     elsif ( -f $file )
843 :     {
844 :     open( $fh, "<$file") || die "input_file_handle could not open '$file'.\n";
845 :     $close = 1;
846 :     }
847 :     else
848 :     {
849 :     die "input_file_handle could not find file '$file'.\n";
850 :     }
851 :     }
852 :     else
853 :     {
854 :     $fh = \*STDIN;
855 :     $close = 0;
856 :     }
857 :    
858 :     return ( $fh, $close );
859 :     }
860 :    
861 :    
862 :     #-------------------------------------------------------------------------------
863 : overbeek 1.1 # Get an output file handle, and boolean on whether to close or not:
864 :     #
865 :     # ( \*FH, $close ) = output_file_handle( $filename );
866 :     # ( \*FH, $close ) = output_file_handle( \*FH );
867 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
868 :     #
869 :     #-------------------------------------------------------------------------------
870 :    
871 :     sub output_file_handle
872 :     {
873 :     my ( $file, $umask ) = @_;
874 :    
875 :     my ( $fh, $close );
876 :     if ( defined $file )
877 :     {
878 :     if ( ref $file eq 'GLOB' )
879 :     {
880 :     $fh = $file;
881 :     $close = 0;
882 :     }
883 :     else
884 :     {
885 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
886 :     chmod 0664, $file; # Seems to work on open file!
887 :     $close = 1;
888 :     }
889 :     }
890 :     else
891 :     {
892 :     $fh = \*STDOUT;
893 :     $close = 0;
894 :     }
895 :    
896 :     return ( $fh, $close );
897 :     }
898 :    
899 :    
900 :     #-------------------------------------------------------------------------------
901 :     #
902 :     # print_blast_as_records( $file, \@queries )
903 :     # print_blast_as_records( \*FH, \@queries )
904 :     # print_blast_as_records( \@queries ) # STDOUT
905 :     #
906 :     # \@queries = read_blast_from_records( $file )
907 :     # \@queries = read_blast_from_records( \*FH )
908 :     # \@queries = read_blast_from_records( ) # STDIN
909 :     #
910 :     # Three output record types:
911 :     #
912 :     # Query \t qid \t qdef \t dlen
913 :     # > \t sid \t sdef \t slen
914 :     # HSP \t ...
915 :     #
916 :     #-------------------------------------------------------------------------------
917 :     sub print_blast_as_records
918 :     {
919 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
920 :     my ( $fh, $close ) = output_file_handle( $file );
921 :    
922 :     my $queries = shift;
923 :     $queries && ref $queries eq 'ARRAY'
924 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
925 :     and return;
926 :    
927 :     foreach my $qdata ( @$queries )
928 :     {
929 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
930 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
931 :     and return;
932 :    
933 :     my $subjcts = pop @$qdata;
934 :     $subjcts && ref $subjcts eq 'ARRAY'
935 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
936 :     and return;
937 :    
938 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
939 :     foreach my $sdata ( @$subjcts )
940 :     {
941 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
942 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
943 :     and return;
944 :    
945 :     my $hsps = pop @$sdata;
946 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
947 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
948 :     and return;
949 :    
950 :     print $fh join( "\t", '>', @$sdata ), "\n";
951 :    
952 :     foreach my $hsp ( @$hsps )
953 :     {
954 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
955 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
956 :     and return;
957 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
958 :     }
959 :     }
960 :     }
961 :    
962 :     close $fh if $close;
963 :     }
964 :    
965 :    
966 :     sub read_blast_from_records
967 :     {
968 :     my ( $file ) = @_;
969 : golsen 1.3
970 :     my ( $fh, $close ) = input_file_handle( $file );
971 :     $fh or print STDERR "read_blast_from_records could not open input file.\n"
972 : overbeek 1.1 and return wantarray ? () : [];
973 : golsen 1.3
974 : overbeek 1.1 my @queries = ();
975 :     my @qdata;
976 :     my @subjects = ();
977 :     my @sdata;
978 :     my @hsps = ();
979 :    
980 :     local $_;
981 :     while ( defined( $_ = <$fh> ) )
982 :     {
983 :     chomp;
984 :     my ( $type, @datum ) = split /\t/;
985 :     if ( $type eq 'Query' )
986 :     {
987 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
988 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
989 :     @qdata = @datum;
990 :     @sdata = ();
991 :     @subjects = ();
992 :     @hsps = ();
993 :     }
994 :     elsif ( $type eq '>' )
995 :     {
996 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
997 :     @sdata = @datum;
998 :     @hsps = ();
999 :     }
1000 :     elsif ( $type eq 'HSP' )
1001 :     {
1002 :     push @hsps, \@datum;
1003 :     }
1004 :     }
1005 :    
1006 :     close $fh if $close;
1007 :    
1008 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1009 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1010 :    
1011 :     wantarray ? @queries : \@queries;
1012 :     }
1013 :    
1014 :    
1015 : golsen 1.3 sub verify_db
1016 :     {
1017 :     my ( $db, $type ) = @_;
1018 :    
1019 :     my @args;
1020 : golsen 1.8 if ( $type =~ m/^p/i )
1021 : golsen 1.3 {
1022 : golsen 1.8 @args = ( "-p", "T", "-i", $db ) if ((! -s "$db.psq") || (-M "$db.psq" > -M $db));
1023 : golsen 1.3 }
1024 :     else
1025 :     {
1026 : golsen 1.8 @args = ( "-p", "F", "-i", $db ) if ((! -s "$db.nsq") || (-M "$db.nsq" > -M $db));
1027 : golsen 1.3 }
1028 : golsen 1.8 @args or return ( -s $db ? 1 : 0 );
1029 :    
1030 : golsen 1.3 #
1031 : golsen 1.8 # Find formatdb appropriate for the excecution environemnt.
1032 : golsen 1.3 #
1033 :    
1034 :     my $prog = SeedAware::executable_for( 'formatdb' );
1035 :     if ( ! $prog )
1036 :     {
1037 :     warn "gjoalignandtree::verify_db: formatdb program not found.\n";
1038 :     return 0;
1039 :     }
1040 :    
1041 :     my $rc = system( $prog, @args );
1042 :    
1043 :     if ( $rc != 0 )
1044 :     {
1045 :     my $cmd = join( ' ', $prog, @args );
1046 :     warn "gjoalignandtree::verify_db: formatdb failed with rc = $rc: $cmd\n";
1047 :     return 0;
1048 :     }
1049 :    
1050 :     return 1;
1051 :     }
1052 :    
1053 :    
1054 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3