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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3