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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3