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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3