[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.4 - (view) (download) (as text)

1 : overbeek 1.1 use FIG;
2 : arodri7 1.3 use Digest::MD5;
3 : overbeek 1.1 my $fig = new FIG;
4 : arodri7 1.2 my $dbf = $fig->db_handle;
5 : overbeek 1.1 use FIGO;
6 :     my $figO = new FIGO;
7 :    
8 :     use strict;
9 :    
10 :     my $sim_thresh = 1.0e-10;
11 : arodri7 1.2 my($usage,$fams,$outD,$proc,$log_file);
12 :     $usage = "usage: generate_partition_worker Families OutDir Process < FIGfams > Partitions";
13 :    
14 :     while ( $ARGV[0] =~ /^-/ )
15 :     {
16 :     $_ = shift @ARGV;
17 :     if ($_ =~ s/^-log//) { $log_file = ($_ || shift @ARGV) }
18 :     else { print STDERR "Bad flag: '$_'\n$usage"; exit 1 }
19 :     }
20 : overbeek 1.1
21 :     (
22 :     ($fams = shift @ARGV) && open(FAMS,"<$fams") &&
23 :     ($outD = shift @ARGV) &&
24 :     ($proc = shift @ARGV)
25 :     )
26 :     || die $usage;
27 :    
28 :     open(OUT,">>$outD/$proc") || die "could not open $outD/$proc";
29 :     my $ofh = select OUT; $| = 1; select $ofh;
30 :     my %done_fam;
31 :    
32 : arodri7 1.3 my(%fam_to_pegs,%peg_to_fams,%peg_to_md5);
33 : overbeek 1.1 while (defined($_ = <FAMS>))
34 :     {
35 :     if ($_ =~ /^(\S+)\t(\S+)/)
36 :     {
37 :     my($fam,$peg) = ($1,$2);
38 : arodri7 1.3 if (((is_fid($peg)) && (! $fig->is_deleted_fid($peg)) && (! $done_fam{$fam})) ||
39 :     (!is_fid($peg)))
40 : overbeek 1.1 {
41 :     push(@{$fam_to_pegs{$fam}},$peg);
42 :     push(@{$peg_to_fams{$peg}},$fam);
43 : arodri7 1.3 #my $md5 = Digest::MD5::md5_hex( uc $seq );
44 :     #$peg_to_md5{$peg}=$md5;
45 : overbeek 1.1 }
46 :     }
47 :     }
48 :     close(FAMS);
49 :     #print STDERR "loaded FIGfamsd\n";
50 :    
51 :     my $peg;
52 :     while (defined($_ = <STDIN>))
53 :     {
54 : arodri7 1.2 if (!$fig){
55 :     $fig = new FIG;
56 :     $dbf = $fig->db_handle;
57 :     }
58 :    
59 : arodri7 1.3 chomp $_;
60 : arodri7 1.2 my $cond = "fam = '".$_."'";
61 :     my $x = $dbf->SQL("select fam from tmp_sync_partitions where $cond");
62 :    
63 : disz 1.4 # if (($_ =~ /(FIG\d+)/) && (! $done_fam{$1}))
64 :     if (($_ =~ /(FIG\d+)/) && (! $done_fam{$1}) && (! scalar @$x > 0))
65 : overbeek 1.1 {
66 :     my $ff = $1;
67 : arodri7 1.2 print STDERR "$proc processing $ff\n";
68 : overbeek 1.1
69 :     my $todo = { $ff => 1 };
70 :     my $loaded_fams = {};
71 : arodri7 1.3 my $state = [$fig,$figO,\%done_fam,$todo,$loaded_fams,\%fam_to_pegs,\%peg_to_fams,$sim_thresh,\%peg_to_md5];
72 : overbeek 1.1 &load_fam($ff,$state);
73 :    
74 :     my @todo;
75 :     while ((@todo = keys(%$todo)) > 0)
76 :     {
77 :     my $fam = $todo[0];
78 :     &process($fam,$state);
79 :     # print STDERR &Dumper(['todo',$todo]);
80 :     }
81 :     }
82 :     }
83 :    
84 :     sub load_fam {
85 :     my($ff,$state) = @_;
86 :    
87 : arodri7 1.3 my($fig,$figO,$done_fam,$todo,$loaded_fams,$fam_to_pegs,$peg_to_fams,$sim_thresh,$peg_to_md5) = @$state;
88 : overbeek 1.1
89 :     my $loaded_fam = $loaded_fams->{$ff};
90 :     if (! $loaded_fam)
91 :     {
92 :     my $pegs = $fam_to_pegs{$ff};
93 :     my @clean_pegs = ();
94 : arodri7 1.3 my $id_flag = 'fid';
95 : overbeek 1.1 foreach my $peg (@$pegs)
96 :     {
97 : arodri7 1.3 my $featureO = new FeatureO($figO,$peg) if (&is_fid($peg));
98 : overbeek 1.1 if ($featureO && (! $fig->possibly_truncated($peg)) && (! $featureO->possible_frameshift))
99 :     {
100 :     push(@clean_pegs,$peg);
101 :     }
102 : arodri7 1.3 #elsif (! is_fid($peg))
103 :     #{
104 :     # # get the md5 sequence
105 :     # my $md5 = $peg_to_md5->{$peg};
106 :     # push(@clean_pegs, "gnl|md5|".$md5);
107 :     # $id_flag = 'non_fid';
108 :     #}
109 : overbeek 1.1 }
110 :    
111 :     my @sims = ();
112 : arodri7 1.3 my @presims = ();
113 :     if ($id_flag eq 'fid')
114 :     {
115 :     @presims = $fig->sims(\@clean_pegs,100000,$sim_thresh,'fig');
116 :     }
117 :     else {
118 :     @presims = $fig->sims(\@clean_pegs,100000,$sim_thresh,'raw');
119 :     }
120 :    
121 :     foreach my $sim (@presims)
122 : overbeek 1.1 {
123 :     my $match1 = $sim->e1 - $sim->b1;
124 :     my $match2 = $sim->e2 - $sim->b2;
125 :     if ((($match1 / $sim->ln1) > 0.7) &&
126 :     (($match2 / $sim->ln2) > 0.7))
127 :     {
128 :     push(@sims,$sim);
129 :     }
130 :     }
131 :    
132 :     my $others_not_in_fam = {};
133 :     my $other_fams = {};
134 :     my $worst_hit = 1000;
135 :     foreach my $sim (sort { ($a->id1 cmp $b->id1) or ($b->bsc <=> $a->bsc) } @sims)
136 :     {
137 :     my $id1 = $sim->id1;
138 :     my $id2 = $sim->id2;
139 :     my $in = $peg_to_fams{$id2};
140 :     my $matched1 = ($sim->e1 + 1) - $sim->b1;
141 :     if (! $in)
142 :     {
143 :     $others_not_in_fam->{$id2} = 1;
144 :     }
145 :     else
146 :     {
147 :     foreach my $fam_in (@$in)
148 :     {
149 :     if ($fam_in eq $ff)
150 :     {
151 :     # print STDERR &Dumper([$worst_hit,$sim->bsc/$matched1]);
152 :     if (((! defined($worst_hit)) || (($sim->bsc / $matched1) < $worst_hit)))
153 :     {
154 :     $worst_hit = sprintf("%5.3f",($sim->bsc / $matched1));
155 :     }
156 :     }
157 :     else
158 :     {
159 :     if (! $other_fams->{$fam_in})
160 :     {
161 :     # print STDERR "$ff - $fam_in $id1 $id2\n";
162 :     }
163 :     $other_fams->{$fam_in} = 1;
164 :     }
165 :     }
166 :     }
167 :     }
168 :     $loaded_fams->{$ff} = [$worst_hit,$other_fams,$others_not_in_fam];
169 :     }
170 :     # print STDERR &Dumper(["loaded $ff",$loaded_fams->{$ff}]);
171 :     }
172 :    
173 :     sub process {
174 :     my($fam,$state) = @_;
175 :    
176 : arodri7 1.3 my($fig,$figO,$done_fam,$todo,$loaded_fams,$fam_to_pegs,$peg_to_fams,$sim_thresh,$peg_to_md5) = @$state;
177 : overbeek 1.1 my $loaded = $loaded_fams->{$fam};
178 :     my($worst_hit,$other_fams,$others_not_in_fam) = @$loaded;
179 :    
180 :     my @contained_in = ();
181 :     my @not_contained_in = ();
182 :     foreach my $famO (keys(%$other_fams))
183 :     {
184 :     &load_fam($famO,$state);
185 :     my $other_fams2 = $loaded_fams->{$famO}->[1];
186 :     my $contained = 1;
187 :     foreach my $famO2 (keys(%$other_fams2))
188 :     {
189 :     if (($famO2 ne $fam) && (! $other_fams->{$famO2}))
190 :     {
191 :     $contained = 0;
192 :     }
193 :     }
194 :    
195 :     if ($contained)
196 :     {
197 :     push(@contained_in,$famO);
198 :     }
199 :     else
200 :     {
201 :     push(@not_contained_in,$famO);
202 :     }
203 :     }
204 :    
205 :     my $set_to_write = [$fam,@contained_in];
206 :     &write_set($set_to_write,$state,\@not_contained_in);
207 :     foreach my $famD (@$set_to_write)
208 :     {
209 :     if ($todo->{$famD}) { delete $todo->{$famD} }
210 :     $done_fam->{$famD} = 1;
211 : arodri7 1.2
212 :     # add to db_handle to show that figfam is done for all other processes
213 :     if (!$fig){
214 :     $fig = new FIG;
215 :     $dbf = $fig->db_handle;
216 :     }
217 :     $dbf->SQL("insert into tmp_sync_partitions (fam) values ('$famD')");
218 : overbeek 1.1 }
219 :    
220 :     foreach my $famO (@not_contained_in)
221 :     {
222 : arodri7 1.2 if (!$fig){
223 :     $fig = new FIG;
224 :     $dbf = $fig->db_handle;
225 :     }
226 :    
227 :     my $cond = "fam = '".$famO."'";
228 :     my $x = $dbf->SQL("select fam from tmp_sync_partitions where $cond");
229 :    
230 :     # if (! $done_fam->{$famO})
231 :     if ( (! $done_fam->{$famO}) && (! scalar @$x > 0) )
232 : overbeek 1.1 {
233 :     $todo->{$famO} = 1;
234 :     }
235 :     }
236 :     }
237 :    
238 :     sub write_set {
239 :     my($set,$state,$not_contained_in) = @_;
240 :    
241 : arodri7 1.3 my($fig,$figO,$done_fam,$todo,$loaded_fams,$fam_to_pegs,$peg_to_fams,$sim_thresh,$peg_to_md5) = @$state;
242 : overbeek 1.1 my %all_pegs;
243 :    
244 :     my @fams_with_worst = ();
245 :     foreach my $fam (@$set)
246 :     {
247 :     my $loaded = $loaded_fams->{$fam};
248 :     my($worst_hit,$other_fams,$others_not_in_fam) = @$loaded;
249 :     push(@fams_with_worst,"$fam,$worst_hit");
250 :     my $pegs = $fam_to_pegs->{$fam};
251 :     foreach my $peg (@$pegs,keys(%$others_not_in_fam))
252 :     {
253 :     $all_pegs{$peg} = 1;
254 :     }
255 :     }
256 :    
257 :     foreach my $fam (@$not_contained_in)
258 :     {
259 :     my $loaded = $loaded_fams->{$fam};
260 :     my($worst_hit,$other_fams,$others_not_in_fam) = @$loaded;
261 :     my $pegs = $fam_to_pegs->{$fam};
262 :     foreach my $peg (@$pegs,keys(%$others_not_in_fam))
263 :     {
264 :     $all_pegs{$peg} = 1;
265 :     }
266 :     }
267 :    
268 :     foreach $peg (sort { &FIG::by_fig_id($a,$b) } keys(%all_pegs))
269 :     {
270 :     print OUT "$peg\n";
271 :     }
272 :     print OUT "//\t",join(" ",@fams_with_worst),"\n";
273 :     }
274 :    
275 : arodri7 1.3 sub is_fid {
276 :     my ($fid) = @_;
277 :    
278 :     if ($fid =~ /^fig\|\d+\.\d+\.peg.\d+/){
279 :     return 1;
280 :     }
281 :     return 0;
282 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3