[Bio] / FigKernelScripts / indirect_functional_coupling_mcauliffe.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/indirect_functional_coupling_mcauliffe.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 ########################################################################
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 : overbeek 1.2 use strict;
20 : overbeek 1.1 use FIG;
21 :     my $fig = new FIG;
22 :     use FIGStatisticalModels;
23 :     my $fs = new FIGStatisticalModels;
24 :    
25 : overbeek 1.2 # $usage = "usage: indirect_functional_coupling < PEGS > PEG-PEG-SC-PEG3-PEG4-BBH13-FC34-BBH42";
26 : overbeek 1.1
27 : overbeek 1.4 while (<STDIN>)
28 : overbeek 1.1 {
29 : overbeek 1.4 chomp;
30 :     my $peg;
31 :     if (/^fig\|/) {$peg=$_}
32 :     else {$peg = $fig->by_alias($_)}
33 :     next unless ($fig->is_real_feature($peg));
34 :     foreach my $tuple (sort { $b->[1] <=> $a->[1] } &indirectly_coupled($fig,$peg))
35 : overbeek 1.1 {
36 : overbeek 1.4 print join("\t",($peg,@$tuple)),"\n";
37 : overbeek 1.1 }
38 :     }
39 :    
40 : overbeek 1.4
41 : overbeek 1.1 sub indirectly_coupled {
42 :     my($fig,$peg1) = @_;
43 :     my(%best_sc);
44 :    
45 :     my $genome1 = &FIG::genome_of($peg1);
46 :    
47 :     foreach my $bbh13 ($fig->bbhs($peg1))
48 :     {
49 : overbeek 1.4 my ($peg3, $bbh_eval13, $bbh_score13) = @$bbh13;
50 :     unless (defined $bbh_score13) {die "Couldn't get a normalized bit score for the BBHs. Are your BBHs upto date?"}
51 :    
52 :     foreach my $fc34 ($fig->coupled_to($peg3))
53 :     {
54 :     my($peg4,$fc_sc34) = @$fc34;
55 : overbeek 1.1 foreach my $bbh42 (grep { &FIG::genome_of($_->[0]) eq $genome1 } $fig->bbhs($peg4))
56 :     {
57 : overbeek 1.4 my($peg2,$bbh_eval24, $bbh_score24) = @$bbh42;
58 : overbeek 1.5 if ($peg1 ne $peg2)
59 : overbeek 1.1 {
60 : overbeek 1.5 unless (defined $bbh_score24) {die "Couldn't get a normalized bit score for the BBHs. Are your BBHs upto date?"}
61 :     my $psc = &score($fig,$peg1,$peg2,$peg3,$peg4,$bbh_score13,$fc_sc34,$bbh_score24);
62 :     my $best = $best_sc{$peg2};
63 :     if ((! defined($best)) || ($best->[0] < $psc))
64 :     {
65 :     $best_sc{$peg2} = [$psc,$peg3,$peg4,$bbh_score13,$fc_sc34,$bbh_score24];
66 :     }
67 : overbeek 1.1 }
68 :     }
69 :     }
70 :     }
71 :     return map { [$_,@{$best_sc{$_}}] } keys(%best_sc);
72 :     }
73 :    
74 :     sub score {
75 : overbeek 1.4 my($fig,$peg1,$peg2,$peg3,$peg4,$sc_bbh13,$fc_sc34,$sc_bbh24) = @_;
76 : overbeek 1.1
77 : overbeek 1.4 unless (defined $sc_bbh13 && $fc_sc34 && $sc_bbh24) {print STDERR "ERROR: $sc_bbh13 && $fc_sc34 && $sc_bbh24\n"; return 0}
78 : overbeek 1.1 my $p1 = &score_bbh($fig,$peg1,$sc_bbh13);
79 : overbeek 1.4 my $p2 = $fs->fc_score2probability($fc_sc34);
80 :     my $p3 = &score_bbh($fig,$peg4,$sc_bbh24);
81 : overbeek 1.1 my $sc = sprintf("%0.2f",$p1 * $p2 * $p3);
82 :     return $sc;
83 :     }
84 :    
85 :     sub score_bbh {
86 : overbeek 1.4 my($fig,$peg,$bsc) = @_;
87 : overbeek 1.1 my($i,$bins);
88 :    
89 : overbeek 1.4 my $func = scalar($fig->function_of($peg));
90 :     my $tc=0;
91 :     $tc = 1 if ($func =~ /aminotransferase|system|component|oxidase|regulator|cytochrome|specific|permease|transcriptional|transport|dehydrogenase/i);
92 : overbeek 1.1
93 : overbeek 1.4 return $fs->fc_bitscore2probability($bsc, $tc);
94 :    
95 : overbeek 1.1 }
96 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3