[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.6, Thu Jan 20 22:05:57 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;  use gjoseqlib;
5  use Data::Dumper;  use Data::Dumper;
6  use strict;  use strict;
# Line 17  Line 18 
18                          parse_probe_format_1                          parse_probe_format_1
19                          parse_probe_format_2                          parse_probe_format_2
20                          parse_probe_format_3                          parse_probe_format_3
21                            parse_probe_format_shew
22                          parse_probe_format_native);                          parse_probe_format_native);
23    
24  =head3 new  =head3 new
# Line 347  Line 349 
349          else          else
350          {          {
351              print STDERR "failed to parse: $_";              print STDERR "failed to parse: $_";
352              die "aborted";              return undef;
353          }          }
354      }      }
355    
# Line 357  Line 359 
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_native  sub parse_probe_format_native
# Line 400  Line 436 
436  {  {
437      my($self, $probes) = @_;      my($self, $probes) = @_;
438    
439      my($probe_suffix) = $probes =~ /(\.[^.]+)$/;      my($probe_suffix) = $probes =~ m,(\.[^/.]+)$,;
440    
441      my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");      my $my_probes = catfile($self->expr_dir, "probes.in$probe_suffix");
442    
# Line 593  Line 629 
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");
741        if (! -f $output)
742        {
743            confess("Output file $output was not generated");
744        }
745    }
746    
747  sub compute_atomic_regulons  sub compute_atomic_regulons
748  {  {
749      my($self, $pearson_cutoff) = @_;      my($self, $pearson_cutoff) = @_;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3