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

Annotation of /FigKernelPackages/ProtFamsLite.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : arodri7 1.1 # -*- perl -*-
2 :     ########################################################################
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     ########################################################################
18 :    
19 :     package ProtFamsLite;
20 :    
21 :     use strict;
22 :     use DB_File;
23 :    
24 :     #use ProtFamLite;
25 :     use Tracer;
26 :    
27 :     use Data::Dumper;
28 :     use Carp;
29 :     use Digest::MD5;
30 : arodri7 1.2 use FIG;
31 : arodri7 1.1
32 :     # This is the constructor. Presumably, $class is 'ProtFamsLite'.
33 :     #
34 :    
35 :     sub new {
36 :     my($class,$fam_data) = @_;
37 :    
38 :     my $protfams = {};
39 :    
40 :     defined($fam_data) || return undef;
41 :     $protfams->{dir} = $fam_data;
42 : arodri7 1.2 $protfams->{root} = $fam_data;
43 :     $protfams->{fig} = new FIG;
44 : arodri7 1.1
45 :     bless $protfams,$class;
46 :     return $protfams;
47 :     }
48 :    
49 :     sub create_tie
50 :     {
51 :     my($self, $hash, $file, $type) = @_;
52 :    
53 :     my $mode = -w $file ? O_RDWR : O_RDONLY;
54 :    
55 :     return tie %$hash, 'DB_File', $file, $mode, 0666, $type;
56 :     }
57 :    
58 :     sub check_db_role {
59 :     my($self) = @_;
60 :     my $fam_data = $self->{dir};
61 :    
62 :     my $db_role = "$fam_data/role.db";
63 :     if (! -s $db_role) { return undef }
64 :    
65 :     my %role_hash;
66 :     my $role_hash_tie = $self->create_tie(\%role_hash, $db_role, $DB_HASH);
67 :     $role_hash_tie || die "tie $db_role failed: $!";
68 :     $self->{role_db} = \%role_hash;
69 :     }
70 :    
71 :     sub check_db_function {
72 :     my($self) = @_;
73 :     my $fam_data = $self->{dir};
74 :    
75 :     my $db_function = "$fam_data/function.db";
76 :     if (! -s $db_function) { return undef }
77 :     my %function_hash;
78 :     my $function_hash_tie = $self->create_tie(\%function_hash, $db_function, $DB_HASH);
79 :     $function_hash_tie || die "tie $db_function failed: $!";
80 :     $self->{function_db} = \%function_hash;
81 :     }
82 :    
83 : arodri7 1.2 sub check_db_family_function {
84 :     my($self) = @_;
85 :     my $fam_data = $self->{dir};
86 :    
87 :     my $db_function = "$fam_data/fam_function.db";
88 :     if (! -s $db_function) { return undef }
89 :     my %function_hash;
90 :     my $function_hash_tie = $self->create_tie(\%function_hash, $db_function, $DB_HASH);
91 :     $function_hash_tie || die "tie $db_function failed: $!";
92 :     $self->{family_function_db} = \%function_hash;
93 :     }
94 :    
95 : arodri7 1.1 sub check_db_prot_to_fams {
96 :     my($self) = @_;
97 :     my $fam_data = $self->{dir};
98 :    
99 :     my $db_prot_to_fams = "$fam_data/protein.db";
100 :     if (! -s $db_prot_to_fams) { return undef }
101 :     my %prot_to_fams_hash;
102 :     my $prot_to_fams_hash_tie = $self->create_tie(\%prot_to_fams_hash, $db_prot_to_fams, $DB_HASH);
103 :     $prot_to_fams_hash_tie || die "tie $db_prot_to_fams failed: $!";
104 :     $self->{prot_to_fams_db} = \%prot_to_fams_hash;
105 :     }
106 :    
107 :     sub check_db_genome_to_fams {
108 :     my($self) = @_;
109 :     my $fam_data = $self->{dir};
110 :    
111 :     my $db_genome_to_fams = "$fam_data/genome.db";
112 :     if (! -s $db_genome_to_fams) { return undef }
113 :     my %genome_to_fams_hash;
114 :     my $genome_to_fams_hash_tie = $self->create_tie(\%genome_to_fams_hash, $db_genome_to_fams, $DB_HASH);
115 :     $genome_to_fams_hash_tie || die "tie $db_genome_to_fams failed: $!";
116 :     $self->{genome_to_fams_db} = \%genome_to_fams_hash;
117 :     }
118 :    
119 :     sub check_db_relevant_prot_data {
120 :     my($self) = @_;
121 :     my $fam_data = $self->{dir};
122 :    
123 :     my $db_relevant_prot_data = "$fam_data/relevant.prot.data.db";
124 :     if (! -s $db_relevant_prot_data) { return undef }
125 :     my %relevant_prot_data_hash;
126 :     my $relevant_prot_data_hash_tie = $self->create_tie(\%relevant_prot_data_hash, $db_relevant_prot_data, $DB_HASH);
127 :     $relevant_prot_data_hash_tie || die "tie $db_relevant_prot_data failed: $!";
128 :     $self->{relevant_prot_data_db} = \%relevant_prot_data_hash;
129 :     }
130 :    
131 :     sub check_db_PDB_connections {
132 :     my($self) = @_;
133 :    
134 :     my $fam_data = $self->{dir};
135 :     my $db_PDB_connections = "$fam_data/PDB.connections.db";
136 :     if (! -s $db_PDB_connections) { return undef }
137 :     my %PDB_hash;
138 :     my $PDB_hash_tie = $self->create_tie(\%PDB_hash, $db_PDB_connections, $DB_HASH);
139 :     $PDB_hash_tie || die "tie $db_PDB_connections failed: $!";
140 :     $self->{PDB_connections_db} = \%PDB_hash;
141 :     }
142 :    
143 :     sub check_db_md5_to_fams {
144 :     my($self) = @_;
145 :    
146 :     my $fam_data = $self->{dir};
147 :     my $db_md5_fams = "$fam_data/md5.protfams.db";
148 :     if (! -s $db_md5_fams) { return undef }
149 :     my %md5_hash;
150 :     my $md5_hash_tie = $self->create_tie(\%md5_hash, $db_md5_fams, $DB_HASH);
151 :     $md5_hash_tie || die "tie $db_md5_fams failed: $!";
152 :     $self->{md5_fams_db} = \%md5_hash;
153 :     }
154 :    
155 : arodri7 1.2 sub check_db_md5_to_prots {
156 :     my($self) = @_;
157 :    
158 :     my $fam_data = $self->{dir};
159 :     my $db_md5_fams = "$fam_data/md5.prots.db";
160 :     if (! -s $db_md5_fams) { return undef }
161 :     my %md5_hash;
162 :     my $md5_hash_tie = $self->create_tie(\%md5_hash, $db_md5_fams, $DB_HASH);
163 :     $md5_hash_tie || die "tie $db_md5_fams failed: $!";
164 :     $self->{md5_prots_db} = \%md5_hash;
165 :     }
166 :    
167 : arodri7 1.1 sub check_db_family_id_map {
168 :     my($self) = @_;
169 :     my $fam_data = $self->{dir};
170 :    
171 :     my $db_map = "$fam_data/family_id_map.db";
172 :     if (! -s $db_map) { return undef }
173 :    
174 :     my %map_hash;
175 :     my $map_hash_tie = $self->create_tie(\%map_hash, $db_map, $DB_HASH);
176 :     $map_hash_tie || die "tie $db_map failed: $!";
177 :     $self->{map_db} = \%map_hash;
178 :     }
179 :    
180 :     sub PDB_connections {
181 :     my($self,$fam,$raw) = @_;
182 :    
183 :     $self->check_db_PDB_connections;
184 :     my $sims = $self->{PDB_connections_db}->{$fam};
185 :     my @sims = map { $_ =~ /pdb\|([0-9a-zA-Z]+)/; [$1,[split(/\t/,$_)]] } split(/\n/,$sims);
186 :     if (! $raw) { @sims = map { $_->[0] } grep { ($_->[1]->[11] > 0.5) && ((($_->[1]->[4] - $_->[1]->[3]) / $_->[1]->[5]) > 0.8) } @sims}
187 :     return \@sims;
188 :     }
189 :    
190 :     sub function_of {
191 :     my($self,$prot,$ignore_comments) = @_;
192 :    
193 :     $self->check_db_relevant_prot_data; # lazily tie DB
194 :    
195 :     my $prot_data;
196 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
197 :     {
198 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
199 :     if ($ignore_comments)
200 :     {
201 :     $func =~ s/\s*\#.*$//;
202 :     }
203 :     return $func;
204 :     }
205 :     return "";
206 :     }
207 :    
208 : arodri7 1.2 sub function_of_family {
209 :     my($self,$fam,$ignore_comments) = @_;
210 :    
211 :     $self->check_db_family_function; # lazily tie DB
212 :    
213 :     my $func;
214 :     if ($func = $self->{family_function_db}->{$fam})
215 :     {
216 :     return $func;
217 :     }
218 :     }
219 :    
220 : arodri7 1.1 sub get_translation {
221 :     my($self,$prot,$ignore_comments) = @_;
222 :    
223 :     $self->check_db_relevant_prot_data; # lazily tie DB
224 :    
225 :     my $prot_data;
226 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
227 :     {
228 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
229 :     return $seq;
230 :     }
231 :     return "";
232 :     }
233 :    
234 : arodri7 1.2 sub get_translation_bulk {
235 :     my($self,$prots,$ignore_comments) = @_;
236 :     my $seqs = {};
237 :    
238 :     $self->check_db_relevant_prot_data; # lazily tie DB
239 :    
240 :     my $prot_data;
241 :     foreach my $prot (@$prots)
242 :     {
243 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
244 :     {
245 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
246 :     $seqs->{$prot} = $seq;
247 :     }
248 :     }
249 :    
250 :     return $seqs;
251 :     }
252 :    
253 : arodri7 1.1
254 :     sub md5_function_of {
255 :     my ($self, $md5, $preferred_db) = @_;
256 :    
257 :     # find if md5 is in a family
258 :    
259 :     # if in_family get the function
260 :    
261 :     # else go to seed to get function
262 :     my (@tuples, $function);
263 :     if (@tuples = $self->{fig}->mapped_prot_ids(qq(gnl|md5|$md5))) {
264 :     my @tmp;
265 :     if (@tmp = grep { $_->[0] =~ m/^fig\|\d+\.\d+\.peg\.\d+$/o } @tuples) {
266 :     #...Pick the function of one of the FIDs
267 :     }
268 :     elsif (@tmp = grep { $_->[0] =~ m/^sp\|\S+$/o } @tuples) {
269 :     #...Pick the function of one of the SP's
270 :     }
271 :     elsif (@tmp = grep { $_->[0] =~ m/^pir\|\S+$/o } @tuples) {
272 :     #...DWIM
273 :     }
274 :     elsif (@tmp = grep { $_->[0] =~ m/^kegg\|\S+$/o } @tuples) {
275 :     }
276 :     elsif (@tmp = grep { $_->[0] =~ m/^gi\|\S+$/o } @tuples) {
277 :     }
278 :     else {
279 :     $function = qq(Hypothetical protein);
280 :     }
281 :     return $function;
282 :     }
283 :     }
284 :    
285 :     sub org_of {
286 :     my($self,$prot) = @_;
287 :    
288 :     $self->check_db_relevant_prot_data; # lazily tie DB
289 :    
290 :     my $prot_data;
291 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
292 :     {
293 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
294 :     return $org;
295 :     }
296 :     return "";
297 :     }
298 :    
299 :     sub path_to {
300 :     my ($self,$family_name) = @_;
301 :    
302 :     $self->check_db_family_id_map; # lazily tie DB
303 :    
304 :     my $internal_id;
305 : arodri7 1.2 if ($internal_id = $self->{map_db}->{$family_name})
306 : arodri7 1.1 {
307 :     #my($internal_id) = split(/\t/,$family_name);
308 :    
309 :     my $group = substr($internal_id, -3);
310 :     $group = (qq(0) x (3-length($group))) . $group;
311 :    
312 :     my $path_to_family = qq($self->{dir}/FAMS/$group/$internal_id);
313 :    
314 :     return $path_to_family;
315 :     }
316 :     return "";
317 :     }
318 :    
319 :     sub seq_of {
320 :     my($self,$prot) = @_;
321 :    
322 :     $self->check_db_relevant_prot_data; # lazily tie DB
323 :    
324 :     my $prot_data;
325 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
326 :     {
327 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
328 :     return $seq;
329 :     }
330 :     return "";
331 :     }
332 :    
333 :     sub aliases_of {
334 :     my($self,$prot) = @_;
335 :    
336 :     $self->check_db_relevant_prot_data; # lazily tie DB
337 :    
338 :     my $prot_data;
339 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
340 :     {
341 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
342 :     return wantarray() ? split(/,/,$aliases) : $aliases;
343 :     }
344 :     return "";
345 :     }
346 :    
347 :     sub families_implementing_role {
348 :     my($self,$role) = @_;
349 :    
350 :     $self->check_db_role; # lazily tie DB
351 :    
352 :     my $fams = $self->{role_db}->{$role};
353 :     return $fams ? split("\t",$fams) : ();
354 :     }
355 :    
356 :     sub families_with_function {
357 :     my($self,$function) = @_;
358 :    
359 :     $self->check_db_function; # lazily tie DB
360 :    
361 :     my $fams = $self->{function_db}->{$function};
362 :     return $fams ? split("\t",$fams) : ();
363 :     }
364 :    
365 :     sub family_containing_prot {
366 :     my($self,$prot) = @_;
367 :     my @fams = $self->families_containing_prot($prot);
368 :     return (@fams > 0) ? $fams[0] : undef;
369 :     }
370 :    
371 :     sub families_containing_prot {
372 :     my($self,$prot) = @_;
373 :    
374 :     $self->check_db_prot_to_fams; # lazily tie DB
375 :     my $fams = $self->{prot_to_fams_db}->{$prot};
376 :     return $fams ? split("\t",$fams) : ();
377 :     }
378 :    
379 :     sub families_in_genome {
380 :     my($self,$genome) = @_;
381 :    
382 :     $self->check_db_genome_to_fams; # lazily tie DB
383 :    
384 :     my $fams = $self->{genome_to_fams_db}->{$genome};
385 :     return $fams ? split("\t",$fams) : ();
386 :     }
387 :    
388 :     sub families_containing_md5 {
389 :     my ($self,$md5) = @_;
390 :    
391 :     $self->check_db_md5_to_fams; # lazily tie DB
392 :    
393 :     my $fams = $self->{md5_fams_db}->{$md5};
394 :     if ($fams){
395 :     return $fams;
396 :     }
397 :     else{
398 :     return undef;
399 :     }
400 :     }
401 :    
402 : arodri7 1.2 sub proteins_containing_md5 {
403 :     my ($self,$md5) = @_;
404 :    
405 :     $self->check_db_md5_to_prots; # lazily tie DB
406 :    
407 :     my $prots = $self->{md5_prots_db}->{$md5};
408 :     if ($prots){
409 :     return $prots;
410 :     }
411 :     else{
412 :     return undef;
413 :     }
414 :     }
415 :    
416 : arodri7 1.1 sub all_families {
417 :     my($self) = @_;
418 :    
419 :     return sort map { chomp; $_ } `cut -f1 $self->{dir}/family.functions`;
420 :     }
421 :    
422 :     sub is_prots_in_family {
423 :     my($self,$seqs) = @_;
424 :     my $id_2_protfam = {};
425 :    
426 :     my $protD = $self->{dir};
427 :     my $protfam_file = "$protD/families.3c";
428 :    
429 :     my $contents;
430 :     sysopen(FAM_PROTS, $protfam_file, 0) or die "could not open file '$protfam_file': $!";
431 :     sysread(FAM_PROTS, $contents, 1000000000);
432 :     close(FAM_PROTS) or die "could not close file '$protfam_file': $!";
433 :     my %protfam_prots = map { $_ =~ /^(\S+)\t(\S+)\t/; $2 => $1 } split("\n", $contents);
434 :    
435 :     foreach (@$seqs){
436 :     if ($protfam_prots{$_}){
437 :     $id_2_protfam->{$_} = $protfam_prots{$_};
438 :     }
439 :     }
440 :     return $id_2_protfam;
441 :     }
442 :    
443 :     sub place_in_family {
444 : arodri7 1.2 my($self,$seq,$debug,$loose,$debug_prot,$nuc,$md5_flag) = @_;
445 : arodri7 1.1 my($protfam,$should_be,$sims);
446 :    
447 :     my $old_sep = $/;
448 :     $/ = "\n";
449 :    
450 :     my $dir = $self->{dir};
451 :     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";
452 :     open(TMP,">$tmpF")
453 :     || die "could not open $tmpF";
454 :     print TMP ">query\n$seq\n";
455 :     close(TMP);
456 :    
457 :     # check if sequence is in any protfam by md5 checksum
458 :     my $md5 = Digest::MD5::md5_hex( uc $seq );
459 :     my $fams = $self->families_containing_md5($md5);
460 : arodri7 1.2 if ((!$md5_flag) && ($fams) && ($fams !~ /\,/)){ # do regular check if the checksum returns more than one protfam
461 :     my $got = new ProtFamLite($self,$fams);
462 : arodri7 1.1 ($should_be, $sims) = $got->should_be_member($seq,$debug,$loose,$debug_prot,$nuc);
463 :     return ($got,$sims);
464 :     }
465 :    
466 :     if ($nuc){
467 :     # open(BLAST,"blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g F | head -n 20 |")
468 :     Open(\*BLAST,"$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T |");
469 :     print STDERR "blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T\n" if ($debug);
470 :     }
471 :     else{
472 :     # open(BLAST,"blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastp | head -n 20 |")
473 :     Open(\*BLAST, "$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastp |");
474 :    
475 :     }
476 :    
477 :     my(%seen);
478 :     my $got = undef;
479 :     my $min_sc = 1;
480 :     my $checked = 0;
481 :     while ((! $got) && ($checked < 30) && (defined($_ = <BLAST>)))
482 :     {
483 :     if ($debug) { print STDERR $_ }
484 :     chop;
485 :     my $sim = [split(/\t/,$_)];
486 :     if ($min_sc > $sim->[10]) { $min_sc = $sim->[10] }
487 :     my $fam_id = ($sim->[1] =~ /(.*)-/) ? $1 : "";
488 :     if (! $seen{$fam_id})
489 :     {
490 :     $seen{$fam_id} = 1;
491 :    
492 :     # skip if the e-value is 1 or higher
493 :     if ($sim->[10] >= 1){
494 :     print STDERR "SIM evalue: " . $sim->[10] . "\n" if ($debug);
495 :     $checked++;
496 :     next;
497 :     }
498 :    
499 : arodri7 1.2 $protfam = new ProtFamLite($self,$fam_id);
500 : arodri7 1.1
501 :     if (not defined($protfam))
502 :     {
503 :     next;
504 :     }
505 :     else
506 :     {
507 :     $checked++;
508 :     if ($debug) { print STDERR "checking family $fam_id ",$protfam->family_function,"\n" }
509 :     ($should_be, $sims) = $protfam->should_be_member($seq,$debug,$loose,$debug_prot,$nuc);
510 :     if ($debug) { print STDERR " should_be=$should_be\n" }
511 :     my $psc;
512 :     if ($should_be && defined($psc = &psc($sims->[0])) && ($psc <= $min_sc))
513 :     {
514 :     if ($debug) { print STDERR "min_sc=$min_sc best=$psc\n" }
515 :     $got = $protfam;
516 :     }
517 :     }
518 :     }
519 :     }
520 :    
521 :     if ($debug)
522 :     {
523 :     while ($_ && (defined($_ = <BLAST>))) ### JUST PRINT OUT REMAINING BLASTS AGAINST REPDB
524 :     {
525 :     print STDERR $_;
526 :     }
527 :     }
528 :     $/ = $old_sep;
529 :     close(BLAST);
530 :    
531 :     unlink($tmpF);
532 :     return $got ? ($got,$sims) : (undef,undef);
533 :     }
534 :    
535 :     =head3
536 :     usage: my $out = $protfams->place_in_family_bulk($seqs);
537 :    
538 :     The genereal overview of this function is that it does an initial blast of all
539 :     input sequences against the repdb with one call, and then does individual blast
540 :     calls for each input sequence to determine to which family it should be placed.
541 :    
542 :     input
543 :     the only difference to place_in_family is that the $seqs parameter is a
544 :     list of sequences to be placed in a family.
545 :     The list can be either a number of sequences in fasta format or an array
546 :     of just the sequences to be placed in a family
547 :    
548 :     i.e. ->>>> [">seq1\nADTGGHHH\n", ">seq2\nRARGHTGK\n", ....]
549 :     or
550 :     ["ADTGGHHH", "RARGHTGK", ....]
551 :    
552 :     output
553 :     return reference to a hash where the keys are the input sequence
554 :     and the values are tuples same as what place in family returns: [$family, $sims]
555 :    
556 :     $family is a hash reference to the Protfam it was placed in. If the sequence
557 :     was not placed in a family, then it will be an empty tuple.
558 :    
559 :     $sims is a reference to the similarities of the Protfam selected
560 :    
561 :    
562 :     =cut
563 :    
564 :     sub place_in_family_bulk {
565 :     my($self,$seqs,$debug,$loose,$debug_prot,$nuc) = @_;
566 :     my($protfam,$should_be,$sims);
567 :     my $out = {};
568 :     my @run_blast = ();
569 :    
570 :     my $dir = $self->{dir};
571 :     my $old_sep = $/;
572 :     $/ = "\n";
573 :     my $seq_list = {};
574 :     my $count = 0;
575 :    
576 :     foreach my $fasta (@$seqs){
577 :     my ($id, $header, $seq);
578 :     if ($fasta =~ m/^>(.*?)\n(.*?)\n/){
579 :     #($header,$seq) = $fasta =~ m/^>(.*?)\n(.*?)\n/;
580 :     $header = $1;
581 :     $seq = $2;
582 :     ($id) = split (/\s+/, $header);
583 :     }
584 :     else{
585 :     $id = "fasta".$count;
586 :     $seq = $fasta;
587 :     $count++;
588 :     $fasta = ">$id\n$seq\n";
589 :     }
590 :    
591 :    
592 :     # check if sequence is in any protein fam by md5 checksum
593 :     my $fams;
594 :     if (!$nuc){
595 :     my $md5 = Digest::MD5::md5_hex( uc $seq );
596 :     $fams = $self->families_containing_md5($md5);
597 :     }
598 :     if (($fams) && ($fams !~ /\,/)){ # do regular check if the checksum returns more than one prtofam
599 : arodri7 1.2 my $got = new ProtFamLite($self, $fams);
600 : arodri7 1.1 ($should_be, $sims) = $got->should_be_member($seq,$debug,$loose,$debug_prot,$nuc);
601 :     $out->{$seq} = [$got, $sims];
602 :     #push (@$out, [$seq, $got, $sims]);
603 :     }
604 :     else{
605 :     $seq_list->{$id} = $seq;
606 :     push (@run_blast, $fasta);
607 :     }
608 :     }
609 :    
610 :     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";
611 :     open(TMP,">$tmpF")
612 :     || die "could not open $tmpF";
613 :     print TMP @run_blast;
614 :     close(TMP);
615 :     my @tmp;
616 :    
617 :     if ($nuc){
618 :     @tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T` || die "could not open blast";
619 :     print STDERR "$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T\n" if ($debug);
620 :     }
621 :     else{
622 :     @tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastp` || die "could not open blast";
623 :     print STDERR "$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T\n" if ($debug);
624 :     }
625 :    
626 :     my $prev_id;
627 :     my $seen = {};
628 :     my $got = {};
629 :     my $min_sc = {};
630 :     my $checked = {};
631 :    
632 :     while ( (defined( $_ = shift @tmp)) )
633 :     {
634 :     my $sim = [split(/\t/,$_)];
635 :    
636 :     next if ($got->{$sim->[0]});
637 :     next if ($checked->{$sim->[0]} > 30);
638 :     if (!$min_sc->{$sim->[0]}){
639 :     $min_sc->{$sim->[0]} = 1;
640 :     }
641 :    
642 :     if ($debug) { print STDERR $_ }
643 :     chop;
644 :     if ($min_sc->{$sim->[0]} > $sim->[10]) { $min_sc->{$sim->[0]} = $sim->[10] }
645 :     my $fam_id = ($sim->[1] =~ /(.*)-/) ? $1 : "";
646 :     if (! $seen->{$sim->[0]}->{$fam_id})
647 :     {
648 :     $seen->{$sim->[0]}->{$fam_id} = 1;
649 :    
650 :     # skip if the e-value is 1 or higher
651 :     if ($sim->[10] >= 1){
652 :     #print STDERR "SIM evalue: " . $sim->[10] . "\n" if ($debug);
653 :     $checked->{$sim->[0]}++;
654 :     next;
655 :     }
656 :    
657 : arodri7 1.2 $protfam = new ProtFamLite($self, $fam_id);
658 : arodri7 1.1
659 :     if (not defined($protfam))
660 :     {
661 :     next;
662 :     }
663 :     else
664 :     {
665 :     $checked->{$sim->[0]}++;
666 :     if ($debug) { print STDERR "checking family $fam_id ",$protfam->family_function,"\n" }
667 :     my ($should_be, $sims) = $protfam->should_be_member($seq_list->{$sim->[0]},$debug,$loose,$debug_prot,$nuc);
668 :     if ($debug) { print STDERR " should_be=$should_be\n" }
669 :     my $psc;
670 :     if ($should_be && defined($psc = &psc($sims->[0])) && ($psc <= $min_sc->{$sim->[0]}))
671 :     {
672 :     if ($debug) { print STDERR "min_sc=$min_sc best=$psc\n" }
673 :     $got->{$sim->[0]} = $protfam;
674 :     #push (@$out, [$seq_list->{$sim->[0]}, $got->{$sim->[0]}, $sims]);
675 :     $out->{$seq_list->{$sim->[0]}} = [$got->{$sim->[0]}, $sims];
676 :     }
677 :     }
678 :     }
679 :    
680 :     # if ($debug)
681 :     # {
682 :     # my $blast_line;
683 :     # while ($blast_line && (defined($blast_line = shift @blast_tmp))) ### JUST PRINT OUT REMAINING BLASTS AGAINST REPDB
684 :     # {
685 :     # print STDERR $blast_line;
686 :     # }
687 :     # }
688 :     #
689 :     # if (!$got){
690 :     # push (@$out, [$prev_id, undef, undef]);
691 :     # }
692 :    
693 :     }
694 :    
695 :     $/ = $old_sep;
696 :     # close(BLAST);
697 :    
698 :     foreach my $id (keys %$seq_list)
699 :     {
700 :     if (!$got->{$id})
701 :     {
702 :     #push (@$out, [$seq_list->{$id}, undef, undef]);
703 :     $out->{$seq_list->{$id}} = [undef, undef];
704 :     }
705 :     }
706 :    
707 :     unlink($tmpF);
708 :     return $out;
709 :     }
710 :    
711 :     sub psc {
712 :     my ($sim) = @_;
713 :     return ($sim->[10] =~ /^e-/) ? "1.0" . $sim->[10] : $sim->[10];
714 :     }
715 :    
716 :    
717 :     =head3
718 :     usage: $protfams->family_functions();
719 :    
720 :     returns a hash of all the functions for all protein fams from the family.functions file
721 :    
722 :     =cut
723 :    
724 :     sub family_functions {
725 :     my($self) = @_;
726 :     my $ffD = $self->{dir};
727 :     my $ff_file = "$ffD/family.functions";
728 :    
729 :     my $contents;
730 :     sysopen(FAM_FUNC, $ff_file, 0) or die "could not open file '$ff_file': $!";
731 :     sysread(FAM_FUNC, $contents, 1000000000);
732 :     close(FAM_FUNC) or die "could not close file '$ff_file': $!";
733 :     my %ff_name = map {split(/\t/)} split("\n", $contents);
734 :    
735 :     return \%ff_name;
736 :     }
737 :    
738 :    
739 :    
740 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3