[Bio] / FigKernelPackages / ExpressionDir.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 package ExpressionDir;
2 :    
3 :     use Data::Dumper;
4 :     use strict;
5 :     use SeedAware;
6 :     use File::Copy;
7 : olson 1.2 use File::Temp 'tempfile';
8 : olson 1.1 use File::Spec::Functions;
9 :     use base 'Class::Accessor';
10 :     use Carp;
11 :     use Fcntl ':seek';
12 :    
13 :     __PACKAGE__->mk_accessors(qw(genome_dir expr_dir error_dir));
14 :    
15 :     =head3 new
16 :    
17 : olson 1.3 my $edir = ExpressionDir->new($genome_source);
18 : olson 1.1
19 :     Create a new ExpressionDir object to be associated with the given genome directory.
20 :     If a subdirectory ExpressionData does not yet exist, one is created.
21 :    
22 :     =cut
23 :    
24 :     sub new
25 :     {
26 :     my($class, $genome_dir) = @_;
27 :    
28 :     my $edir = catfile($genome_dir, "ExpressionData");
29 :     if (! -d $edir)
30 :     {
31 :     mkdir($edir);
32 :     }
33 :     my $errdir = catfile($genome_dir, 'errors');
34 :     if (! -d $errdir)
35 :     {
36 :     mkdir($errdir);
37 :     }
38 :    
39 :    
40 :     my $self = {
41 :     genome_dir => $genome_dir,
42 :     expr_dir => $edir,
43 :     error_dir => $errdir,
44 :     };
45 :     return bless $self, $class;
46 :     }
47 :    
48 :     sub parse_probe_format_1
49 :     {
50 :     my($self, $in_file, $out_file) = @_;
51 :    
52 :     my($fh);
53 :    
54 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
55 :     my $l = <$fh>;
56 :     chomp $l;
57 :     $l =~ s/\r//;
58 :     my @hdrs = split(/\t/, $l);
59 :     my %hdrs;
60 :     $hdrs{$hdrs[$_]} = $_ for 0..$#hdrs;
61 :    
62 :     my $x_col = $hdrs{"Probe X"};
63 :     my $y_col = $hdrs{"Probe Y"};
64 :     my $seq_col = $hdrs{"Probe Sequence"};
65 :     if (!(defined($x_col) && defined($y_col) && defined($seq_col)))
66 :     {
67 :     close($fh);
68 :     return undef;
69 :     }
70 :    
71 :     my $out;
72 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
73 :    
74 :     while (<$fh>)
75 :     {
76 :     chomp;
77 :     s/\r//g;
78 :     my @flds = split(/\t/,$_);
79 :     my($x,$y,$seq);
80 :     $x = $flds[$x_col];
81 :     $y = $flds[$y_col];
82 :     $seq = $flds[$seq_col];
83 :     my $id = "$x\_$y";
84 :     print $out "$id\t$seq\n";
85 :     }
86 :     close($fh);
87 :     close($out);
88 :     return 1;
89 :     }
90 :    
91 :     sub parse_probe_format_2
92 :     {
93 :     my($self, $in_file, $out_file) = @_;
94 :    
95 :     my($fh);
96 :    
97 :     local $/ = "\n>";
98 :    
99 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
100 :     my $l = <$fh>;
101 :     chomp $l;
102 :     $l =~ s/\r//;
103 :    
104 :     if ($l !~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
105 :     {
106 :     close($fh);
107 :     return undef;
108 :     }
109 :     seek($fh, 0, SEEK_SET);
110 :    
111 :     my $out;
112 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
113 :    
114 :     while (<$fh>)
115 :     {
116 :     chomp;
117 :    
118 :     if ($_ =~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
119 :     {
120 :     if (length($3) < 15)
121 :     {
122 :     close($fh);
123 :     confess "Bad length at line $. of $in_file";
124 :     }
125 :     print $out "$1\_$2\t$3\n";
126 :     }
127 :     else
128 :     {
129 :     confess "Bad input at line $. of $in_file";
130 :     }
131 :     }
132 :     close($out);
133 :     close($fh);
134 :     return 1;
135 :     }
136 :    
137 :     #
138 :     # Our "native" format, used for passing through pre-parsed data.
139 :     #
140 :     sub parse_probe_format_3
141 :     {
142 :     my($self, $in_file, $out_file) = @_;
143 :    
144 :     my($fh);
145 :    
146 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
147 :     my $l = <$fh>;
148 :     chomp $l;
149 :     $l =~ s/\r//;
150 :    
151 :     if ($l !~ /^\d+_\d+\t[ACGT]+$/)
152 :     {
153 :     close($fh);
154 :     return undef;
155 :     }
156 :     seek($fh, 0, SEEK_SET);
157 :    
158 :     my $out;
159 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
160 :    
161 :     while (<$fh>)
162 :     {
163 :     if ($_ =~ /^\d+_\d+\t[ACGT]+$/)
164 :     {
165 :     print $out $_;
166 :     }
167 :     else
168 :     {
169 :     confess "Bad input at line $. of $in_file";
170 :     }
171 :     }
172 :     close($out);
173 :     close($fh);
174 :     return 1;
175 :     }
176 :    
177 :     sub compute_probe_to_peg
178 :     {
179 :     my($self, $probes) = @_;
180 :    
181 :     my $my_probes = catfile($self->expr_dir, "probes.in");
182 :    
183 :     copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
184 :    
185 :     my $probes_fasta = catfile($self->expr_dir, "probes.fasta");
186 :    
187 :     #
188 :     # Attempt to translate probe file.
189 :     #
190 :     my $success;
191 :     for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))
192 :     {
193 :     if ($self->$meth($my_probes, $probes_fasta))
194 :     {
195 :     print STDERR "Translated $probes to $probes_fasta using $meth\n";
196 :     $success = 1;
197 :     last;
198 :     }
199 :     else
200 :     {
201 :     print STDERR "Failed to translate $probes to $probes_fasta using $meth\n";
202 :     }
203 :     }
204 :     if (!$success)
205 :     {
206 :     confess "Could not translate $probes\n";
207 :     }
208 :    
209 :     my $peg_probe_table = catfile($self->expr_dir, 'peg.probe.table');
210 :     my $probe_occ_table = catfile($self->expr_dir, 'probe.occ.table');
211 :    
212 : olson 1.2 my $feature_dir = catfile($self->genome_dir, "Features");
213 :     my @tbls;
214 :     for my $ftype (qw(peg rna))
215 :     {
216 :     my $tfile = catfile($feature_dir, $ftype, 'tbl');
217 :     if (-f $tfile)
218 :     {
219 :     push(@tbls, $tfile);
220 :     }
221 :     }
222 :     if (@tbls == 0)
223 :     {
224 :     confess "Could not find any tbl files in $feature_dir";
225 :     }
226 :    
227 : olson 1.1 system_with_redirect([executable_for("make_probes_to_genes"),
228 :     $probes_fasta,
229 :     catfile($self->genome_dir, 'contigs'),
230 : olson 1.2 $tbls[0],
231 : olson 1.1 $peg_probe_table,
232 : olson 1.2 $probe_occ_table,
233 :     @tbls[1..$#tbls],
234 :     ],
235 : olson 1.1 { stderr => catfile($self->expr_dir, 'problems') });
236 :    
237 :    
238 :     system_with_redirect([executable_for("remove_multiple_occurring_probes"),
239 :     $peg_probe_table,
240 :     $probe_occ_table],
241 :     { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
242 :    
243 :     $self->make_missing_probes($peg_probe_table, $probes_fasta,
244 :     catfile($self->expr_dir, 'probe.no.match'));
245 :     }
246 : olson 1.2
247 : olson 1.1 sub make_missing_probes
248 :     {
249 :     my($self, $probe_table, $probes, $output) = @_;
250 :     open(MATCH,"<", $probe_table) or die "Cannot open $probe_table: $!";
251 :     open(PROBES,"<", $probes) or die "Cannot open $probes: $!";
252 :     open(OUTPUT, ">", $output) or die "Cannot open $output: $!";
253 :     my %locations;
254 :     while(<MATCH>)
255 :     {
256 :     chomp;
257 :     my($peg,$loc)=split "\t";
258 :     $locations{$loc} = $peg;
259 :     }
260 :    
261 :     while(<PROBES>)
262 :     {
263 :     chomp;
264 :     my($loc,$seq) = split "\t";
265 :     print OUTPUT $loc, "\n" if ! exists $locations{$loc};
266 :     }
267 :     close(MATCH);
268 :     close(PROBES);
269 :     close(OUTPUT);
270 :     }
271 :    
272 : olson 1.2 #
273 :     # we don't copy the experiment files in here because
274 :     # they may be very large. This may change.
275 :     #
276 :     # We do copy the cdf.
277 :     #
278 :     sub compute_rma_normalized
279 :     {
280 :     my($self, $cdf_file, $expt_dir) = @_;
281 :    
282 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
283 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
284 :    
285 :     #
286 :     # We need to build the R library for this cdf.
287 :     #
288 :     my($fh, $tempfile) = tempfile();
289 :     #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
290 :    
291 :     my $cdf_path = $self->expr_dir;
292 :     my $libdir = catfile($self->expr_dir, "r_lib");
293 :     -d $libdir or mkdir $libdir;
294 :     my $pkgdir = catfile($self->expr_dir, "r_pkg");
295 :     -d $pkgdir or mkdir $pkgdir;
296 :    
297 :     print $fh "library(makecdfenv);\n";
298 :     print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
299 :     close($fh);
300 :     system("cat", $tempfile);
301 :     system("Rscript", $tempfile);
302 :     system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
303 :    
304 :     local($ENV{R_LIBS}) = $libdir;
305 :     my @cmd = (executable_for("RunRMA"),
306 :     "data",
307 :     catfile($self->expr_dir, "peg.probe.table"),
308 :     catfile($self->expr_dir, "probe.no.match"),
309 :     $expt_dir,
310 :     $self->expr_dir);
311 :    
312 :     my $rc = system_with_redirect(\@cmd);
313 :     # my $rc = system_with_redirect(\@cmd, { stderr => "rma.stderr" });
314 :     if ($rc != 0)
315 :     {
316 :     confess("Error running RMA: rc=$rc cmd=@cmd");
317 :     }
318 :     my $output = catfile($self->expr_dir, "rma_normalized.tab");
319 :     if (! -f $output)
320 :     {
321 :     confess("Output file $output was not generated");
322 :     }
323 :     }
324 :    
325 :     sub compute_atomic_regulons
326 :     {
327 :     my($self) = @_;
328 :    
329 :     my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
330 :     my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
331 :     my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
332 :     }
333 :    
334 : olson 1.3 sub get_experiment_names
335 :     {
336 :     my($self) = @_;
337 :     my $f = catfile($self->expr_dir, "experiment.names");
338 :     my $fh;
339 :     open($fh, "<", $f) or confess "Could not open $f: $!";
340 :     my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
341 :     close($fh);
342 :     return @out;
343 :     }
344 :    
345 :     sub get_experiment_on_off_calls
346 :     {
347 :     my($self, $expt) = @_;
348 :    
349 :     my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
350 :     my $fh;
351 :     open($fh, "<", $f) or confess "Could not open $f: $!";
352 :     my $names = <$fh>;
353 :     chomp $names;
354 :     my @names = split(/\t/, $names);
355 :     my $idx = 0;
356 :     my $expt_idx;
357 :     foreach my $n (@names)
358 :     {
359 :     if ($n eq $expt)
360 :     {
361 :     $expt_idx = $idx;
362 :     last;
363 :     }
364 :     $idx++;
365 :     }
366 :     if (!defined($expt_idx))
367 :     {
368 :     confess("Could not find experiment $expt in $f");
369 :     }
370 :    
371 :     my $calls = {};
372 :     while (<$fh>)
373 :     {
374 :     chomp;
375 :     my($peg, @calls) = split(/\t/);
376 :     #
377 :     # +1 because the row[0] element is the peg, and our index is
378 :     # zero-based.
379 :     #
380 :     $calls->{$peg} = $calls[$expt_idx + 1];
381 :     }
382 :    
383 :     close($fh);
384 :     return($calls);
385 :    
386 :     }
387 :    
388 :     =head3 save_model_gene_activity
389 :    
390 :     $e->save_model_gene_activity($data)
391 :    
392 :     Save the results of a modeling run for a given experiment.
393 :    
394 :     $data is of the form { experiment_id => $data_hash }
395 :    
396 :     =cut
397 :    
398 :     sub save_model_gene_activity
399 :     {
400 :     my($self, $data) = @_;
401 :     }
402 : olson 1.2
403 : olson 1.1 1;
404 : olson 1.2
405 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3