[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.5, Fri Jan 14 20:06:51 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    our @probe_parsers = qw(parse_probe_format_1lq
16                            parse_probe_format_1
17                            parse_probe_format_2
18                            parse_probe_format_3);
19    
20  =head3 new  =head3 new
21    
22      my $edir = ExpressionDir->new($genome_dir);      my $edir = ExpressionDir->new($expr_dir);
23    
24    Create a new ExpressionDir object from an existing expression dir.
25    
 Create a new ExpressionDir object to be associated with the given genome directory.  
 If a subdirectory ExpressionData does not yet exist, one is created.  
26    
27  =cut  =cut
28    
29  sub new  sub new
30  {  {
31      my($class, $genome_dir) = @_;      my($class, $expr_dir) = @_;
32    
33        my $gfile = catfile($expr_dir, "GENOME_ID");
34        open(GF, "<", $gfile) or die "Cannot open $expr_dir/GENOME_ID: $!";
35        my $genome_id = <GF>;
36        chomp $genome_id;
37        close(GF);
38    
39        my $self = {
40            genome_dir => catfile($expr_dir, $genome_id),
41            genome_id => $genome_id,
42            expr_dir => $expr_dir,
43        };
44        return bless $self, $class;
45    }
46    
47    =head3 create
48    
49        my $edir = ExpressionDir->create($expr_dir, $genome_id, $genome_src)
50    
51    Create a new expression directory from the given genome id and genome data source.
52    The data source is either a FIG or FIGV object, or a SAPserver object
53    that points at a server from which the data can be extracted.
54    
55    =cut
56    
57      my $edir = catfile($genome_dir, "ExpressionData");  sub create
     if (! -d $edir)  
58      {      {
59          mkdir($edir);      my($class, $expr_dir, $genome_id, $genome_src) = @_;
60    
61        if (! -d $expr_dir)
62        {
63            mkdir($expr_dir);
64      }      }
65      my $errdir = catfile($genome_dir, 'errors');      my $genome_dir = catfile($expr_dir, $genome_id);
66      if (! -d $errdir)      if (! -d $genome_dir)
67      {      {
68          mkdir($errdir);          mkdir($genome_dir);
69      }      }
70    
71        open(GF, ">", catfile($expr_dir, "GENOME_ID"));
72        print GF "$genome_id\n";
73        close(GF);
74    
75      my $self = {      my $self = {
76          genome_dir => $genome_dir,          genome_dir => $genome_dir,
77          expr_dir => $edir,          genome_id => $genome_id,
78          error_dir => $errdir,          expr_dir => $expr_dir,
79      };      };
80      return bless $self, $class;      bless $self, $class;
81    
82        if (ref($genome_src) =~ /^FIG/)
83        {
84            $self->create_from_fig($genome_src);
85        }
86        elsif (ref($genome_src =~ /^SAP/))
87        {
88            $self->create_from_sap($genome_src);
89        }
90        else
91        {
92            confess "Unknown genome source\n";
93        }
94        return $self;
95    }
96    
97    sub create_from_fig
98    {
99        my($self, $fig) = @_;
100    
101        my $gdir = $fig->organism_directory($self->genome_id);
102    
103        if (! -d $gdir)
104        {
105            confess "Genome directory $gdir not found";
106        }
107    
108        copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));
109        mkdir(catfile($self->genome_dir, "Features"));
110        my @pegs;
111        for my $ftype (qw(peg rna))
112        {
113            my $ofdir = catfile($gdir, "Features", $ftype);
114            my $nfdir = catfile($self->genome_dir, "Features", $ftype);
115    
116            if (open(OT, "<", "$ofdir/tbl"))
117            {
118                mkdir($nfdir);
119                open(NT, ">", "$nfdir/tbl") or confess "Cannot write $nfdir/tbl: $!";
120                while (<OT>)
121                {
122                    my($id) = /^(\S+)\t/;
123                    if (!$fig->is_deleted_fid($id))
124                    {
125                        print NT $_;
126                        push(@pegs, $id) if $ftype eq 'peg';
127                    }
128                }
129                close(OT);
130                close(NT);
131            }
132            copy("$ofdir/fasta", "$nfdir/fasta");
133        }
134    
135        open(SS, ">", catfile($self->genome_dir, "subsystem.data"));
136        my %phash = $fig->subsystems_for_pegs_complete(\@pegs);
137        for my $peg (keys %phash)
138        {
139            my $list = $phash{$peg};
140            next unless $list;
141    
142            for my $ent (@$list)
143            {
144                print SS join("\t", $peg, @$ent), "\n";
145            }
146        }
147        close(SS);
148    
149        open(AF, ">", catfile($self->genome_dir, "assigned_functions"));
150        my $fns = $fig->function_of_bulk(\@pegs);
151        for my $peg (@pegs)
152        {
153            print AF join("\t", $peg, $fns->{$peg}), "\n";
154        }
155        close(AF);
156    }
157    
158    sub create_from_sap
159    {
160        my($self, $sap) = @_;
161        confess "create_from_sap not yet implemented";
162    }
163    
164    sub parse_probe_format_1lq
165    {
166        my($self, $in_file, $out_file) = @_;
167    
168        my($fh);
169    
170        if ($in_file !~ /\.1lq$/)
171        {
172            return undef;
173        }
174    
175        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
176    
177        my $out;
178        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
179    
180        # Skip 3 header lines.
181        $_ = <$fh> for 1..3;
182    
183        while (defined($_ = <$fh>))
184        {
185            if ($_ =~ /(\d+)\s+(\d+)\s+([ACGT]+)\s+(-?\d+)\s/)
186            {
187                if (length($3) < 15)
188                {
189                    close($fh);
190                    close($out);
191                    confess "Bad length at line $. of $in_file";
192                    return undef;
193                }
194                next if ($4 =~ /\d+3$/); #mismatch probe
195                my($x,$y,$seq) = ($1,$2,$3);
196                $seq = scalar reverse $seq;
197                print $out "$x\_$y\t$seq\n";
198            }
199            else
200            {
201                #
202                # We expect some lines not to match.
203                #
204            }
205        }
206    
207        close($fh);
208        close($out);
209        return 1;
210  }  }
211    
212  sub parse_probe_format_1  sub parse_probe_format_1
# Line 178  Line 342 
342  {  {
343      my($self, $probes) = @_;      my($self, $probes) = @_;
344    
345      my $my_probes = catfile($self->expr_dir, "probes.in");      my($probe_suffix) = $probes =~ /(\.[^.]+)$/;
346    
347        my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");
348    
349      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
350    
351      my $probes_fasta = catfile($self->expr_dir, "probes.fasta");      my $probes_fasta = catfile($self->expr_dir, "probes");
352    
353      #      #
354      # Attempt to translate probe file.      # Attempt to translate probe file.
355      #      #
356      my $success;      my $success;
357      for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))      for my $meth (@probe_parsers)
358      {      {
359          if ($self->$meth($my_probes, $probes_fasta))          if ($self->$meth($my_probes, $probes_fasta))
360          {          {
# Line 224  Line 390 
390          confess "Could not find any tbl files in $feature_dir";          confess "Could not find any tbl files in $feature_dir";
391      }      }
392    
393      system_with_redirect([executable_for("make_probes_to_genes"),      $self->run([executable_for("make_probes_to_genes"),
394                            $probes_fasta,                            $probes_fasta,
395                            catfile($self->genome_dir, 'contigs'),                            catfile($self->genome_dir, 'contigs'),
396                            $tbls[0],                            $tbls[0],
# Line 235  Line 401 
401                           { stderr => catfile($self->expr_dir, 'problems') });                           { stderr => catfile($self->expr_dir, 'problems') });
402    
403    
404      system_with_redirect([executable_for("remove_multiple_occurring_probes"),      $self->run([executable_for("remove_multiple_occurring_probes"),
405                            $peg_probe_table,                            $peg_probe_table,
406                            $probe_occ_table],                            $probe_occ_table],
407                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
# Line 297  Line 463 
463      print $fh "library(makecdfenv);\n";      print $fh "library(makecdfenv);\n";
464      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);
465      close($fh);      close($fh);
     system("cat", $tempfile);  
466      system("Rscript", $tempfile);      system("Rscript", $tempfile);
467      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
468    
469      local($ENV{R_LIBS}) = $libdir;      local($ENV{R_LIBS}) = $libdir;
470      my @cmd = (executable_for("RunRMA"),      $self->run([executable_for("RunRMA"),
471                 "data",                 "data",
472                 catfile($self->expr_dir, "peg.probe.table"),                 catfile($self->expr_dir, "peg.probe.table"),
473                 catfile($self->expr_dir, "probe.no.match"),                 catfile($self->expr_dir, "probe.no.match"),
474                 $expt_dir,                 $expt_dir,
475                 $self->expr_dir);                  $self->expr_dir]);
476    
     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");  
     }  
477      my $output = catfile($self->expr_dir, "rma_normalized.tab");      my $output = catfile($self->expr_dir, "rma_normalized.tab");
478      if (! -f $output)      if (! -f $output)
479      {      {
# Line 324  Line 483 
483    
484  sub compute_atomic_regulons  sub compute_atomic_regulons
485  {  {
486      my($self) = @_;      my($self, $pearson_cutoff) = @_;
487    
488        $pearson_cutoff ||= 0.7;
489    
490      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
491      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
492      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
493        my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
494        my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
495    
496        $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
497               { stdout => $coreg_clusters });
498        $self->run([executable_for("make_coreg_conjectures_based_on_subsys"), $self->expr_dir],
499               { stdout => $coreg_subsys });
500    
501        $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
502               { stdout => $merged_clusters });
503        $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
504    
505        if (-s $probes_always_on == 0)
506        {
507            confess "No always-on probes were found";
508        }
509    
510        $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
511               { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
512    
513        $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
514               { stdout => catfile($self->expr_dir, "split.clusters") });
515    
516        $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
517    }
518    
519    sub run
520    {
521        my($self, $cmd, $redirect) = @_;
522    
523        print "Run @$cmd\n";
524        my $rc = system_with_redirect($cmd, $redirect);
525        if ($rc != 0)
526        {
527            confess "Command failed: @$cmd\n";
528        }
529    }
530    
531    sub get_experiment_names
532    {
533        my($self) = @_;
534        my $f = catfile($self->expr_dir, "experiment.names");
535        my $fh;
536        open($fh, "<", $f) or confess "Could not open $f: $!";
537        my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
538        close($fh);
539        return @out;
540    }
541    
542    sub get_experiment_on_off_calls
543    {
544        my($self, $expt) = @_;
545    
546        my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
547        my $fh;
548        open($fh, "<", $f) or confess "Could not open $f: $!";
549        my $names = <$fh>;
550        chomp $names;
551        my @names = split(/\t/, $names);
552        my $idx = 0;
553        my $expt_idx;
554        foreach my $n (@names)
555        {
556            if ($n eq $expt)
557            {
558                $expt_idx = $idx;
559                last;
560            }
561            $idx++;
562        }
563        if (!defined($expt_idx))
564        {
565            confess("Could not find experiment $expt in $f");
566        }
567    
568        my $calls = {};
569        while (<$fh>)
570        {
571            chomp;
572            my($peg, @calls) = split(/\t/);
573            #
574            # +1 because the row[0] element is the peg, and our index is
575            # zero-based.
576            #
577            $calls->{$peg} = $calls[$expt_idx + 1];
578        }
579    
580        close($fh);
581        return($calls);
582    
583    }
584    
585    =head3 save_model_gene_activity
586    
587        $e->save_model_gene_activity($data)
588    
589    Save the results of a modeling run for a given experiment.
590    
591    $data is of the form { experiment_id => $data_hash }
592    
593    =cut
594    
595    sub save_model_gene_activity
596    {
597        my($self, $data) = @_;
598    }
599    
600    sub all_features
601    {
602        my($self, $type) = @_;
603    
604        my @ftypes;
605        my $fdir = catfile($self->genome_dir, "Features");
606        if (defined($type))
607        {
608            @ftypes = ($type);
609        }
610        else
611        {
612            opendir(D, $fdir);
613            @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
614            closedir(D);
615        }
616        my @out;
617        for my $ftype (@ftypes)
618        {
619            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
620            {
621                push(@out, map { /^(\S+)/; $1 } <TBL>);
622                close(TBL);
623            }
624        }
625        return @out;
626    }
627    
628    sub fid_locations
629    {
630        my($self, $fids) = @_;
631    
632        my %fids;
633        $fids{$_}++ for @$fids;
634    
635        my $genome_id = $self->genome_id;
636    
637        my $fdir = catfile($self->genome_dir, "Features");
638        opendir(D, $fdir);
639        my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
640        closedir(D);
641    
642        my $out = {};
643        for my $ftype (@ftypes)
644        {
645            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
646            {
647                while (<TBL>)
648                {
649                    my($id, $locs) = /^(\S+)\t(\S+)\t/;
650    
651                    if ($fids{$id})
652                    {
653                        $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
654                    }
655                }
656                close(TBL);
657            }
658  }  }
659        return $out;
660    }
661    
662    sub ids_in_subsystems
663    {
664        my($self) = @_;
665    
666        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
667        my $res = {};
668        while (<SS>)
669        {
670            chomp;
671            my($peg, $ss, $role, $variant) = split(/\t/);
672            $ss =~ s/\s+/_/g;
673            push(@{$res->{$ss}->{$role}}, $peg);
674        }
675        close(SS);
676        return $res;
677    }
678    
679    sub ids_to_subsystems
680    {
681        my($self, $ids) = @_;
682    
683        my %ids;
684        $ids{$_} = 1 for @$ids;
685    
686        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
687        my $res = {};
688        while (<SS>)
689        {
690            chomp;
691            my($peg, $ss, $role, $variant) = split(/\t/);
692            if ($ids{$peg})
693            {
694                push(@{$res->{$peg}}, $ss);
695            }
696        }
697        close(SS);
698        return $res;
699    }
700    
701    sub ids_to_functions
702    {
703        my($self, $ids) = @_;
704        open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
705        my %ids;
706        $ids{$_} = 1 for @$ids;
707        my $res = {};
708    
709        while (<AF>)
710        {
711            chomp;
712            my($id, $fn) = split(/\t/);
713            $res->{$id} = $fn if $ids{$id};
714        }
715        close(AF);
716        return $res;
717    }
718    
719  1;  1;
720    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3