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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.20 - (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 : golsen 1.17 local $_;
92 :     while ( defined( $_ = shift @seqs0 ) )
93 : overbeek 1.1 {
94 :     push @seqs, $_;
95 :     push @seqs, pop @seqs0 if @seqs0;
96 :     }
97 :     @seqs = reverse @seqs;
98 :     }
99 :     elsif ( $order =~ /^close/i )
100 :     {
101 :     my $pref_len;
102 :     if ( defined $opts->{ length } ) # caller supplies preferred length?
103 :     {
104 :     $opts->{ order } = 'closest_to_length';
105 :     $pref_len = $opts->{ length } + 0;
106 :     }
107 :     else # preferred length is median
108 :     {
109 :     $opts->{ order } = 'closest_to_median';
110 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
111 :     my $i = int( @seqs0 / 2 );
112 :     $pref_len = ( @seqs0 % 2 ) ? length( $seqs0[$i]->[2] )
113 :     : ( length( $seqs0[$i-1]->[2] ) + length( $seqs0[$i]->[2] ) ) / 2;
114 :     }
115 :    
116 :     @seqs = map { $_->[0] }
117 :     sort { $a->[1] <=> $b->[1] } # closest to preferred first
118 :     map { [ $_, abs( length( $_->[2] ) - $pref_len ) ] } # dist from pref?
119 :     @$seqs;
120 :     }
121 :     else
122 :     {
123 :     $opts->{ order } = 'increasing';
124 :     @seqs = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
125 :     }
126 :    
127 :     wantarray ? @seqs : \@seqs;
128 :     }
129 :    
130 :    
131 :     #-------------------------------------------------------------------------------
132 :     # Write representative seqeunce groups to a file
133 :     #
134 :     # write_representative_groups( \%rep, $file )
135 :     #
136 :     #-------------------------------------------------------------------------------
137 :    
138 :     sub write_representative_groups
139 :     {
140 :     my ( $repH, $file ) = @_;
141 :     $repH && (ref $repH eq 'HASH') && keys %$repH
142 :     or print STDERR "write_representative_groups called with invalid rep hash.\n"
143 :     and return;
144 :    
145 :     my ( $fh, $close ) = &output_file_handle( $file );
146 :    
147 :     # Order keys from largest to smallest set, then alphabetical
148 :    
149 :     my @keys = sort { @{ $repH->{$b} } <=> @{ $repH->{$a} } || lc $a cmp lc $b } keys %$repH;
150 :    
151 :     foreach ( @keys ) { print $fh join( "\t", ( $_, @{ $repH->{ $_ } } ) ), "\n" }
152 :    
153 :     close $fh if $close;
154 :     }
155 :    
156 :    
157 :     #-------------------------------------------------------------------------------
158 :     #
159 :     # @align = trim_align_to_median_ends( \@align, \%opts )
160 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
161 :     #
162 :     # Options:
163 :     #
164 : fangfang 1.5 # begin => bool # Trim start (specifically)
165 :     # end => bool # Trim end (specifically)
166 :     # fract_cov => fract # Fraction of sequences to be covered (D: 0.75)
167 : overbeek 1.1 #
168 :     #-------------------------------------------------------------------------------
169 :     sub trim_align_to_median_ends
170 :     {
171 :     my ( $align, $opts ) = @_;
172 : golsen 1.8
173 : overbeek 1.1 $align && ref $align eq 'ARRAY' && @$align
174 :     or print STDERR "trim_align_to_median_ends called with invalid sequence list.\n"
175 :     and return wantarray ? () : [];
176 :    
177 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
178 :     my $tr_beg = $opts->{ begin } || $opts->{ beg } || $opts->{ start } ? 1 : 0;
179 :     my $tr_end = $opts->{ end } ? 1 : 0;
180 :     $tr_beg = $tr_end = 1 if ! ( $tr_beg || $tr_end );
181 : fangfang 1.5 my $frac = $opts->{ fract_cov } || 0.75;
182 : overbeek 1.1
183 : fangfang 1.5 my( @ngap1, @ngap2);
184 : overbeek 1.1 foreach my $seq ( @$align )
185 :     {
186 :     my( $b, $e ) = $seq->[2] =~ /^(-*).*[^-](-*)$/;
187 : fangfang 1.5 push @ngap1, length( $b || '' );
188 :     push @ngap2, length( $e || '' );
189 : overbeek 1.1 }
190 :    
191 : fangfang 1.5 @ngap1 = sort { $a <=> $b } @ngap1;
192 :     @ngap2 = sort { $a <=> $b } @ngap2;
193 : overbeek 1.1
194 : fangfang 1.5 my $ngap1 = $tr_beg ? $ngap1[ int( @ngap1 * $frac ) ] : 0;
195 :     my $ngap2 = $tr_end ? $ngap2[ int( @ngap2 * $frac ) ] : 0;
196 : overbeek 1.1
197 :     my $ori_len = length( $align->[0]->[2] );
198 : fangfang 1.5 my $new_len = $ori_len - ( $ngap1 + $ngap2 );
199 :     my @align2 = map { [ @$_[0,1], substr( $_->[2], $ngap1, $new_len ) ] }
200 : overbeek 1.1 @$align;
201 :    
202 :     wantarray ? @align2 : \@align2;
203 :     }
204 :    
205 :    
206 :     #-------------------------------------------------------------------------------
207 :     #
208 : golsen 1.3 # print_alignment_as_pseudoclustal( $align, \%opts )
209 :     #
210 : overbeek 1.1 # Options:
211 :     #
212 :     # file => $filename # supply a file name to open and write
213 :     # file => \*FH # supply an open file handle (D = STDOUT)
214 :     # line => $linelen # residues per line (D = 60)
215 :     # lower => $bool # all lower case sequence
216 :     # upper => $bool # all upper case sequence
217 :     #
218 :     #-------------------------------------------------------------------------------
219 :    
220 :     sub print_alignment_as_pseudoclustal
221 :     {
222 :     my ( $align, $opts ) = @_;
223 :     $align && ref $align eq 'ARRAY' && @$align
224 :     or print STDERR "print_alignment_as_pseudoclustal called with invalid sequence list.\n"
225 :     and return wantarray ? () : [];
226 :    
227 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
228 :     my $line_len = $opts->{ line } || 60;
229 :     my $case = $opts->{ upper } ? 1 : $opts->{ lower } ? -1 : 0;
230 :    
231 :     my ( $fh, $close ) = &output_file_handle( $opts->{ file } );
232 :    
233 :     my $namelen = 0;
234 :     foreach ( @$align ) { $namelen = length $_->[0] if $namelen < length $_->[0] }
235 :     my $fmt = "%-${namelen}s %s\n";
236 :    
237 :     my $id;
238 :     my @lines = map { $id = $_->[0]; [ map { sprintf $fmt, $id, $_ }
239 : golsen 1.3 map { $case < 0 ? lc $_ : $case > 0 ? uc $_ : $_ } # map sequence only
240 : overbeek 1.1 $_->[2] =~ m/.{1,$line_len}/g
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 : golsen 1.20 #
293 :     # 1. Minimum terminal gaps
294 : overbeek 1.1 # 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 : golsen 1.20 my ( $r0 ) = map { $_->[0] } # sequence entry
305 :     sort { $a->[1] <=> $b->[1] || $b->[2] <=> $a->[2] } # min terminal gaps, max aas
306 :     map { my $tgap = ( $_->[2] =~ /^(-+)/ ? length( $1 ) : 0 )
307 :     + ( $_->[2] =~ /(-+)$/ ? length( $1 ) : 0 );
308 :     my $naa = $_->[2] =~ tr/ACDEFGHIKLMNPQRSTVWYacdefghiklmnpqrstvwy//;
309 :     [ $_, $tgap, $naa ]
310 :     }
311 :     @$align;
312 : overbeek 1.1
313 : golsen 1.20 my $rep = [ @$r0 ]; # Make a copy
314 :     $rep->[2] =~ s/[^A-Za-z]+//g; # Compress to letters
315 : overbeek 1.1
316 : golsen 1.20 $rep;
317 : overbeek 1.1 }
318 :    
319 :    
320 :     #-------------------------------------------------------------------------------
321 : golsen 1.20 # Use a profile to extract a region from a sequence. By default, the hsps
322 :     # can be chained together, walking both directions from the highest scoring
323 :     # hsp. Only one region will be extracted per database sequence.
324 : overbeek 1.1 #
325 : golsen 1.3 # \@seq = extract_with_psiblast( $db, $profile, \%opts )
326 :     # ( \@seq, \%report ) = extract_with_psiblast( $db, $profile, \%opts )
327 :     #
328 :     # $db can be $db_file_name or \@db_seqs
329 :     # $profile can be $profile_file_name or \@profile_seqs
330 :     #
331 :     # If supplied as file, profile is pseudoclustal
332 :     #
333 :     # Report records:
334 :     #
335 : fangfang 1.13 # [ $sid, $score, $e_value, $slen, $status,
336 :     # $frac_id, $frac_pos,
337 : golsen 1.3 # $q_uncov_n_term, $q_uncov_c_term,
338 :     # $s_uncov_n_term, $s_uncov_c_term ]
339 :     #
340 :     # Options:
341 :     #
342 : golsen 1.20 # best_hsp => $bool # only report the best hsp (no merging)
343 : golsen 1.3 # e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
344 : golsen 1.8 # max_e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
345 : golsen 1.20 # max_gap => $max_gap_size # maximum region unmatched in both query and subject (D = 100)
346 :     # max_offset => $max_offset # maximum asymmetry in query and subject gap size (D = 1000)
347 :     # # (i.e., very large insertions are allowed)
348 :     # max_q_uncov => $aa_per_end # maximum unmatched query, per end (D = 20)
349 :     # max_q_uncov_c => $aa_per_end # maximum unmatched query, C-term (D = 20)
350 :     # max_q_uncov_n => $aa_per_end # maximum unmatched query, N-term (D = 20)
351 : golsen 1.3 # min_ident => $frac_ident # minimum fraction identity (D = 0.15)
352 :     # min_positive => $frac_positive # minimum fraction positive scoring (D = 0.20)
353 : fangfang 1.11 # min_frac_cov => $frac_cov # minimum fraction coverage of query and subject sequence (D = 0.20)
354 : fangfang 1.10 # min_q_cov => $frac_cov # minimum fraction coverage of query sequence (D = 0.50)
355 : fangfang 1.11 # min_s_cov => $frac_cov # minimum fraction coverage of subject sequence (D = 0.01)
356 : golsen 1.8 # n_result => $max_results # maxiumn restuls returned (D = 1000)
357 : golsen 1.16 # n_cpu => $n_thread # number of blastpgp threads (D = 2)
358 : golsen 1.8 # n_thread => $n_thread # number of blastpgp threads (D = 2)
359 : golsen 1.20 # one_hsp => $bool # only report the best hsp (no merging)
360 : golsen 1.3 # query => $q_file_name # query sequence file (D = most complete)
361 :     # query => \@q_seq_entry # query sequence (D = most complete)
362 : golsen 1.20 # query_used => \@q_seq_entry # output of the query sequence used
363 : golsen 1.3 # stderr => $file # blastpgp stderr (D = /dev/stderr)
364 : golsen 1.20 # strict => $bool # only report matches that are unique
365 :     # # and continuous (one hsp)
366 : golsen 1.3 #
367 : golsen 1.20 # If supplied, the query must be identical to a profile sequence, but with no gaps.
368 : golsen 1.3 #-------------------------------------------------------------------------------
369 :    
370 :     sub extract_with_psiblast
371 :     {
372 : fangfang 1.6 my ( $seqs, $profile, $opts ) = @_;
373 : golsen 1.3
374 :     $opts && ref $opts eq 'HASH' or $opts = {};
375 :    
376 : golsen 1.20 # These 3 options set parameters for blastpgp:
377 :    
378 :     $opts->{ e_value } ||= $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
379 :     $opts->{ n_result } ||= $opts->{ nresult } || 1000;
380 :     $opts->{ n_thread } ||= $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 4;
381 :    
382 :     # These options set values used in this routine:
383 :    
384 : fangfang 1.11 my $max_q_uncov_c = $opts->{ max_q_uncov_c } || $opts->{ max_q_uncov } || 20;
385 :     my $max_q_uncov_n = $opts->{ max_q_uncov_n } || $opts->{ max_q_uncov } || 20;
386 :     my $min_q_cov = $opts->{ min_q_cov } || $opts->{ min_frac_cov } || 0.50;
387 :     my $min_s_cov = $opts->{ min_s_cov } || $opts->{ min_frac_cov } || 0.01;
388 : golsen 1.20 my $min_ident = $opts->{ min_ident } || $opts->{ min_identity } || 0.15;
389 :     my $min_pos = $opts->{ min_positive } || 0.20;
390 :     my $strict = $opts->{ strict } || 0;
391 :     my $best_hsp = $opts->{ best_hsp } || $opts->{ one_hsp } || 0;
392 : golsen 1.3
393 : fangfang 1.6 my $blast = blastpgp( $seqs, $profile, $opts );
394 : golsen 1.20 $blast && @$blast
395 :     or return wantarray ? ( [], {} ) : [];
396 : golsen 1.3
397 :     my ( $qid, $qdef, $qlen, $qhits ) = @{ $blast->[0] };
398 :    
399 :     my @trimmed;
400 : golsen 1.20 my %seqs;
401 : golsen 1.3 my %report;
402 :    
403 :     foreach my $sdata ( @$qhits )
404 :     {
405 :     my ( $sid, $sdef, $slen, $hsps ) = @$sdata;
406 :    
407 : golsen 1.20 my $hsp0;
408 :     if ( $hsp0 = $strict ? one_real_hsp( $hsps ) : $best_hsp ? $hsps->[0] : pseudo_hsp( $hsps, $opts ) )
409 : golsen 1.3 {
410 :     # [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
411 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
412 :    
413 : golsen 1.20 my ( $scr, $exp, $nmat, $nid, $npos, $q1, $q2, $s1, $s2, $sseq ) = @$hsp0[ 0, 1, 4, 5, 6, 9, 10, 12, 13, 14 ];
414 : golsen 1.3
415 :     my $status;
416 : fangfang 1.10 if ( $q1-1 > $max_q_uncov_n ) { $status = 'missing start' }
417 :     elsif ( $qlen-$q2 > $max_q_uncov_c ) { $status = 'missing end' }
418 :     elsif ( $nid / $nmat < $min_ident ) { $status = 'low identity' }
419 :     elsif ( $npos / $nmat < $min_pos ) { $status = 'low positives' }
420 :     elsif ( ($q2-$q1+1)/$qlen < $min_q_cov ) { $status = 'low coverage' }
421 :     elsif ( ($s2-$s1+1)/$slen < $min_s_cov ) { $status = 'long subject' }
422 : golsen 1.3 else
423 :     {
424 : golsen 1.20 if ( $sseq )
425 :     {
426 :     $sseq =~ s/-+//g;
427 :     }
428 :     else
429 :     {
430 :     if ( ! keys %seqs )
431 :     {
432 :     if ( ref( $seqs ) eq 'ARRAY' )
433 :     {
434 :     %seqs = map { $_->[0] => $_ } @$seqs;
435 :     }
436 :     elsif ( -f $seqs )
437 :     {
438 :     %seqs = map { $_->[0] => $_ } gjoseqlib::read_fasta( $seqs );
439 :     }
440 :     else
441 :     {
442 :     print STDERR "extract_with_psiblast could not locate sequences.\n";
443 :     return wantarray ? ( [], {} ) : [];
444 :     }
445 :     }
446 :     $sseq = substr( $seqs{$sid}->[2] || '', $s1, $s2-$s1+1 );
447 :     }
448 : golsen 1.3 $sdef .= " ($s1-$s2/$slen)" if ( ( $s1 > 1 ) || ( $s2 < $slen ) );
449 :     push @trimmed, [ $sid, $sdef, $sseq ];
450 :     $status = 'included';
451 :     }
452 :    
453 : fangfang 1.13 my $frac_id = sprintf("%.3f", $nid/$nmat);
454 :     my $frac_pos = sprintf("%.3f", $npos/$nmat);
455 :    
456 :     $report{ $sid } = [ $sid, $scr, $exp, $slen, $status, $frac_id, $frac_pos, $q1-1, $qlen-$q2, $s1-1, $slen-$s2 ];
457 : golsen 1.3 }
458 :     }
459 :    
460 :     wantarray ? ( \@trimmed, \%report ) : \@trimmed;
461 :     }
462 :    
463 :    
464 :     #-------------------------------------------------------------------------------
465 : golsen 1.20 # Allow fragmentary matches inside of the highest-scoring hsp (or extending
466 :     # at most 10 amino acids beyond):
467 : golsen 1.3 #
468 : golsen 1.20 # $hsp = one_real_hsp( \@hsps )
469 : golsen 1.3 #
470 :     #-------------------------------------------------------------------------------
471 :    
472 :     sub one_real_hsp
473 :     {
474 :     my ( $hsps ) = @_;
475 : golsen 1.20 return undef if ! ( $hsps && ( ref( $hsps ) eq 'ARRAY' ) && @$hsps );
476 :     return $hsps->[0] if @$hsps == 1;
477 :    
478 :     my $hsp0 = $hsps->[0];
479 :     my $q1_min = $hsp0->[ 9] - 10;
480 :     my $q2_max = $hsp0->[10] + 10;
481 :     foreach my $hsp ( @$hsps[ 1 .. (@$hsps-1) ] )
482 :     {
483 :     return undef if ( $hsp->[ 9] < $q1_min ) || ( $hsp->[10] > $q2_max );
484 :     }
485 :    
486 :     $hsp0;
487 :     }
488 :    
489 :    
490 :     #-------------------------------------------------------------------------------
491 :     # Assemble a pseudo HSP from discontinuous BLAST HSPs:
492 :     #
493 :     # $hsp = pseudo_hsp( \@hsps, \%options )
494 :     #
495 :     #-------------------------------------------------------------------------------
496 :     sub pseudo_hsp
497 :     {
498 :     my ( $hsps, $options ) = @_;
499 :     return undef if ! ( $hsps && ( ref( $hsps ) eq 'ARRAY' ) && @$hsps );
500 :     return $hsps->[0] if @$hsps == 1;
501 :     $options ||= {};
502 :     my $max_gap = $options->{ max_gap } || 100; # Maximum unmatched region
503 :     my $max_off = $options->{ max_off } || $options->{ max_offset } || 1000; # Maximum offset
504 : golsen 1.3
505 : golsen 1.20 my @hsps = sort { s_mid( $a ) <=> s_mid( $b ) } @$hsps;
506 :     my $i0 = 0;
507 :     my $scr = 0;
508 :     my $s;
509 :     for ( my $i = 0; $i < @hsps; $i++ )
510 : golsen 1.3 {
511 : golsen 1.20 if ( ( $s = $hsps[$i]->[0] ) > $scr ) { $i0 = $i; $scr = $s }
512 : golsen 1.3 }
513 :    
514 : golsen 1.20 my $hsp0 = [ @{$hsps[$i0]} ];
515 :    
516 :     for ( my $i = $i0 - 1; $i >= 0; $i-- )
517 :     {
518 :     my ( $q1, $q2, $s1, $s2 ) = @$hsp0[ 9, 10, 12, 13 ];
519 :     my $hsp = $hsps[$i];
520 :     next unless $hsp->[9] < $q1 && $hsp->[12] < $s1;
521 :     my $q_gap = $q1 - $hsp->[10] - 1;
522 :     my $s_gap = $s1 - $hsp->[13] - 1;
523 :     last if ( $q_gap > $max_gap ) && ( $s_gap > $max_gap );
524 :     last if ( abs( $q_gap - $s_gap ) > $max_off );
525 :    
526 :     # It is not worth trying to maintain the subject sequence
527 :    
528 :     $hsp0->[14] = '';
529 :    
530 :     # If there is overlap, prorate the scores
531 :    
532 :     my $factor = 1;
533 :     my $gap = $q_gap < $s_gap ? $q_gap : $s_gap;
534 :     if ( $q_gap < 0 || $s_gap < 0 )
535 :     {
536 :     if ( $q_gap < $s_gap )
537 :     {
538 :     my $len = $hsp->[10] - $hsp->[9] + 1;
539 :     $factor = ( $len + $q_gap ) / $len;
540 :     }
541 :     else
542 :     {
543 :     my $len = $hsp->[13] - $hsp->[12] + 1;
544 :     $factor = ( $len + $s_gap ) / $len;
545 :     }
546 :     }
547 :    
548 :     $hsp0->[ 0] += $factor * $hsp->[0];
549 :     $hsp0->[ 2] += $factor * $hsp->[2];
550 :     $hsp0->[ 3] += $factor * $hsp->[3];
551 :     $hsp0->[ 4] += $factor * $hsp->[4];
552 :     $hsp0->[ 9] = $hsp->[ 9];
553 :     $hsp0->[12] = $hsp->[12];
554 :     }
555 :    
556 :     for ( my $i = $i0 + 1; $i < @hsps; $i++ )
557 :     {
558 :     my ( $q1, $q2, $s1, $s2 ) = @$hsp0[ 9, 10, 12, 13 ];
559 :     my $hsp = $hsps[$i];
560 :     next unless $hsp->[10] > $q2 && $hsp->[13] > $s2;
561 :     my $q_gap = $hsp->[ 9] - $q2 - 1;
562 :     my $s_gap = $hsp->[12] - $s2 - 1;
563 :     last if ( $q_gap > $max_gap ) && ( $s_gap > $max_gap );
564 :     last if ( abs( $q_gap - $s_gap ) > $max_off );
565 :    
566 :     # It is not worth trying to maintain the subject sequence
567 :    
568 :     $hsp0->[14] = '';
569 :    
570 :     # If there is overlap, prorate the scores
571 :    
572 :     my $factor = 1;
573 :     my $gap = $q_gap < $s_gap ? $q_gap : $s_gap;
574 :     if ( $q_gap < 0 || $s_gap < 0 )
575 :     {
576 :     if ( $q_gap < $s_gap )
577 :     {
578 :     my $len = $hsp->[10] - $hsp->[9] + 1;
579 :     $factor = ( $len + $q_gap ) / $len;
580 :     }
581 :     else
582 :     {
583 :     my $len = $hsp->[13] - $hsp->[12] + 1;
584 :     $factor = ( $len + $s_gap ) / $len;
585 :     }
586 :     }
587 :    
588 :     my $bitscr = $factor * $hsp->[0];
589 :     $hsp0->[ 0] += $bitscr;
590 :     $hsp0->[ 1] *= 0.5 ** $bitscr;
591 :     $hsp0->[ 2] += $factor * $hsp->[2];
592 :     $hsp0->[ 3] += $factor * $hsp->[3];
593 :     $hsp0->[ 4] += $factor * $hsp->[4];
594 :     $hsp0->[10] = $hsp->[10];
595 :     $hsp0->[13] = $hsp->[13];
596 :     }
597 :    
598 :     $hsp0->[0] = int( $hsp0->[0] );
599 :     $hsp0->[2] = int( $hsp0->[2] );
600 :     $hsp0->[3] = int( $hsp0->[3] );
601 :     $hsp0->[4] = int( $hsp0->[4] );
602 :    
603 :     $hsp0;
604 : golsen 1.3 }
605 :    
606 :    
607 : golsen 1.20 sub s_mid { $_[0]->[12] + $_[0]->[13] }
608 :    
609 :    
610 : golsen 1.3 #-------------------------------------------------------------------------------
611 : golsen 1.20 # Do a PSI-BLAST search:
612 : golsen 1.3 #
613 : golsen 1.8 # $structured_blast = blastpgp( $dbfile, $profilefile, \%options )
614 :     # $structured_blast = blastpgp( \@dbseq, $profilefile, \%options )
615 :     # $structured_blast = blastpgp( $dbfile, \@profileseq, \%options )
616 :     # $structured_blast = blastpgp( \@dbseq, \@profileseq, \%options )
617 : overbeek 1.1 #
618 :     # Required:
619 :     #
620 : golsen 1.3 # $db_file or \@db_seq
621 :     # $profile_file or \@profile_seq
622 :     #
623 :     # Options:
624 :     #
625 : golsen 1.8 # e_value => $max_e_value # maximum E-value of matches (D = 0.01)
626 :     # max_e_value => $max_e_value # maximum E-value of matches (D = 0.01)
627 :     # n_result => $max_results # maximum matches returned (D = 1000)
628 : golsen 1.20 # n_thread => $n_thread # number of blastpgp threads (D = 4)
629 : golsen 1.8 # stderr => $file # place to send blast stderr (D = /dev/null)
630 :     # query => $q_file_name # most complete
631 :     # query => \@q_seq_entry # most complete
632 : fangfang 1.14 # query_used => \@q_seq_entry # output
633 : overbeek 1.1 #
634 :     #-------------------------------------------------------------------------------
635 :    
636 :     sub blastpgp
637 :     {
638 :     my ( $db, $profile, $opts ) = @_;
639 :     $opts ||= {};
640 : golsen 1.8 my $tmp = SeedAware::location_of_tmp( $opts );
641 : overbeek 1.1
642 :     my ( $dbfile, $rm_db );
643 :     if ( defined $db && ref $db )
644 :     {
645 :     ref $db eq 'ARRAY'
646 :     && @$db
647 :     || print STDERR "blastpgp requires one or more database sequences.\n"
648 :     && return undef;
649 : golsen 1.20 my $dbfh;
650 :     ( $dbfh, $dbfile ) = SeedAware::open_tmp_file( "blastpgp_db", '', $tmp );
651 :     gjoseqlib::write_fasta( $dbfh, $db );
652 :     close( $dbfh );
653 : overbeek 1.1 $rm_db = 1;
654 :     }
655 :     elsif ( defined $db && -f $db )
656 :     {
657 :     $dbfile = $db;
658 :     }
659 :     else
660 :     {
661 :     die "blastpgp requires database.";
662 :     }
663 : golsen 1.3 verify_db( $dbfile, 'P' ); # protein
664 : overbeek 1.1
665 :     my ( $proffile, $rm_profile );
666 :     if ( defined $profile && ref $profile )
667 :     {
668 :     ref $profile eq 'ARRAY'
669 :     && @$profile
670 : golsen 1.20 or print STDERR "blastpgp requires one or more profile sequences.\n"
671 :     and return undef;
672 :     my $proffh;
673 :     ( $proffh, $proffile ) = SeedAware::open_tmp_file( "blastpgp_profile", '', $tmp );
674 :     # Uppercase is critical to the behavior of blastpgp
675 :     &print_alignment_as_pseudoclustal( $profile, { file => $proffh, upper => 1 } );
676 :     close( $proffh );
677 : overbeek 1.1 $rm_profile = 1;
678 :     }
679 :     elsif ( defined $profile && -f $profile )
680 :     {
681 :     $proffile = $profile;
682 :     }
683 :     else
684 :     {
685 : golsen 1.20 print STDERR "blastpgp requires a profile."
686 :     and return undef;
687 : overbeek 1.1 }
688 :    
689 :     my ( $qfile, $rm_query );
690 :     my $query = $opts->{ query };
691 :     if ( defined $query && ref $query )
692 :     {
693 : golsen 1.3 ref $query eq 'ARRAY' && @$query == 3
694 :     or print STDERR "blastpgp invalid query sequence.\n"
695 :     and return undef;
696 :    
697 : golsen 1.20 my $qfh;
698 :     ( $qfh, $qfile ) = SeedAware::open_tmp_file( "blastpgp_query", '', $tmp );
699 :     gjoseqlib::write_fasta( $qfh, [$query] );
700 :     close( $qfh );
701 : overbeek 1.1 $rm_query = 1;
702 :     }
703 :     elsif ( defined $query && -f $query )
704 :     {
705 :     $qfile = $query;
706 : golsen 1.16 ( $query ) = gjoseqlib::read_fasta( $qfile );
707 : golsen 1.20 $query
708 :     or print STDERR "blastpgp failed to get query from '$qfile'."
709 :     and return undef;
710 : overbeek 1.1 }
711 :     elsif ( $profile ) # Build it from profile
712 :     {
713 : golsen 1.20 $query = &representative_for_profile( $profile )
714 :     or print STDERR "blastpgp failed to get query from the profile."
715 :     and return undef;
716 :     my $qfh;
717 :     ( $qfh, $qfile ) = SeedAware::open_tmp_file( "blastpgp_query", '', $tmp );
718 :     gjoseqlib::write_fasta( $qfh, [$query] );
719 :     close( $qfh );
720 : overbeek 1.1 $rm_query = 1;
721 :     }
722 :     else
723 : golsen 1.3 {
724 : golsen 1.20 print STDERR "blastpgp requires a query."
725 :     and return undef;
726 : overbeek 1.1 }
727 :    
728 : fangfang 1.15 $opts->{ query_used } = $query;
729 : fangfang 1.13
730 : golsen 1.20 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
731 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 4;
732 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
733 : fangfang 1.13
734 : golsen 1.9 my $blastpgp = SeedAware::executable_for( 'blastpgp' )
735 :     or print STDERR "Could not find executable for program 'blastpgp'.\n"
736 :     and return undef;
737 :    
738 :     my @cmd = ( $blastpgp,
739 : golsen 1.8 '-j', 1, # function is profile alignment
740 :     '-B', $proffile, # location of profile
741 : golsen 1.3 '-d', $dbfile,
742 : golsen 1.8 '-i', $qfile,
743 :     '-e', $e_val,
744 : golsen 1.3 '-b', $nkeep,
745 :     '-v', $nkeep,
746 : golsen 1.8 '-t', 1 # issues warning if not set
747 : golsen 1.3 );
748 : golsen 1.8 push @cmd, ( '-a', $n_cpu ) if $n_cpu;
749 : golsen 1.3
750 : golsen 1.8 my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
751 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
752 : golsen 1.3 or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
753 :     and return undef;
754 :     my $out = &gjoparseblast::structured_blast_output( $blastfh, 1 ); # selfmatches okay
755 :     close $blastfh;
756 : overbeek 1.1
757 :     if ( $rm_db )
758 :     {
759 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
760 :     unlink @files if @files;
761 :     }
762 :     unlink $proffile if $rm_profile;
763 :     unlink $qfile if $rm_query;
764 :    
765 :     return $out;
766 :     }
767 :    
768 :    
769 :     #-------------------------------------------------------------------------------
770 : golsen 1.8 # Do psiblast against tranlated genomic DNA
771 :     #
772 :     # $structured_blast = blastpgpn( $dbfile, $profilefile, \%options )
773 :     # $structured_blast = blastpgpn( \@dbseq, $profilefile, \%options )
774 :     # $structured_blast = blastpgpn( $dbfile, \@profileseq, \%options )
775 :     # $structured_blast = blastpgpn( \@dbseq, \@profileseq, \%options )
776 :     #
777 :     # Required:
778 :     #
779 :     # $db_file or \@db_seq
780 :     # $profile_file or \@profile_seq
781 :     #
782 :     # Options:
783 :     #
784 :     # aa_db_file => $trans_file # put translation db here
785 :     # e_value => $max_e_value # D = 0.01
786 :     # max_e_val => $max_e_value # D = 0.01
787 :     # n_cpu => $n_cpu # synonym of n_thread
788 :     # n_result => $max_seq # D = 1000
789 : golsen 1.20 # n_thread => $n_cpu # D = 4
790 : golsen 1.8 # query => $q_file_name # most complete
791 :     # query => \@q_seq_entry # most complete
792 : fangfang 1.14 # query_used => \@q_seq_entry # output
793 : golsen 1.8 # stderr => $file # place to send blast stderr (D = /dev/null)
794 :     # tmp => $temp_dir # location for temp files
795 :     #
796 :     # depricated alternatives:
797 :     #
798 :     # prot_file => $trans_file # put translation db here
799 :     #-------------------------------------------------------------------------------
800 :    
801 :     sub blastpgpn
802 :     {
803 :     my ( $ndb, $profile, $opts ) = @_;
804 :     $opts ||= {};
805 :     my $tmp = SeedAware::location_of_tmp( $opts );
806 :    
807 :     $opts->{ aa_db_file } ||= $opts->{ prot_file } if defined $opts->{ prot_file };
808 : golsen 1.20 my $dbfile = $opts->{ aa_db_file };
809 :     my $rm_db = ! ( $dbfile && -f $dbfile );
810 : golsen 1.8 if ( defined $dbfile && -f $dbfile && -s $dbfile )
811 :     {
812 :     # The tranaslated sequence database exists
813 :     }
814 :     elsif ( defined $ndb )
815 :     {
816 :     if ( ref $ndb eq 'ARRAY' && @$ndb )
817 :     {
818 :     ref $ndb eq 'ARRAY'
819 :     && @$ndb
820 :     || print STDERR "Bad sequence reference passed to blastpgpn.\n"
821 :     && return undef;
822 : golsen 1.20 my $dbfh;
823 :     if ( $dbfile )
824 :     {
825 :     open( $dbfh, '>', $dbfile );
826 :     }
827 :     else
828 :     {
829 :     ( $dbfh, $dbfile ) = SeedAware::open_tmp_file( "blastpgpn_db", '', $tmp );
830 :     }
831 :     $dbfh or print STDERR 'Could not open $dbfile.'
832 :     and return undef;
833 :     foreach ( @$ndb )
834 :     {
835 :     gjoseqlib::write_alignment_as_fasta( $dbfh, six_translations( $_ ) );
836 :     }
837 :     close( $dbfh );
838 : golsen 1.8 }
839 :     elsif ( -f $ndb && -s $ndb )
840 :     {
841 : golsen 1.20 my $dbfh;
842 :     open( $dbfh, '>', $dbfile )
843 : golsen 1.8 or print STDERR "Could not open protein database file '$dbfile'.\n"
844 :     and return undef;
845 :     my $entry;
846 :     while ( $entry = gjoseqlib::next_fasta_entry( $ndb ) )
847 :     {
848 : golsen 1.20 gjoseqlib::write_alignment_as_fasta( $dbfh, six_translations( $entry ) );
849 : golsen 1.8 }
850 : golsen 1.20 close( $dbfh );
851 : golsen 1.8 }
852 :     else
853 :     {
854 : golsen 1.20 print STDERR "blastpgpn requires a sequence database."
855 :     and return undef;
856 : golsen 1.8 }
857 :     }
858 : golsen 1.20 else
859 :     {
860 :     die "blastpgpn requires a sequence database.";
861 :     }
862 : golsen 1.8 verify_db( $dbfile, 'P' ); # protein
863 :    
864 :     my ( $proffile, $rm_profile );
865 :     if ( defined $profile && ref $profile )
866 :     {
867 :     ref $profile eq 'ARRAY'
868 :     && @$profile
869 :     || print STDERR "blastpgpn requires one or more profile sequences.\n"
870 :     && return undef;
871 : golsen 1.20 my $proffh;
872 :     ( $proffh, $proffile ) = SeedAware::open_tmp_file( "blastpgpn_profile", '', $tmp );
873 :     &print_alignment_as_pseudoclustal( $proffh, { file => $proffile, upper => 1 } );
874 :     close( $proffh );
875 : golsen 1.8 $rm_profile = 1;
876 :     }
877 :     elsif ( defined $profile && -f $profile )
878 :     {
879 :     $proffile = $profile;
880 :     }
881 :     else
882 :     {
883 : golsen 1.20 print STDERR "blastpgpn requires profile."
884 :     and return undef;
885 : golsen 1.8 }
886 :    
887 :     my ( $qfile, $rm_query );
888 :     my $query = $opts->{ query };
889 :     if ( defined $query && ref $query )
890 :     {
891 :     ref $query eq 'ARRAY' && @$query == 3
892 :     or print STDERR "blastpgpn invalid query sequence.\n"
893 :     and return undef;
894 :    
895 : golsen 1.20 my $qfh;
896 :     ( $qfh, $qfile ) = SeedAware::open_tmp_file( "blastpgpn_query", '', $tmp );
897 :     gjoseqlib::write_fasta( $qfh, [$query] );
898 :     close( $qfh );
899 : golsen 1.8 $rm_query = 1;
900 :     }
901 :     elsif ( defined $query && -f $query )
902 :     {
903 :     $qfile = $query;
904 : golsen 1.16 ( $query ) = gjoseqlib::read_fasta( $qfile );
905 : golsen 1.8 }
906 :     elsif ( $profile ) # Build it from profile
907 :     {
908 :     $query = &representative_for_profile( $profile );
909 : golsen 1.20 my $qfh;
910 :     ( $qfh, $qfile ) = SeedAware::open_tmp_file( "blastpgpn_query", '', $tmp );
911 :     gjoseqlib::write_fasta( $qfh, [$query] );
912 :     close( $qfh );
913 : golsen 1.8 $rm_query = 1;
914 :     }
915 :     else
916 :     {
917 : golsen 1.20 print STDERR "blastpgpn requires a query."
918 :     and return undef;
919 : golsen 1.8 }
920 :    
921 : fangfang 1.15 $opts->{ query_used } = $query;
922 : fangfang 1.13
923 : golsen 1.20 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
924 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 4;
925 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
926 : golsen 1.8
927 : golsen 1.9 my $blastpgp = SeedAware::executable_for( 'blastpgp' )
928 :     or print STDERR "Could not find executable for program 'blastpgp'.\n"
929 :     and return undef;
930 :    
931 :     my @cmd = ( $blastpgp,
932 : golsen 1.8 '-j', 1, # function is profile alignment
933 :     '-B', $proffile, # location of protein profile
934 :     '-d', $dbfile, # protein database
935 :     '-i', $qfile, # one protein from the profile
936 :     '-e', $e_val,
937 :     '-b', $nkeep,
938 :     '-v', $nkeep,
939 :     '-t', 1, # issues warning if not set
940 :     );
941 :     push @cmd, ( '-a', $n_cpu ) if $n_cpu;
942 :    
943 :     my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
944 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
945 :     or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
946 :     and return undef;
947 :     my $out = &gjoparseblast::structured_blast_output( $blastfh );
948 :     close $blastfh;
949 :    
950 :     if ( $rm_db )
951 :     {
952 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
953 :     unlink @files if @files;
954 :     }
955 :     unlink $proffile if $rm_profile;
956 :     unlink $qfile if $rm_query;
957 :    
958 :     # Fix the blastp output to look like tblastn output
959 :    
960 :     foreach my $qdata ( @$out )
961 :     {
962 :     my %sdata; # There are now multiple sequences per subject DNA
963 :     foreach my $sdata ( @{ $qdata->[3] } )
964 :     {
965 :     my ( $sid, $sdef, undef, $hsps ) = @$sdata;
966 :     my $fr;
967 :     ( $sid, $fr ) = $sid =~ m/^(.*)\.([-+]\d)$/;
968 :     my ( $b, $e, $slen ) = $sdef =~ m/(\d+)-(\d+)\/(\d+)$/;
969 :     $sdef =~ s/ ?\S+$//;
970 :     my $sd2 = $sdata{ $sid } ||= [ $sid, $sdef, $slen, [] ];
971 :     foreach ( @$hsps ) { adjust_hsp( $_, $fr, $b ); push @{$sd2->[3]}, $_ }
972 :     }
973 :    
974 :     # Order the hsps for each subject
975 :    
976 :     foreach my $sid ( keys %sdata )
977 :     {
978 :     my $hsps = $sdata{ $sid }->[3];
979 :     @$hsps = sort { $b->[0] <=> $a->[0] } @$hsps;
980 :     }
981 :    
982 :     # Order the subjects for the query
983 :    
984 :     @{$qdata->[3]} = map { $_->[0] } # remove score
985 :     sort { $b->[1] <=> $a->[1] } # sort by score
986 :     map { [ $_, $_->[3]->[0]->[0] ] } # tag with top score
987 :     map { $sdata{ $_ } } # subject data
988 :     keys %sdata; # subject ids
989 :     }
990 :    
991 :     return $out;
992 :     }
993 :    
994 :    
995 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
996 :     # When search is versus six frame translation, there is a need to adjust
997 :     # the frame and location information in the hsp back to the DNA coordinates.
998 :     #
999 :     # adjust_hsp( $hsp, $frame, $begin )
1000 :     #
1001 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
1002 :     # scr Eval nseg Eval naln nid npos ngap fr q1 q2 qseq s1 s2 sseq
1003 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1004 :     sub adjust_hsp
1005 :     {
1006 :     my ( $hsp, $fr, $b ) = @_;
1007 :     $hsp->[8] = $fr;
1008 :     if ( $fr > 0 )
1009 :     {
1010 :     $hsp->[12] = $b + 3 * ( $hsp->[12] - 1 );
1011 :     $hsp->[13] = $b + 3 * ( $hsp->[13] - 1 ) + 2;
1012 :     }
1013 :     else
1014 :     {
1015 :     $hsp->[12] = $b - 3 * ( $hsp->[12] - 1 );
1016 :     $hsp->[13] = $b - 3 * ( $hsp->[13] - 1 ) - 2;
1017 :     }
1018 :     }
1019 :    
1020 :    
1021 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1022 :     # Do a six frame translation for use by blastpgpn. This modifications of
1023 :     # the identifiers and definitions are essential to the interpretation of
1024 :     # the blast results. The program 'translate_fasta_6' produces the same
1025 :     # output format, and is much faster.
1026 :     #
1027 :     # @translations = six_translations( $nucleotide_entry )
1028 :     #
1029 :     # The ids are modified by adding ".frame" (+1, +2, +3, -1, -2, -3).
1030 :     # The definition is mofidified by adding " begin-end/of_length".
1031 :     # NCBI frame numbers reverse strand translation frames from the end of the
1032 :     # sequence (i.e., the beginning of the complement of the strand).
1033 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1034 :    
1035 :     sub six_translations
1036 :     {
1037 :     my ( $id, $def, $seq ) = map { defined($_) ? $_ : '' } @{ $_[0] };
1038 :     my $l = length( $seq );
1039 :    
1040 :     return () if $l < 15;
1041 :    
1042 :     # fr beg end
1043 :     my @intervals = ( [ '+1', 1, $l - ( $l % 3 ) ],
1044 :     [ '+2', 2, $l - ( ($l-1) % 3 ) ],
1045 :     [ '+3', 3, $l - ( ($l-2) % 3 ) ],
1046 :     [ '-1', $l, 1 + ( $l % 3 ) ],
1047 :     [ '-2', $l-1, 1 + ( ($l-1) % 3 ) ],
1048 :     [ '-3', $l-2, 1 + ( ($l-2) % 3 ) ]
1049 :     );
1050 :     my ( $fr, $b, $e );
1051 :    
1052 :     map { ( $fr, $b, $e ) = @$_;
1053 :     [ "$id.$fr", "$def $b-$e/$l", gjoseqlib::translate_seq( gjoseqlib::DNA_subseq( \$seq, $b, $e ) ) ]
1054 :     } @intervals;
1055 :     }
1056 :    
1057 :    
1058 :     #-------------------------------------------------------------------------------
1059 : golsen 1.3 # Get an input file handle, and boolean on whether to close or not:
1060 :     #
1061 :     # ( \*FH, $close ) = input_file_handle( $filename );
1062 :     # ( \*FH, $close ) = input_file_handle( \*FH );
1063 :     # ( \*FH, $close ) = input_file_handle( ); # D = STDIN
1064 :     #
1065 :     #-------------------------------------------------------------------------------
1066 :    
1067 :     sub input_file_handle
1068 :     {
1069 :     my ( $file, $umask ) = @_;
1070 :    
1071 :     my ( $fh, $close );
1072 :     if ( defined $file )
1073 :     {
1074 :     if ( ref $file eq 'GLOB' )
1075 :     {
1076 :     $fh = $file;
1077 :     $close = 0;
1078 :     }
1079 :     elsif ( -f $file )
1080 :     {
1081 :     open( $fh, "<$file") || die "input_file_handle could not open '$file'.\n";
1082 :     $close = 1;
1083 :     }
1084 :     else
1085 :     {
1086 :     die "input_file_handle could not find file '$file'.\n";
1087 :     }
1088 :     }
1089 :     else
1090 :     {
1091 :     $fh = \*STDIN;
1092 :     $close = 0;
1093 :     }
1094 :    
1095 :     return ( $fh, $close );
1096 :     }
1097 :    
1098 :    
1099 :     #-------------------------------------------------------------------------------
1100 : overbeek 1.1 # Get an output file handle, and boolean on whether to close or not:
1101 :     #
1102 :     # ( \*FH, $close ) = output_file_handle( $filename );
1103 :     # ( \*FH, $close ) = output_file_handle( \*FH );
1104 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
1105 :     #
1106 :     #-------------------------------------------------------------------------------
1107 :    
1108 :     sub output_file_handle
1109 :     {
1110 :     my ( $file, $umask ) = @_;
1111 :    
1112 :     my ( $fh, $close );
1113 :     if ( defined $file )
1114 :     {
1115 :     if ( ref $file eq 'GLOB' )
1116 :     {
1117 :     $fh = $file;
1118 :     $close = 0;
1119 :     }
1120 :     else
1121 :     {
1122 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
1123 :     chmod 0664, $file; # Seems to work on open file!
1124 :     $close = 1;
1125 :     }
1126 :     }
1127 :     else
1128 :     {
1129 :     $fh = \*STDOUT;
1130 :     $close = 0;
1131 :     }
1132 :    
1133 :     return ( $fh, $close );
1134 :     }
1135 :    
1136 :    
1137 :     #-------------------------------------------------------------------------------
1138 :     #
1139 :     # print_blast_as_records( $file, \@queries )
1140 :     # print_blast_as_records( \*FH, \@queries )
1141 :     # print_blast_as_records( \@queries ) # STDOUT
1142 :     #
1143 :     # \@queries = read_blast_from_records( $file )
1144 :     # \@queries = read_blast_from_records( \*FH )
1145 :     # \@queries = read_blast_from_records( ) # STDIN
1146 :     #
1147 :     # Three output record types:
1148 :     #
1149 :     # Query \t qid \t qdef \t dlen
1150 :     # > \t sid \t sdef \t slen
1151 :     # HSP \t ...
1152 :     #
1153 :     #-------------------------------------------------------------------------------
1154 :     sub print_blast_as_records
1155 :     {
1156 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
1157 :     my ( $fh, $close ) = output_file_handle( $file );
1158 :    
1159 :     my $queries = shift;
1160 :     $queries && ref $queries eq 'ARRAY'
1161 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
1162 :     and return;
1163 :    
1164 :     foreach my $qdata ( @$queries )
1165 :     {
1166 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
1167 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
1168 :     and return;
1169 :    
1170 :     my $subjcts = pop @$qdata;
1171 :     $subjcts && ref $subjcts eq 'ARRAY'
1172 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
1173 :     and return;
1174 :    
1175 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
1176 :     foreach my $sdata ( @$subjcts )
1177 :     {
1178 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
1179 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
1180 :     and return;
1181 :    
1182 :     my $hsps = pop @$sdata;
1183 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
1184 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
1185 :     and return;
1186 :    
1187 :     print $fh join( "\t", '>', @$sdata ), "\n";
1188 :    
1189 :     foreach my $hsp ( @$hsps )
1190 :     {
1191 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
1192 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
1193 :     and return;
1194 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
1195 :     }
1196 :     }
1197 :     }
1198 :    
1199 :     close $fh if $close;
1200 :     }
1201 :    
1202 :    
1203 :     sub read_blast_from_records
1204 :     {
1205 :     my ( $file ) = @_;
1206 : golsen 1.3
1207 :     my ( $fh, $close ) = input_file_handle( $file );
1208 :     $fh or print STDERR "read_blast_from_records could not open input file.\n"
1209 : overbeek 1.1 and return wantarray ? () : [];
1210 : golsen 1.3
1211 : overbeek 1.1 my @queries = ();
1212 :     my @qdata;
1213 :     my @subjects = ();
1214 :     my @sdata;
1215 :     my @hsps = ();
1216 :    
1217 :     local $_;
1218 :     while ( defined( $_ = <$fh> ) )
1219 :     {
1220 :     chomp;
1221 :     my ( $type, @datum ) = split /\t/;
1222 :     if ( $type eq 'Query' )
1223 :     {
1224 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1225 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1226 :     @qdata = @datum;
1227 :     @sdata = ();
1228 :     @subjects = ();
1229 :     @hsps = ();
1230 :     }
1231 :     elsif ( $type eq '>' )
1232 :     {
1233 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1234 :     @sdata = @datum;
1235 :     @hsps = ();
1236 :     }
1237 :     elsif ( $type eq 'HSP' )
1238 :     {
1239 :     push @hsps, \@datum;
1240 :     }
1241 :     }
1242 :    
1243 :     close $fh if $close;
1244 :    
1245 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1246 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1247 :    
1248 :     wantarray ? @queries : \@queries;
1249 :     }
1250 :    
1251 :    
1252 : golsen 1.3 sub verify_db
1253 :     {
1254 :     my ( $db, $type ) = @_;
1255 : fangfang 1.18
1256 : golsen 1.3 my @args;
1257 : golsen 1.8 if ( $type =~ m/^p/i )
1258 : golsen 1.3 {
1259 : fangfang 1.18 @args = ( "-p", "T", "-i", $db ) unless ((-s "$db.psq") && (-M "$db.psq" <= -M $db)) || ((-s "$db.00.psq") && (-M "$db.00.psq" <= -M $db));
1260 : golsen 1.3 }
1261 :     else
1262 :     {
1263 : fangfang 1.18 @args = ( "-p", "F", "-i", $db ) unless ((-s "$db.nsq") && (-M "$db.nsq" <= -M $db)) || ((-s "$db.00.nsq") && (-M "$db.00.nsq" <= -M $db));
1264 : golsen 1.3 }
1265 : fangfang 1.18
1266 : golsen 1.8 @args or return ( -s $db ? 1 : 0 );
1267 :    
1268 : fangfang 1.18
1269 : golsen 1.3 #
1270 : golsen 1.8 # Find formatdb appropriate for the excecution environemnt.
1271 : golsen 1.3 #
1272 :    
1273 :     my $prog = SeedAware::executable_for( 'formatdb' );
1274 :     if ( ! $prog )
1275 :     {
1276 :     warn "gjoalignandtree::verify_db: formatdb program not found.\n";
1277 :     return 0;
1278 :     }
1279 :    
1280 :     my $rc = system( $prog, @args );
1281 :    
1282 :     if ( $rc != 0 )
1283 :     {
1284 :     my $cmd = join( ' ', $prog, @args );
1285 :     warn "gjoalignandtree::verify_db: formatdb failed with rc = $rc: $cmd\n";
1286 :     return 0;
1287 :     }
1288 :    
1289 :     return 1;
1290 :     }
1291 :    
1292 :    
1293 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3