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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (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 :    
219 :     my ( $rep ) = map { $->[0] } # sequence entry
220 :     sort { $b->[1] <=> $a->[1] } # max nongaps
221 :     map { [ $_, tr/-//c ] } # count nongaps
222 :     @cand;
223 :    
224 :     $rep = [ @$rep ]; # Make a copy
225 :     $rep->[2] =~ s/-+//g; # Compress sequence
226 :    
227 :     return $rep;
228 :     }
229 :    
230 :    
231 :     #-------------------------------------------------------------------------------
232 :     #
233 :     # $structured_blast = blastpgp( $db, $profile, $options )
234 :     #
235 :     # Required:
236 :     # db_file or db_seq
237 :     # profile_file or profile_seq
238 :     #
239 :     # Optional:
240 :     # query => q_file
241 :     # query => q_seq
242 :     #
243 :     #-------------------------------------------------------------------------------
244 :    
245 :     sub blastpgp
246 :     {
247 :     my ( $db, $profile, $opts ) = @_;
248 :     $opts ||= {};
249 :    
250 :     my ( $dbfile, $rm_db );
251 :     if ( defined $db && ref $db )
252 :     {
253 :     ref $db eq 'ARRAY'
254 :     && @$db
255 :     || print STDERR "blastpgp requires one or more database sequences.\n"
256 :     && return undef;
257 :     $dbfile = "blastpgp_db_$$";
258 :     gjoseqlib::print_alignment_as_fasta( $dbfile, $db );
259 :     $rm_db = 1;
260 :     }
261 :     elsif ( defined $db && -f $db )
262 :     {
263 :     $dbfile = $db;
264 :     }
265 :     else
266 :     {
267 :     die "blastpgp requires database.";
268 :     }
269 :     SEED_utils::verify_db( $db, 'T' ); # protein
270 :    
271 :     my ( $proffile, $rm_profile );
272 :     if ( defined $profile && ref $profile )
273 :     {
274 :     ref $profile eq 'ARRAY'
275 :     && @$profile
276 :     || print STDERR "blastpgp requires one or more profile sequences.\n"
277 :     && return undef;
278 :     $proffile = "blastpgp_profile_$$";
279 :     &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
280 :     $rm_profile = 1;
281 :     }
282 :     elsif ( defined $profile && -f $profile )
283 :     {
284 :     $proffile = $profile;
285 :     }
286 :     else
287 :     {
288 :     die "blastpgp requires profile.";
289 :     }
290 :    
291 :     my ( $qfile, $rm_query );
292 :     my $query = $opts->{ query };
293 :     if ( defined $query && ref $query )
294 :     {
295 :     ref $query eq 'ARRAY'
296 :     && @$query
297 :     || print STDERR "blastpgp invalid query sequence.\n"
298 :     && return undef;
299 :     $qfile = "blastpgp_query_$$";
300 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
301 :     $rm_query = 1;
302 :     }
303 :     elsif ( defined $query && -f $query )
304 :     {
305 :     $qfile = $query;
306 :     }
307 :     elsif ( $profile ) # Build it from profile
308 :     {
309 :     $query = &representative_for_profile( $profile );
310 :     $qfile = "blastpgp_query_$$";
311 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
312 :     $rm_query = 1;
313 :     }
314 :     else
315 :     {;
316 :     die "blastpgp requires database.";
317 :     }
318 :    
319 :     my $nkeep = $opts->{ nresult } || 5000;
320 :     my $cmd = "blastpgp -B '$proffile' -j 1 -d '$dbfile -b $nkeep -v $nkeep -i '$qfile'";
321 :    
322 :     open BLAST, "$cmd |" or print STDERR "Failed to open: '$cmd'.\n"
323 :     and return undef;
324 :     my $out = &gjoparseblast::structured_blast_output( \*BLAST, 1 ); # selfmatches okay
325 :     close BLAST;
326 :    
327 :     if ( $rm_db )
328 :     {
329 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
330 :     unlink @files if @files;
331 :     }
332 :     unlink $proffile if $rm_profile;
333 :     unlink $qfile if $rm_query;
334 :    
335 :     return $out;
336 :     }
337 :    
338 :    
339 :     #-------------------------------------------------------------------------------
340 :     # Get an output file handle, and boolean on whether to close or not:
341 :     #
342 :     # ( \*FH, $close ) = output_file_handle( $filename );
343 :     # ( \*FH, $close ) = output_file_handle( \*FH );
344 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
345 :     #
346 :     #-------------------------------------------------------------------------------
347 :    
348 :     sub output_file_handle
349 :     {
350 :     my ( $file, $umask ) = @_;
351 :    
352 :     my ( $fh, $close );
353 :     if ( defined $file )
354 :     {
355 :     if ( ref $file eq 'GLOB' )
356 :     {
357 :     $fh = $file;
358 :     $close = 0;
359 :     }
360 :     else
361 :     {
362 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
363 :     chmod 0664, $file; # Seems to work on open file!
364 :     $close = 1;
365 :     }
366 :     }
367 :     else
368 :     {
369 :     $fh = \*STDOUT;
370 :     $close = 0;
371 :     }
372 :    
373 :     return ( $fh, $close );
374 :     }
375 :    
376 :    
377 :     #-------------------------------------------------------------------------------
378 :     #
379 :     # print_blast_as_records( $file, \@queries )
380 :     # print_blast_as_records( \*FH, \@queries )
381 :     # print_blast_as_records( \@queries ) # STDOUT
382 :     #
383 :     # \@queries = read_blast_from_records( $file )
384 :     # \@queries = read_blast_from_records( \*FH )
385 :     # \@queries = read_blast_from_records( ) # STDIN
386 :     #
387 :     # Three output record types:
388 :     #
389 :     # Query \t qid \t qdef \t dlen
390 :     # > \t sid \t sdef \t slen
391 :     # HSP \t ...
392 :     #
393 :     #-------------------------------------------------------------------------------
394 :     sub print_blast_as_records
395 :     {
396 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
397 :     my ( $fh, $close ) = output_file_handle( $file );
398 :    
399 :     my $queries = shift;
400 :     $queries && ref $queries eq 'ARRAY'
401 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
402 :     and return;
403 :    
404 :     foreach my $qdata ( @$queries )
405 :     {
406 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
407 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
408 :     and return;
409 :    
410 :     my $subjcts = pop @$qdata;
411 :     $subjcts && ref $subjcts eq 'ARRAY'
412 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
413 :     and return;
414 :    
415 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
416 :     foreach my $sdata ( @$subjcts )
417 :     {
418 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
419 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
420 :     and return;
421 :    
422 :     my $hsps = pop @$sdata;
423 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
424 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
425 :     and return;
426 :    
427 :     print $fh join( "\t", '>', @$sdata ), "\n";
428 :    
429 :     foreach my $hsp ( @$hsps )
430 :     {
431 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
432 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
433 :     and return;
434 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
435 :     }
436 :     }
437 :     }
438 :    
439 :     close $fh if $close;
440 :     }
441 :    
442 :    
443 :     sub read_blast_from_records
444 :     {
445 :     my ( $file ) = @_;
446 :     my ( $fh, $close );
447 :     if ( $file && ref $file eq 'GLOB' )
448 :     {
449 :     $fh = $file;
450 :     }
451 :     elsif ( defined $file )
452 :     {
453 :     -f $file && open( $fh, "<$file" )
454 :     or print STDERR "read_blast_from_records could not open '$file'.\n"
455 :     and return wantarray ? () : [];
456 :     $close = 1;
457 :     }
458 :     else
459 :     {
460 :     $fh = \@STDIN;
461 :     }
462 :     my @queries = ();
463 :     my @qdata;
464 :     my @subjects = ();
465 :     my @sdata;
466 :     my @hsps = ();
467 :    
468 :     local $_;
469 :     while ( defined( $_ = <$fh> ) )
470 :     {
471 :     chomp;
472 :     my ( $type, @datum ) = split /\t/;
473 :     if ( $type eq 'Query' )
474 :     {
475 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
476 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
477 :     @qdata = @datum;
478 :     @sdata = ();
479 :     @subjects = ();
480 :     @hsps = ();
481 :     }
482 :     elsif ( $type eq '>' )
483 :     {
484 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
485 :     @sdata = @datum;
486 :     @hsps = ();
487 :     }
488 :     elsif ( $type eq 'HSP' )
489 :     {
490 :     push @hsps, \@datum;
491 :     }
492 :     }
493 :    
494 :     close $fh if $close;
495 :    
496 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
497 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
498 :    
499 :     wantarray ? @queries : \@queries;
500 :     }
501 :    
502 :    
503 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3