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

Annotation of /FigKernelScripts/partition_prots.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use DB_File;
2 :     use strict;
3 :    
4 :     use FIG;
5 :     my $fig = new FIG;
6 :    
7 :     my $max_sz;
8 : overbeek 1.2 my $usage = "usage: partition_prots MaxPartSz";
9 : overbeek 1.1 (
10 :     ($max_sz = shift @ARGV)
11 :     )
12 :     || die $usage;
13 :    
14 :    
15 :     my %prot_hash;
16 :     my %seen_hash;
17 : overbeek 1.3 &build_hashes(\%prot_hash,\%seen_hash,$max_sz);
18 : overbeek 1.1
19 :     my @all = sort keys(%prot_hash);
20 :     my $n = @all;
21 :     print STDERR "$n proteins to process\n";
22 :    
23 : overbeek 1.4 my $got_sims = 0;
24 :     my $no_sims = 0;
25 : overbeek 1.1 my $i;
26 :     for ($i=0; ($i < @all); $i++)
27 :     {
28 :     if (! $seen_hash{$all[$i]})
29 :     {
30 : overbeek 1.2 print STDERR "processing $all[$i], $i of $#all\n";
31 : overbeek 1.4 &process($fig,$all[$i],\%seen_hash,\%prot_hash,$max_sz,\$n,\$got_sims,\$no_sims);
32 : overbeek 1.1 }
33 : overbeek 1.3 else
34 :     {
35 :     print STDERR "seen $all[$i]\n";
36 :     }
37 : overbeek 1.1 }
38 :    
39 :     sub process {
40 : overbeek 1.4 my($fig,$prot,$seen,$prot_hash,$max_sz,$nP,$gotP,$nosimsP) = @_;
41 : overbeek 1.1
42 :     my %in;
43 :     $in{$prot} = 1;
44 :     my %closest;
45 :     $closest{$prot} = 0;
46 :    
47 :     my $sz = 1;
48 :     my $peg = $prot;
49 :     while (($sz < $max_sz) && defined($peg))
50 :     {
51 :     $$nP--;
52 :     delete $closest{$peg};
53 :     my @sims = $fig->sims($peg,10000,1,'raw');
54 : overbeek 1.4 if (@sims > 0)
55 :     {
56 :     $$got_sims++;
57 :     }
58 :     else
59 :     {
60 :     $$no_sims++;
61 :     if ($$no_sims > (10 * $$got_sims))
62 :     {
63 :     die "sims do not seem to be working properly";
64 :     }
65 :     }
66 :    
67 : overbeek 1.1 foreach my $sim (@sims)
68 :     {
69 :     my $id2 = $sim->id2;
70 :     if ($prot_hash->{$id2})
71 :     {
72 :     if (! $in{$id2})
73 :     {
74 :     $in{$id2} = 1;
75 :     $sz++;
76 :     $closest{$id2} = $sim->bsc;
77 :     }
78 :     elsif (defined($closest{$id2}) && ($closest{$id2} < $sim->bsc))
79 :     {
80 :     $closest{$id2} = $sim->bsc;
81 :     }
82 :     }
83 :     }
84 :     $seen->{$peg} = 1;
85 :     print "IN\t$prot\t$peg\n";
86 : overbeek 1.2 $peg = &the_closest(\%closest,$seen);
87 : overbeek 1.1 }
88 :     print STDERR "$$nP left\n";
89 :     foreach $_ (sort keys(%in))
90 :     {
91 :     print join("\t",('SET',$prot,$_)),"\n";
92 :     }
93 :     }
94 :    
95 : overbeek 1.2 sub the_closest {
96 :     my($closest,$seen) = @_;
97 : overbeek 1.1
98 :     my($peg,$x,$peg1,$x1);
99 :     while (($peg1,$x1) = each %$closest)
100 :     {
101 : overbeek 1.3 if ((! $seen->{$peg1}) && ((! defined($peg1)) || ($x1 > $x)))
102 : overbeek 1.1 {
103 :     $peg = $peg1;
104 :     $x = $x1;
105 :     }
106 :     }
107 :     return $peg;
108 :     }
109 :    
110 :     sub build_hashes {
111 : overbeek 1.3 my($prot_hash,$seen_hash,$sz) = @_;
112 : overbeek 1.1
113 :     my($prot_hash_tie,$seen_hash_tie);
114 : overbeek 1.3 # $prot_hash_tie = tie %$prot_hash, 'DB_File',"partition_prot_hash_$sz.db",O_RDWR,0666,$DB_BTREE;
115 :     # $seen_hash_tie = tie %$seen_hash, 'DB_File',"seen_hash_$sz.db",O_RDWR,0666,$DB_BTREE;
116 : overbeek 1.2 # return;
117 : overbeek 1.1
118 : overbeek 1.3 $prot_hash_tie = tie %$prot_hash, 'DB_File',"partition_prot_hash_$sz.db",O_CREAT,0666,$DB_BTREE;
119 : overbeek 1.1 $prot_hash_tie || die "tie failed";
120 :    
121 : overbeek 1.3 $seen_hash_tie = tie %$seen_hash, 'DB_File',"seen_hash_$sz.db",O_CREAT,0666,$DB_BTREE;
122 : overbeek 1.1 $seen_hash_tie || die "tie failed";
123 :    
124 :     open(SYMS,"<$FIG_Config::global/peg.synonyms")
125 :     || die "could not open peg.synonyms";
126 :     while (defined($_ = <SYMS>))
127 :     {
128 :     chop;
129 :     my($head,$rest) = split(/\t/,$_);
130 :     if ($rest =~ /fig\|/)
131 :     {
132 :     my @fig = map { ($_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/) ? $1 : () } split(/;/,$rest);
133 :     @fig = grep { $fig->is_complete(&FIG::genome_of($_)) } @fig;
134 :     if (@fig > 0)
135 :     {
136 :     $head =~ s/,\d+$//;
137 : overbeek 1.3 $prot_hash->{$head} = 1;
138 : overbeek 1.1 foreach my $peg (@fig)
139 :     {
140 : overbeek 1.3 $seen_hash->{$peg} = 1;
141 : overbeek 1.1 }
142 :     }
143 :     }
144 :     }
145 :     foreach my $genome ($fig->genomes('complete'))
146 :     {
147 :     foreach my $peg ($fig->all_features($genome,'peg'))
148 :     {
149 : overbeek 1.3 if (! $seen_hash->{$peg})
150 : overbeek 1.1 {
151 : overbeek 1.3 $prot_hash->{$peg} = 1;
152 : overbeek 1.1 }
153 :     }
154 :     }
155 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3