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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package gjoalignandtree;
2 :    
3 :     #
4 :    
5 :     use strict;
6 :     use gjoseqlib;
7 :    
8 :     #-------------------------------------------------------------------------------
9 :     # Order a set of sequences:
10 :     #
11 :     # @seqs = order_sequences_by_length( \@seqs, \%opts )
12 :     # \@seqs = order_sequences_by_length( \@seqs, \%opts )
13 :     #
14 :     # Options:
15 :     #
16 :     # order => key # increasing (D), decreasing, median_outward,
17 :     # # closest_to_median, closest_to_length
18 :     # length => prefered_length
19 :     #
20 :     #-------------------------------------------------------------------------------
21 :    
22 :     sub order_sequences_by_length
23 :     {
24 :     my ( $seqs, $opts ) = @_;
25 :     $seqs && ref $seqs eq 'ARRAY' && @$seqs
26 :     or print STDERR "order_sequences_by_length called with invalid sequence list.\n"
27 :     and return wantarray ? () : [];
28 :    
29 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
30 :     my $order = $opts->{ order } || 'increasing';
31 :    
32 :     my @seqs = ();
33 :     if ( $order =~ /^dec/i )
34 :     {
35 :     $opts->{ order } = 'decreaseing';
36 :     @seqs = sort { length( $b->[2] ) <=> length( $a->[2] ) } @$seqs;
37 :     }
38 :     elsif ( $order =~ /^med/i )
39 :     {
40 :     $opts->{ order } = 'median_outward';
41 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
42 :     while ( $_ = shift @seqs0 )
43 :     {
44 :     push @seqs, $_;
45 :     push @seqs, pop @seqs0 if @seqs0;
46 :     }
47 :     @seqs = reverse @seqs;
48 :     }
49 :     elsif ( $order =~ /^close/i )
50 :     {
51 :     my $pref_len;
52 :     if ( defined $opts->{ length } ) # caller supplies preferred length?
53 :     {
54 :     $opts->{ order } = 'closest_to_length';
55 :     $pref_len = $opts->{ length } + 0;
56 :     }
57 :     else # preferred length is median
58 :     {
59 :     $opts->{ order } = 'closest_to_median';
60 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
61 :     my $i = int( @seqs0 / 2 );
62 :     $pref_len = ( @seqs0 % 2 ) ? length( $seqs0[$i]->[2] )
63 :     : ( length( $seqs0[$i-1]->[2] ) + length( $seqs0[$i]->[2] ) ) / 2;
64 :     }
65 :    
66 :     @seqs = map { $_->[0] }
67 :     sort { $a->[1] <=> $b->[1] } # closest to preferred first
68 :     map { [ $_, abs( length( $_->[2] ) - $pref_len ) ] } # dist from pref?
69 :     @$seqs;
70 :     }
71 :     else
72 :     {
73 :     $opts->{ order } = 'increasing';
74 :     @seqs = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
75 :     }
76 :    
77 :     wantarray ? @seqs : \@seqs;
78 :     }
79 :    
80 :    
81 :     #-------------------------------------------------------------------------------
82 :     # Write representative seqeunce groups to a file
83 :     #
84 :     # write_representative_groups( \%rep, $file )
85 :     #
86 :     #-------------------------------------------------------------------------------
87 :    
88 :     sub write_representative_groups
89 :     {
90 :     my ( $repH, $file ) = @_;
91 :     $repH && (ref $repH eq 'HASH') && keys %$repH
92 :     or print STDERR "write_representative_groups called with invalid rep hash.\n"
93 :     and return;
94 :    
95 :     my ( $fh, $close ) = &output_file_handle( $file );
96 :    
97 :     # Order keys from largest to smallest set, then alphabetical
98 :    
99 :     my @keys = sort { @{ $repH->{$b} } <=> @{ $repH->{$a} } || lc $a cmp lc $b } keys %$repH;
100 :    
101 :     foreach ( @keys ) { print $fh join( "\t", ( $_, @{ $repH->{ $_ } } ) ), "\n" }
102 :    
103 :     close $fh if $close;
104 :     }
105 :    
106 :    
107 :     #-------------------------------------------------------------------------------
108 :     #
109 :     # @align = trim_align_to_median_ends( \@align, \%opts )
110 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
111 :     #
112 :     # Options:
113 :     #
114 :     # begin => bool # Trim start (specifically)
115 :     # end => bool # Trim end (specifically)
116 :     #
117 :     #-------------------------------------------------------------------------------
118 :     sub trim_align_to_median_ends
119 :     {
120 :     my ( $align, $opts ) = @_;
121 :     $align && ref $align eq 'ARRAY' && @$align
122 :     or print STDERR "trim_align_to_median_ends called with invalid sequence list.\n"
123 :     and return wantarray ? () : [];
124 :    
125 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
126 :     my $tr_beg = $opts->{ begin } || $opts->{ beg } || $opts->{ start } ? 1 : 0;
127 :     my $tr_end = $opts->{ end } ? 1 : 0;
128 :     $tr_beg = $tr_end = 1 if ! ( $tr_beg || $tr_end );
129 :    
130 :     my( @pos1, @pos2, $i );
131 :     foreach my $seq ( @$align )
132 :     {
133 :     my( $b, $e ) = $seq->[2] =~ /^(-*).*[^-](-*)$/;
134 :     push @pos1, length( $b || '' );
135 :     push @pos2, length( $e || '' );
136 :     }
137 :    
138 :     @pos1 = sort { $a <=> $b } @pos1;
139 :     @pos2 = sort { $a <=> $b } @pos2;
140 :    
141 :     my $pos1 = $tr_beg ? $pos1[ int( @pos1/2 ) ] : 0;
142 :     my $pos2 = $tr_end ? $pos2[ int( @pos2/2 ) ] : 0;
143 :    
144 :     my $ori_len = length( $align->[0]->[2] );
145 :     my $new_len = $ori_len - ( $pos1 + $pos2 );
146 :     my @align2 = map { [ @$_[0,1], substr( $_->[2], $pos1, $new_len ) ] }
147 :     @$align;
148 :    
149 :     wantarray ? @align2 : \@align2;
150 :     }
151 :    
152 :    
153 :     #-------------------------------------------------------------------------------
154 :     #
155 :     # Options:
156 :     #
157 :     # file => $filename # supply a file name to open and write
158 :     # file => \*FH # supply an open file handle (D = STDOUT)
159 :     # line => $linelen # residues per line (D = 60)
160 :     # lower => $bool # all lower case sequence
161 :     # upper => $bool # all upper case sequence
162 :     #
163 :     #-------------------------------------------------------------------------------
164 :    
165 :     sub print_alignment_as_pseudoclustal
166 :     {
167 :     my ( $align, $opts ) = @_;
168 :     $align && ref $align eq 'ARRAY' && @$align
169 :     or print STDERR "print_alignment_as_pseudoclustal called with invalid sequence list.\n"
170 :     and return wantarray ? () : [];
171 :    
172 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
173 :     my $line_len = $opts->{ line } || 60;
174 :     my $case = $opts->{ upper } ? 1 : $opts->{ lower } ? -1 : 0;
175 :    
176 :     my ( $fh, $close ) = &output_file_handle( $opts->{ file } );
177 :    
178 :     my $namelen = 0;
179 :     foreach ( @$align ) { $namelen = length $_->[0] if $namelen < length $_->[0] }
180 :     my $fmt = "%-${namelen}s %s\n";
181 :    
182 :     my $id;
183 :     my @lines = map { $id = $_->[0]; [ map { sprintf $fmt, $id, $_ }
184 :     map { $case < 0 ? lc $_ : $case > 0 ? up $_ : $_ } # map sequence only
185 :     $_->[2] =~ m/.{1,$line_len}/g
186 :     ] }
187 :    
188 :     @$align;
189 :    
190 :     my $ngroup = @{ $lines[0] };
191 :     for ( my $i = 0; $i < $ngroup; $i++ )
192 :     {
193 :     foreach ( @lines ) { print $fh $_->[$i] if $_->[$i] }
194 :     print $fh "\n";
195 :     }
196 :    
197 :     close $fh if $close;
198 :     }
199 :    
200 :    
201 :     #-------------------------------------------------------------------------------
202 :     # The profile 'query' sequence:
203 :     # 1. No terminal gaps, if possible
204 :     # 1.1 Otherwise, no gaps at start
205 :     # 2. Longest sequence passing above
206 :     #
207 :     # $prof_rep = representative_for_profile( $align )
208 :     #-------------------------------------------------------------------------------
209 :     sub representative_for_profile
210 :     {
211 :     my ( $align ) = @_;
212 :     $align && ref $align eq 'ARRAY' && @$align
213 :     or die "representative_for_profile called with invalid sequence list.\n";
214 :    
215 :     my @cand;
216 :     @cand = grep { $_->[2] =~ /^[^-].*[^-]$/ } @$align; # No end gaps
217 :     @cand = grep { $_->[2] =~ /^[^-]/ } @$align if ! @cand; # No start gaps
218 : fangfang 1.2 @cand = @$align if ! @cand;
219 : overbeek 1.1
220 :     my ( $rep ) = map { $->[0] } # sequence entry
221 :     sort { $b->[1] <=> $a->[1] } # max nongaps
222 :     map { [ $_, tr/-//c ] } # count nongaps
223 :     @cand;
224 :    
225 :     $rep = [ @$rep ]; # Make a copy
226 :     $rep->[2] =~ s/-+//g; # Compress sequence
227 :    
228 :     return $rep;
229 :     }
230 :    
231 :    
232 :     #-------------------------------------------------------------------------------
233 :     #
234 :     # $structured_blast = blastpgp( $db, $profile, $options )
235 :     #
236 :     # Required:
237 :     # db_file or db_seq
238 :     # profile_file or profile_seq
239 :     #
240 :     # Optional:
241 :     # query => q_file
242 :     # query => q_seq
243 :     #
244 :     #-------------------------------------------------------------------------------
245 :    
246 :     sub blastpgp
247 :     {
248 :     my ( $db, $profile, $opts ) = @_;
249 :     $opts ||= {};
250 :    
251 :     my ( $dbfile, $rm_db );
252 :     if ( defined $db && ref $db )
253 :     {
254 :     ref $db eq 'ARRAY'
255 :     && @$db
256 :     || print STDERR "blastpgp requires one or more database sequences.\n"
257 :     && return undef;
258 :     $dbfile = "blastpgp_db_$$";
259 :     gjoseqlib::print_alignment_as_fasta( $dbfile, $db );
260 :     $rm_db = 1;
261 :     }
262 :     elsif ( defined $db && -f $db )
263 :     {
264 :     $dbfile = $db;
265 :     }
266 :     else
267 :     {
268 :     die "blastpgp requires database.";
269 :     }
270 :     SEED_utils::verify_db( $db, 'T' ); # protein
271 :    
272 :     my ( $proffile, $rm_profile );
273 :     if ( defined $profile && ref $profile )
274 :     {
275 :     ref $profile eq 'ARRAY'
276 :     && @$profile
277 :     || print STDERR "blastpgp requires one or more profile sequences.\n"
278 :     && return undef;
279 :     $proffile = "blastpgp_profile_$$";
280 :     &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
281 :     $rm_profile = 1;
282 :     }
283 :     elsif ( defined $profile && -f $profile )
284 :     {
285 :     $proffile = $profile;
286 :     }
287 :     else
288 :     {
289 :     die "blastpgp requires profile.";
290 :     }
291 :    
292 :     my ( $qfile, $rm_query );
293 :     my $query = $opts->{ query };
294 :     if ( defined $query && ref $query )
295 :     {
296 :     ref $query eq 'ARRAY'
297 :     && @$query
298 :     || print STDERR "blastpgp invalid query sequence.\n"
299 :     && return undef;
300 :     $qfile = "blastpgp_query_$$";
301 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
302 :     $rm_query = 1;
303 :     }
304 :     elsif ( defined $query && -f $query )
305 :     {
306 :     $qfile = $query;
307 :     }
308 :     elsif ( $profile ) # Build it from profile
309 :     {
310 :     $query = &representative_for_profile( $profile );
311 :     $qfile = "blastpgp_query_$$";
312 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
313 :     $rm_query = 1;
314 :     }
315 :     else
316 :     {;
317 :     die "blastpgp requires database.";
318 :     }
319 :    
320 :     my $nkeep = $opts->{ nresult } || 5000;
321 :     my $cmd = "blastpgp -B '$proffile' -j 1 -d '$dbfile -b $nkeep -v $nkeep -i '$qfile'";
322 :    
323 :     open BLAST, "$cmd |" or print STDERR "Failed to open: '$cmd'.\n"
324 :     and return undef;
325 :     my $out = &gjoparseblast::structured_blast_output( \*BLAST, 1 ); # selfmatches okay
326 :     close BLAST;
327 :    
328 :     if ( $rm_db )
329 :     {
330 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
331 :     unlink @files if @files;
332 :     }
333 :     unlink $proffile if $rm_profile;
334 :     unlink $qfile if $rm_query;
335 :    
336 :     return $out;
337 :     }
338 :    
339 :    
340 :     #-------------------------------------------------------------------------------
341 :     # Get an output file handle, and boolean on whether to close or not:
342 :     #
343 :     # ( \*FH, $close ) = output_file_handle( $filename );
344 :     # ( \*FH, $close ) = output_file_handle( \*FH );
345 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
346 :     #
347 :     #-------------------------------------------------------------------------------
348 :    
349 :     sub output_file_handle
350 :     {
351 :     my ( $file, $umask ) = @_;
352 :    
353 :     my ( $fh, $close );
354 :     if ( defined $file )
355 :     {
356 :     if ( ref $file eq 'GLOB' )
357 :     {
358 :     $fh = $file;
359 :     $close = 0;
360 :     }
361 :     else
362 :     {
363 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
364 :     chmod 0664, $file; # Seems to work on open file!
365 :     $close = 1;
366 :     }
367 :     }
368 :     else
369 :     {
370 :     $fh = \*STDOUT;
371 :     $close = 0;
372 :     }
373 :    
374 :     return ( $fh, $close );
375 :     }
376 :    
377 :    
378 :     #-------------------------------------------------------------------------------
379 :     #
380 :     # print_blast_as_records( $file, \@queries )
381 :     # print_blast_as_records( \*FH, \@queries )
382 :     # print_blast_as_records( \@queries ) # STDOUT
383 :     #
384 :     # \@queries = read_blast_from_records( $file )
385 :     # \@queries = read_blast_from_records( \*FH )
386 :     # \@queries = read_blast_from_records( ) # STDIN
387 :     #
388 :     # Three output record types:
389 :     #
390 :     # Query \t qid \t qdef \t dlen
391 :     # > \t sid \t sdef \t slen
392 :     # HSP \t ...
393 :     #
394 :     #-------------------------------------------------------------------------------
395 :     sub print_blast_as_records
396 :     {
397 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
398 :     my ( $fh, $close ) = output_file_handle( $file );
399 :    
400 :     my $queries = shift;
401 :     $queries && ref $queries eq 'ARRAY'
402 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
403 :     and return;
404 :    
405 :     foreach my $qdata ( @$queries )
406 :     {
407 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
408 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
409 :     and return;
410 :    
411 :     my $subjcts = pop @$qdata;
412 :     $subjcts && ref $subjcts eq 'ARRAY'
413 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
414 :     and return;
415 :    
416 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
417 :     foreach my $sdata ( @$subjcts )
418 :     {
419 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
420 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
421 :     and return;
422 :    
423 :     my $hsps = pop @$sdata;
424 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
425 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
426 :     and return;
427 :    
428 :     print $fh join( "\t", '>', @$sdata ), "\n";
429 :    
430 :     foreach my $hsp ( @$hsps )
431 :     {
432 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
433 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
434 :     and return;
435 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
436 :     }
437 :     }
438 :     }
439 :    
440 :     close $fh if $close;
441 :     }
442 :    
443 :    
444 :     sub read_blast_from_records
445 :     {
446 :     my ( $file ) = @_;
447 :     my ( $fh, $close );
448 :     if ( $file && ref $file eq 'GLOB' )
449 :     {
450 :     $fh = $file;
451 :     }
452 :     elsif ( defined $file )
453 :     {
454 :     -f $file && open( $fh, "<$file" )
455 :     or print STDERR "read_blast_from_records could not open '$file'.\n"
456 :     and return wantarray ? () : [];
457 :     $close = 1;
458 :     }
459 :     else
460 :     {
461 :     $fh = \@STDIN;
462 :     }
463 :     my @queries = ();
464 :     my @qdata;
465 :     my @subjects = ();
466 :     my @sdata;
467 :     my @hsps = ();
468 :    
469 :     local $_;
470 :     while ( defined( $_ = <$fh> ) )
471 :     {
472 :     chomp;
473 :     my ( $type, @datum ) = split /\t/;
474 :     if ( $type eq 'Query' )
475 :     {
476 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
477 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
478 :     @qdata = @datum;
479 :     @sdata = ();
480 :     @subjects = ();
481 :     @hsps = ();
482 :     }
483 :     elsif ( $type eq '>' )
484 :     {
485 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
486 :     @sdata = @datum;
487 :     @hsps = ();
488 :     }
489 :     elsif ( $type eq 'HSP' )
490 :     {
491 :     push @hsps, \@datum;
492 :     }
493 :     }
494 :    
495 :     close $fh if $close;
496 :    
497 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
498 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
499 :    
500 :     wantarray ? @queries : \@queries;
501 :     }
502 :    
503 :    
504 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3