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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.9 - (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 : olson 1.8 $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
602 : olson 1.6 }
603 :    
604 : olson 1.7 sub compute_rma_normalized_from_locus_tags
605 :     {
606 :     my($self, $cdf_file, $tag_file, $expt_dir) = @_;
607 :    
608 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
609 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
610 :    
611 :     my $my_tags = catfile($self->expr_dir, "expr.tags");
612 :     copy($tag_file, $my_tags) or confess "Cannot copy $tag_file to $my_tags: $!";
613 :    
614 :     #
615 :     # Create the sif2peg mapping.
616 :     #
617 :    
618 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
619 :    
620 :     {
621 :     my %tags;
622 :     #
623 :     # Read the tbl file for the organism and create map from locus tag -> peg.
624 :     #
625 :     my $tbl = catfile($self->genome_dir, "Features", "peg", "tbl");
626 :     if (my $fh = FileHandle->new($tbl, "<"))
627 :     {
628 :     while (<$fh>)
629 :     {
630 :     chomp;
631 :     my($fid, $loc, @aliases) = split(/\t/);
632 :     for my $a (@aliases)
633 :     {
634 :     $tags{$a} = $fid;
635 :     if ($a =~ /^(\S+)\|(.*)/)
636 :     {
637 :     $tags{$2} = $fid;
638 :     }
639 :     }
640 :     }
641 :     close($fh);
642 :     }
643 :     else
644 :     {
645 :     confess "Could not open tbl file $tbl: $!";
646 :     }
647 :    
648 :     my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
649 :     my $in = FileHandle->new($my_tags, "<") or confess "Cannot open $my_tags: $!";
650 :    
651 :     # Per Matt DeJongh: Note that in some cases multiple chip ids
652 :     # map to the same locus tag; in this case ignore the chip id
653 :     # that has "_s_" in it.
654 :    
655 :     my %locus_map;
656 :     while (<$in>)
657 :     {
658 :     chomp;
659 :     my($chip_id, $tag) = split(/\t/);
660 :     next if $tag eq '';
661 :    
662 :     if (exists($locus_map{$tag}))
663 :     {
664 :     print "Dup: chip_id=$chip_id tag=$tag old=$locus_map{$tag}\n";
665 : olson 1.8 next if $chip_id =~ /_s_/ ;
666 :     next if $chip_id =~ /^Rick/ && $self->genome_id eq '452659.3';
667 : olson 1.7 }
668 :     $locus_map{$tag} = $chip_id;
669 :     }
670 :     close($in);
671 :    
672 :     for my $locus_tag (sort keys %locus_map)
673 :     {
674 :     my $chip_id = $locus_map{$locus_tag};
675 :     my $peg = $tags{$locus_tag};
676 :     if ($peg ne '')
677 :     {
678 :     print $out "$chip_id\t$peg\n";
679 :     }
680 :     }
681 :     close($out);
682 :     }
683 :    
684 :     if (! -s $sif2peg)
685 :     {
686 :     confess "No probe to peg mappings were found\n";
687 :     }
688 : olson 1.8
689 :     $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
690 :     }
691 :    
692 :     sub compute_rma_normalized_from_pegidcorr
693 :     {
694 :     my($self, $cdf_file, $corr_file, $expt_dir) = @_;
695 :    
696 :     my $my_cdf = catfile($self->expr_dir, "expr.cdf");
697 :     copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
698 :    
699 :     my $my_corr = catfile($self->expr_dir, "peg.id.corr");
700 :     copy($corr_file, $my_corr) or confess "Cannot copy $corr_file to $my_corr: $!";
701 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
702 :     #
703 :     # Create the sif2peg mapping.
704 :     #
705 :    
706 :     my $sif2peg = catfile($self->expr_dir, "sif2peg");
707 :    
708 :     #
709 :     # The peg.id.corr table is of the form peg \t chip-id \t something-else
710 :     # Just rewrite into chip-id \t peg.
711 :     #
712 :    
713 :     my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
714 :     my $in = FileHandle->new($my_corr, "<") or confess "Cannot open $my_corr: $!";
715 :     while (<$in>)
716 :     {
717 :     chomp;
718 :     my($peg, $chip_id, undef) = split(/\t/);
719 :    
720 :     print $out "$chip_id\t$peg\n";
721 :     }
722 :     close($in);
723 :     close($out);
724 :    
725 :     if (! -s $sif2peg)
726 :     {
727 :     confess "No probe to peg mappings were found\n";
728 :     }
729 :     $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
730 :     }
731 :    
732 :     sub compute_rma_normalized_using_sif2peg
733 :     {
734 :     my($self, $sif2peg, $expt_dir) = @_;
735 : olson 1.7
736 :     #
737 :     # We need to build the R library for this cdf.
738 :     #
739 :     my($fh, $tempfile) = tempfile();
740 :     #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
741 :    
742 :     my $cdf_path = $self->expr_dir;
743 :     my $libdir = catfile($self->expr_dir, "r_lib");
744 :     -d $libdir or mkdir $libdir;
745 :     my $pkgdir = catfile($self->expr_dir, "r_pkg");
746 :     -d $pkgdir or mkdir $pkgdir;
747 :    
748 :     print $fh "library(makecdfenv);\n";
749 :     print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
750 :     close($fh);
751 :     system("Rscript", $tempfile);
752 :     system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
753 :    
754 :     local($ENV{R_LIBS}) = $libdir;
755 :     $self->run([executable_for("RunRMA_SIF_format"),
756 :     "data",
757 :     $expt_dir,
758 :     $self->expr_dir]);
759 :    
760 :     my $output = catfile($self->expr_dir, "rma_normalized.tab");
761 :     if (! -f $output)
762 :     {
763 :     confess("Output file $output was not generated");
764 :     }
765 :     }
766 :    
767 : olson 1.2 sub compute_atomic_regulons
768 :     {
769 : olson 1.4 my($self, $pearson_cutoff) = @_;
770 :    
771 :     $pearson_cutoff ||= 0.7;
772 : olson 1.2
773 :     my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
774 :     my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
775 :     my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
776 : olson 1.4 my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
777 :     my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
778 :    
779 :     $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
780 :     { stdout => $coreg_clusters });
781 :     $self->run([executable_for("make_coreg_conjectures_based_on_subsys"), $self->expr_dir],
782 :     { stdout => $coreg_subsys });
783 :    
784 :     $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
785 :     { stdout => $merged_clusters });
786 :     $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
787 :    
788 : olson 1.5 if (-s $probes_always_on == 0)
789 :     {
790 :     confess "No always-on probes were found";
791 :     }
792 :    
793 : olson 1.4 $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
794 :     { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
795 :    
796 :     $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
797 :     { stdout => catfile($self->expr_dir, "split.clusters") });
798 :    
799 : olson 1.5 $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
800 : olson 1.4 }
801 :    
802 :     sub run
803 :     {
804 :     my($self, $cmd, $redirect) = @_;
805 :    
806 :     print "Run @$cmd\n";
807 :     my $rc = system_with_redirect($cmd, $redirect);
808 :     if ($rc != 0)
809 :     {
810 :     confess "Command failed: @$cmd\n";
811 :     }
812 : olson 1.2 }
813 :    
814 : olson 1.3 sub get_experiment_names
815 :     {
816 :     my($self) = @_;
817 :     my $f = catfile($self->expr_dir, "experiment.names");
818 :     my $fh;
819 :     open($fh, "<", $f) or confess "Could not open $f: $!";
820 :     my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
821 :     close($fh);
822 :     return @out;
823 :     }
824 :    
825 :     sub get_experiment_on_off_calls
826 :     {
827 :     my($self, $expt) = @_;
828 :    
829 :     my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
830 :     my $fh;
831 :     open($fh, "<", $f) or confess "Could not open $f: $!";
832 :     my $names = <$fh>;
833 :     chomp $names;
834 :     my @names = split(/\t/, $names);
835 :     my $idx = 0;
836 :     my $expt_idx;
837 :     foreach my $n (@names)
838 :     {
839 :     if ($n eq $expt)
840 :     {
841 :     $expt_idx = $idx;
842 :     last;
843 :     }
844 :     $idx++;
845 :     }
846 :     if (!defined($expt_idx))
847 :     {
848 :     confess("Could not find experiment $expt in $f");
849 :     }
850 :    
851 :     my $calls = {};
852 :     while (<$fh>)
853 :     {
854 :     chomp;
855 :     my($peg, @calls) = split(/\t/);
856 :     #
857 :     # +1 because the row[0] element is the peg, and our index is
858 :     # zero-based.
859 :     #
860 :     $calls->{$peg} = $calls[$expt_idx + 1];
861 :     }
862 :    
863 :     close($fh);
864 :     return($calls);
865 :    
866 :     }
867 :    
868 :     =head3 save_model_gene_activity
869 :    
870 :     $e->save_model_gene_activity($data)
871 :    
872 :     Save the results of a modeling run for a given experiment.
873 :    
874 :     $data is of the form { experiment_id => $data_hash }
875 :    
876 :     =cut
877 :    
878 :     sub save_model_gene_activity
879 :     {
880 :     my($self, $data) = @_;
881 :     }
882 : olson 1.2
883 : olson 1.4 sub all_features
884 :     {
885 :     my($self, $type) = @_;
886 :    
887 :     my @ftypes;
888 :     my $fdir = catfile($self->genome_dir, "Features");
889 :     if (defined($type))
890 :     {
891 :     @ftypes = ($type);
892 :     }
893 :     else
894 :     {
895 :     opendir(D, $fdir);
896 :     @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
897 :     closedir(D);
898 :     }
899 :     my @out;
900 :     for my $ftype (@ftypes)
901 :     {
902 :     if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
903 :     {
904 :     push(@out, map { /^(\S+)/; $1 } <TBL>);
905 :     close(TBL);
906 :     }
907 :     }
908 :     return @out;
909 :     }
910 :    
911 :     sub fid_locations
912 :     {
913 :     my($self, $fids) = @_;
914 :    
915 :     my %fids;
916 :     $fids{$_}++ for @$fids;
917 :    
918 :     my $genome_id = $self->genome_id;
919 :    
920 :     my $fdir = catfile($self->genome_dir, "Features");
921 :     opendir(D, $fdir);
922 :     my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
923 :     closedir(D);
924 :    
925 :     my $out = {};
926 :     for my $ftype (@ftypes)
927 :     {
928 :     if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
929 :     {
930 :     while (<TBL>)
931 :     {
932 :     my($id, $locs) = /^(\S+)\t(\S+)\t/;
933 :    
934 :     if ($fids{$id})
935 :     {
936 :     $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
937 :     }
938 :     }
939 :     close(TBL);
940 :     }
941 :     }
942 :     return $out;
943 :     }
944 :    
945 :     sub ids_in_subsystems
946 :     {
947 :     my($self) = @_;
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 :     $ss =~ s/\s+/_/g;
956 :     push(@{$res->{$ss}->{$role}}, $peg);
957 :     }
958 :     close(SS);
959 :     return $res;
960 :     }
961 :    
962 :     sub ids_to_subsystems
963 :     {
964 :     my($self, $ids) = @_;
965 :    
966 :     my %ids;
967 :     $ids{$_} = 1 for @$ids;
968 :    
969 :     open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
970 :     my $res = {};
971 :     while (<SS>)
972 :     {
973 :     chomp;
974 :     my($peg, $ss, $role, $variant) = split(/\t/);
975 :     if ($ids{$peg})
976 :     {
977 :     push(@{$res->{$peg}}, $ss);
978 :     }
979 :     }
980 :     close(SS);
981 :     return $res;
982 :     }
983 :    
984 :     sub ids_to_functions
985 :     {
986 :     my($self, $ids) = @_;
987 :     open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
988 :     my %ids;
989 :     $ids{$_} = 1 for @$ids;
990 :     my $res = {};
991 :    
992 :     while (<AF>)
993 :     {
994 :     chomp;
995 :     my($id, $fn) = split(/\t/);
996 :     $res->{$id} = $fn if $ids{$id};
997 :     }
998 :     close(AF);
999 :     return $res;
1000 :     }
1001 :    
1002 : overbeek 1.9 sub best_pearson_corr {
1003 :     my($self,$pegs1,$cutoff) = @_;
1004 :    
1005 :     my @pegs2 = $self->all_features('peg');
1006 :     my $handle = $self->get_pc_hash_strip($pegs1,\@pegs2);
1007 :    
1008 :     my %ok;
1009 :     my $i;
1010 :     for ($i=0; ($i < @$pegs1); $i++)
1011 :     {
1012 :     foreach my $peg2 ( @pegs2 )
1013 :     {
1014 :     my $pc = &pearson_corr($handle,$pegs1->[$i],$peg2);
1015 :     if (abs($pc >= $cutoff))
1016 :     {
1017 :     $ok{$pegs1->[$i]} -> {$peg2} = $pc;
1018 :     }
1019 :     }
1020 :     }
1021 :     return \%ok;
1022 :     }
1023 :    
1024 :     sub pearson_corr {
1025 :     my($hash,$peg1,$peg2) = @_;
1026 :     my $v = $hash->{$peg1}->{$peg2};
1027 :     return defined($v) ? sprintf("%0.3f",$v) : " ";
1028 :     }
1029 :    
1030 :     sub get_pc_hash_strip {
1031 :     my($self,$pegs1,$pegs2) = @_;
1032 :     my $corrH = $self->get_corr;
1033 :     my $hash = &compute_pc_strip($pegs1,$pegs2,$corrH);
1034 :     return $hash;
1035 :     }
1036 :    
1037 :     sub get_corr {
1038 :     my($self) = @_;
1039 :    
1040 :     my $dir = $self->expr_dir;
1041 :     my $rawF = "$dir/rma_normalized.tab";
1042 :     my %gene_to_values;
1043 :     open(RAW,"<$rawF") || die "could not open $rawF";
1044 :     while (<RAW>)
1045 :     {
1046 :     chomp;
1047 :     my ($gene_id, @gxp_values) = split("\t");
1048 :     $gene_to_values{$gene_id} = \@gxp_values;
1049 :     }
1050 :     close(RAW);
1051 :     return \%gene_to_values;
1052 :     }
1053 :    
1054 :     sub compute_pc_strip {
1055 :     my ($pegs1,$pegs2, $gxp_hash) = @_;
1056 :     my %values = ();
1057 :    
1058 :     for (my $i = 0; $i < @$pegs1; $i++)
1059 :     {
1060 :     my $stat = Statistics::Descriptive::Full->new();
1061 :     $stat->add_data(@{$gxp_hash->{$pegs1->[$i]}});
1062 :    
1063 :     foreach my $peg2 (@$pegs2)
1064 :     {
1065 :     if ($pegs1->[$i] ne $peg2)
1066 :     {
1067 :     my ($q, $m, $r, $err) = $stat->least_squares_fit(@{$gxp_hash->{$peg2}});
1068 :     $values{$pegs1->[$i]}->{$peg2} = $r;
1069 :     }
1070 :     }
1071 :     }
1072 :    
1073 :     return \%values;
1074 :     }
1075 :    
1076 :    
1077 : olson 1.1 1;
1078 : olson 1.2

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3