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

Annotation of /FigKernelScripts/partition_prots.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (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 :     my $i;
24 :     for ($i=0; ($i < @all); $i++)
25 :     {
26 :     if (! $seen_hash{$all[$i]})
27 :     {
28 : overbeek 1.2 print STDERR "processing $all[$i], $i of $#all\n";
29 : overbeek 1.1 &process($fig,$all[$i],\%seen_hash,\%prot_hash,$max_sz,\$n);
30 :     }
31 : overbeek 1.3 else
32 :     {
33 :     print STDERR "seen $all[$i]\n";
34 :     }
35 : overbeek 1.1 }
36 :    
37 :     sub process {
38 : overbeek 1.2 my($fig,$prot,$seen,$prot_hash,$max_sz,$nP) = @_;
39 : overbeek 1.1
40 :     my %in;
41 :     $in{$prot} = 1;
42 :     my %closest;
43 :     $closest{$prot} = 0;
44 :    
45 :     my $sz = 1;
46 :     my $peg = $prot;
47 :     while (($sz < $max_sz) && defined($peg))
48 :     {
49 :     $$nP--;
50 :     delete $closest{$peg};
51 :     my @sims = $fig->sims($peg,10000,1,'raw');
52 :     foreach my $sim (@sims)
53 :     {
54 :     my $id2 = $sim->id2;
55 :     if ($prot_hash->{$id2})
56 :     {
57 :     if (! $in{$id2})
58 :     {
59 :     $in{$id2} = 1;
60 :     $sz++;
61 :     $closest{$id2} = $sim->bsc;
62 :     }
63 :     elsif (defined($closest{$id2}) && ($closest{$id2} < $sim->bsc))
64 :     {
65 :     $closest{$id2} = $sim->bsc;
66 :     }
67 :     }
68 :     }
69 :     $seen->{$peg} = 1;
70 :     print "IN\t$prot\t$peg\n";
71 : overbeek 1.2 $peg = &the_closest(\%closest,$seen);
72 : overbeek 1.1 }
73 :     print STDERR "$$nP left\n";
74 :     foreach $_ (sort keys(%in))
75 :     {
76 :     print join("\t",('SET',$prot,$_)),"\n";
77 :     }
78 :     }
79 :    
80 : overbeek 1.2 sub the_closest {
81 :     my($closest,$seen) = @_;
82 : overbeek 1.1
83 :     my($peg,$x,$peg1,$x1);
84 :     while (($peg1,$x1) = each %$closest)
85 :     {
86 : overbeek 1.3 if ((! $seen->{$peg1}) && ((! defined($peg1)) || ($x1 > $x)))
87 : overbeek 1.1 {
88 :     $peg = $peg1;
89 :     $x = $x1;
90 :     }
91 :     }
92 :     return $peg;
93 :     }
94 :    
95 :     sub build_hashes {
96 : overbeek 1.3 my($prot_hash,$seen_hash,$sz) = @_;
97 : overbeek 1.1
98 :     my($prot_hash_tie,$seen_hash_tie);
99 : overbeek 1.3 # $prot_hash_tie = tie %$prot_hash, 'DB_File',"partition_prot_hash_$sz.db",O_RDWR,0666,$DB_BTREE;
100 :     # $seen_hash_tie = tie %$seen_hash, 'DB_File',"seen_hash_$sz.db",O_RDWR,0666,$DB_BTREE;
101 : overbeek 1.2 # return;
102 : overbeek 1.1
103 : overbeek 1.3 $prot_hash_tie = tie %$prot_hash, 'DB_File',"partition_prot_hash_$sz.db",O_CREAT,0666,$DB_BTREE;
104 : overbeek 1.1 $prot_hash_tie || die "tie failed";
105 :    
106 : overbeek 1.3 $seen_hash_tie = tie %$seen_hash, 'DB_File',"seen_hash_$sz.db",O_CREAT,0666,$DB_BTREE;
107 : overbeek 1.1 $seen_hash_tie || die "tie failed";
108 :    
109 :     open(SYMS,"<$FIG_Config::global/peg.synonyms")
110 :     || die "could not open peg.synonyms";
111 :     while (defined($_ = <SYMS>))
112 :     {
113 :     chop;
114 :     my($head,$rest) = split(/\t/,$_);
115 :     if ($rest =~ /fig\|/)
116 :     {
117 :     my @fig = map { ($_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/) ? $1 : () } split(/;/,$rest);
118 :     @fig = grep { $fig->is_complete(&FIG::genome_of($_)) } @fig;
119 :     if (@fig > 0)
120 :     {
121 :     $head =~ s/,\d+$//;
122 : overbeek 1.3 $prot_hash->{$head} = 1;
123 : overbeek 1.1 foreach my $peg (@fig)
124 :     {
125 : overbeek 1.3 $seen_hash->{$peg} = 1;
126 : overbeek 1.1 }
127 :     }
128 :     }
129 :     }
130 :     foreach my $genome ($fig->genomes('complete'))
131 :     {
132 :     foreach my $peg ($fig->all_features($genome,'peg'))
133 :     {
134 : overbeek 1.3 if (! $seen_hash->{$peg})
135 : overbeek 1.1 {
136 : overbeek 1.3 $prot_hash->{$peg} = 1;
137 : overbeek 1.1 }
138 :     }
139 :     }
140 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3