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

Annotation of /FigKernelScripts/generate_partitions_worker.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use FIG;
2 :     my $fig = new FIG;
3 :     use FIGO;
4 :     my $figO = new FIGO;
5 :    
6 :     use strict;
7 :    
8 :     my $sim_thresh = 1.0e-10;
9 :     my($usage,$fams,$outD,$doneF,$proc);
10 :     $usage = "usage: generate_partition_worker Families OutDir DoneF Process < FIGfams > Partitions";
11 :    
12 :     (
13 :     ($fams = shift @ARGV) && open(FAMS,"<$fams") &&
14 :     ($outD = shift @ARGV) &&
15 :     ($doneF = shift @ARGV) &&
16 :     ($proc = shift @ARGV)
17 :     )
18 :     || die $usage;
19 :    
20 :     open(OUT,">>$outD/$proc") || die "could not open $outD/$proc";
21 :     my $ofh = select OUT; $| = 1; select $ofh;
22 :    
23 :     my %done_fam;
24 :     foreach $_ (`cat $doneF`)
25 :     {
26 :     chop;
27 :     $done_fam{$_} = 1;
28 :     }
29 :    
30 :     my(%fam_to_pegs,%peg_to_fams);
31 :     while (defined($_ = <FAMS>))
32 :     {
33 :     if ($_ =~ /^(\S+)\t(\S+)/)
34 :     {
35 :     my($fam,$peg) = ($1,$2);
36 :     if ((! $fig->is_deleted_fid($peg)) && (! $done_fam{$fam}))
37 :     {
38 :     push(@{$fam_to_pegs{$fam}},$peg);
39 :     push(@{$peg_to_fams{$peg}},$fam);
40 :     }
41 :     }
42 :     }
43 :     close(FAMS);
44 :     #print STDERR "loaded FIGfamsd\n";
45 :    
46 :     my $peg;
47 :     while (defined($_ = <STDIN>))
48 :     {
49 :     if (($_ =~ /(FIG\d{6})/) && (! $done_fam{$1}))
50 :     {
51 :     print STDERR "$proc processing $1\n";
52 :     my $ff = $1;
53 :    
54 :     my $todo = { $ff => 1 };
55 :     my $loaded_fams = {};
56 :     my $state = [$fig,$figO,\%done_fam,$todo,$loaded_fams,\%fam_to_pegs,\%peg_to_fams,$sim_thresh];
57 :     &load_fam($ff,$state);
58 :    
59 :     my @todo;
60 :     while ((@todo = keys(%$todo)) > 0)
61 :     {
62 :     my $fam = $todo[0];
63 :     &process($fam,$state);
64 :     # print STDERR &Dumper(['todo',$todo]);
65 :     }
66 :     }
67 :     }
68 :    
69 :     sub load_fam {
70 :     my($ff,$state) = @_;
71 :    
72 :     my($fig,$figO,$done_fam,$todo,$loaded_fams,$fam_to_pegs,$peg_to_fams,$sim_thresh) = @$state;
73 :    
74 :     my $loaded_fam = $loaded_fams->{$ff};
75 :     if (! $loaded_fam)
76 :     {
77 :     my $pegs = $fam_to_pegs{$ff};
78 :     my @clean_pegs = ();
79 :     foreach my $peg (@$pegs)
80 :     {
81 :     my $featureO = new FeatureO($figO,$peg);
82 :     if ($featureO && (! $fig->possibly_truncated($peg)) && (! $featureO->possible_frameshift))
83 :     {
84 :     push(@clean_pegs,$peg);
85 :     }
86 :     }
87 :    
88 :     my @sims = ();
89 :     foreach my $sim ($fig->sims(\@clean_pegs,100000,$sim_thresh,'fig'))
90 :     {
91 :     my $match1 = $sim->e1 - $sim->b1;
92 :     my $match2 = $sim->e2 - $sim->b2;
93 :     if ((($match1 / $sim->ln1) > 0.7) &&
94 :     (($match2 / $sim->ln2) > 0.7))
95 :     {
96 :     push(@sims,$sim);
97 :     }
98 :     }
99 :    
100 :     my $others_not_in_fam = {};
101 :     my $other_fams = {};
102 :     my $worst_hit = 1000;
103 :     foreach my $sim (sort { ($a->id1 cmp $b->id1) or ($b->bsc <=> $a->bsc) } @sims)
104 :     {
105 :     my $id1 = $sim->id1;
106 :     my $id2 = $sim->id2;
107 :     my $in = $peg_to_fams{$id2};
108 :     my $matched1 = ($sim->e1 + 1) - $sim->b1;
109 :     if (! $in)
110 :     {
111 :     $others_not_in_fam->{$id2} = 1;
112 :     }
113 :     else
114 :     {
115 :     foreach my $fam_in (@$in)
116 :     {
117 :     if ($fam_in eq $ff)
118 :     {
119 :     # print STDERR &Dumper([$worst_hit,$sim->bsc/$matched1]);
120 :     if (((! defined($worst_hit)) || (($sim->bsc / $matched1) < $worst_hit)))
121 :     {
122 :     $worst_hit = sprintf("%5.3f",($sim->bsc / $matched1));
123 :     }
124 :     }
125 :     else
126 :     {
127 :     if (! $other_fams->{$fam_in})
128 :     {
129 :     # print STDERR "$ff - $fam_in $id1 $id2\n";
130 :     }
131 :     $other_fams->{$fam_in} = 1;
132 :     }
133 :     }
134 :     }
135 :     }
136 :     $loaded_fams->{$ff} = [$worst_hit,$other_fams,$others_not_in_fam];
137 :     }
138 :     # print STDERR &Dumper(["loaded $ff",$loaded_fams->{$ff}]);
139 :     }
140 :    
141 :     sub process {
142 :     my($fam,$state) = @_;
143 :    
144 :     my($fig,$figO,$done_fam,$todo,$loaded_fams,$fam_to_pegs,$peg_to_fams,$sim_thresh) = @$state;
145 :     my $loaded = $loaded_fams->{$fam};
146 :     my($worst_hit,$other_fams,$others_not_in_fam) = @$loaded;
147 :    
148 :     my @contained_in = ();
149 :     my @not_contained_in = ();
150 :     foreach my $famO (keys(%$other_fams))
151 :     {
152 :     &load_fam($famO,$state);
153 :     my $other_fams2 = $loaded_fams->{$famO}->[1];
154 :     my $contained = 1;
155 :     foreach my $famO2 (keys(%$other_fams2))
156 :     {
157 :     if (($famO2 ne $fam) && (! $other_fams->{$famO2}))
158 :     {
159 :     $contained = 0;
160 :     }
161 :     }
162 :    
163 :     if ($contained)
164 :     {
165 :     push(@contained_in,$famO);
166 :     }
167 :     else
168 :     {
169 :     push(@not_contained_in,$famO);
170 :     }
171 :     }
172 :    
173 :     my $set_to_write = [$fam,@contained_in];
174 :     &write_set($set_to_write,$state,\@not_contained_in);
175 :     foreach my $famD (@$set_to_write)
176 :     {
177 :     if ($todo->{$famD}) { delete $todo->{$famD} }
178 :     $done_fam->{$famD} = 1;
179 :     }
180 :    
181 :     foreach my $famO (@not_contained_in)
182 :     {
183 :     if (! $done_fam->{$famO})
184 :     {
185 :     $todo->{$famO} = 1;
186 :     }
187 :     }
188 :     }
189 :    
190 :     sub write_set {
191 :     my($set,$state,$not_contained_in) = @_;
192 :    
193 :     my($fig,$figO,$done_fam,$todo,$loaded_fams,$fam_to_pegs,$peg_to_fams,$sim_thresh) = @$state;
194 :     my %all_pegs;
195 :    
196 :     my @fams_with_worst = ();
197 :     foreach my $fam (@$set)
198 :     {
199 :     my $loaded = $loaded_fams->{$fam};
200 :     my($worst_hit,$other_fams,$others_not_in_fam) = @$loaded;
201 :     push(@fams_with_worst,"$fam,$worst_hit");
202 :     my $pegs = $fam_to_pegs->{$fam};
203 :     foreach my $peg (@$pegs,keys(%$others_not_in_fam))
204 :     {
205 :     $all_pegs{$peg} = 1;
206 :     }
207 :     }
208 :    
209 :     foreach my $fam (@$not_contained_in)
210 :     {
211 :     my $loaded = $loaded_fams->{$fam};
212 :     my($worst_hit,$other_fams,$others_not_in_fam) = @$loaded;
213 :     my $pegs = $fam_to_pegs->{$fam};
214 :     foreach my $peg (@$pegs,keys(%$others_not_in_fam))
215 :     {
216 :     $all_pegs{$peg} = 1;
217 :     }
218 :     }
219 :    
220 :     foreach $peg (sort { &FIG::by_fig_id($a,$b) } keys(%all_pegs))
221 :     {
222 :     print OUT "$peg\n";
223 :     }
224 :     print OUT "//\t",join(" ",@fams_with_worst),"\n";
225 :     }
226 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3