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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (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 : olson 1.4 __PACKAGE__->mk_accessors(qw(genome_dir expr_dir genome_id));
14 : olson 1.1
15 :     =head3 new
16 :    
17 : olson 1.4 my $edir = ExpressionDir->new($expr_dir);
18 :    
19 :     Create a new ExpressionDir object from an existing expression dir.
20 : olson 1.1
21 :    
22 :     =cut
23 :    
24 :     sub new
25 :     {
26 : olson 1.4 my($class, $expr_dir) = @_;
27 :    
28 :     my $gfile = catfile($expr_dir, "GENOME_ID");
29 :     open(GF, "<", $gfile) or die "Cannot open $expr_dir/GENOME_ID: $!";
30 :     my $genome_id = <GF>;
31 :     chomp $genome_id;
32 :     close(GF);
33 :    
34 :     my $self = {
35 :     genome_dir => catfile($expr_dir, $genome_id),
36 :     genome_id => $genome_id,
37 :     expr_dir => $expr_dir,
38 :     };
39 :     return bless $self, $class;
40 :     }
41 :    
42 :     =head3 create
43 : olson 1.1
44 : olson 1.4 my $edir = ExpressionDir->create($expr_dir, $genome_id, $genome_src)
45 :    
46 :     Create a new expression directory from the given genome id and genome data source.
47 :     The data source is either a FIG or FIGV object, or a SAPserver object
48 :     that points at a server from which the data can be extracted.
49 :    
50 :     =cut
51 :    
52 :     sub create
53 :     {
54 :     my($class, $expr_dir, $genome_id, $genome_src) = @_;
55 :    
56 :     if (! -d $expr_dir)
57 : olson 1.1 {
58 : olson 1.4 mkdir($expr_dir);
59 : olson 1.1 }
60 : olson 1.4 my $genome_dir = catfile($expr_dir, $genome_id);
61 :     if (! -d $genome_dir)
62 : olson 1.1 {
63 : olson 1.4 mkdir($genome_dir);
64 : olson 1.1 }
65 : olson 1.4
66 :     open(GF, ">", catfile($expr_dir, "GENOME_ID"));
67 :     print GF "$genome_id\n";
68 :     close(GF);
69 : olson 1.1
70 :     my $self = {
71 :     genome_dir => $genome_dir,
72 : olson 1.4 genome_id => $genome_id,
73 :     expr_dir => $expr_dir,
74 : olson 1.1 };
75 : olson 1.4 bless $self, $class;
76 :    
77 :     if (ref($genome_src) =~ /^FIG/)
78 :     {
79 :     $self->create_from_fig($genome_src);
80 :     }
81 :     elsif (ref($genome_src =~ /^SAP/))
82 :     {
83 :     $self->create_from_sap($genome_src);
84 :     }
85 :     else
86 :     {
87 :     confess "Unknown genome source\n";
88 :     }
89 :     return $self;
90 :     }
91 :    
92 :     sub create_from_fig
93 :     {
94 :     my($self, $fig) = @_;
95 :    
96 :     my $gdir = $fig->organism_directory($self->genome_id);
97 :     copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));
98 :     mkdir(catfile($self->genome_dir, "Features"));
99 :     my @pegs;
100 :     for my $ftype (qw(peg rna))
101 :     {
102 :     my $ofdir = catfile($gdir, "Features", $ftype);
103 :     my $nfdir = catfile($self->genome_dir, "Features", $ftype);
104 :    
105 :     if (open(OT, "<", "$ofdir/tbl"))
106 :     {
107 :     mkdir($nfdir);
108 :     open(NT, ">", "$nfdir/tbl") or confess "Cannot write $nfdir/tbl: $!";
109 :     while (<OT>)
110 :     {
111 :     my($id) = /^(\S+)\t/;
112 :     if (!$fig->is_deleted_fid($id))
113 :     {
114 :     print NT $_;
115 :     push(@pegs, $id) if $ftype eq 'peg';
116 :     }
117 :     }
118 :     close(OT);
119 :     close(NT);
120 :     }
121 :     copy("$ofdir/fasta", "$nfdir/fasta");
122 :     }
123 :    
124 :     open(SS, ">", catfile($self->genome_dir, "subsystem.data"));
125 :     my %phash = $fig->subsystems_for_pegs_complete(\@pegs);
126 :     for my $peg (keys %phash)
127 :     {
128 :     my $list = $phash{$peg};
129 :     next unless $list;
130 :    
131 :     for my $ent (@$list)
132 :     {
133 :     print SS join("\t", $peg, @$ent), "\n";
134 :     }
135 :     }
136 :     close(SS);
137 :    
138 :     open(AF, ">", catfile($self->genome_dir, "assigned_functions"));
139 :     my $fns = $fig->function_of_bulk(\@pegs);
140 :     for my $peg (@pegs)
141 :     {
142 :     print AF join("\t", $peg, $fns->{$peg}), "\n";
143 :     }
144 :     close(AF);
145 :     }
146 :    
147 :     sub create_from_sap
148 :     {
149 :     my($self, $sap) = @_;
150 :     confess "create_from_sap not yet implemented";
151 : olson 1.1 }
152 :    
153 :     sub parse_probe_format_1
154 :     {
155 :     my($self, $in_file, $out_file) = @_;
156 :    
157 :     my($fh);
158 :    
159 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
160 :     my $l = <$fh>;
161 :     chomp $l;
162 :     $l =~ s/\r//;
163 :     my @hdrs = split(/\t/, $l);
164 :     my %hdrs;
165 :     $hdrs{$hdrs[$_]} = $_ for 0..$#hdrs;
166 :    
167 :     my $x_col = $hdrs{"Probe X"};
168 :     my $y_col = $hdrs{"Probe Y"};
169 :     my $seq_col = $hdrs{"Probe Sequence"};
170 :     if (!(defined($x_col) && defined($y_col) && defined($seq_col)))
171 :     {
172 :     close($fh);
173 :     return undef;
174 :     }
175 :    
176 :     my $out;
177 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
178 :    
179 :     while (<$fh>)
180 :     {
181 :     chomp;
182 :     s/\r//g;
183 :     my @flds = split(/\t/,$_);
184 :     my($x,$y,$seq);
185 :     $x = $flds[$x_col];
186 :     $y = $flds[$y_col];
187 :     $seq = $flds[$seq_col];
188 :     my $id = "$x\_$y";
189 :     print $out "$id\t$seq\n";
190 :     }
191 :     close($fh);
192 :     close($out);
193 :     return 1;
194 :     }
195 :    
196 :     sub parse_probe_format_2
197 :     {
198 :     my($self, $in_file, $out_file) = @_;
199 :    
200 :     my($fh);
201 :    
202 :     local $/ = "\n>";
203 :    
204 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
205 :     my $l = <$fh>;
206 :     chomp $l;
207 :     $l =~ s/\r//;
208 :    
209 :     if ($l !~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
210 :     {
211 :     close($fh);
212 :     return undef;
213 :     }
214 :     seek($fh, 0, SEEK_SET);
215 :    
216 :     my $out;
217 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
218 :    
219 :     while (<$fh>)
220 :     {
221 :     chomp;
222 :    
223 :     if ($_ =~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
224 :     {
225 :     if (length($3) < 15)
226 :     {
227 :     close($fh);
228 :     confess "Bad length at line $. of $in_file";
229 :     }
230 :     print $out "$1\_$2\t$3\n";
231 :     }
232 :     else
233 :     {
234 :     confess "Bad input at line $. of $in_file";
235 :     }
236 :     }
237 :     close($out);
238 :     close($fh);
239 :     return 1;
240 :     }
241 :    
242 :     #
243 :     # Our "native" format, used for passing through pre-parsed data.
244 :     #
245 :     sub parse_probe_format_3
246 :     {
247 :     my($self, $in_file, $out_file) = @_;
248 :    
249 :     my($fh);
250 :    
251 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
252 :     my $l = <$fh>;
253 :     chomp $l;
254 :     $l =~ s/\r//;
255 :    
256 :     if ($l !~ /^\d+_\d+\t[ACGT]+$/)
257 :     {
258 :     close($fh);
259 :     return undef;
260 :     }
261 :     seek($fh, 0, SEEK_SET);
262 :    
263 :     my $out;
264 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
265 :    
266 :     while (<$fh>)
267 :     {
268 :     if ($_ =~ /^\d+_\d+\t[ACGT]+$/)
269 :     {
270 :     print $out $_;
271 :     }
272 :     else
273 :     {
274 :     confess "Bad input at line $. of $in_file";
275 :     }
276 :     }
277 :     close($out);
278 :     close($fh);
279 :     return 1;
280 :     }
281 :    
282 :     sub compute_probe_to_peg
283 :     {
284 :     my($self, $probes) = @_;
285 :    
286 :     my $my_probes = catfile($self->expr_dir, "probes.in");
287 :    
288 :     copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
289 :    
290 :     my $probes_fasta = catfile($self->expr_dir, "probes.fasta");
291 :    
292 :     #
293 :     # Attempt to translate probe file.
294 :     #
295 :     my $success;
296 :     for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))
297 :     {
298 :     if ($self->$meth($my_probes, $probes_fasta))
299 :     {
300 :     print STDERR "Translated $probes to $probes_fasta using $meth\n";
301 :     $success = 1;
302 :     last;
303 :     }
304 :     else
305 :     {
306 :     print STDERR "Failed to translate $probes to $probes_fasta using $meth\n";
307 :     }
308 :     }
309 :     if (!$success)
310 :     {
311 :     confess "Could not translate $probes\n";
312 :     }
313 :    
314 :     my $peg_probe_table = catfile($self->expr_dir, 'peg.probe.table');
315 :     my $probe_occ_table = catfile($self->expr_dir, 'probe.occ.table');
316 :    
317 : olson 1.2 my $feature_dir = catfile($self->genome_dir, "Features");
318 :     my @tbls;
319 :     for my $ftype (qw(peg rna))
320 :     {
321 :     my $tfile = catfile($feature_dir, $ftype, 'tbl');
322 :     if (-f $tfile)
323 :     {
324 :     push(@tbls, $tfile);
325 :     }
326 :     }
327 :     if (@tbls == 0)
328 :     {
329 :     confess "Could not find any tbl files in $feature_dir";
330 :     }
331 :    
332 : olson 1.4 $self->run([executable_for("make_probes_to_genes"),
333 :     $probes_fasta,
334 :     catfile($self->genome_dir, 'contigs'),
335 :     $tbls[0],
336 :     $peg_probe_table,
337 :     $probe_occ_table,
338 :     @tbls[1..$#tbls],
339 :     ],
340 :     { stderr => catfile($self->expr_dir, 'problems') });
341 : olson 1.1
342 :    
343 : olson 1.4 $self->run([executable_for("remove_multiple_occurring_probes"),
344 :     $peg_probe_table,
345 :     $probe_occ_table],
346 :     { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
347 : olson 1.1
348 :     $self->make_missing_probes($peg_probe_table, $probes_fasta,
349 :     catfile($self->expr_dir, 'probe.no.match'));
350 :     }
351 : olson 1.2
352 : olson 1.1 sub make_missing_probes
353 :     {
354 :     my($self, $probe_table, $probes, $output) = @_;
355 :     open(MATCH,"<", $probe_table) or die "Cannot open $probe_table: $!";
356 :     open(PROBES,"<", $probes) or die "Cannot open $probes: $!";
357 :     open(OUTPUT, ">", $output) or die "Cannot open $output: $!";
358 :     my %locations;
359 :     while(<MATCH>)
360 :     {
361 :     chomp;
362 :     my($peg,$loc)=split "\t";
363 :     $locations{$loc} = $peg;
364 :     }
365 :    
366 :     while(<PROBES>)
367 :     {
368 :     chomp;
369 :     my($loc,$seq) = split "\t";
370 :     print OUTPUT $loc, "\n" if ! exists $locations{$loc};
371 :     }
372 :     close(MATCH);
373 :     close(PROBES);
374 :     close(OUTPUT);
375 :     }
376 :    
377 : olson 1.2 #
378 :     # we don't copy the experiment files in here because
379 :     # they may be very large. This may change.
380 :     #
381 :     # We do copy the cdf.
382 :     #
383 :     sub compute_rma_normalized
384 :     {
385 :     my($self, $cdf_file, $expt_dir) = @_;
386 :    
387 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
388 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
389 :    
390 :     #
391 :     # We need to build the R library for this cdf.
392 :     #
393 :     my($fh, $tempfile) = tempfile();
394 :     #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
395 :    
396 :     my $cdf_path = $self->expr_dir;
397 :     my $libdir = catfile($self->expr_dir, "r_lib");
398 :     -d $libdir or mkdir $libdir;
399 :     my $pkgdir = catfile($self->expr_dir, "r_pkg");
400 :     -d $pkgdir or mkdir $pkgdir;
401 :    
402 :     print $fh "library(makecdfenv);\n";
403 :     print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
404 :     close($fh);
405 :     system("Rscript", $tempfile);
406 :     system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
407 :    
408 :     local($ENV{R_LIBS}) = $libdir;
409 : olson 1.4 $self->run([executable_for("RunRMA"),
410 :     "data",
411 :     catfile($self->expr_dir, "peg.probe.table"),
412 :     catfile($self->expr_dir, "probe.no.match"),
413 :     $expt_dir,
414 :     $self->expr_dir]);
415 :    
416 : olson 1.2 my $output = catfile($self->expr_dir, "rma_normalized.tab");
417 :     if (! -f $output)
418 :     {
419 :     confess("Output file $output was not generated");
420 :     }
421 :     }
422 :    
423 :     sub compute_atomic_regulons
424 :     {
425 : olson 1.4 my($self, $pearson_cutoff) = @_;
426 :    
427 :     $pearson_cutoff ||= 0.7;
428 : olson 1.2
429 :     my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
430 :     my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
431 :     my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
432 : olson 1.4 my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
433 :     my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
434 :    
435 :    
436 :     $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
437 :     { stdout => $coreg_clusters });
438 :     $self->run([executable_for("make_coreg_conjectures_based_on_subsys"), $self->expr_dir],
439 :     { stdout => $coreg_subsys });
440 :    
441 :     $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
442 :     { stdout => $merged_clusters });
443 :     $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
444 :    
445 :     $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
446 :     { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
447 :    
448 :     $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
449 :     { stdout => catfile($self->expr_dir, "split.clusters") });
450 :    
451 :     }
452 :    
453 :     sub run
454 :     {
455 :     my($self, $cmd, $redirect) = @_;
456 :    
457 :     print "Run @$cmd\n";
458 :     my $rc = system_with_redirect($cmd, $redirect);
459 :     if ($rc != 0)
460 :     {
461 :     confess "Command failed: @$cmd\n";
462 :     }
463 : olson 1.2 }
464 :    
465 : olson 1.3 sub get_experiment_names
466 :     {
467 :     my($self) = @_;
468 :     my $f = catfile($self->expr_dir, "experiment.names");
469 :     my $fh;
470 :     open($fh, "<", $f) or confess "Could not open $f: $!";
471 :     my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
472 :     close($fh);
473 :     return @out;
474 :     }
475 :    
476 :     sub get_experiment_on_off_calls
477 :     {
478 :     my($self, $expt) = @_;
479 :    
480 :     my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
481 :     my $fh;
482 :     open($fh, "<", $f) or confess "Could not open $f: $!";
483 :     my $names = <$fh>;
484 :     chomp $names;
485 :     my @names = split(/\t/, $names);
486 :     my $idx = 0;
487 :     my $expt_idx;
488 :     foreach my $n (@names)
489 :     {
490 :     if ($n eq $expt)
491 :     {
492 :     $expt_idx = $idx;
493 :     last;
494 :     }
495 :     $idx++;
496 :     }
497 :     if (!defined($expt_idx))
498 :     {
499 :     confess("Could not find experiment $expt in $f");
500 :     }
501 :    
502 :     my $calls = {};
503 :     while (<$fh>)
504 :     {
505 :     chomp;
506 :     my($peg, @calls) = split(/\t/);
507 :     #
508 :     # +1 because the row[0] element is the peg, and our index is
509 :     # zero-based.
510 :     #
511 :     $calls->{$peg} = $calls[$expt_idx + 1];
512 :     }
513 :    
514 :     close($fh);
515 :     return($calls);
516 :    
517 :     }
518 :    
519 :     =head3 save_model_gene_activity
520 :    
521 :     $e->save_model_gene_activity($data)
522 :    
523 :     Save the results of a modeling run for a given experiment.
524 :    
525 :     $data is of the form { experiment_id => $data_hash }
526 :    
527 :     =cut
528 :    
529 :     sub save_model_gene_activity
530 :     {
531 :     my($self, $data) = @_;
532 :     }
533 : olson 1.2
534 : olson 1.4 sub all_features
535 :     {
536 :     my($self, $type) = @_;
537 :    
538 :     my @ftypes;
539 :     my $fdir = catfile($self->genome_dir, "Features");
540 :     if (defined($type))
541 :     {
542 :     @ftypes = ($type);
543 :     }
544 :     else
545 :     {
546 :     opendir(D, $fdir);
547 :     @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
548 :     closedir(D);
549 :     }
550 :     my @out;
551 :     for my $ftype (@ftypes)
552 :     {
553 :     if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
554 :     {
555 :     push(@out, map { /^(\S+)/; $1 } <TBL>);
556 :     close(TBL);
557 :     }
558 :     }
559 :     return @out;
560 :     }
561 :    
562 :     sub fid_locations
563 :     {
564 :     my($self, $fids) = @_;
565 :    
566 :     my %fids;
567 :     $fids{$_}++ for @$fids;
568 :    
569 :     my $genome_id = $self->genome_id;
570 :    
571 :     my $fdir = catfile($self->genome_dir, "Features");
572 :     opendir(D, $fdir);
573 :     my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
574 :     closedir(D);
575 :    
576 :     my $out = {};
577 :     for my $ftype (@ftypes)
578 :     {
579 :     if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
580 :     {
581 :     while (<TBL>)
582 :     {
583 :     my($id, $locs) = /^(\S+)\t(\S+)\t/;
584 :    
585 :     if ($fids{$id})
586 :     {
587 :     $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
588 :     }
589 :     }
590 :     close(TBL);
591 :     }
592 :     }
593 :     return $out;
594 :     }
595 :    
596 :     sub ids_in_subsystems
597 :     {
598 :     my($self) = @_;
599 :    
600 :     open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
601 :     my $res = {};
602 :     while (<SS>)
603 :     {
604 :     chomp;
605 :     my($peg, $ss, $role, $variant) = split(/\t/);
606 :     $ss =~ s/\s+/_/g;
607 :     push(@{$res->{$ss}->{$role}}, $peg);
608 :     }
609 :     close(SS);
610 :     return $res;
611 :     }
612 :    
613 :     sub ids_to_subsystems
614 :     {
615 :     my($self, $ids) = @_;
616 :    
617 :     my %ids;
618 :     $ids{$_} = 1 for @$ids;
619 :    
620 :     open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
621 :     my $res = {};
622 :     while (<SS>)
623 :     {
624 :     chomp;
625 :     my($peg, $ss, $role, $variant) = split(/\t/);
626 :     if ($ids{$peg})
627 :     {
628 :     push(@{$res->{$peg}}, $ss);
629 :     }
630 :     }
631 :     close(SS);
632 :     return $res;
633 :     }
634 :    
635 :     sub ids_to_functions
636 :     {
637 :     my($self, $ids) = @_;
638 :     open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
639 :     my %ids;
640 :     $ids{$_} = 1 for @$ids;
641 :     my $res = {};
642 :    
643 :     while (<AF>)
644 :     {
645 :     chomp;
646 :     my($id, $fn) = split(/\t/);
647 :     $res->{$id} = $fn if $ids{$id};
648 :     }
649 :     close(AF);
650 :     return $res;
651 :     }
652 :    
653 : olson 1.1 1;
654 : olson 1.2
655 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3