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

Annotation of /FigKernelScripts/partition_prots.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3