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

Annotation of /FigKernelScripts/partition_prots.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3