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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11 - (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 : overbeek 1.11 my $genome_ss_dir = catfile($gdir, "Subsystems");
146 :     if (-d $genome_ss_dir && -s "$genome_ss_dir/bindings" && -s "$genome_ss_dir/subsystems")
147 :     {
148 :     my $my_ss_dir = $self->genome_dir . "/Subsystems";
149 :     mkdir($my_ss_dir) or die "Cannot mkdir $my_ss_dir: $!";
150 :     for my $f (qw(bindings subsystems))
151 :     {
152 :     copy("$genome_ss_dir/$f", "$my_ss_dir/$f") or die "Cannot copy $genome_ss_dir/$f to $my_ss_dir/$f: $!";
153 :     }
154 :     }
155 :    
156 : olson 1.4 open(SS, ">", catfile($self->genome_dir, "subsystem.data"));
157 :     my %phash = $fig->subsystems_for_pegs_complete(\@pegs);
158 :     for my $peg (keys %phash)
159 :     {
160 :     my $list = $phash{$peg};
161 : olson 1.6 next unless $list;
162 : olson 1.4
163 :     for my $ent (@$list)
164 :     {
165 :     print SS join("\t", $peg, @$ent), "\n";
166 :     }
167 :     }
168 :     close(SS);
169 :    
170 :     open(AF, ">", catfile($self->genome_dir, "assigned_functions"));
171 :     my $fns = $fig->function_of_bulk(\@pegs);
172 :     for my $peg (@pegs)
173 :     {
174 :     print AF join("\t", $peg, $fns->{$peg}), "\n";
175 :     }
176 :     close(AF);
177 : olson 1.6
178 :     open(PD, ">", catfile($self->genome_dir, "peg_dna.fasta"));
179 :     for my $peg (@pegs)
180 :     {
181 :     my $dna = $fig->dna_seq($self->genome_id, split(/,/, $locs{$peg}));
182 :     if ($dna eq '')
183 :     {
184 :     die "no dna for ", $self->genome_id, " $peg $locs{$peg}\n";
185 :     }
186 :     print_alignment_as_fasta(\*PD, [$peg, undef, $dna]);
187 :     }
188 :     close(PD);
189 : olson 1.4 }
190 :    
191 :     sub create_from_sap
192 :     {
193 :     my($self, $sap) = @_;
194 :     confess "create_from_sap not yet implemented";
195 : olson 1.1 }
196 :    
197 : olson 1.5 sub parse_probe_format_1lq
198 :     {
199 :     my($self, $in_file, $out_file) = @_;
200 :    
201 :     my($fh);
202 :    
203 :     if ($in_file !~ /\.1lq$/)
204 :     {
205 :     return undef;
206 :     }
207 :    
208 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
209 :    
210 :     my $out;
211 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
212 :    
213 :     # Skip 3 header lines.
214 :     $_ = <$fh> for 1..3;
215 :    
216 :     while (defined($_ = <$fh>))
217 :     {
218 :     if ($_ =~ /(\d+)\s+(\d+)\s+([ACGT]+)\s+(-?\d+)\s/)
219 :     {
220 :     if (length($3) < 15)
221 :     {
222 :     close($fh);
223 :     close($out);
224 :     confess "Bad length at line $. of $in_file";
225 :     return undef;
226 :     }
227 :     next if ($4 =~ /\d+3$/); #mismatch probe
228 :     my($x,$y,$seq) = ($1,$2,$3);
229 :     $seq = scalar reverse $seq;
230 :     print $out "$x\_$y\t$seq\n";
231 :     }
232 :     else
233 :     {
234 :     #
235 :     # We expect some lines not to match.
236 :     #
237 :     }
238 :     }
239 :    
240 :     close($fh);
241 :     close($out);
242 :     return 1;
243 :     }
244 :    
245 : olson 1.1 sub parse_probe_format_1
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 :     my @hdrs = split(/\t/, $l);
256 :     my %hdrs;
257 :     $hdrs{$hdrs[$_]} = $_ for 0..$#hdrs;
258 :    
259 :     my $x_col = $hdrs{"Probe X"};
260 :     my $y_col = $hdrs{"Probe Y"};
261 :     my $seq_col = $hdrs{"Probe Sequence"};
262 :     if (!(defined($x_col) && defined($y_col) && defined($seq_col)))
263 :     {
264 :     close($fh);
265 :     return undef;
266 :     }
267 :    
268 :     my $out;
269 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
270 :    
271 :     while (<$fh>)
272 :     {
273 :     chomp;
274 :     s/\r//g;
275 :     my @flds = split(/\t/,$_);
276 :     my($x,$y,$seq);
277 :     $x = $flds[$x_col];
278 :     $y = $flds[$y_col];
279 :     $seq = $flds[$seq_col];
280 :     my $id = "$x\_$y";
281 :     print $out "$id\t$seq\n";
282 :     }
283 :     close($fh);
284 :     close($out);
285 :     return 1;
286 :     }
287 :    
288 :     sub parse_probe_format_2
289 :     {
290 :     my($self, $in_file, $out_file) = @_;
291 :    
292 :     my($fh);
293 :    
294 :     local $/ = "\n>";
295 :    
296 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
297 :     my $l = <$fh>;
298 :     chomp $l;
299 :     $l =~ s/\r//;
300 :    
301 :     if ($l !~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
302 :     {
303 :     close($fh);
304 :     return undef;
305 :     }
306 :     seek($fh, 0, SEEK_SET);
307 :    
308 :     my $out;
309 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
310 :    
311 :     while (<$fh>)
312 :     {
313 :     chomp;
314 :    
315 :     if ($_ =~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
316 :     {
317 :     if (length($3) < 15)
318 :     {
319 :     close($fh);
320 :     confess "Bad length at line $. of $in_file";
321 :     }
322 :     print $out "$1\_$2\t$3\n";
323 :     }
324 :     else
325 :     {
326 :     confess "Bad input at line $. of $in_file";
327 :     }
328 :     }
329 :     close($out);
330 :     close($fh);
331 :     return 1;
332 :     }
333 :    
334 : olson 1.6 sub parse_probe_format_3
335 :     {
336 :     my($self, $in_file, $out_file) = @_;
337 :    
338 :     my($fh);
339 :    
340 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
341 :     my $out;
342 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
343 :    
344 :     local $/ = "\n>";
345 :    
346 :     while (defined($_ = <$fh>))
347 :     {
348 :     if ($_ =~ /^>?\S+\s+(\d+)\s+(\d+)[^\n]+\n(\S+)/s)
349 :     {
350 :     my $x = $1;
351 :     my $y = $2;
352 :     my $seq = $3;
353 :     if ($seq ne "!")
354 :     {
355 :     if ($seq !~ /^[ACGT]+$/) { die $_ }
356 :     if (length($seq) < 15) { print STDERR "BAD: $_"; die "Failed" }
357 :     print $out "$x\_$y\t$seq\n";
358 :     }
359 :     }
360 :     else
361 :     {
362 :     print STDERR "failed to parse: $_";
363 : olson 1.7 return undef;
364 : olson 1.6 }
365 :     }
366 :    
367 :     close($out);
368 :     close($fh);
369 :     return 1;
370 :     }
371 :    
372 : olson 1.1 #
373 : olson 1.7 # This one showed up in the shewanella data.
374 :     #
375 :     sub parse_probe_format_shew
376 :     {
377 :     my($self, $in_file, $out_file) = @_;
378 :    
379 :     my($fh);
380 :    
381 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
382 :     my $l = <$fh>;
383 :     chomp $l;
384 :     $l =~ s/\r//;
385 :    
386 :     if ($l !~ /x\ty\tprobe_type\tsequence/)
387 :     {
388 :     close($fh);
389 :     return undef;
390 :     }
391 :    
392 :     my $out;
393 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
394 :    
395 :     while (<$fh>)
396 :     {
397 :     if ($_ =~ /^(\d+)\t(\d+)\tPM\t+([ACGT]+)/)
398 :     {
399 :     print $out "$1\_$2\t$3\n";
400 :     }
401 :     }
402 :     close($out);
403 :     close($fh);
404 :     return 1;
405 :     }
406 :     #
407 : olson 1.1 # Our "native" format, used for passing through pre-parsed data.
408 :     #
409 : olson 1.6 sub parse_probe_format_native
410 : olson 1.1 {
411 :     my($self, $in_file, $out_file) = @_;
412 :    
413 :     my($fh);
414 :    
415 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
416 :     my $l = <$fh>;
417 :     chomp $l;
418 :     $l =~ s/\r//;
419 :    
420 :     if ($l !~ /^\d+_\d+\t[ACGT]+$/)
421 :     {
422 :     close($fh);
423 :     return undef;
424 :     }
425 :     seek($fh, 0, SEEK_SET);
426 :    
427 :     my $out;
428 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
429 :    
430 :     while (<$fh>)
431 :     {
432 :     if ($_ =~ /^\d+_\d+\t[ACGT]+$/)
433 :     {
434 :     print $out $_;
435 :     }
436 :     else
437 :     {
438 :     confess "Bad input at line $. of $in_file";
439 :     }
440 :     }
441 :     close($out);
442 :     close($fh);
443 :     return 1;
444 :     }
445 :    
446 :     sub compute_probe_to_peg
447 :     {
448 :     my($self, $probes) = @_;
449 :    
450 : olson 1.7 my($probe_suffix) = $probes =~ m,(\.[^/.]+)$,;
451 : olson 1.5
452 :     my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");
453 : olson 1.1
454 :     copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
455 :    
456 : olson 1.5 my $probes_fasta = catfile($self->expr_dir, "probes");
457 : olson 1.1
458 :     #
459 :     # Attempt to translate probe file.
460 :     #
461 :     my $success;
462 : olson 1.5 for my $meth (@probe_parsers)
463 : olson 1.1 {
464 :     if ($self->$meth($my_probes, $probes_fasta))
465 :     {
466 :     print STDERR "Translated $probes to $probes_fasta using $meth\n";
467 :     $success = 1;
468 :     last;
469 :     }
470 :     else
471 :     {
472 :     print STDERR "Failed to translate $probes to $probes_fasta using $meth\n";
473 :     }
474 :     }
475 :     if (!$success)
476 :     {
477 :     confess "Could not translate $probes\n";
478 :     }
479 :    
480 :     my $peg_probe_table = catfile($self->expr_dir, 'peg.probe.table');
481 :     my $probe_occ_table = catfile($self->expr_dir, 'probe.occ.table');
482 :    
483 : olson 1.2 my $feature_dir = catfile($self->genome_dir, "Features");
484 :     my @tbls;
485 :     for my $ftype (qw(peg rna))
486 :     {
487 :     my $tfile = catfile($feature_dir, $ftype, 'tbl');
488 :     if (-f $tfile)
489 :     {
490 :     push(@tbls, $tfile);
491 :     }
492 :     }
493 :     if (@tbls == 0)
494 :     {
495 :     confess "Could not find any tbl files in $feature_dir";
496 :     }
497 :    
498 : olson 1.4 $self->run([executable_for("make_probes_to_genes"),
499 :     $probes_fasta,
500 :     catfile($self->genome_dir, 'contigs'),
501 :     $tbls[0],
502 :     $peg_probe_table,
503 :     $probe_occ_table,
504 :     @tbls[1..$#tbls],
505 :     ],
506 :     { stderr => catfile($self->expr_dir, 'problems') });
507 : olson 1.1
508 :    
509 : olson 1.4 $self->run([executable_for("remove_multiple_occurring_probes"),
510 :     $peg_probe_table,
511 : olson 1.10 ],
512 : olson 1.4 { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
513 : olson 1.1
514 :     $self->make_missing_probes($peg_probe_table, $probes_fasta,
515 :     catfile($self->expr_dir, 'probe.no.match'));
516 : olson 1.10 $self->make_missing_probes(catfile($self->expr_dir, 'peg.probe.table.no.multiple'), $probes_fasta,
517 :     catfile($self->expr_dir, 'probe.no.multiple.no.match'));
518 : olson 1.1 }
519 : olson 1.2
520 : olson 1.1 sub make_missing_probes
521 :     {
522 :     my($self, $probe_table, $probes, $output) = @_;
523 :     open(MATCH,"<", $probe_table) or die "Cannot open $probe_table: $!";
524 :     open(PROBES,"<", $probes) or die "Cannot open $probes: $!";
525 :     open(OUTPUT, ">", $output) or die "Cannot open $output: $!";
526 :     my %locations;
527 :     while(<MATCH>)
528 :     {
529 :     chomp;
530 :     my($peg,$loc)=split "\t";
531 :     $locations{$loc} = $peg;
532 :     }
533 :    
534 :     while(<PROBES>)
535 :     {
536 :     chomp;
537 :     my($loc,$seq) = split "\t";
538 :     print OUTPUT $loc, "\n" if ! exists $locations{$loc};
539 :     }
540 :     close(MATCH);
541 :     close(PROBES);
542 :     close(OUTPUT);
543 :     }
544 :    
545 : olson 1.2 #
546 :     # we don't copy the experiment files in here because
547 :     # they may be very large. This may change.
548 :     #
549 :     # We do copy the cdf.
550 :     #
551 :     sub compute_rma_normalized
552 :     {
553 :     my($self, $cdf_file, $expt_dir) = @_;
554 :    
555 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
556 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
557 :    
558 :     #
559 :     # We need to build the R library for this cdf.
560 :     #
561 :     my($fh, $tempfile) = tempfile();
562 :     #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
563 :    
564 :     my $cdf_path = $self->expr_dir;
565 :     my $libdir = catfile($self->expr_dir, "r_lib");
566 :     -d $libdir or mkdir $libdir;
567 :     my $pkgdir = catfile($self->expr_dir, "r_pkg");
568 :     -d $pkgdir or mkdir $pkgdir;
569 :    
570 :     print $fh "library(makecdfenv);\n";
571 :     print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
572 :     close($fh);
573 :     system("Rscript", $tempfile);
574 :     system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
575 :    
576 :     local($ENV{R_LIBS}) = $libdir;
577 : olson 1.4 $self->run([executable_for("RunRMA"),
578 :     "data",
579 : olson 1.10 catfile($self->expr_dir, "peg.probe.table.no.multiple"),
580 :     catfile($self->expr_dir, "probe.no.multiple.no.match"),
581 : olson 1.4 $expt_dir,
582 :     $self->expr_dir]);
583 :    
584 : olson 1.2 my $output = catfile($self->expr_dir, "rma_normalized.tab");
585 :     if (! -f $output)
586 :     {
587 :     confess("Output file $output was not generated");
588 :     }
589 :     }
590 :    
591 : olson 1.6 sub compute_rma_normalized_from_sif
592 :     {
593 :     my($self, $cdf_file, $sif_file, $expt_dir) = @_;
594 :    
595 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
596 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
597 :    
598 :     my $my_sif = catfile($self->expr_dir, "expr.sif");
599 :     copy($sif_file, $my_sif) or confess "Cannot copy $sif_file to $my_sif: $!";
600 :    
601 :     #
602 :     # Create the sif2peg mapping.
603 :     #
604 :    
605 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
606 :     if (! -f $sif2peg)
607 :     {
608 :     $self->run([executable_for('map_sif_to_pegs'),
609 :     $my_sif,
610 :     catfile($self->genome_dir, "peg_dna.fasta")],
611 :     { stdout => $sif2peg });
612 :     }
613 :    
614 : olson 1.8 $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
615 : olson 1.6 }
616 :    
617 : olson 1.7 sub compute_rma_normalized_from_locus_tags
618 :     {
619 :     my($self, $cdf_file, $tag_file, $expt_dir) = @_;
620 :    
621 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
622 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
623 :    
624 :     my $my_tags = catfile($self->expr_dir, "expr.tags");
625 :     copy($tag_file, $my_tags) or confess "Cannot copy $tag_file to $my_tags: $!";
626 :    
627 :     #
628 :     # Create the sif2peg mapping.
629 :     #
630 :    
631 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
632 :    
633 :     {
634 :     my %tags;
635 :     #
636 :     # Read the tbl file for the organism and create map from locus tag -> peg.
637 :     #
638 :     my $tbl = catfile($self->genome_dir, "Features", "peg", "tbl");
639 :     if (my $fh = FileHandle->new($tbl, "<"))
640 :     {
641 :     while (<$fh>)
642 :     {
643 :     chomp;
644 :     my($fid, $loc, @aliases) = split(/\t/);
645 :     for my $a (@aliases)
646 :     {
647 :     $tags{$a} = $fid;
648 :     if ($a =~ /^(\S+)\|(.*)/)
649 :     {
650 :     $tags{$2} = $fid;
651 :     }
652 :     }
653 :     }
654 :     close($fh);
655 :     }
656 :     else
657 :     {
658 :     confess "Could not open tbl file $tbl: $!";
659 :     }
660 :    
661 :     my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
662 :     my $in = FileHandle->new($my_tags, "<") or confess "Cannot open $my_tags: $!";
663 :    
664 :     # Per Matt DeJongh: Note that in some cases multiple chip ids
665 :     # map to the same locus tag; in this case ignore the chip id
666 :     # that has "_s_" in it.
667 :    
668 :     my %locus_map;
669 :     while (<$in>)
670 :     {
671 :     chomp;
672 :     my($chip_id, $tag) = split(/\t/);
673 :     next if $tag eq '';
674 :    
675 :     if (exists($locus_map{$tag}))
676 :     {
677 :     print "Dup: chip_id=$chip_id tag=$tag old=$locus_map{$tag}\n";
678 : olson 1.8 next if $chip_id =~ /_s_/ ;
679 :     next if $chip_id =~ /^Rick/ && $self->genome_id eq '452659.3';
680 : olson 1.7 }
681 :     $locus_map{$tag} = $chip_id;
682 :     }
683 :     close($in);
684 :    
685 :     for my $locus_tag (sort keys %locus_map)
686 :     {
687 :     my $chip_id = $locus_map{$locus_tag};
688 :     my $peg = $tags{$locus_tag};
689 :     if ($peg ne '')
690 :     {
691 :     print $out "$chip_id\t$peg\n";
692 :     }
693 :     }
694 :     close($out);
695 :     }
696 :    
697 :     if (! -s $sif2peg)
698 :     {
699 :     confess "No probe to peg mappings were found\n";
700 :     }
701 : olson 1.8
702 :     $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
703 :     }
704 :    
705 :     sub compute_rma_normalized_from_pegidcorr
706 :     {
707 :     my($self, $cdf_file, $corr_file, $expt_dir) = @_;
708 :    
709 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
710 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
711 :    
712 :     my $my_corr = catfile($self->expr_dir, "peg.id.corr");
713 :     copy($corr_file, $my_corr) or confess "Cannot copy $corr_file to $my_corr: $!";
714 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
715 :     #
716 :     # Create the sif2peg mapping.
717 :     #
718 :    
719 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
720 :    
721 :     #
722 :     # The peg.id.corr table is of the form peg \t chip-id \t something-else
723 :     # Just rewrite into chip-id \t peg.
724 :     #
725 :    
726 :     my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
727 :     my $in = FileHandle->new($my_corr, "<") or confess "Cannot open $my_corr: $!";
728 :     while (<$in>)
729 :     {
730 :     chomp;
731 :     my($peg, $chip_id, undef) = split(/\t/);
732 :    
733 :     print $out "$chip_id\t$peg\n";
734 :     }
735 :     close($in);
736 :     close($out);
737 :    
738 :     if (! -s $sif2peg)
739 :     {
740 :     confess "No probe to peg mappings were found\n";
741 :     }
742 :     $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
743 :     }
744 :    
745 :     sub compute_rma_normalized_using_sif2peg
746 :     {
747 :     my($self, $sif2peg, $expt_dir) = @_;
748 : olson 1.7
749 :     #
750 :     # We need to build the R library for this cdf.
751 :     #
752 :     my($fh, $tempfile) = tempfile();
753 :     #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
754 :    
755 :     my $cdf_path = $self->expr_dir;
756 :     my $libdir = catfile($self->expr_dir, "r_lib");
757 :     -d $libdir or mkdir $libdir;
758 :     my $pkgdir = catfile($self->expr_dir, "r_pkg");
759 :     -d $pkgdir or mkdir $pkgdir;
760 :    
761 :     print $fh "library(makecdfenv);\n";
762 :     print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
763 :     close($fh);
764 :     system("Rscript", $tempfile);
765 :     system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
766 :    
767 :     local($ENV{R_LIBS}) = $libdir;
768 :     $self->run([executable_for("RunRMA_SIF_format"),
769 :     "data",
770 :     $expt_dir,
771 :     $self->expr_dir]);
772 :    
773 :     my $output = catfile($self->expr_dir, "rma_normalized.tab");
774 :     if (! -f $output)
775 :     {
776 :     confess("Output file $output was not generated");
777 :     }
778 :     }
779 :    
780 : olson 1.2 sub compute_atomic_regulons
781 :     {
782 : olson 1.4 my($self, $pearson_cutoff) = @_;
783 :    
784 :     $pearson_cutoff ||= 0.7;
785 : olson 1.2
786 :     my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
787 :     my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
788 :     my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
789 : olson 1.4 my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
790 :     my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
791 :    
792 : overbeek 1.11
793 : olson 1.4 $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
794 :     { stdout => $coreg_clusters });
795 : overbeek 1.11
796 :     my $genome_ss_dir = $self->genome_dir . "/Subsystems";
797 :     $self->run([executable_for("make_coreg_conjectures_based_on_subsys"),
798 :     $self->expr_dir,
799 :     (-d $genome_ss_dir ? $genome_ss_dir : ()),
800 :     ],
801 : olson 1.4 { stdout => $coreg_subsys });
802 :    
803 :     $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
804 :     { stdout => $merged_clusters });
805 :     $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
806 :    
807 : olson 1.10 if (-s $pegs_always_on == 0)
808 : olson 1.5 {
809 : olson 1.10 confess "No always-on pegs were found";
810 : olson 1.5 }
811 :    
812 : olson 1.4 $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
813 :     { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
814 :    
815 :     $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
816 :     { stdout => catfile($self->expr_dir, "split.clusters") });
817 :    
818 : olson 1.5 $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
819 : olson 1.4 }
820 :    
821 :     sub run
822 :     {
823 :     my($self, $cmd, $redirect) = @_;
824 :    
825 :     print "Run @$cmd\n";
826 :     my $rc = system_with_redirect($cmd, $redirect);
827 :     if ($rc != 0)
828 :     {
829 :     confess "Command failed: @$cmd\n";
830 :     }
831 : olson 1.2 }
832 :    
833 : olson 1.3 sub get_experiment_names
834 :     {
835 :     my($self) = @_;
836 :     my $f = catfile($self->expr_dir, "experiment.names");
837 :     my $fh;
838 :     open($fh, "<", $f) or confess "Could not open $f: $!";
839 :     my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
840 :     close($fh);
841 :     return @out;
842 :     }
843 :    
844 :     sub get_experiment_on_off_calls
845 :     {
846 :     my($self, $expt) = @_;
847 :    
848 :     my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
849 :     my $fh;
850 :     open($fh, "<", $f) or confess "Could not open $f: $!";
851 :     my $names = <$fh>;
852 :     chomp $names;
853 :     my @names = split(/\t/, $names);
854 :     my $idx = 0;
855 :     my $expt_idx;
856 :     foreach my $n (@names)
857 :     {
858 :     if ($n eq $expt)
859 :     {
860 :     $expt_idx = $idx;
861 :     last;
862 :     }
863 :     $idx++;
864 :     }
865 :     if (!defined($expt_idx))
866 :     {
867 :     confess("Could not find experiment $expt in $f");
868 :     }
869 :    
870 :     my $calls = {};
871 :     while (<$fh>)
872 :     {
873 :     chomp;
874 :     my($peg, @calls) = split(/\t/);
875 :     #
876 :     # +1 because the row[0] element is the peg, and our index is
877 :     # zero-based.
878 :     #
879 :     $calls->{$peg} = $calls[$expt_idx + 1];
880 :     }
881 :    
882 :     close($fh);
883 :     return($calls);
884 :    
885 :     }
886 :    
887 :     =head3 save_model_gene_activity
888 :    
889 :     $e->save_model_gene_activity($data)
890 :    
891 :     Save the results of a modeling run for a given experiment.
892 :    
893 :     $data is of the form { experiment_id => $data_hash }
894 :    
895 :     =cut
896 :    
897 :     sub save_model_gene_activity
898 :     {
899 :     my($self, $data) = @_;
900 :     }
901 : olson 1.2
902 : olson 1.4 sub all_features
903 :     {
904 :     my($self, $type) = @_;
905 :    
906 :     my @ftypes;
907 :     my $fdir = catfile($self->genome_dir, "Features");
908 :     if (defined($type))
909 :     {
910 :     @ftypes = ($type);
911 :     }
912 :     else
913 :     {
914 :     opendir(D, $fdir);
915 :     @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
916 :     closedir(D);
917 :     }
918 :     my @out;
919 :     for my $ftype (@ftypes)
920 :     {
921 :     if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
922 :     {
923 :     push(@out, map { /^(\S+)/; $1 } <TBL>);
924 :     close(TBL);
925 :     }
926 :     }
927 :     return @out;
928 :     }
929 :    
930 :     sub fid_locations
931 :     {
932 :     my($self, $fids) = @_;
933 :    
934 :     my %fids;
935 :     $fids{$_}++ for @$fids;
936 :    
937 :     my $genome_id = $self->genome_id;
938 :    
939 :     my $fdir = catfile($self->genome_dir, "Features");
940 :     opendir(D, $fdir);
941 :     my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
942 :     closedir(D);
943 :    
944 :     my $out = {};
945 :     for my $ftype (@ftypes)
946 :     {
947 :     if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
948 :     {
949 :     while (<TBL>)
950 :     {
951 :     my($id, $locs) = /^(\S+)\t(\S+)\t/;
952 :    
953 :     if ($fids{$id})
954 :     {
955 :     $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
956 :     }
957 :     }
958 :     close(TBL);
959 :     }
960 :     }
961 :     return $out;
962 :     }
963 :    
964 :     sub ids_in_subsystems
965 :     {
966 :     my($self) = @_;
967 :    
968 : olson 1.10 my $dir = $self->genome_dir;
969 :     my $fh;
970 :     if (!open($fh, "<", "$dir/Subsystems/bindings"))
971 :     {
972 :     warn "No bindings file, falling back to old method\n";
973 :     return $self->ids_in_subsystems_old();
974 :     }
975 :    
976 :     my $res;
977 :     while (<$fh>)
978 :     {
979 :     chomp;
980 :     my($ss, $role, $fid) = split(/\t/);
981 :     $ss =~ s/\s+/_/g;
982 :     push(@{$res->{$ss}->{$role}}, $fid);
983 :     }
984 :     close($fh);
985 :     return $res;
986 :     }
987 :    
988 :     sub ids_to_subsystems
989 :     {
990 :     my($self, $ids) = @_;
991 :    
992 :     my $dir = $self->genome_dir;
993 :     my $fh;
994 :     if (!open($fh, "<", "$dir/Subsystems/bindings"))
995 :     {
996 :     warn "No bindings file, falling back to old method\n";
997 :     return $self->ids_to_subsystems_old($ids);
998 :     }
999 :    
1000 :     my %ids;
1001 :     $ids{$_} = 1 for @$ids;
1002 :    
1003 :     my $res = {};
1004 :     while (<$fh>)
1005 :     {
1006 :     chomp;
1007 :     my($ss, $role, $fid) = split(/\t/);
1008 :     if ($ids{$fid})
1009 :     {
1010 :     push(@{$res->{$fid}}, $ss);
1011 :     }
1012 :     }
1013 :     close(SS);
1014 :    
1015 :     return $res;
1016 :     }
1017 :    
1018 :     sub ids_in_subsystems_old
1019 :     {
1020 :     my($self) = @_;
1021 :    
1022 : olson 1.4 open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
1023 :     my $res = {};
1024 :     while (<SS>)
1025 :     {
1026 :     chomp;
1027 :     my($peg, $ss, $role, $variant) = split(/\t/);
1028 :     $ss =~ s/\s+/_/g;
1029 :     push(@{$res->{$ss}->{$role}}, $peg);
1030 :     }
1031 :     close(SS);
1032 :     return $res;
1033 :     }
1034 :    
1035 : olson 1.10 sub ids_to_subsystems_old
1036 : olson 1.4 {
1037 :     my($self, $ids) = @_;
1038 :    
1039 :     my %ids;
1040 :     $ids{$_} = 1 for @$ids;
1041 :    
1042 :     open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
1043 :     my $res = {};
1044 :     while (<SS>)
1045 :     {
1046 :     chomp;
1047 :     my($peg, $ss, $role, $variant) = split(/\t/);
1048 :     if ($ids{$peg})
1049 :     {
1050 :     push(@{$res->{$peg}}, $ss);
1051 :     }
1052 :     }
1053 :     close(SS);
1054 :     return $res;
1055 :     }
1056 :    
1057 :     sub ids_to_functions
1058 :     {
1059 :     my($self, $ids) = @_;
1060 :     open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
1061 :     my %ids;
1062 :     $ids{$_} = 1 for @$ids;
1063 :     my $res = {};
1064 :    
1065 :     while (<AF>)
1066 :     {
1067 :     chomp;
1068 :     my($id, $fn) = split(/\t/);
1069 :     $res->{$id} = $fn if $ids{$id};
1070 :     }
1071 :     close(AF);
1072 :     return $res;
1073 :     }
1074 :    
1075 : overbeek 1.9 sub best_pearson_corr {
1076 :     my($self,$pegs1,$cutoff) = @_;
1077 :    
1078 :     my @pegs2 = $self->all_features('peg');
1079 :     my $handle = $self->get_pc_hash_strip($pegs1,\@pegs2);
1080 :    
1081 :     my %ok;
1082 :     my $i;
1083 :     for ($i=0; ($i < @$pegs1); $i++)
1084 :     {
1085 :     foreach my $peg2 ( @pegs2 )
1086 :     {
1087 :     my $pc = &pearson_corr($handle,$pegs1->[$i],$peg2);
1088 :     if (abs($pc >= $cutoff))
1089 :     {
1090 :     $ok{$pegs1->[$i]} -> {$peg2} = $pc;
1091 :     }
1092 :     }
1093 :     }
1094 :     return \%ok;
1095 :     }
1096 :    
1097 :     sub pearson_corr {
1098 :     my($hash,$peg1,$peg2) = @_;
1099 :     my $v = $hash->{$peg1}->{$peg2};
1100 :     return defined($v) ? sprintf("%0.3f",$v) : " ";
1101 :     }
1102 :    
1103 :     sub get_pc_hash_strip {
1104 :     my($self,$pegs1,$pegs2) = @_;
1105 :     my $corrH = $self->get_corr;
1106 :     my $hash = &compute_pc_strip($pegs1,$pegs2,$corrH);
1107 :     return $hash;
1108 :     }
1109 :    
1110 :     sub get_corr {
1111 :     my($self) = @_;
1112 :    
1113 :     my $dir = $self->expr_dir;
1114 :     my $rawF = "$dir/rma_normalized.tab";
1115 :     my %gene_to_values;
1116 :     open(RAW,"<$rawF") || die "could not open $rawF";
1117 :     while (<RAW>)
1118 :     {
1119 :     chomp;
1120 :     my ($gene_id, @gxp_values) = split("\t");
1121 :     $gene_to_values{$gene_id} = \@gxp_values;
1122 :     }
1123 :     close(RAW);
1124 :     return \%gene_to_values;
1125 :     }
1126 :    
1127 :     sub compute_pc_strip {
1128 :     my ($pegs1,$pegs2, $gxp_hash) = @_;
1129 :     my %values = ();
1130 :    
1131 :     for (my $i = 0; $i < @$pegs1; $i++)
1132 :     {
1133 :     my $stat = Statistics::Descriptive::Full->new();
1134 :     $stat->add_data(@{$gxp_hash->{$pegs1->[$i]}});
1135 :    
1136 :     foreach my $peg2 (@$pegs2)
1137 :     {
1138 :     if ($pegs1->[$i] ne $peg2)
1139 :     {
1140 :     my ($q, $m, $r, $err) = $stat->least_squares_fit(@{$gxp_hash->{$peg2}});
1141 :     $values{$pegs1->[$i]}->{$peg2} = $r;
1142 :     }
1143 :     }
1144 :     }
1145 :    
1146 :     return \%values;
1147 :     }
1148 :    
1149 :    
1150 : olson 1.1 1;
1151 : olson 1.2

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3