[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.6, Thu Jan 20 22:05:57 2011 UTC
# Line 1  Line 1 
1  package ExpressionDir;  package ExpressionDir;
2    
3    use gjoseqlib;
4  use Data::Dumper;  use Data::Dumper;
5  use strict;  use strict;
6  use SeedAware;  use SeedAware;
# Line 10  Line 11 
11  use Carp;  use Carp;
12  use Fcntl ':seek';  use Fcntl ':seek';
13    
14  __PACKAGE__->mk_accessors(qw(genome_dir expr_dir error_dir));  __PACKAGE__->mk_accessors(qw(genome_dir expr_dir genome_id));
15    
16    our @probe_parsers = qw(parse_probe_format_1lq
17                            parse_probe_format_1
18                            parse_probe_format_2
19                            parse_probe_format_3
20                            parse_probe_format_native);
21    
22  =head3 new  =head3 new
23    
24      my $edir = ExpressionDir->new($genome_dir);      my $edir = ExpressionDir->new($expr_dir);
25    
26    Create a new ExpressionDir object from an existing expression dir.
27    
 Create a new ExpressionDir object to be associated with the given genome directory.  
 If a subdirectory ExpressionData does not yet exist, one is created.  
28    
29  =cut  =cut
30    
31  sub new  sub new
32  {  {
33      my($class, $genome_dir) = @_;      my($class, $expr_dir) = @_;
34    
35        my $gfile = catfile($expr_dir, "GENOME_ID");
36        open(GF, "<", $gfile) or die "Cannot open $expr_dir/GENOME_ID: $!";
37        my $genome_id = <GF>;
38        chomp $genome_id;
39        close(GF);
40    
41        my $self = {
42            genome_dir => catfile($expr_dir, $genome_id),
43            genome_id => $genome_id,
44            expr_dir => $expr_dir,
45        };
46        return bless $self, $class;
47    }
48    
49    =head3 create
50    
51      my $edir = catfile($genome_dir, "ExpressionData");      my $edir = ExpressionDir->create($expr_dir, $genome_id, $genome_src)
52      if (! -d $edir)  
53    Create a new expression directory from the given genome id and genome data source.
54    The data source is either a FIG or FIGV object, or a SAPserver object
55    that points at a server from which the data can be extracted.
56    
57    =cut
58    
59    sub create
60      {      {
61          mkdir($edir);      my($class, $expr_dir, $genome_id, $genome_src) = @_;
62    
63        if (! -d $expr_dir)
64        {
65            mkdir($expr_dir);
66      }      }
67      my $errdir = catfile($genome_dir, 'errors');      my $genome_dir = catfile($expr_dir, $genome_id);
68      if (! -d $errdir)      if (! -d $genome_dir)
69      {      {
70          mkdir($errdir);          mkdir($genome_dir);
71      }      }
72    
73        open(GF, ">", catfile($expr_dir, "GENOME_ID"));
74        print GF "$genome_id\n";
75        close(GF);
76    
77      my $self = {      my $self = {
78          genome_dir => $genome_dir,          genome_dir => $genome_dir,
79          expr_dir => $edir,          genome_id => $genome_id,
80          error_dir => $errdir,          expr_dir => $expr_dir,
81      };      };
82      return bless $self, $class;      bless $self, $class;
83    
84        if (ref($genome_src) =~ /^FIG/)
85        {
86            $self->create_from_fig($genome_src);
87        }
88        elsif (ref($genome_src =~ /^SAP/))
89        {
90            $self->create_from_sap($genome_src);
91        }
92        else
93        {
94            confess "Unknown genome source\n";
95        }
96        return $self;
97    }
98    
99    sub create_from_fig
100    {
101        my($self, $fig) = @_;
102    
103        my $gdir = $fig->organism_directory($self->genome_id);
104    
105        if (! -d $gdir)
106        {
107            confess "Genome directory $gdir not found";
108        }
109    
110        copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));
111        mkdir(catfile($self->genome_dir, "Features"));
112        my @pegs;
113        my %locs;
114        for my $ftype (qw(peg rna))
115        {
116            my $ofdir = catfile($gdir, "Features", $ftype);
117            my $nfdir = catfile($self->genome_dir, "Features", $ftype);
118    
119            if (open(OT, "<", "$ofdir/tbl"))
120            {
121                mkdir($nfdir);
122                open(NT, ">", "$nfdir/tbl") or confess "Cannot write $nfdir/tbl: $!";
123                while (<OT>)
124                {
125                    my($id) = /^(\S+)\t(\S+)/;
126                    if (!$fig->is_deleted_fid($id))
127                    {
128                        print NT $_;
129                        if ($ftype eq 'peg')
130                        {
131                            push(@pegs, $id);
132                            $locs{$id} = $2;
133                        }
134    
135                    }
136                }
137                close(OT);
138                close(NT);
139            }
140            copy("$ofdir/fasta", "$nfdir/fasta");
141        }
142    
143        open(SS, ">", catfile($self->genome_dir, "subsystem.data"));
144        my %phash = $fig->subsystems_for_pegs_complete(\@pegs);
145        for my $peg (keys %phash)
146        {
147            my $list = $phash{$peg};
148            next unless $list;
149    
150            for my $ent (@$list)
151            {
152                print SS join("\t", $peg, @$ent), "\n";
153            }
154        }
155        close(SS);
156    
157        open(AF, ">", catfile($self->genome_dir, "assigned_functions"));
158        my $fns = $fig->function_of_bulk(\@pegs);
159        for my $peg (@pegs)
160        {
161            print AF join("\t", $peg, $fns->{$peg}), "\n";
162        }
163        close(AF);
164    
165        open(PD, ">", catfile($self->genome_dir, "peg_dna.fasta"));
166        for my $peg (@pegs)
167        {
168            my $dna = $fig->dna_seq($self->genome_id, split(/,/, $locs{$peg}));
169            if ($dna eq '')
170            {
171                die "no dna for ", $self->genome_id, " $peg $locs{$peg}\n";
172            }
173            print_alignment_as_fasta(\*PD, [$peg, undef, $dna]);
174        }
175        close(PD);
176    }
177    
178    sub create_from_sap
179    {
180        my($self, $sap) = @_;
181        confess "create_from_sap not yet implemented";
182    }
183    
184    sub parse_probe_format_1lq
185    {
186        my($self, $in_file, $out_file) = @_;
187    
188        my($fh);
189    
190        if ($in_file !~ /\.1lq$/)
191        {
192            return undef;
193        }
194    
195        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
196    
197        my $out;
198        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
199    
200        # Skip 3 header lines.
201        $_ = <$fh> for 1..3;
202    
203        while (defined($_ = <$fh>))
204        {
205            if ($_ =~ /(\d+)\s+(\d+)\s+([ACGT]+)\s+(-?\d+)\s/)
206            {
207                if (length($3) < 15)
208                {
209                    close($fh);
210                    close($out);
211                    confess "Bad length at line $. of $in_file";
212                    return undef;
213                }
214                next if ($4 =~ /\d+3$/); #mismatch probe
215                my($x,$y,$seq) = ($1,$2,$3);
216                $seq = scalar reverse $seq;
217                print $out "$x\_$y\t$seq\n";
218            }
219            else
220            {
221                #
222                # We expect some lines not to match.
223                #
224            }
225        }
226    
227        close($fh);
228        close($out);
229        return 1;
230  }  }
231    
232  sub parse_probe_format_1  sub parse_probe_format_1
# Line 134  Line 318 
318      return 1;      return 1;
319  }  }
320    
321    sub parse_probe_format_3
322    {
323        my($self, $in_file, $out_file) = @_;
324    
325        my($fh);
326    
327        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
328        my $out;
329        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
330    
331        local $/ = "\n>";
332    
333        while (defined($_ = <$fh>))
334        {
335            if ($_ =~ /^>?\S+\s+(\d+)\s+(\d+)[^\n]+\n(\S+)/s)
336            {
337                my $x = $1;
338                my $y = $2;
339                my $seq = $3;
340                if ($seq ne "!")
341                {
342                    if ($seq !~ /^[ACGT]+$/) { die $_ }
343                    if (length($seq) < 15) { print STDERR "BAD: $_"; die "Failed" }
344                    print $out "$x\_$y\t$seq\n";
345                }
346            }
347            else
348            {
349                print STDERR "failed to parse: $_";
350                die "aborted";
351            }
352        }
353    
354        close($out);
355        close($fh);
356        return 1;
357    }
358    
359  #  #
360  # Our "native" format, used for passing through pre-parsed data.  # Our "native" format, used for passing through pre-parsed data.
361  #  #
362  sub parse_probe_format_3  sub parse_probe_format_native
363  {  {
364      my($self, $in_file, $out_file) = @_;      my($self, $in_file, $out_file) = @_;
365    
# Line 178  Line 400 
400  {  {
401      my($self, $probes) = @_;      my($self, $probes) = @_;
402    
403      my $my_probes = catfile($self->expr_dir, "probes.in");      my($probe_suffix) = $probes =~ /(\.[^.]+)$/;
404    
405        my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");
406    
407      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
408    
409      my $probes_fasta = catfile($self->expr_dir, "probes.fasta");      my $probes_fasta = catfile($self->expr_dir, "probes");
410    
411      #      #
412      # Attempt to translate probe file.      # Attempt to translate probe file.
413      #      #
414      my $success;      my $success;
415      for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))      for my $meth (@probe_parsers)
416      {      {
417          if ($self->$meth($my_probes, $probes_fasta))          if ($self->$meth($my_probes, $probes_fasta))
418          {          {
# Line 224  Line 448 
448          confess "Could not find any tbl files in $feature_dir";          confess "Could not find any tbl files in $feature_dir";
449      }      }
450    
451      system_with_redirect([executable_for("make_probes_to_genes"),      $self->run([executable_for("make_probes_to_genes"),
452                            $probes_fasta,                            $probes_fasta,
453                            catfile($self->genome_dir, 'contigs'),                            catfile($self->genome_dir, 'contigs'),
454                            $tbls[0],                            $tbls[0],
# Line 235  Line 459 
459                           { stderr => catfile($self->expr_dir, 'problems') });                           { stderr => catfile($self->expr_dir, 'problems') });
460    
461    
462      system_with_redirect([executable_for("remove_multiple_occurring_probes"),      $self->run([executable_for("remove_multiple_occurring_probes"),
463                            $peg_probe_table,                            $peg_probe_table,
464                            $probe_occ_table],                            $probe_occ_table],
465                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
# Line 297  Line 521 
521      print $fh "library(makecdfenv);\n";      print $fh "library(makecdfenv);\n";
522      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);
523      close($fh);      close($fh);
     system("cat", $tempfile);  
524      system("Rscript", $tempfile);      system("Rscript", $tempfile);
525      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
526    
527      local($ENV{R_LIBS}) = $libdir;      local($ENV{R_LIBS}) = $libdir;
528      my @cmd = (executable_for("RunRMA"),      $self->run([executable_for("RunRMA"),
529                 "data",                 "data",
530                 catfile($self->expr_dir, "peg.probe.table"),                 catfile($self->expr_dir, "peg.probe.table"),
531                 catfile($self->expr_dir, "probe.no.match"),                 catfile($self->expr_dir, "probe.no.match"),
532                 $expt_dir,                 $expt_dir,
533                 $self->expr_dir);                  $self->expr_dir]);
534    
535      my $rc = system_with_redirect(\@cmd);      my $output = catfile($self->expr_dir, "rma_normalized.tab");
536  #    my $rc = system_with_redirect(\@cmd, { stderr => "rma.stderr" });      if (! -f $output)
537      if ($rc != 0)      {
538            confess("Output file $output was not generated");
539        }
540    }
541    
542    sub compute_rma_normalized_from_sif
543    {
544        my($self, $cdf_file, $sif_file, $expt_dir) = @_;
545    
546        my $my_cdf = catfile($self->expr_dir, "expr.cdf");
547        copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
548    
549        my $my_sif = catfile($self->expr_dir, "expr.sif");
550        copy($sif_file, $my_sif) or confess "Cannot copy $sif_file to $my_sif: $!";
551    
552        #
553        # Create the sif2peg mapping.
554        #
555    
556        my $sif2peg = catfile($self->expr_dir, "sif2peg");
557        if (! -f $sif2peg)
558      {      {
559          confess("Error running RMA: rc=$rc cmd=@cmd");          $self->run([executable_for('map_sif_to_pegs'),
560                        $my_sif,
561                        catfile($self->genome_dir, "peg_dna.fasta")],
562                   { stdout => $sif2peg });
563      }      }
564    
565        #
566        # We need to build the R library for this cdf.
567        #
568        my($fh, $tempfile) = tempfile();
569    #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
570    
571        my $cdf_path = $self->expr_dir;
572        my $libdir = catfile($self->expr_dir, "r_lib");
573        -d $libdir or mkdir $libdir;
574        my $pkgdir = catfile($self->expr_dir, "r_pkg");
575        -d $pkgdir or mkdir $pkgdir;
576    
577        print $fh "library(makecdfenv);\n";
578        print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
579        close($fh);
580        system("Rscript", $tempfile);
581        system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
582    
583        local($ENV{R_LIBS}) = $libdir;
584        $self->run([executable_for("RunRMA_SIF_format"),
585                    "data",
586                    $expt_dir,
587                    $self->expr_dir]);
588    
589      my $output = catfile($self->expr_dir, "rma_normalized.tab");      my $output = catfile($self->expr_dir, "rma_normalized.tab");
590      if (! -f $output)      if (! -f $output)
591      {      {
# Line 324  Line 595 
595    
596  sub compute_atomic_regulons  sub compute_atomic_regulons
597  {  {
598      my($self) = @_;      my($self, $pearson_cutoff) = @_;
599    
600        $pearson_cutoff ||= 0.7;
601    
602      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
603      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
604      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
605        my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
606        my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
607    
608        $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
609               { stdout => $coreg_clusters });
610        $self->run([executable_for("make_coreg_conjectures_based_on_subsys"), $self->expr_dir],
611               { stdout => $coreg_subsys });
612    
613        $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
614               { stdout => $merged_clusters });
615        $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
616    
617        if (-s $probes_always_on == 0)
618        {
619            confess "No always-on probes were found";
620        }
621    
622        $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
623               { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
624    
625        $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
626               { stdout => catfile($self->expr_dir, "split.clusters") });
627    
628        $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
629    }
630    
631    sub run
632    {
633        my($self, $cmd, $redirect) = @_;
634    
635        print "Run @$cmd\n";
636        my $rc = system_with_redirect($cmd, $redirect);
637        if ($rc != 0)
638        {
639            confess "Command failed: @$cmd\n";
640        }
641    }
642    
643    sub get_experiment_names
644    {
645        my($self) = @_;
646        my $f = catfile($self->expr_dir, "experiment.names");
647        my $fh;
648        open($fh, "<", $f) or confess "Could not open $f: $!";
649        my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
650        close($fh);
651        return @out;
652    }
653    
654    sub get_experiment_on_off_calls
655    {
656        my($self, $expt) = @_;
657    
658        my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
659        my $fh;
660        open($fh, "<", $f) or confess "Could not open $f: $!";
661        my $names = <$fh>;
662        chomp $names;
663        my @names = split(/\t/, $names);
664        my $idx = 0;
665        my $expt_idx;
666        foreach my $n (@names)
667        {
668            if ($n eq $expt)
669            {
670                $expt_idx = $idx;
671                last;
672            }
673            $idx++;
674        }
675        if (!defined($expt_idx))
676        {
677            confess("Could not find experiment $expt in $f");
678        }
679    
680        my $calls = {};
681        while (<$fh>)
682        {
683            chomp;
684            my($peg, @calls) = split(/\t/);
685            #
686            # +1 because the row[0] element is the peg, and our index is
687            # zero-based.
688            #
689            $calls->{$peg} = $calls[$expt_idx + 1];
690        }
691    
692        close($fh);
693        return($calls);
694    
695    }
696    
697    =head3 save_model_gene_activity
698    
699        $e->save_model_gene_activity($data)
700    
701    Save the results of a modeling run for a given experiment.
702    
703    $data is of the form { experiment_id => $data_hash }
704    
705    =cut
706    
707    sub save_model_gene_activity
708    {
709        my($self, $data) = @_;
710    }
711    
712    sub all_features
713    {
714        my($self, $type) = @_;
715    
716        my @ftypes;
717        my $fdir = catfile($self->genome_dir, "Features");
718        if (defined($type))
719        {
720            @ftypes = ($type);
721        }
722        else
723        {
724            opendir(D, $fdir);
725            @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
726            closedir(D);
727        }
728        my @out;
729        for my $ftype (@ftypes)
730        {
731            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
732            {
733                push(@out, map { /^(\S+)/; $1 } <TBL>);
734                close(TBL);
735            }
736        }
737        return @out;
738    }
739    
740    sub fid_locations
741    {
742        my($self, $fids) = @_;
743    
744        my %fids;
745        $fids{$_}++ for @$fids;
746    
747        my $genome_id = $self->genome_id;
748    
749        my $fdir = catfile($self->genome_dir, "Features");
750        opendir(D, $fdir);
751        my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
752        closedir(D);
753    
754        my $out = {};
755        for my $ftype (@ftypes)
756        {
757            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
758            {
759                while (<TBL>)
760                {
761                    my($id, $locs) = /^(\S+)\t(\S+)\t/;
762    
763                    if ($fids{$id})
764                    {
765                        $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
766                    }
767                }
768                close(TBL);
769            }
770        }
771        return $out;
772    }
773    
774    sub ids_in_subsystems
775    {
776        my($self) = @_;
777    
778        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
779        my $res = {};
780        while (<SS>)
781        {
782            chomp;
783            my($peg, $ss, $role, $variant) = split(/\t/);
784            $ss =~ s/\s+/_/g;
785            push(@{$res->{$ss}->{$role}}, $peg);
786  }  }
787        close(SS);
788        return $res;
789    }
790    
791    sub ids_to_subsystems
792    {
793        my($self, $ids) = @_;
794    
795        my %ids;
796        $ids{$_} = 1 for @$ids;
797    
798        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
799        my $res = {};
800        while (<SS>)
801        {
802            chomp;
803            my($peg, $ss, $role, $variant) = split(/\t/);
804            if ($ids{$peg})
805            {
806                push(@{$res->{$peg}}, $ss);
807            }
808        }
809        close(SS);
810        return $res;
811    }
812    
813    sub ids_to_functions
814    {
815        my($self, $ids) = @_;
816        open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
817        my %ids;
818        $ids{$_} = 1 for @$ids;
819        my $res = {};
820    
821        while (<AF>)
822        {
823            chomp;
824            my($id, $fn) = split(/\t/);
825            $res->{$id} = $fn if $ids{$id};
826        }
827        close(AF);
828        return $res;
829    }
830    
831  1;  1;
832    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3