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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3