[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.7, Mon Jan 24 18:16:35 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)
     if ($rc != 0)  
573      {      {
574          confess("Error running RMA: rc=$rc cmd=@cmd");          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        #
602        # We need to build the R library for this cdf.
603        #
604        my($fh, $tempfile) = tempfile();
605    #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
606    
607        my $cdf_path = $self->expr_dir;
608        my $libdir = catfile($self->expr_dir, "r_lib");
609        -d $libdir or mkdir $libdir;
610        my $pkgdir = catfile($self->expr_dir, "r_pkg");
611        -d $pkgdir or mkdir $pkgdir;
612    
613        print $fh "library(makecdfenv);\n";
614        print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
615        close($fh);
616        system("Rscript", $tempfile);
617        system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
618    
619        local($ENV{R_LIBS}) = $libdir;
620        $self->run([executable_for("RunRMA_SIF_format"),
621                    "data",
622                    $expt_dir,
623                    $self->expr_dir]);
624    
625        my $output = catfile($self->expr_dir, "rma_normalized.tab");
626        if (! -f $output)
627        {
628            confess("Output file $output was not generated");
629        }
630    }
631    
632    sub compute_rma_normalized_from_locus_tags
633    {
634        my($self, $cdf_file, $tag_file, $expt_dir) = @_;
635    
636        my $my_cdf = catfile($self->expr_dir, "expr.cdf");
637        copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
638    
639        my $my_tags = catfile($self->expr_dir, "expr.tags");
640        copy($tag_file, $my_tags) or confess "Cannot copy $tag_file to $my_tags: $!";
641    
642        #
643        # Create the sif2peg mapping.
644        #
645    
646        my $sif2peg = catfile($self->expr_dir, "sif2peg");
647    
648        {
649            my %tags;
650            #
651            # Read the tbl file for the organism and create map from locus tag -> peg.
652            #
653            my $tbl = catfile($self->genome_dir, "Features", "peg", "tbl");
654            if (my $fh = FileHandle->new($tbl, "<"))
655            {
656                while (<$fh>)
657                {
658                    chomp;
659                    my($fid, $loc, @aliases) = split(/\t/);
660                    for my $a (@aliases)
661                    {
662                        $tags{$a} = $fid;
663                        if ($a =~ /^(\S+)\|(.*)/)
664                        {
665                            $tags{$2} = $fid;
666                        }
667                    }
668                }
669                close($fh);
670            }
671            else
672            {
673                confess "Could not open tbl file $tbl: $!";
674            }
675    
676            my $out = FileHandle->new($sif2peg, ">") or confess "Cannot open $sif2peg for writing: $!";
677            my $in = FileHandle->new($my_tags, "<") or confess "Cannot open $my_tags: $!";
678    
679            # Per Matt DeJongh: Note that in some cases multiple chip ids
680            # map to the same locus tag; in this case ignore the chip id
681            # that has "_s_" in it.
682    
683            my %locus_map;
684            while (<$in>)
685            {
686                chomp;
687                my($chip_id, $tag) = split(/\t/);
688                next if $tag eq '';
689    
690                if (exists($locus_map{$tag}))
691                {
692                    print "Dup: chip_id=$chip_id tag=$tag  old=$locus_map{$tag}\n";
693                    next if $chip_id =~ /_s_/;
694                }
695                $locus_map{$tag} = $chip_id;
696            }
697            close($in);
698    
699            for my $locus_tag (sort keys %locus_map)
700            {
701                my $chip_id = $locus_map{$locus_tag};
702                my $peg = $tags{$locus_tag};
703                if ($peg ne '')
704                {
705                    print $out "$chip_id\t$peg\n";
706                }
707            }
708            close($out);
709        }
710    
711        if (! -s $sif2peg)
712        {
713            confess "No probe to peg mappings were found\n";
714        }
715    
716        #
717        # We need to build the R library for this cdf.
718        #
719        my($fh, $tempfile) = tempfile();
720    #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
721    
722        my $cdf_path = $self->expr_dir;
723        my $libdir = catfile($self->expr_dir, "r_lib");
724        -d $libdir or mkdir $libdir;
725        my $pkgdir = catfile($self->expr_dir, "r_pkg");
726        -d $pkgdir or mkdir $pkgdir;
727    
728        print $fh "library(makecdfenv);\n";
729        print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
730        close($fh);
731        system("Rscript", $tempfile);
732        system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
733    
734        local($ENV{R_LIBS}) = $libdir;
735        $self->run([executable_for("RunRMA_SIF_format"),
736                    "data",
737                    $expt_dir,
738                    $self->expr_dir]);
739    
740      my $output = catfile($self->expr_dir, "rma_normalized.tab");      my $output = catfile($self->expr_dir, "rma_normalized.tab");
741      if (! -f $output)      if (! -f $output)
742      {      {
# Line 324  Line 746 
746    
747  sub compute_atomic_regulons  sub compute_atomic_regulons
748  {  {
749      my($self) = @_;      my($self, $pearson_cutoff) = @_;
750    
751        $pearson_cutoff ||= 0.7;
752    
753      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");      my $coreg_clusters = catfile($self->expr_dir, "coregulated.clusters");
754      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");      my $coreg_subsys = catfile($self->expr_dir, "coregulated.subsys");
755      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");      my $merged_clusters = catfile($self->expr_dir, "merged.clusters");
756        my $probes_always_on = catfile($self->expr_dir, "probes.always.on");
757        my $pegs_always_on = catfile($self->expr_dir, "pegs.always.on");
758    
759        $self->run([executable_for("call_coregulated_clusters_on_chromosome"), $self->expr_dir],
760               { stdout => $coreg_clusters });
761        $self->run([executable_for("make_coreg_conjectures_based_on_subsys"), $self->expr_dir],
762               { stdout => $coreg_subsys });
763    
764        $self->run([executable_for("filter_and_merge_gene_sets"), $self->expr_dir, $coreg_clusters, $coreg_subsys],
765               { stdout => $merged_clusters });
766        $self->run([executable_for("get_ON_probes"), $self->expr_dir, $probes_always_on, $pegs_always_on]);
767    
768        if (-s $probes_always_on == 0)
769        {
770            confess "No always-on probes were found";
771  }  }
772    
773        $self->run([executable_for("Pipeline"), $pegs_always_on, $merged_clusters, $self->expr_dir],
774               { stdout => catfile($self->expr_dir, "comments.by.Pipeline.R") });
775    
776        $self->run([executable_for("SplitGeneSets"), $merged_clusters, $pearson_cutoff, $self->expr_dir],
777               { stdout => catfile($self->expr_dir, "split.clusters") });
778    
779        $self->run([executable_for("compute_atomic_regulons_for_dir"), $self->expr_dir]);
780    }
781    
782    sub run
783    {
784        my($self, $cmd, $redirect) = @_;
785    
786        print "Run @$cmd\n";
787        my $rc = system_with_redirect($cmd, $redirect);
788        if ($rc != 0)
789        {
790            confess "Command failed: @$cmd\n";
791        }
792    }
793    
794    sub get_experiment_names
795    {
796        my($self) = @_;
797        my $f = catfile($self->expr_dir, "experiment.names");
798        my $fh;
799        open($fh, "<", $f) or confess "Could not open $f: $!";
800        my @out = map { chomp; my($num, $name) = split(/\t/); $name } <$fh>;
801        close($fh);
802        return @out;
803    }
804    
805    sub get_experiment_on_off_calls
806    {
807        my($self, $expt) = @_;
808    
809        my $f= catfile($self->expr_dir, "final_on_off_calls.txt");
810        my $fh;
811        open($fh, "<", $f) or confess "Could not open $f: $!";
812        my $names = <$fh>;
813        chomp $names;
814        my @names = split(/\t/, $names);
815        my $idx = 0;
816        my $expt_idx;
817        foreach my $n (@names)
818        {
819            if ($n eq $expt)
820            {
821                $expt_idx = $idx;
822                last;
823            }
824            $idx++;
825        }
826        if (!defined($expt_idx))
827        {
828            confess("Could not find experiment $expt in $f");
829        }
830    
831        my $calls = {};
832        while (<$fh>)
833        {
834            chomp;
835            my($peg, @calls) = split(/\t/);
836            #
837            # +1 because the row[0] element is the peg, and our index is
838            # zero-based.
839            #
840            $calls->{$peg} = $calls[$expt_idx + 1];
841        }
842    
843        close($fh);
844        return($calls);
845    
846    }
847    
848    =head3 save_model_gene_activity
849    
850        $e->save_model_gene_activity($data)
851    
852    Save the results of a modeling run for a given experiment.
853    
854    $data is of the form { experiment_id => $data_hash }
855    
856    =cut
857    
858    sub save_model_gene_activity
859    {
860        my($self, $data) = @_;
861    }
862    
863    sub all_features
864    {
865        my($self, $type) = @_;
866    
867        my @ftypes;
868        my $fdir = catfile($self->genome_dir, "Features");
869        if (defined($type))
870        {
871            @ftypes = ($type);
872        }
873        else
874        {
875            opendir(D, $fdir);
876            @ftypes = grep { -f catfile($fdir, $_) && /^\./ } readdir(D);
877            closedir(D);
878        }
879        my @out;
880        for my $ftype (@ftypes)
881        {
882            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
883            {
884                push(@out, map { /^(\S+)/; $1 } <TBL>);
885                close(TBL);
886            }
887        }
888        return @out;
889    }
890    
891    sub fid_locations
892    {
893        my($self, $fids) = @_;
894    
895        my %fids;
896        $fids{$_}++ for @$fids;
897    
898        my $genome_id = $self->genome_id;
899    
900        my $fdir = catfile($self->genome_dir, "Features");
901        opendir(D, $fdir);
902        my @ftypes = grep { -d catfile($fdir, $_) && ! /^\./ } readdir(D);
903        closedir(D);
904    
905        my $out = {};
906        for my $ftype (@ftypes)
907        {
908            if (open(TBL, "<", catfile($fdir, $ftype, "tbl")))
909            {
910                while (<TBL>)
911                {
912                    my($id, $locs) = /^(\S+)\t(\S+)\t/;
913    
914                    if ($fids{$id})
915                    {
916                        $out->{$id} = "$genome_id:" . SeedUtils::boundary_loc($locs);
917                    }
918                }
919                close(TBL);
920            }
921        }
922        return $out;
923    }
924    
925    sub ids_in_subsystems
926    {
927        my($self) = @_;
928    
929        open(SS, "<", catfile($self->genome_dir, "subsystem.data"));
930        my $res = {};
931        while (<SS>)
932        {
933            chomp;
934            my($peg, $ss, $role, $variant) = split(/\t/);
935            $ss =~ s/\s+/_/g;
936            push(@{$res->{$ss}->{$role}}, $peg);
937        }
938        close(SS);
939        return $res;
940    }
941    
942    sub ids_to_subsystems
943    {
944        my($self, $ids) = @_;
945    
946        my %ids;
947        $ids{$_} = 1 for @$ids;
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            if ($ids{$peg})
956            {
957                push(@{$res->{$peg}}, $ss);
958            }
959        }
960        close(SS);
961        return $res;
962    }
963    
964    sub ids_to_functions
965    {
966        my($self, $ids) = @_;
967        open(AF, "<", catfile($self->genome_dir, "assigned_functions"));
968        my %ids;
969        $ids{$_} = 1 for @$ids;
970        my $res = {};
971    
972        while (<AF>)
973        {
974            chomp;
975            my($id, $fn) = split(/\t/);
976            $res->{$id} = $fn if $ids{$id};
977        }
978        close(AF);
979        return $res;
980    }
981    
982  1;  1;
983    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3