[Bio] / FigKernelPackages / SampleDir.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/SampleDir.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Wed Oct 20 17:20:07 2010 UTC revision 1.2, Tue Oct 26 15:09:35 2010 UTC
# Line 1  Line 1 
1  #  
2  # Utilities for dealing with samples, function summaries,  # Utilities for dealing with samples, function summaries,
3  # OTU summaries, and the vectors that represent them.  # OTU summaries, and the vectors that represent them.
4  #  #
# Line 32  Line 32 
32    
33  use strict;  use strict;
34  use SeedUtils;  use SeedUtils;
35    use gjoseqlib;
36  use ANNOserver;  use ANNOserver;
37  use Data::Dumper;  use Data::Dumper;
38  use YAML::Any qw(Dump Load DumpFile LoadFile);  use YAML::Any qw(Dump Load DumpFile LoadFile);
39    use Time::HiRes 'gettimeofday';
40    
41  use base qw(Class::Accessor);  use base qw(Class::Accessor);
42    
# Line 46  Line 48 
48    
49  sub create  sub create
50  {  {
51      my($class, $dir, $sample_file) = @_;      my($class, $dir, $sample_file, $name, $desc) = @_;
52    
53      if (-d $dir)      if (-d $dir)
54      {      {
# Line 61  Line 63 
63      mkdir($dir) or die "SampleDir::create: mkdir $dir failed: $!";      mkdir($dir) or die "SampleDir::create: mkdir $dir failed: $!";
64      my $seq_type = SeedUtils::validate_fasta_file($sample_file, "$dir/sample.fa");      my $seq_type = SeedUtils::validate_fasta_file($sample_file, "$dir/sample.fa");
65    
66        my $fh;
67        if (open($fh, ">", "$dir/NAME"))
68        {
69            print $fh "$name\n";
70        }
71        close($fh);
72        if (open($fh, ">", "$dir/DESCRIPTION"))
73        {
74            print $fh $desc;
75        }
76        close($fh);
77    
78      my $obj = $class->new($dir);      my $obj = $class->new($dir);
79      $obj->sequence_type($seq_type);      $obj->sequence_type($seq_type);
80    
# Line 77  Line 91 
91          args => \%args,          args => \%args,
92      };      };
93    
94      return bless $self, $class;      bless $self, $class;
95    
96        $self->load_stats();
97    
98        return $self;
99    }
100    
101    sub compute_and_save_stats
102    {
103        my($self) = @_;
104    
105        my $fh;
106        if (!open($fh, "<", "$self->{dir}/sample.fa"))
107        {
108            warn "Cannot open $self->{dir}/sample.fa: $!";
109            return;
110        }
111    
112        my $n = 0;
113        my $tot = 0;
114        my $num_at = 0;
115        my $num_gc = 0;
116        my($max, $min, @lens);
117        while (my($id, $com, $seq) = gjoseqlib::read_next_fasta_seq($fh))
118        {
119            my $len = length($seq);
120            $tot += $len;
121            $n++;
122            push(@lens, $len);
123            $max = $len if !defined($max) || $len > $max;
124            $min = $len if !defined($min) || $len < $min;
125            $num_gc += ($seq =~ tr/gcGC//);
126            $num_at += ($seq =~ tr/atAT//);
127        }
128        my $avg = $tot / $n;
129        my $median = (sort { $a <=> $b } @lens)[int($n / 2)];
130        my $data = {
131            total_size => $tot,
132            count => $n,
133            min => $min,
134            max => $max,
135            mean => $avg,
136            median => $median,
137            gc_content => 100 * (($num_gc + 1) / ($num_gc + $num_at + 2)),
138        };
139    
140        $self->{statistics} = $data;
141    
142        DumpFile("$self->{dir}/statistics.yml", $data);
143    }
144    
145    sub load_stats
146    {
147        my($self) = @_;
148        my $file = "$self->{dir}/statistics.yml";
149    
150        if (-s $file)
151        {
152            $self->{statistics} = LoadFile($file);
153        }
154        else
155        {
156            $self->compute_and_save_stats();
157        }
158    }
159    
160    sub get_statistics
161    {
162        my($self) = @_;
163        return $self->{statistics};
164  }  }
165    
166  sub name  sub name
# Line 91  Line 174 
174  sub description  sub description
175  {  {
176      my($self) = @_;      my($self) = @_;
177      my $n = &SeedUtils::file_head("$self->{dir}/DESCRIPTION",1 );      my $n = &SeedUtils::file_read("$self->{dir}/DESCRIPTION");
     chomp $n;  
178      return $n;      return $n;
179  }  }
180    
# Line 166  Line 248 
248      my $fn_sum_fh;      my $fn_sum_fh;
249      open($fn_sum_fh, ">", "$analysis_dir/sample.fn.sum") or die "Cannot write $analysis_dir/sample.fn.sum: $!";      open($fn_sum_fh, ">", "$analysis_dir/sample.fn.sum") or die "Cannot write $analysis_dir/sample.fn.sum: $!";
250    
251        my $start_time = gettimeofday;
252      my $rh = $self->anno->assign_functions_to_dna(-input => $fh, %params);      my $rh = $self->anno->assign_functions_to_dna(-input => $fh, %params);
253      my $totF = 0;      my $totF = 0;
254      my $totO = 0;      my $totO = 0;
# Line 189  Line 272 
272              $totF++;              $totF++;
273          }          }
274      }      }
275        my $end_time = gettimeofday;
276    
277      close($details_fh);      close($details_fh);
278    
# Line 210  Line 294 
294      }      }
295      close($otu_sum_fh);      close($otu_sum_fh);
296    
297        #
298        # Summary data.
299        #
300        my $sum = {
301            hits_with_function => $totF,
302            hits_with_otu => $totO,
303            distinct_functions => scalar(keys %fn_summary),
304            distinct_otus => scalar(keys %otu_summary),
305            elapsed_time => ($end_time - $start_time),
306        };
307        DumpFile("$analysis_dir/summary.yml", $sum);
308    
309        return ($ds, $analysis_num);
310    }
311    
312    sub get_analysis_list
313    {
314        my($self) = @_;
315    
316        my $out = {};
317        if (opendir(my $dh1, $self->dir))
318        {
319            while (my $f1 = readdir($dh1))
320            {
321                my $p1 = "$self->{dir}/$f1";
322                if (opendir(my $dh2, $p1))
323                {
324                    while (my $f2 = readdir($dh2))
325                    {
326                        my $p2 = "$p1/$f2";
327                        if ($f2 =~ /^\d+$/ && -d $p2)
328                        {
329                            push(@{$out->{$f1}}, $f2);
330                        }
331                    }
332                    closedir($dh2);
333                }
334            }
335            closedir($dh1);
336        }
337        return $out;
338    }
339    
340    sub get_analysis
341    {
342        my($self, $which, $n) = @_;
343        my $dir = "$self->{dir}/$which/$n";
344        if (-d $dir)
345        {
346            return SampleAnalysis->new($self, $which, $n, $dir);
347        }
348        else
349        {
350            return undef;
351        }
352    }
353    
354    package SampleAnalysis;
355    use strict;
356    use YAML::Any qw(Dump Load DumpFile LoadFile);
357    
358    use base 'Class::Accessor';
359    
360    __PACKAGE__->mk_accessors(qw(sample dataset index dir));
361    
362    sub new
363    {
364        my($class, $sample, $dataset, $index, $dir) = @_;
365    
366        my $params = eval { LoadFile("$dir/params.yml"); };
367        my $summary = eval { LoadFile("$dir/summary.yml"); };
368    
369        my $self = {
370            sample => $sample,
371            dataset => $dataset,
372            index => $index,
373            dir => $dir,
374            params => ($params || {}),
375            summary => ($summary || {}),
376        };
377        return bless $self, $class;
378    }
379    
380    sub get_parameters
381    {
382        my($self) = @_;
383        return $self->{params};
384    }
385    
386    sub get_summary
387    {
388        my($self) = @_;
389        return $self->{summary};
390    }
391    
392    sub get_function_file
393    {
394        my($self) = @_;
395        return "$self->{dir}/sample.fn.sum";
396    }
397    
398    sub get_otu_file
399    {
400        my($self) = @_;
401        return "$self->{dir}/sample.otu.sum";
402    }
403    
404    sub get_all_hits_file
405    {
406        my($self) = @_;
407        return "$self->{dir}/sample.out";
408    
409        my $grid = Wx::Grid->new($self, -1);
410    }
411    
412    sub get_function_vector
413    {
414        my($self, $basis) = @_;
415    
416        return $basis->create_vector($self->get_function_file);
417    }
418    
419    sub get_otu_vector
420    {
421        my($self, $basis) = @_;
422    
423        return $basis->create_vector($self->get_otu_file);
424    }
425    
426    package VectorBasis;
427    use strict;
428    use PDL;
429    
430    =head1 DESCRIPTION
431    
432    A vector basis is used to define the indexes in the function or OTU vectors
433    to be used for each function or OTU string. The definition of these is
434    kept on the SEED servers and is retrieved via the get_vector_basis_sets
435    call. The resulting data file contains tab-separated pairs
436    (index, name). The indexes here are 1-based, so we convert them to be
437    0-based.
438    
439    =cut
440    
441    sub new
442    {
443        my($class, $file) = @_;
444    
445        my $fh;
446        if (!open($fh, "<", $file))
447        {
448            die "VectorBasis::new: cannot open $file: $!";
449        }
450    
451        my $by_idx = [];
452        my $by_name = {};
453        my $n = 0;
454        while (<$fh>)
455        {
456            chomp;
457            my($idx, $str) = split(/\t/);
458            $idx--;
459            $by_idx->[$idx] = $str;
460            $by_name->{$str} = $idx;
461            $n++;
462        }
463    
464        my $self = {
465            file => $file,
466            by_index => $by_idx,
467            by_name => $by_name,
468            length => $n,
469        };
470    
471        return bless $self, $class;
472    }
473    
474    =head2 create_vector
475    
476    Create a PDL vector from the given file. It is to be of the format of the
477    function and OTU summary files - tab separated triples (count, fraction, name).
478    
479    =cut
480    
481    sub create_vector
482    {
483        my($self, $file) = @_;
484        my $fh;
485        if (!open($fh, "<", $file))
486        {
487            die "VectorBasis::create_vector: cannot open $file: $!";
488        }
489    
490        my $vec = zeroes($self->{length});
491    
492        while (<$fh>)
493        {
494            chomp;
495            my($count, $frac, $str) = split(/\t/);
496            my $idx = $self->{by_name}->{$str};
497            if (defined($idx))
498            {
499                $vec->index($idx) += $count;
500            }
501        }
502        return $vec;
503  }  }
504    
505    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3