[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.5, Fri Jan 14 20:06:51 2011 UTC revision 1.6, Thu Jan 20 22:05:57 2011 UTC
# Line 1  Line 1 
1  package ExpressionDir;  package ExpressionDir;
2    
3    use gjoseqlib;
4  use Data::Dumper;  use Data::Dumper;
5  use strict;  use strict;
6  use SeedAware;  use SeedAware;
# Line 15  Line 16 
16  our @probe_parsers = qw(parse_probe_format_1lq  our @probe_parsers = qw(parse_probe_format_1lq
17                          parse_probe_format_1                          parse_probe_format_1
18                          parse_probe_format_2                          parse_probe_format_2
19                          parse_probe_format_3);                          parse_probe_format_3
20                            parse_probe_format_native);
21    
22  =head3 new  =head3 new
23    
# Line 108  Line 110 
110      copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));      copy(catfile($gdir, "contigs"), catfile($self->genome_dir, "contigs"));
111      mkdir(catfile($self->genome_dir, "Features"));      mkdir(catfile($self->genome_dir, "Features"));
112      my @pegs;      my @pegs;
113        my %locs;
114      for my $ftype (qw(peg rna))      for my $ftype (qw(peg rna))
115      {      {
116          my $ofdir = catfile($gdir, "Features", $ftype);          my $ofdir = catfile($gdir, "Features", $ftype);
# Line 119  Line 122 
122              open(NT, ">", "$nfdir/tbl") or confess "Cannot write $nfdir/tbl: $!";              open(NT, ">", "$nfdir/tbl") or confess "Cannot write $nfdir/tbl: $!";
123              while (<OT>)              while (<OT>)
124              {              {
125                  my($id) = /^(\S+)\t/;                  my($id) = /^(\S+)\t(\S+)/;
126                  if (!$fig->is_deleted_fid($id))                  if (!$fig->is_deleted_fid($id))
127                  {                  {
128                      print NT $_;                      print NT $_;
129                      push(@pegs, $id) if $ftype eq 'peg';                      if ($ftype eq 'peg')
130                        {
131                            push(@pegs, $id);
132                            $locs{$id} = $2;
133                        }
134    
135                  }                  }
136              }              }
137              close(OT);              close(OT);
# Line 153  Line 161 
161          print AF join("\t", $peg, $fns->{$peg}), "\n";          print AF join("\t", $peg, $fns->{$peg}), "\n";
162      }      }
163      close(AF);      close(AF);
164    
165        open(PD, ">", catfile($self->genome_dir, "peg_dna.fasta"));
166        for my $peg (@pegs)
167        {
168            my $dna = $fig->dna_seq($self->genome_id, split(/,/, $locs{$peg}));
169            if ($dna eq '')
170            {
171                die "no dna for ", $self->genome_id, " $peg $locs{$peg}\n";
172            }
173            print_alignment_as_fasta(\*PD, [$peg, undef, $dna]);
174        }
175        close(PD);
176  }  }
177    
178  sub create_from_sap  sub create_from_sap
# Line 298  Line 318 
318      return 1;      return 1;
319  }  }
320    
321    sub parse_probe_format_3
322    {
323        my($self, $in_file, $out_file) = @_;
324    
325        my($fh);
326    
327        open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
328        my $out;
329        open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
330    
331        local $/ = "\n>";
332    
333        while (defined($_ = <$fh>))
334        {
335            if ($_ =~ /^>?\S+\s+(\d+)\s+(\d+)[^\n]+\n(\S+)/s)
336            {
337                my $x = $1;
338                my $y = $2;
339                my $seq = $3;
340                if ($seq ne "!")
341                {
342                    if ($seq !~ /^[ACGT]+$/) { die $_ }
343                    if (length($seq) < 15) { print STDERR "BAD: $_"; die "Failed" }
344                    print $out "$x\_$y\t$seq\n";
345                }
346            }
347            else
348            {
349                print STDERR "failed to parse: $_";
350                die "aborted";
351            }
352        }
353    
354        close($out);
355        close($fh);
356        return 1;
357    }
358    
359  #  #
360  # Our "native" format, used for passing through pre-parsed data.  # Our "native" format, used for passing through pre-parsed data.
361  #  #
362  sub parse_probe_format_3  sub parse_probe_format_native
363  {  {
364      my($self, $in_file, $out_file) = @_;      my($self, $in_file, $out_file) = @_;
365    
# Line 481  Line 539 
539      }      }
540  }  }
541    
542    sub compute_rma_normalized_from_sif
543    {
544        my($self, $cdf_file, $sif_file, $expt_dir) = @_;
545    
546        my $my_cdf = catfile($self->expr_dir, "expr.cdf");
547        copy($cdf_file, $my_cdf) or confess "Cannot copy $cdf_file to $my_cdf: $!";
548    
549        my $my_sif = catfile($self->expr_dir, "expr.sif");
550        copy($sif_file, $my_sif) or confess "Cannot copy $sif_file to $my_sif: $!";
551    
552        #
553        # Create the sif2peg mapping.
554        #
555    
556        my $sif2peg = catfile($self->expr_dir, "sif2peg");
557        if (! -f $sif2peg)
558        {
559            $self->run([executable_for('map_sif_to_pegs'),
560                        $my_sif,
561                        catfile($self->genome_dir, "peg_dna.fasta")],
562                   { stdout => $sif2peg });
563        }
564    
565        #
566        # We need to build the R library for this cdf.
567        #
568        my($fh, $tempfile) = tempfile();
569    #m = make.cdf.package("S_aureus.cdf", cdf.path="..",packagename="foo",package.path="/tmp")
570    
571        my $cdf_path = $self->expr_dir;
572        my $libdir = catfile($self->expr_dir, "r_lib");
573        -d $libdir or mkdir $libdir;
574        my $pkgdir = catfile($self->expr_dir, "r_pkg");
575        -d $pkgdir or mkdir $pkgdir;
576    
577        print $fh "library(makecdfenv);\n";
578        print $fh qq(make.cdf.package("expr.cdf", cdf.path="$cdf_path", packagename="datacdf", package.path="$pkgdir", species="genome name");\n);
579        close($fh);
580        system("Rscript", $tempfile);
581        system("R", "CMD", "INSTALL", "-l", $libdir, "$pkgdir/datacdf");
582    
583        local($ENV{R_LIBS}) = $libdir;
584        $self->run([executable_for("RunRMA_SIF_format"),
585                    "data",
586                    $expt_dir,
587                    $self->expr_dir]);
588    
589        my $output = catfile($self->expr_dir, "rma_normalized.tab");
590        if (! -f $output)
591        {
592            confess("Output file $output was not generated");
593        }
594    }
595    
596  sub compute_atomic_regulons  sub compute_atomic_regulons
597  {  {
598      my($self, $pearson_cutoff) = @_;      my($self, $pearson_cutoff) = @_;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3