[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.8, Mon Jan 24 19:42:26 2011 UTC
# Line 1  Line 1 
1  package ExpressionDir;  package ExpressionDir;
2    
3    use FileHandle;
4    use gjoseqlib;
5  use Data::Dumper;  use Data::Dumper;
6  use strict;  use strict;
7  use SeedAware;  use SeedAware;
# Line 10  Line 12 
12  use Carp;  use Carp;
13  use Fcntl ':seek';  use Fcntl ':seek';
14    
15  __PACKAGE__->mk_accessors(qw(genome_dir expr_dir error_dir));  __PACKAGE__->mk_accessors(qw(genome_dir expr_dir genome_id));
16    
17    our @probe_parsers = qw(parse_probe_format_1lq
18                            parse_probe_format_1
19                            parse_probe_format_2
20                            parse_probe_format_3
21                            parse_probe_format_shew
22                            parse_probe_format_native);
23    
24  =head3 new  =head3 new
25    
26      my $edir = ExpressionDir->new($genome_dir);      my $edir = ExpressionDir->new($expr_dir);
27    
28    Create a new ExpressionDir object from an existing expression dir.
29    
 Create a new ExpressionDir object to be associated with the given genome directory.  
 If a subdirectory ExpressionData does not yet exist, one is created.  
30    
31  =cut  =cut
32    
33  sub new  sub new
34  {  {
35      my($class, $genome_dir) = @_;      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    
53        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      my $edir = catfile($genome_dir, "ExpressionData");      if (! -d $expr_dir)
     if (! -d $edir)  
66      {      {
67          mkdir($edir);          mkdir($expr_dir);
68      }      }
69      my $errdir = catfile($genome_dir, 'errors');      my $genome_dir = catfile($expr_dir, $genome_id);
70      if (! -d $errdir)      if (! -d $genome_dir)
71      {      {
72          mkdir($errdir);          mkdir($genome_dir);
73      }      }
74    
75        open(GF, ">", catfile($expr_dir, "GENOME_ID"));
76        print GF "$genome_id\n";
77        close(GF);
78    
79      my $self = {      my $self = {
80          genome_dir => $genome_dir,          genome_dir => $genome_dir,
81          expr_dir => $edir,          genome_id => $genome_id,
82          error_dir => $errdir,          expr_dir => $expr_dir,
83      };      };
84      return bless $self, $class;      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    
107        if (! -d $gdir)
108        {
109            confess "Genome directory $gdir not found";
110        }
111    
112        copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));
113        mkdir(catfile($self->genome_dir, "Features"));
114        my @pegs;
115        my %locs;
116        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                    my($id) = /^(\S+)\t(\S+)/;
128                    if (!$fig->is_deleted_fid($id))
129                    {
130                        print NT $_;
131                        if ($ftype eq 'peg')
132                        {
133                            push(@pegs, $id);
134                            $locs{$id} = $2;
135                        }
136    
137                    }
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            next unless $list;
151    
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    
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    }
179    
180    sub create_from_sap
181    {
182        my($self, $sap) = @_;
183        confess "create_from_sap not yet implemented";
184    }
185    
186    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  sub parse_probe_format_1  sub parse_probe_format_1
# Line 134  Line 320 
320      return 1;      return 1;
321  }  }
322    
323    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                return undef;
353            }
354        }
355    
356        close($out);
357        close($fh);
358        return 1;
359    }
360    
361    #
362    # 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  # Our "native" format, used for passing through pre-parsed data.  # Our "native" format, used for passing through pre-parsed data.
397  #  #
398  sub parse_probe_format_3  sub parse_probe_format_native
399  {  {
400      my($self, $in_file, $out_file) = @_;      my($self, $in_file, $out_file) = @_;
401    
# Line 178  Line 436 
436  {  {
437      my($self, $probes) = @_;      my($self, $probes) = @_;
438    
439      my $my_probes = catfile($self->expr_dir, "probes.in");      my($probe_suffix) = $probes =~ m,(\.[^/.]+)$,;
440    
441        my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");
442    
443      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
444    
445      my $probes_fasta = catfile($self->expr_dir, "probes.fasta");      my $probes_fasta = catfile($self->expr_dir, "probes");
446    
447      #      #
448      # Attempt to translate probe file.      # Attempt to translate probe file.
449      #      #
450      my $success;      my $success;
451      for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))      for my $meth (@probe_parsers)
452      {      {
453          if ($self->$meth($my_probes, $probes_fasta))          if ($self->$meth($my_probes, $probes_fasta))
454          {          {
# Line 224  Line 484 
484          confess "Could not find any tbl files in $feature_dir";          confess "Could not find any tbl files in $feature_dir";
485      }      }
486    
487      system_with_redirect([executable_for("make_probes_to_genes"),      $self->run([executable_for("make_probes_to_genes"),
488                            $probes_fasta,                            $probes_fasta,
489                            catfile($self->genome_dir, 'contigs'),                            catfile($self->genome_dir, 'contigs'),
490                            $tbls[0],                            $tbls[0],
# Line 235  Line 495 
495                           { stderr => catfile($self->expr_dir, 'problems') });                           { stderr => catfile($self->expr_dir, 'problems') });
496    
497    
498      system_with_redirect([executable_for("remove_multiple_occurring_probes"),      $self->run([executable_for("remove_multiple_occurring_probes"),
499                            $peg_probe_table,                            $peg_probe_table,
500                            $probe_occ_table],                            $probe_occ_table],
501                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
# Line 297  Line 557 
557      print $fh "library(makecdfenv);\n";      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);      print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
559      close($fh);      close($fh);
     system("cat", $tempfile);  
560      system("Rscript", $tempfile);      system("Rscript", $tempfile);
561      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
562    
563      local($ENV{R_LIBS}) = $libdir;      local($ENV{R_LIBS}) = $libdir;
564      my @cmd = (executable_for("RunRMA"),      $self->run([executable_for("RunRMA"),
565                 "data",                 "data",
566                 catfile($self->expr_dir, "peg.probe.table"),                 catfile($self->expr_dir, "peg.probe.table"),
567                 catfile($self->expr_dir, "probe.no.match"),                 catfile($self->expr_dir, "probe.no.match"),
568                 $expt_dir,                 $expt_dir,
569                 $self->expr_dir);                  $self->expr_dir]);
570    
571      my $rc = system_with_redirect(\@cmd);      my $output = catfile($self->expr_dir, "rma_normalized.tab");
572  #    my $rc = system_with_redirect(\@cmd, { stderr => "rma.stderr" });      if (! -f $output)
573      if ($rc != 0)      {
574            confess("Output file $output was not generated");
575        }
576    }
577    
578    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        $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
602    }
603    
604    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          confess("Error running RMA: rc=$rc cmd=@cmd");                  print "Dup: chip_id=$chip_id tag=$tag  old=$locus_map{$tag}\n";
665                    next if $chip_id =~ /_s_/ ;
666                    next if $chip_id =~ /^Rick/ && $self->genome_id eq '452659.3';
667      }      }
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    
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    
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");      my $output = catfile($self->expr_dir, "rma_normalized.tab");
761      if (! -f $output)      if (! -f $output)
762      {      {
# Line 324  Line 766 
766    
767  sub compute_atomic_regulons  sub compute_atomic_regulons
768  {  {
769      my($self) = @_;      my($self, $pearson_cutoff) = @_;
770    
771        $pearson_cutoff ||= 0.7;
772    
773      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
774      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
775      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
776        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        if (-s $probes_always_on == 0)
789        {
790            confess "No always-on probes were found";
791  }  }
792    
793        $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        $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
800    }
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    }
813    
814    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    
883    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  1;  1;
1003    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3