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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3