[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.11, Tue Apr 19 19:55:47 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        my $genome_ss_dir = catfile($gdir, "Subsystems");
146        if (-d $genome_ss_dir && -s "$genome_ss_dir/bindings" && -s "$genome_ss_dir/subsystems")
147        {
148            my $my_ss_dir = $self->genome_dir . "/Subsystems";
149            mkdir($my_ss_dir) or die "Cannot mkdir $my_ss_dir: $!";
150            for my $f (qw(bindings subsystems))
151            {
152                copy("$genome_ss_dir/$f", "$my_ss_dir/$f") or die "Cannot copy $genome_ss_dir/$f to $my_ss_dir/$f: $!";
153            }
154        }
155    
156        open(SS, ">", catfile($self->genome_dir, "subsystem.data"));
157        my %phash = $fig->subsystems_for_pegs_complete(\@pegs);
158        for my $peg (keys %phash)
159        {
160            my $list = $phash{$peg};
161            next unless $list;
162    
163            for my $ent (@$list)
164            {
165                print SS join("\t", $peg, @$ent), "\n";
166            }
167        }
168        close(SS);
169    
170        open(AF, ">", catfile($self->genome_dir, "assigned_functions"));
171        my $fns = $fig->function_of_bulk(\@pegs);
172        for my $peg (@pegs)
173        {
174            print AF join("\t", $peg, $fns->{$peg}), "\n";
175        }
176        close(AF);
177    
178        open(PD, ">", catfile($self->genome_dir, "peg_dna.fasta"));
179        for my $peg (@pegs)
180        {
181            my $dna = $fig->dna_seq($self->genome_id, split(/,/, $locs{$peg}));
182            if ($dna eq '')
183            {
184                die "no dna for ", $self->genome_id, " $peg $locs{$peg}\n";
185            }
186            print_alignment_as_fasta(\*PD, [$peg, undef, $dna]);
187        }
188        close(PD);
189    }
190    
191    sub create_from_sap
192    {
193        my($self, $sap) = @_;
194        confess "create_from_sap not yet implemented";
195    }
196    
197    sub parse_probe_format_1lq
198    {
199        my($self, $in_file, $out_file) = @_;
200    
201        my($fh);
202    
203        if ($in_file !~ /\.1lq$/)
204        {
205            return undef;
206        }
207    
208        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
209    
210        my $out;
211        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
212    
213        # Skip 3 header lines.
214        $_ = <$fh> for 1..3;
215    
216        while (defined($_ = <$fh>))
217        {
218            if ($_ =~ /(\d+)\s+(\d+)\s+([ACGT]+)\s+(-?\d+)\s/)
219            {
220                if (length($3) < 15)
221                {
222                    close($fh);
223                    close($out);
224                    confess "Bad length at line $. of $in_file";
225                    return undef;
226                }
227                next if ($4 =~ /\d+3$/); #mismatch probe
228                my($x,$y,$seq) = ($1,$2,$3);
229                $seq = scalar reverse $seq;
230                print $out "$x\_$y\t$seq\n";
231            }
232            else
233            {
234                #
235                # We expect some lines not to match.
236                #
237            }
238        }
239    
240        close($fh);
241        close($out);
242        return 1;
243  }  }
244    
245  sub parse_probe_format_1  sub parse_probe_format_1
# Line 134  Line 331 
331      return 1;      return 1;
332  }  }
333    
334    sub parse_probe_format_3
335    {
336        my($self, $in_file, $out_file) = @_;
337    
338        my($fh);
339    
340        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
341        my $out;
342        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
343    
344        local $/ = "\n>";
345    
346        while (defined($_ = <$fh>))
347        {
348            if ($_ =~ /^>?\S+\s+(\d+)\s+(\d+)[^\n]+\n(\S+)/s)
349            {
350                my $x = $1;
351                my $y = $2;
352                my $seq = $3;
353                if ($seq ne "!")
354                {
355                    if ($seq !~ /^[ACGT]+$/) { die $_ }
356                    if (length($seq) < 15) { print STDERR "BAD: $_"; die "Failed" }
357                    print $out "$x\_$y\t$seq\n";
358                }
359            }
360            else
361            {
362                print STDERR "failed to parse: $_";
363                return undef;
364            }
365        }
366    
367        close($out);
368        close($fh);
369        return 1;
370    }
371    
372    #
373    # This one showed up in the shewanella data.
374    #
375    sub parse_probe_format_shew
376    {
377        my($self, $in_file, $out_file) = @_;
378    
379        my($fh);
380    
381        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
382        my $l = <$fh>;
383        chomp $l;
384        $l =~ s/\r//;
385    
386        if ($l !~ /x\ty\tprobe_type\tsequence/)
387        {
388            close($fh);
389            return undef;
390        }
391    
392        my $out;
393        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
394    
395        while (<$fh>)
396        {
397            if ($_ =~ /^(\d+)\t(\d+)\tPM\t+([ACGT]+)/)
398            {
399                print $out "$1\_$2\t$3\n";
400            }
401        }
402        close($out);
403        close($fh);
404        return 1;
405    }
406  #  #
407  # Our "native" format, used for passing through pre-parsed data.  # Our "native" format, used for passing through pre-parsed data.
408  #  #
409  sub parse_probe_format_3  sub parse_probe_format_native
410  {  {
411      my($self, $in_file, $out_file) = @_;      my($self, $in_file, $out_file) = @_;
412    
# Line 178  Line 447 
447  {  {
448      my($self, $probes) = @_;      my($self, $probes) = @_;
449    
450      my $my_probes = catfile($self->expr_dir, "probes.in");      my($probe_suffix) = $probes =~ m,(\.[^/.]+)$,;
451    
452        my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");
453    
454      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";      copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
455    
456      my $probes_fasta = catfile($self->expr_dir, "probes.fasta");      my $probes_fasta = catfile($self->expr_dir, "probes");
457    
458      #      #
459      # Attempt to translate probe file.      # Attempt to translate probe file.
460      #      #
461      my $success;      my $success;
462      for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))      for my $meth (@probe_parsers)
463      {      {
464          if ($self->$meth($my_probes, $probes_fasta))          if ($self->$meth($my_probes, $probes_fasta))
465          {          {
# Line 224  Line 495 
495          confess "Could not find any tbl files in $feature_dir";          confess "Could not find any tbl files in $feature_dir";
496      }      }
497    
498      system_with_redirect([executable_for("make_probes_to_genes"),      $self->run([executable_for("make_probes_to_genes"),
499                            $probes_fasta,                            $probes_fasta,
500                            catfile($self->genome_dir, 'contigs'),                            catfile($self->genome_dir, 'contigs'),
501                            $tbls[0],                            $tbls[0],
# Line 235  Line 506 
506                           { stderr => catfile($self->expr_dir, 'problems') });                           { stderr => catfile($self->expr_dir, 'problems') });
507    
508    
509      system_with_redirect([executable_for("remove_multiple_occurring_probes"),      $self->run([executable_for("remove_multiple_occurring_probes"),
510                            $peg_probe_table,                            $peg_probe_table,
511                            $probe_occ_table],                  ],
512                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );                            { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
513    
514      $self->make_missing_probes($peg_probe_table, $probes_fasta,      $self->make_missing_probes($peg_probe_table, $probes_fasta,
515                                 catfile($self->expr_dir, 'probe.no.match'));                                 catfile($self->expr_dir, 'probe.no.match'));
516        $self->make_missing_probes(catfile($self->expr_dir, 'peg.probe.table.no.multiple'), $probes_fasta,
517                                   catfile($self->expr_dir, 'probe.no.multiple.no.match'));
518  }  }
519    
520  sub make_missing_probes  sub make_missing_probes
# Line 297  Line 570 
570      print $fh "library(makecdfenv);\n";      print $fh "library(makecdfenv);\n";
571      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);
572      close($fh);      close($fh);
     system("cat", $tempfile);  
573      system("Rscript", $tempfile);      system("Rscript", $tempfile);
574      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");      system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
575    
576      local($ENV{R_LIBS}) = $libdir;      local($ENV{R_LIBS}) = $libdir;
577      my @cmd = (executable_for("RunRMA"),      $self->run([executable_for("RunRMA"),
578                 "data",                 "data",
579                 catfile($self->expr_dir, "peg.probe.table"),                  catfile($self->expr_dir, "peg.probe.table.no.multiple"),
580                 catfile($self->expr_dir, "probe.no.match"),                  catfile($self->expr_dir, "probe.no.multiple.no.match"),
581                 $expt_dir,                 $expt_dir,
582                 $self->expr_dir);                  $self->expr_dir]);
583    
584      my $rc = system_with_redirect(\@cmd);      my $output = catfile($self->expr_dir, "rma_normalized.tab");
585  #    my $rc = system_with_redirect(\@cmd, { stderr => "rma.stderr" });      if (! -f $output)
586      if ($rc != 0)      {
587            confess("Output file $output was not generated");
588        }
589    }
590    
591    sub compute_rma_normalized_from_sif
592    {
593        my($self, $cdf_file, $sif_file, $expt_dir) = @_;
594    
595        my $my_cdf = catfile($self->expr_dir, "expr.cdf");
596        copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
597    
598        my $my_sif = catfile($self->expr_dir, "expr.sif");
599        copy($sif_file, $my_sif) or confess "Cannot copy $sif_file to $my_sif: $!";
600    
601        #
602        # Create the sif2peg mapping.
603        #
604    
605        my $sif2peg = catfile($self->expr_dir, "sif2peg");
606        if (! -f $sif2peg)
607        {
608            $self->run([executable_for('map_sif_to_pegs'),
609                        $my_sif,
610                        catfile($self->genome_dir, "peg_dna.fasta")],
611                   { stdout => $sif2peg });
612        }
613    
614        $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
615    }
616    
617    sub compute_rma_normalized_from_locus_tags
618    {
619        my($self, $cdf_file, $tag_file, $expt_dir) = @_;
620    
621        my $my_cdf = catfile($self->expr_dir, "expr.cdf");
622        copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
623    
624        my $my_tags = catfile($self->expr_dir, "expr.tags");
625        copy($tag_file, $my_tags) or confess "Cannot copy $tag_file to $my_tags: $!";
626    
627        #
628        # Create the sif2peg mapping.
629        #
630    
631        my $sif2peg = catfile($self->expr_dir, "sif2peg");
632    
633        {
634            my %tags;
635            #
636            # Read the tbl file for the organism and create map from locus tag -> peg.
637            #
638            my $tbl = catfile($self->genome_dir, "Features", "peg", "tbl");
639            if (my $fh = FileHandle->new($tbl, "<"))
640            {
641                while (<$fh>)
642                {
643                    chomp;
644                    my($fid, $loc, @aliases) = split(/\t/);
645                    for my $a (@aliases)
646                    {
647                        $tags{$a} = $fid;
648                        if ($a =~ /^(\S+)\|(.*)/)
649                        {
650                            $tags{$2} = $fid;
651                        }
652                    }
653                }
654                close($fh);
655            }
656            else
657      {      {
658          confess("Error running RMA: rc=$rc cmd=@cmd");              confess "Could not open tbl file $tbl: $!";
659            }
660    
661            my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
662            my $in = FileHandle->new($my_tags, "<") or confess "Cannot open $my_tags: $!";
663    
664            # Per Matt DeJongh: Note that in some cases multiple chip ids
665            # map to the same locus tag; in this case ignore the chip id
666            # that has "_s_" in it.
667    
668            my %locus_map;
669            while (<$in>)
670            {
671                chomp;
672                my($chip_id, $tag) = split(/\t/);
673                next if $tag eq '';
674    
675                if (exists($locus_map{$tag}))
676                {
677                    print "Dup: chip_id=$chip_id tag=$tag  old=$locus_map{$tag}\n";
678                    next if $chip_id =~ /_s_/ ;
679                    next if $chip_id =~ /^Rick/ && $self->genome_id eq '452659.3';
680                }
681                $locus_map{$tag} = $chip_id;
682            }
683            close($in);
684    
685            for my $locus_tag (sort keys %locus_map)
686            {
687                my $chip_id = $locus_map{$locus_tag};
688                my $peg = $tags{$locus_tag};
689                if ($peg ne '')
690                {
691                    print $out "$chip_id\t$peg\n";
692                }
693      }      }
694            close($out);
695        }
696    
697        if (! -s $sif2peg)
698        {
699            confess "No probe to peg mappings were found\n";
700        }
701    
702        $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
703    }
704    
705    sub compute_rma_normalized_from_pegidcorr
706    {
707        my($self, $cdf_file, $corr_file, $expt_dir) = @_;
708    
709        my $my_cdf = catfile($self->expr_dir, "expr.cdf");
710        copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
711    
712        my $my_corr = catfile($self->expr_dir, "peg.id.corr");
713        copy($corr_file, $my_corr) or confess "Cannot copy $corr_file to $my_corr: $!";
714        my $sif2peg = catfile($self->expr_dir, "sif2peg");
715        #
716        # Create the sif2peg mapping.
717        #
718    
719        my $sif2peg = catfile($self->expr_dir, "sif2peg");
720    
721        #
722        # The peg.id.corr table is of the form peg \t chip-id \t something-else
723        # Just rewrite into chip-id \t peg.
724        #
725    
726        my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
727        my $in = FileHandle->new($my_corr, "<") or confess "Cannot open $my_corr: $!";
728        while (<$in>)
729        {
730            chomp;
731            my($peg, $chip_id, undef) = split(/\t/);
732    
733            print $out "$chip_id\t$peg\n";
734        }
735        close($in);
736        close($out);
737    
738        if (! -s $sif2peg)
739        {
740            confess "No probe to peg mappings were found\n";
741        }
742        $self->compute_rma_normalized_using_sif2peg($sif2peg, $expt_dir);
743    }
744    
745    sub compute_rma_normalized_using_sif2peg
746    {
747        my($self, $sif2peg, $expt_dir) = @_;
748    
749        #
750        # We need to build the R library for this cdf.
751        #
752        my($fh, $tempfile) = tempfile();
753    #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
754    
755        my $cdf_path = $self->expr_dir;
756        my $libdir = catfile($self->expr_dir, "r_lib");
757        -d $libdir or mkdir $libdir;
758        my $pkgdir = catfile($self->expr_dir, "r_pkg");
759        -d $pkgdir or mkdir $pkgdir;
760    
761        print $fh "library(makecdfenv);\n";
762        print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
763        close($fh);
764        system("Rscript", $tempfile);
765        system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
766    
767        local($ENV{R_LIBS}) = $libdir;
768        $self->run([executable_for("RunRMA_SIF_format"),
769                    "data",
770                    $expt_dir,
771                    $self->expr_dir]);
772    
773      my $output = catfile($self->expr_dir, "rma_normalized.tab");      my $output = catfile($self->expr_dir, "rma_normalized.tab");
774      if (! -f $output)      if (! -f $output)
775      {      {
# Line 324  Line 779 
779    
780  sub compute_atomic_regulons  sub compute_atomic_regulons
781  {  {
782      my($self) = @_;      my($self, $pearson_cutoff) = @_;
783    
784        $pearson_cutoff ||= 0.7;
785    
786      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
787      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
788      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
789        my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
790        my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
791    
792    
793        $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
794               { stdout => $coreg_clusters });
795    
796        my $genome_ss_dir = $self->genome_dir . "/Subsystems";
797        $self->run([executable_for("make_coreg_conjectures_based_on_subsys"),
798                    $self->expr_dir,
799                    (-d $genome_ss_dir ? $genome_ss_dir : ()),
800                    ],
801               { stdout => $coreg_subsys });
802    
803        $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
804               { stdout => $merged_clusters });
805        $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
806    
807        if (-s $pegs_always_on == 0)
808        {
809            confess "No always-on pegs were found";
810  }  }
811    
812        $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
813               { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
814    
815  1;      $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
816               { stdout => catfile($self->expr_dir, "split.clusters") });
817    
818        $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
819    }
820    
821    sub run
822    {
823        my($self, $cmd, $redirect) = @_;
824    
825        print "Run @$cmd\n";
826        my $rc = system_with_redirect($cmd, $redirect);
827        if ($rc != 0)
828        {
829            confess "Command failed: @$cmd\n";
830        }
831    }
832    
833    sub get_experiment_names
834    {
835        my($self) = @_;
836        my $f = catfile($self->expr_dir, "experiment.names");
837        my $fh;
838        open($fh, "<", $f) or confess "Could not open $f: $!";
839        my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
840        close($fh);
841        return @out;
842    }
843    
844    sub get_experiment_on_off_calls
845    {
846        my($self, $expt) = @_;
847    
848        my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
849        my $fh;
850        open($fh, "<", $f) or confess "Could not open $f: $!";
851        my $names = <$fh>;
852        chomp $names;
853        my @names = split(/\t/, $names);
854        my $idx = 0;
855        my $expt_idx;
856        foreach my $n (@names)
857        {
858            if ($n eq $expt)
859            {
860                $expt_idx = $idx;
861                last;
862            }
863            $idx++;
864        }
865        if (!defined($expt_idx))
866        {
867            confess("Could not find experiment $expt in $f");
868        }
869    
870        my $calls = {};
871        while (<$fh>)
872        {
873            chomp;
874            my($peg, @calls) = split(/\t/);
875            #
876            # +1 because the row[0] element is the peg, and our index is
877            # zero-based.
878            #
879            $calls->{$peg} = $calls[$expt_idx + 1];
880        }
881    
882        close($fh);
883        return($calls);
884    
885    }
886    
887    =head3 save_model_gene_activity
888    
889        $e->save_model_gene_activity($data)
890    
891    Save the results of a modeling run for a given experiment.
892    
893    $data is of the form { experiment_id => $data_hash }
894    
895    =cut
896    
897    sub save_model_gene_activity
898    {
899        my($self, $data) = @_;
900    }
901    
902    sub all_features
903    {
904        my($self, $type) = @_;
905    
906        my @ftypes;
907        my $fdir = catfile($self->genome_dir, "Features");
908        if (defined($type))
909        {
910            @ftypes = ($type);
911        }
912        else
913        {
914            opendir(D, $fdir);
915            @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
916            closedir(D);
917        }
918        my @out;
919        for my $ftype (@ftypes)
920        {
921            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
922            {
923                push(@out, map { /^(\S+)/; $1 } <TBL>);
924                close(TBL);
925            }
926        }
927        return @out;
928    }
929    
930    sub fid_locations
931    {
932        my($self, $fids) = @_;
933    
934        my %fids;
935        $fids{$_}++ for @$fids;
936    
937        my $genome_id = $self->genome_id;
938    
939        my $fdir = catfile($self->genome_dir, "Features");
940        opendir(D, $fdir);
941        my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
942        closedir(D);
943    
944        my $out = {};
945        for my $ftype (@ftypes)
946        {
947            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
948            {
949                while (<TBL>)
950                {
951                    my($id, $locs) = /^(\S+)\t(\S+)\t/;
952    
953                    if ($fids{$id})
954                    {
955                        $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
956                    }
957                }
958                close(TBL);
959            }
960        }
961        return $out;
962    }
963    
964    sub ids_in_subsystems
965    {
966        my($self) = @_;
967    
968        my $dir = $self->genome_dir;
969        my $fh;
970        if (!open($fh, "<", "$dir/Subsystems/bindings"))
971        {
972            warn "No bindings file, falling back to old method\n";
973            return $self->ids_in_subsystems_old();
974        }
975    
976        my $res;
977        while (<$fh>)
978        {
979            chomp;
980            my($ss, $role, $fid) = split(/\t/);
981            $ss =~ s/\s+/_/g;
982            push(@{$res->{$ss}->{$role}}, $fid);
983        }
984        close($fh);
985        return $res;
986    }
987    
988    sub ids_to_subsystems
989    {
990        my($self, $ids) = @_;
991    
992        my $dir = $self->genome_dir;
993        my $fh;
994        if (!open($fh, "<", "$dir/Subsystems/bindings"))
995        {
996            warn "No bindings file, falling back to old method\n";
997            return $self->ids_to_subsystems_old($ids);
998        }
999    
1000        my %ids;
1001        $ids{$_} = 1 for @$ids;
1002    
1003        my $res = {};
1004        while (<$fh>)
1005        {
1006            chomp;
1007            my($ss, $role, $fid) = split(/\t/);
1008            if ($ids{$fid})
1009            {
1010                push(@{$res->{$fid}}, $ss);
1011            }
1012        }
1013        close(SS);
1014    
1015        return $res;
1016    }
1017    
1018    sub ids_in_subsystems_old
1019    {
1020        my($self) = @_;
1021    
1022        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
1023        my $res = {};
1024        while (<SS>)
1025        {
1026            chomp;
1027            my($peg, $ss, $role, $variant) = split(/\t/);
1028            $ss =~ s/\s+/_/g;
1029            push(@{$res->{$ss}->{$role}}, $peg);
1030        }
1031        close(SS);
1032        return $res;
1033    }
1034    
1035    sub ids_to_subsystems_old
1036    {
1037        my($self, $ids) = @_;
1038    
1039        my %ids;
1040        $ids{$_} = 1 for @$ids;
1041    
1042        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
1043        my $res = {};
1044        while (<SS>)
1045        {
1046            chomp;
1047            my($peg, $ss, $role, $variant) = split(/\t/);
1048            if ($ids{$peg})
1049            {
1050                push(@{$res->{$peg}}, $ss);
1051            }
1052        }
1053        close(SS);
1054        return $res;
1055    }
1056    
1057    sub ids_to_functions
1058    {
1059        my($self, $ids) = @_;
1060        open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
1061        my %ids;
1062        $ids{$_} = 1 for @$ids;
1063        my $res = {};
1064    
1065        while (<AF>)
1066        {
1067            chomp;
1068            my($id, $fn) = split(/\t/);
1069            $res->{$id} = $fn if $ids{$id};
1070        }
1071        close(AF);
1072        return $res;
1073    }
1074    
1075    sub best_pearson_corr {
1076        my($self,$pegs1,$cutoff) = @_;
1077    
1078        my @pegs2 = $self->all_features('peg');
1079        my $handle = $self->get_pc_hash_strip($pegs1,\@pegs2);
1080    
1081        my %ok;
1082        my $i;
1083        for ($i=0; ($i < @$pegs1); $i++)
1084        {
1085            foreach my $peg2 ( @pegs2 )
1086            {
1087                my $pc = &pearson_corr($handle,$pegs1->[$i],$peg2);
1088                if (abs($pc >= $cutoff))
1089                {
1090                    $ok{$pegs1->[$i]} -> {$peg2} = $pc;
1091                }
1092            }
1093        }
1094        return \%ok;
1095    }
1096    
1097    sub pearson_corr {
1098        my($hash,$peg1,$peg2) = @_;
1099        my $v = $hash->{$peg1}->{$peg2};
1100        return defined($v) ? sprintf("%0.3f",$v) : " ";
1101    }
1102    
1103    sub get_pc_hash_strip {
1104        my($self,$pegs1,$pegs2) = @_;
1105        my $corrH = $self->get_corr;
1106        my $hash  = &compute_pc_strip($pegs1,$pegs2,$corrH);
1107        return $hash;
1108    }
1109    
1110    sub get_corr {
1111        my($self) = @_;
1112    
1113        my $dir           = $self->expr_dir;
1114        my $rawF          = "$dir/rma_normalized.tab";
1115        my %gene_to_values;
1116        open(RAW,"<$rawF") || die "could not open $rawF";
1117        while (<RAW>)
1118        {
1119            chomp;
1120            my ($gene_id, @gxp_values) = split("\t");
1121            $gene_to_values{$gene_id} = \@gxp_values;
1122        }
1123        close(RAW);
1124        return \%gene_to_values;
1125    }
1126    
1127    sub compute_pc_strip {
1128        my ($pegs1,$pegs2, $gxp_hash) = @_;
1129        my %values = ();
1130    
1131        for (my $i = 0; $i < @$pegs1; $i++)
1132        {
1133            my $stat = Statistics::Descriptive::Full->new();
1134            $stat->add_data(@{$gxp_hash->{$pegs1->[$i]}});
1135    
1136            foreach my $peg2 (@$pegs2)
1137            {
1138                if ($pegs1->[$i] ne $peg2)
1139                {
1140                    my ($q, $m, $r, $err) = $stat->least_squares_fit(@{$gxp_hash->{$peg2}});
1141                    $values{$pegs1->[$i]}->{$peg2} = $r;
1142                }
1143            }
1144        }
1145    
1146        return \%values;
1147    }
1148    
1149    
1150    1;
1151    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3