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

Diff of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2, Thu Jan 6 22:57:32 2011 UTC revision 1.4, Fri Jan 7 22:20:55 2011 UTC
# Line 10  Line 10 
10  use Carp;  use Carp;
11  use Fcntl ':seek';  use Fcntl ':seek';
12    
13  __PACKAGE__->mk_accessors(qw(genome_dir expr_dir error_dir));  __PACKAGE__->mk_accessors(qw(genome_dir expr_dir genome_id));
14    
15  =head3 new  =head3 new
16    
17      my $edir = ExpressionDir->new($genome_dir);      my $edir = ExpressionDir->new($expr_dir);
18    
19    Create a new ExpressionDir object from an existing expression dir.
20    
 Create a new ExpressionDir object to be associated with the given genome directory.  
 If a subdirectory ExpressionData does not yet exist, one is created.  
21    
22  =cut  =cut
23    
24  sub new  sub new
25  {  {
26      my($class, $genome_dir) = @_;      my($class, $expr_dir) = @_;
27    
28        my $gfile = catfile($expr_dir, "GENOME_ID");
29        open(GF, "<", $gfile) or die "Cannot open $expr_dir/GENOME_ID: $!";
30        my $genome_id = <GF>;
31        chomp $genome_id;
32        close(GF);
33    
34        my $self = {
35            genome_dir => catfile($expr_dir, $genome_id),
36            genome_id => $genome_id,
37            expr_dir => $expr_dir,
38        };
39        return bless $self, $class;
40    }
41    
42    =head3 create
43    
44        my $edir = ExpressionDir->create($expr_dir, $genome_id, $genome_src)
45    
46    Create a new expression directory from the given genome id and genome data source.
47    The data source is either a FIG or FIGV object, or a SAPserver object
48    that points at a server from which the data can be extracted.
49    
50    =cut
51    
52    sub create
53    {
54        my($class, $expr_dir, $genome_id, $genome_src) = @_;
55    
56      my $edir = catfile($genome_dir, "ExpressionData");      if (! -d $expr_dir)
     if (! -d $edir)  
57      {      {
58          mkdir($edir);          mkdir($expr_dir);
59      }      }
60      my $errdir = catfile($genome_dir, 'errors');      my $genome_dir = catfile($expr_dir, $genome_id);
61      if (! -d $errdir)      if (! -d $genome_dir)
62      {      {
63          mkdir($errdir);          mkdir($genome_dir);
64      }      }
65    
66        open(GF, ">", catfile($expr_dir, "GENOME_ID"));
67        print GF "$genome_id\n";
68        close(GF);
69    
70      my $self = {      my $self = {
71          genome_dir => $genome_dir,          genome_dir => $genome_dir,
72          expr_dir => $edir,          genome_id => $genome_id,
73          error_dir => $errdir,          expr_dir => $expr_dir,
74      };      };
75      return bless $self, $class;      bless $self, $class;
76    
77        if (ref($genome_src) =~ /^FIG/)
78        {
79            $self->create_from_fig($genome_src);
80        }
81        elsif (ref($genome_src =~ /^SAP/))
82        {
83            $self->create_from_sap($genome_src);
84        }
85        else
86        {
87            confess "Unknown genome source\n";
88        }
89        return $self;
90    }
91    
92    sub create_from_fig
93    {
94        my($self, $fig) = @_;
95    
96        my $gdir = $fig->organism_directory($self->genome_id);
97        copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));
98        mkdir(catfile($self->genome_dir, "Features"));
99        my @pegs;
100        for my $ftype (qw(peg rna))
101        {
102            my $ofdir = catfile($gdir, "Features", $ftype);
103            my $nfdir = catfile($self->genome_dir, "Features", $ftype);
104    
105            if (open(OT, "<", "$ofdir/tbl"))
106            {
107                mkdir($nfdir);
108                open(NT, ">", "$nfdir/tbl") or confess "Cannot write $nfdir/tbl: $!";
109                while (<OT>)
110                {
111                    my($id) = /^(\S+)\t/;
112                    if (!$fig->is_deleted_fid($id))
113                    {
114                        print NT $_;
115                        push(@pegs, $id) if $ftype eq 'peg';
116                    }
117                }
118                close(OT);
119                close(NT);
120            }
121            copy("$ofdir/fasta", "$nfdir/fasta");
122        }
123    
124        open(SS, ">", catfile($self->genome_dir, "subsystem.data"));
125        my %phash = $fig->subsystems_for_pegs_complete(\@pegs);
126        for my $peg (keys %phash)
127        {
128            my $list = $phash{$peg};
129            next unless $list;
130    
131            for my $ent (@$list)
132            {
133                print SS join("\t", $peg, @$ent), "\n";
134            }
135        }
136        close(SS);
137    
138        open(AF, ">", catfile($self->genome_dir, "assigned_functions"));
139        my $fns = $fig->function_of_bulk(\@pegs);
140        for my $peg (@pegs)
141        {
142            print AF join("\t", $peg, $fns->{$peg}), "\n";
143        }
144        close(AF);
145    }
146    
147    sub create_from_sap
148    {
149        my($self, $sap) = @_;
150        confess "create_from_sap not yet implemented";
151  }  }
152    
153  sub parse_probe_format_1  sub parse_probe_format_1
# Line 224  Line 329 
329          confess "Could not find any tbl files in $feature_dir";          confess "Could not find any tbl files in $feature_dir";
330      }      }
331    
332      system_with_redirect([executable_for("make_probes_to_genes"),      $self->run([executable_for("make_probes_to_genes"),
333                            $probes_fasta,                            $probes_fasta,
334                            catfile($self->genome_dir, 'contigs'),                            catfile($self->genome_dir, 'contigs'),
335                            $tbls[0],                            $tbls[0],
# Line 235  Line 340 
340                           { stderr => catfile($self->expr_dir, 'problems') });                           { stderr => catfile($self->expr_dir, 'problems') });
341    
342    
343      system_with_redirect([executable_for("remove_multiple_occurring_probes"),      $self->run([executable_for("remove_multiple_occurring_probes"),
344                            $peg_probe_table,                            $peg_probe_table,
345                            $probe_occ_table],                            $probe_occ_table],
346                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
# Line 297  Line 402 
402      print $fh "library(makecdfenv);\n";      print $fh "library(makecdfenv);\n";
403      print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);      print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
404      close($fh);      close($fh);
     system("cat", $tempfile);  
405      system("Rscript", $tempfile);      system("Rscript", $tempfile);
406      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
407    
408      local($ENV{R_LIBS}) = $libdir;      local($ENV{R_LIBS}) = $libdir;
409      my @cmd = (executable_for("RunRMA"),      $self->run([executable_for("RunRMA"),
410                 "data",                 "data",
411                 catfile($self->expr_dir, "peg.probe.table"),                 catfile($self->expr_dir, "peg.probe.table"),
412                 catfile($self->expr_dir, "probe.no.match"),                 catfile($self->expr_dir, "probe.no.match"),
413                 $expt_dir,                 $expt_dir,
414                 $self->expr_dir);                  $self->expr_dir]);
415    
     my $rc = system_with_redirect(\@cmd);  
 #    my $rc = system_with_redirect(\@cmd, { stderr => "rma.stderr" });  
     if ($rc != 0)  
     {  
         confess("Error running RMA: rc=$rc cmd=@cmd");  
     }  
416      my $output = catfile($self->expr_dir, "rma_normalized.tab");      my $output = catfile($self->expr_dir, "rma_normalized.tab");
417      if (! -f $output)      if (! -f $output)
418      {      {
# Line 324  Line 422 
422    
423  sub compute_atomic_regulons  sub compute_atomic_regulons
424  {  {
425      my($self) = @_;      my($self, $pearson_cutoff) = @_;
426    
427        $pearson_cutoff ||= 0.7;
428    
429      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
430      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
431      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
432        my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
433        my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
434    
435    
436        $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
437               { stdout => $coreg_clusters });
438        $self->run([executable_for("make_coreg_conjectures_based_on_subsys"), $self->expr_dir],
439               { stdout => $coreg_subsys });
440    
441        $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
442               { stdout => $merged_clusters });
443        $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
444    
445        $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
446               { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
447    
448        $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
449               { stdout => catfile($self->expr_dir, "split.clusters") });
450    
451    }
452    
453    sub run
454    {
455        my($self, $cmd, $redirect) = @_;
456    
457        print "Run @$cmd\n";
458        my $rc = system_with_redirect($cmd, $redirect);
459        if ($rc != 0)
460        {
461            confess "Command failed: @$cmd\n";
462        }
463  }  }
464    
465    sub get_experiment_names
466    {
467        my($self) = @_;
468        my $f = catfile($self->expr_dir, "experiment.names");
469        my $fh;
470        open($fh, "<", $f) or confess "Could not open $f: $!";
471        my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
472        close($fh);
473        return @out;
474    }
475    
476    sub get_experiment_on_off_calls
477    {
478        my($self, $expt) = @_;
479    
480        my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
481        my $fh;
482        open($fh, "<", $f) or confess "Could not open $f: $!";
483        my $names = <$fh>;
484        chomp $names;
485        my @names = split(/\t/, $names);
486        my $idx = 0;
487        my $expt_idx;
488        foreach my $n (@names)
489        {
490            if ($n eq $expt)
491            {
492                $expt_idx = $idx;
493                last;
494            }
495            $idx++;
496        }
497        if (!defined($expt_idx))
498        {
499            confess("Could not find experiment $expt in $f");
500        }
501    
502        my $calls = {};
503        while (<$fh>)
504        {
505            chomp;
506            my($peg, @calls) = split(/\t/);
507            #
508            # +1 because the row[0] element is the peg, and our index is
509            # zero-based.
510            #
511            $calls->{$peg} = $calls[$expt_idx + 1];
512        }
513    
514        close($fh);
515        return($calls);
516    
517    }
518    
519    =head3 save_model_gene_activity
520    
521        $e->save_model_gene_activity($data)
522    
523    Save the results of a modeling run for a given experiment.
524    
525    $data is of the form { experiment_id => $data_hash }
526    
527    =cut
528    
529    sub save_model_gene_activity
530    {
531        my($self, $data) = @_;
532    }
533    
534    sub all_features
535    {
536        my($self, $type) = @_;
537    
538        my @ftypes;
539        my $fdir = catfile($self->genome_dir, "Features");
540        if (defined($type))
541        {
542            @ftypes = ($type);
543        }
544        else
545        {
546            opendir(D, $fdir);
547            @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
548            closedir(D);
549        }
550        my @out;
551        for my $ftype (@ftypes)
552        {
553            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
554            {
555                push(@out, map { /^(\S+)/; $1 } <TBL>);
556                close(TBL);
557            }
558        }
559        return @out;
560    }
561    
562    sub fid_locations
563    {
564        my($self, $fids) = @_;
565    
566        my %fids;
567        $fids{$_}++ for @$fids;
568    
569        my $genome_id = $self->genome_id;
570    
571        my $fdir = catfile($self->genome_dir, "Features");
572        opendir(D, $fdir);
573        my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
574        closedir(D);
575    
576        my $out = {};
577        for my $ftype (@ftypes)
578        {
579            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
580            {
581                while (<TBL>)
582                {
583                    my($id, $locs) = /^(\S+)\t(\S+)\t/;
584    
585                    if ($fids{$id})
586                    {
587                        $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
588                    }
589                }
590                close(TBL);
591            }
592        }
593        return $out;
594    }
595    
596    sub ids_in_subsystems
597    {
598        my($self) = @_;
599    
600        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
601        my $res = {};
602        while (<SS>)
603        {
604            chomp;
605            my($peg, $ss, $role, $variant) = split(/\t/);
606            $ss =~ s/\s+/_/g;
607            push(@{$res->{$ss}->{$role}}, $peg);
608        }
609        close(SS);
610        return $res;
611    }
612    
613    sub ids_to_subsystems
614    {
615        my($self, $ids) = @_;
616    
617        my %ids;
618        $ids{$_} = 1 for @$ids;
619    
620        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
621        my $res = {};
622        while (<SS>)
623        {
624            chomp;
625            my($peg, $ss, $role, $variant) = split(/\t/);
626            if ($ids{$peg})
627            {
628                push(@{$res->{$peg}}, $ss);
629            }
630        }
631        close(SS);
632        return $res;
633    }
634    
635    sub ids_to_functions
636    {
637        my($self, $ids) = @_;
638        open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
639        my %ids;
640        $ids{$_} = 1 for @$ids;
641        my $res = {};
642    
643        while (<AF>)
644        {
645            chomp;
646            my($id, $fn) = split(/\t/);
647            $res->{$id} = $fn if $ids{$id};
648        }
649        close(AF);
650        return $res;
651    }
652    
653  1;  1;
654    

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3