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

Annotation of /FigKernelScripts/FFB2_get_oligos2.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 # -*- perl -*-
2 :    
3 :     use strict;
4 :     use Data::Dumper;
5 :     use Carp;
6 :     use Cache::Memcached::Fast;
7 :     use IPC::Run qw(start finish reap_nb run);
8 :     use Proc::ParallelLoop;
9 :     use IPC::Open3;
10 :     use POSIX ":sys_wait_h";
11 :     use FileHandle;
12 :    
13 :     my $usage = "usage: FFB2_get_oligos families.2c family.functions function.index outdir [memcache-host memcache-port]";
14 :    
15 : olson 1.4 my $min_sz = 7;
16 :     my $max_sz = 12;
17 : olson 1.1 my($fam2cF,$famFuncF,$famFuncIdx, $out_dir);
18 :    
19 :     (
20 :     ($fam2cF = shift @ARGV) && open(FAM2C,"<$fam2cF") &&
21 :     ($famFuncF = shift @ARGV) &&
22 :     ($famFuncIdx = shift @ARGV) &&
23 :     ($out_dir = shift @ARGV)
24 :     )
25 :     || die $usage;
26 :    
27 :     open(FFIDX, ">", $famFuncIdx) or die "Cannot write $famFuncIdx:$ !";
28 :    
29 :     -d $out_dir || mkdir $out_dir || die "Cannot mkdir $out_dir: $!";
30 :    
31 :     my $mchost = shift;
32 :     my $mcport = shift;
33 :     my $mc;
34 :    
35 :     if ($mchost && $mcport)
36 :     {
37 :     $mc = new Cache::Memcached::Fast({ servers => ["$mchost:$mcport"] } );
38 :     }
39 :    
40 :     #
41 :     # We write $chunksize oligos to a sort and later merge the output files into a pipe to usable_motifs.
42 :     #
43 :    
44 :     my $chunksize = 1_000_000;
45 :     my @cmd = ("sort", "-S", "100M");
46 :    
47 :     my $tmp_dir = "$FIG_Config::temp/ff_bundle.$$";
48 :     -d $tmp_dir || mkdir $tmp_dir || die "Cannot mkdir $tmp_dir: $!";
49 :    
50 :     use FIG;
51 :     my $fig = new FIG;
52 : disz 1.2 my %fam_functions = map { $_ =~ /^(FIG\d+)\t(\S.*\S)/; $1 => $2 } `cat $famFuncF`;
53 : olson 1.1 my $nextF = 0;
54 :     my $x = <FAM2C>;
55 :     my %function_index;
56 :     while ($x && ($x =~ /^(\S+)/))
57 :     {
58 :     my $fam = $1;
59 :     my $fam_function = $fam_functions{$fam};
60 :     my $fam_functionI = &fam_function_index(\%function_index,\$nextF,$fam_function);
61 :     while ($x && ($x =~ /^(\S+)\t(\S+)/) && ($1 eq $fam))
62 :     {
63 :     my $peg = $2;
64 :     if (my $seq = &get_translation($fig, $peg))
65 :     {
66 :     my $i;
67 :     my $ln = length($seq);
68 :     for($i=0; ($i <= ($ln - $min_sz)); $i++)
69 :     {
70 :     my $oligo = uc substr($seq,$i,($i < ($ln-$max_sz)) ? $max_sz : ($ln-$i));
71 :     if ($oligo =~ /^([ACDEFGHIKLMNPQRSTVWY])[ACDEFGHIKLMNPQRSTVWY]*$/)
72 :     {
73 : olson 1.3 my $offset = $ln - $i;
74 :     bundle_write($1, $oligo . "\t" . $fam_functionI . "\t" . $fam . "\tOFF" . $offset . "\n");
75 : olson 1.1 }
76 :     }
77 :     }
78 :     $x = <FAM2C>;
79 :     }
80 :     }
81 :     close(FAM2C);
82 :     close FFIDX;
83 :    
84 :     bundle_finish();
85 :    
86 :     my %n_written;
87 :     my %file_idx;
88 :     my %files;
89 :     my %fh;
90 :     my @pending;
91 :    
92 :     sub bundle_write
93 :     {
94 :     my($char, $str) = @_;
95 :    
96 :     my $ent = $fh{$char};
97 :    
98 :     if ($n_written{$char} >= $chunksize)
99 :     {
100 :     bundle_close($char);
101 :     undef $ent;
102 :     bundle_check();
103 :     $n_written{$char} = 0;
104 :     }
105 :     if (!$ent)
106 :     {
107 :     $ent = bundle_open($char);
108 :     }
109 :     my $fh = $ent->{fh};
110 :     print $fh $str;
111 :     $n_written{$char}++;
112 :     }
113 :    
114 :     sub bundle_open
115 :     {
116 :     my($char) = @_;
117 :    
118 :     my $idx = $file_idx{$char} + 0;
119 :     my $file = sprintf("$tmp_dir/bundle.$char.%05d", $idx);
120 :     print "Write to $file\n";
121 :     $file_idx{$char} = $idx + 1;
122 :    
123 :     push @{$files{$char}}, $file;
124 :    
125 :     my($rpipe, $wpipe);
126 :     pipe($rpipe, $wpipe);
127 :    
128 :     my $pid = fork;
129 :     if ($pid == 0)
130 :     {
131 :     open(STDIN, "<&", $rpipe) or die "Cannot dup stdin: $!";
132 :     close($rpipe);
133 :     close($wpipe);
134 :     open(STDOUT, ">", $file) or die "Cannot write $file: $!";
135 :     exec(@cmd);
136 :     die "exec failed: $!";
137 :     }
138 :    
139 :     close($rpipe);
140 :     my $ent = { fh => $wpipe, file => $file, pid => $pid };
141 :     $fh{$char} = $ent;
142 :     return $ent;
143 :     }
144 :    
145 :     sub bundle_close
146 :     {
147 :     my($char) = @_;
148 :     my $ent = $fh{$char};
149 :     return unless $ent;
150 :     $ent->{fh}->close;
151 :     delete $fh{$char};
152 :     push @pending, $ent;
153 :     }
154 :    
155 :     sub bundle_check
156 :     {
157 :     my @np;
158 :     for my $ent (@pending)
159 :     {
160 :     my $r = waitpid($ent->{pid}, WNOHANG);
161 :     if ($r)
162 :     {
163 :     print "Wait $ent->{pid} returns $r err=$?\n";
164 :     }
165 :     else
166 :     {
167 :     push(@np, $ent);
168 :     }
169 :     }
170 :     @pending = @np;
171 :     }
172 :    
173 :     sub bundle_finish
174 :     {
175 :     for my $char (keys %fh)
176 :     {
177 :     bundle_close($char);
178 :     }
179 :    
180 :     for my $ent (@pending)
181 :     {
182 :     print "Wait for $ent->{pid}\n";
183 :     my $r = waitpid($ent->{pid}, 0);
184 :     if ($r)
185 :     {
186 :     print "Wait $ent->{pid} returns $r err=$?\n";
187 :     }
188 :     }
189 :    
190 :     # pareach [ keys %files ], sub {
191 :     # my $char = shift;
192 :     for my $char (keys %files)
193 :     {
194 :     my $out = "$out_dir/all.$char";
195 :    
196 :     my @files = @{$files{$char}};
197 :     my $r = run ["sort", "-m", @files], ">", $out;
198 :     print "run for @files returns $r\n";
199 :     # my $cmd = "ls -l @files; sort -m @files > $out";
200 :     # my $rc = system($cmd);
201 :     # print "rc=$rc: $cmd\n";
202 :     }
203 :     }
204 :    
205 :     sub get_translation
206 :     {
207 :     my($fig, $peg) = @_;
208 :    
209 :     my $tr;
210 :     if ($mc)
211 :     {
212 :     $tr = $mc->get("s:$peg");
213 :     }
214 :     if (!$tr)
215 :     {
216 :     $tr = $fig->get_translation($peg);
217 :     if ($tr && $mc)
218 :     {
219 :     $mc->set("s:$peg", $tr);
220 :     }
221 :     }
222 :     return $tr;
223 :     }
224 :    
225 :     sub fam_function_index {
226 :     my($fam_functions,$nextFP,$fam_function) = @_;
227 :    
228 :     defined($fam_function) || confess "bad";
229 :     my $x = $fam_functions->{$fam_function};
230 :     if (! defined($x))
231 :     {
232 :     $x = $fam_functions->{$fam_function} = $$nextFP;
233 :     print FFIDX "$x\t$fam_function\n";
234 :     $$nextFP++;
235 :     }
236 :     return $x;
237 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3