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

Annotation of /FigKernelScripts/ex_get_subsystem_based_estimates.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use SeedEnv;
2 :     use strict;
3 :     use Carp;
4 :     use Data::Dumper;
5 :     use Getopt::Long;
6 :    
7 :     my $usage = "usage: ex_get_subsystem_based_estimates -d Data\n";
8 :     my $dataD;
9 :     my $rc = GetOptions('d=s' => \$dataD);
10 :    
11 :     if ((! $rc) || (! -d $dataD))
12 :     {
13 :     print STDERR $usage; exit ;
14 :     }
15 :    
16 :     open(GENOME,"<$dataD/genome") || die "$dataD/genome does not exist";
17 :     my $genome = <GENOME>;
18 :     chomp $genome;
19 :     ($genome =~ /^(\d+\.\d+)/) || die "genome: $genome is invalid";
20 :     $genome = $1;
21 :    
22 :     my $pccsF = "$dataD/locus.pccs";
23 :     my $aliasF = "$dataD/aliases";
24 :     open(ALIAS,"<$aliasF") || die "could not open $aliasF";
25 :     my %to_peg;
26 :     while (defined($_ = <ALIAS>))
27 :     {
28 :     chomp;
29 :     my @tmp = split(/\t/,$_);
30 :     my @fig = grep { $_ =~ /^fig\|/ } @tmp;
31 : overbeek 1.2 if (@fig == 1)
32 : overbeek 1.1 {
33 : overbeek 1.2 foreach my $alias (grep { $_ ne $fig[0] } @tmp)
34 : overbeek 1.1 {
35 :     $to_peg{$alias} = $fig[0];
36 :     }
37 :     }
38 :     }
39 :     close(ALIAS);
40 :    
41 :     my $corrH = {};
42 :     open(PCC,"<$dataD/locus.pccs") || die "could not open $dataD/locus.pccs";
43 :     while (defined($_ = <PCC>))
44 :     {
45 :     if (($_ =~ /^(\S+)\t(\S+)\t(\S+)/) &&
46 :     (my $peg1 = $to_peg{$1}) &&
47 :     (my $peg2 = $to_peg{$2}))
48 :     {
49 :     $corrH->{$peg1}->{$peg2} = $corrH->{$peg2}->{$peg1} = $3;
50 :     }
51 :     }
52 :     close(PCC);
53 :    
54 :     my $sapO = SAPserver->new;
55 :    
56 :     my $genomeH = $sapO->genomes_to_subsystems( -ids => [$genome] );
57 :     my @subs = map { ($_->[1] =~ /^\*?(0|-1)$/) ? () : $_->[0] } @{$genomeH->{$genome}};
58 :    
59 :     my $subH = $sapO->ids_in_subsystems( -subsystems => \@subs,
60 :     -genome => $genome);
61 :     my @subs = sort keys %$subH;
62 :    
63 :     my %bad;
64 :    
65 :     foreach my $sub (@subs)
66 :     {
67 :     my %pegs;
68 :     my @pegs;
69 :    
70 :     my $sub_entry = $subH->{$sub};
71 :     @pegs = ();
72 :     foreach my $role (keys(%$sub_entry))
73 :     {
74 :     my $pegs = $sub_entry->{$role};
75 :     foreach $_ (@$pegs) { $pegs{$_} = 1 }
76 :     }
77 :     @pegs = sort { &SeedUtils::by_fig_id($a,$b) } keys(%pegs);
78 :    
79 :     my @sets = grep { @$_ > 1 } split_on_pc(\@pegs,$corrH);
80 :    
81 :     if (@sets > ((@pegs + 2) / 3))
82 :     {
83 :     $bad{$sub} = 1;
84 :     # print STDERR &Dumper([$sub,\@sets]);
85 :     }
86 :     else
87 :     {
88 :     foreach my $set (@sets)
89 :     {
90 :     if (@$set > 1)
91 :     {
92 :     print join(",",@$set),"\tInSubsystem:$sub\n";
93 :     }
94 :     }
95 :     }
96 :     }
97 :     foreach $_ (keys(%bad))
98 :     {
99 :     print STDERR "bad subsystem\t$_\n";
100 :     }
101 :    
102 :     sub split_on_pc {
103 :     my($pegs,$corrH) = @_;
104 :    
105 :     my @sets = ();
106 :     my %used;
107 :     my $i;
108 :     for ($i=0; ($i < (@$pegs - 1)); $i++)
109 :     {
110 :     if (! $used{$pegs->[$i]})
111 :     {
112 :     my @poss = ($pegs->[$i]);
113 :     my $j;
114 :     for ($j=$i+1; ($j < @$pegs); $j++)
115 :     {
116 :     if (&corr($pegs->[$j],\@poss,$corrH))
117 :     {
118 :     push(@poss,$pegs->[$j]);
119 :     $used{$pegs->[$j]} = 1;
120 :     }
121 :     }
122 :     push(@sets,\@poss);
123 :     }
124 :     }
125 :     return @sets;
126 :     }
127 :    
128 :     sub corr {
129 :     my($peg1,$cluster,$corrH) = @_;
130 :    
131 :     my $sum = 0;
132 :     foreach my $peg2 (@$cluster)
133 :     {
134 :     my $v = $corrH->{$peg1}->{$peg2};
135 :     if ((! defined($v)) || ($v < 0.4)) { return 0 }
136 :     $sum += $v;
137 :     }
138 :     return (($sum / @$cluster) >= 0.7);
139 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3