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

Annotation of /FigKernelPackages/SampleDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (view) (download) (as text)

1 : olson 1.2
2 : olson 1.3 #
3 :     # This is a SAS component.
4 :     #
5 :    
6 :    
7 : olson 1.1 # Utilities for dealing with samples, function summaries,
8 :     # OTU summaries, and the vectors that represent them.
9 :     #
10 :    
11 :     #
12 :     # A sample is stored in a directory.
13 :     #
14 :     # The raw DNA for the sample is in sample.fa
15 :     # Metadata about the sample is stored in a set of flat files:
16 :     #
17 :     # NAME
18 :     # DESCRIPTION
19 :     #
20 :     # Each analysis is stored in a directory
21 :     #
22 :     # dataset-name/idx
23 :     #
24 :     # where dataset-name is the name of the Kmer dataset used for the analysis,
25 :     # and idx is a numeric index. Each analysis with a different set of
26 :     # parameters is stored in its own directory.
27 :     #
28 :     # The processing parameters are stored in a YAML file params.yml.
29 :     # The parameters file contains a serialized hash with the following keys:
30 :     # kmer: kmer size
31 :     # max_gap: max-gap parameter
32 :     # min_size: min-size parameter
33 :     # dataset: Kmer dataset identifier
34 :     #
35 :    
36 :     package SampleDir;
37 :    
38 :     use strict;
39 :     use SeedUtils;
40 : olson 1.2 use gjoseqlib;
41 : olson 1.1 use ANNOserver;
42 :     use Data::Dumper;
43 :     use YAML::Any qw(Dump Load DumpFile LoadFile);
44 : olson 1.2 use Time::HiRes 'gettimeofday';
45 : olson 1.1
46 :     use base qw(Class::Accessor);
47 :    
48 :     __PACKAGE__->mk_accessors(qw(dir sequence_type anno));
49 :    
50 :     our %kmer_defaults = (-kmer => 8,
51 :     -minHits => 3,
52 :     -maxGap => 600);
53 :    
54 :     sub create
55 :     {
56 : olson 1.2 my($class, $dir, $sample_file, $name, $desc) = @_;
57 : olson 1.1
58 :     if (-d $dir)
59 :     {
60 :     die "SampleDir::create: directorty $dir already exists\n";
61 :     }
62 :    
63 :     if (! -f $sample_file)
64 :     {
65 :     die "SampleDir::create: sample file $sample_file does not exist\n";
66 :     }
67 :    
68 :     mkdir($dir) or die "SampleDir::create: mkdir $dir failed: $!";
69 :     my $seq_type = SeedUtils::validate_fasta_file($sample_file, "$dir/sample.fa");
70 :    
71 : olson 1.2 my $fh;
72 :     if (open($fh, ">", "$dir/NAME"))
73 :     {
74 :     print $fh "$name\n";
75 :     }
76 :     close($fh);
77 :     if (open($fh, ">", "$dir/DESCRIPTION"))
78 :     {
79 :     print $fh $desc;
80 :     }
81 :     close($fh);
82 :    
83 : olson 1.1 my $obj = $class->new($dir);
84 :     $obj->sequence_type($seq_type);
85 :    
86 :     return $obj;
87 :     }
88 :    
89 :     sub new
90 :     {
91 :     my($class, $dir, %args ) = @_;
92 :    
93 :     my $self = {
94 :     dir => $dir,
95 :     anno => ANNOserver->new(exists($args{-url}) ? (-url => $args{-url}) : ()),
96 :     args => \%args,
97 :     };
98 :    
99 : olson 1.2 bless $self, $class;
100 :    
101 :     $self->load_stats();
102 :    
103 :     return $self;
104 :     }
105 :    
106 :     sub compute_and_save_stats
107 :     {
108 :     my($self) = @_;
109 :    
110 :     my $fh;
111 :     if (!open($fh, "<", "$self->{dir}/sample.fa"))
112 :     {
113 :     warn "Cannot open $self->{dir}/sample.fa: $!";
114 :     return;
115 :     }
116 :    
117 :     my $n = 0;
118 :     my $tot = 0;
119 :     my $num_at = 0;
120 :     my $num_gc = 0;
121 :     my($max, $min, @lens);
122 :     while (my($id, $com, $seq) = gjoseqlib::read_next_fasta_seq($fh))
123 :     {
124 :     my $len = length($seq);
125 :     $tot += $len;
126 :     $n++;
127 :     push(@lens, $len);
128 :     $max = $len if !defined($max) || $len > $max;
129 :     $min = $len if !defined($min) || $len < $min;
130 :     $num_gc += ($seq =~ tr/gcGC//);
131 :     $num_at += ($seq =~ tr/atAT//);
132 :     }
133 :     my $avg = $tot / $n;
134 :     my $median = (sort { $a <=> $b } @lens)[int($n / 2)];
135 :     my $data = {
136 :     total_size => $tot,
137 :     count => $n,
138 :     min => $min,
139 :     max => $max,
140 :     mean => $avg,
141 :     median => $median,
142 :     gc_content => 100 * (($num_gc + 1) / ($num_gc + $num_at + 2)),
143 :     };
144 :    
145 :     $self->{statistics} = $data;
146 :    
147 :     DumpFile("$self->{dir}/statistics.yml", $data);
148 :     }
149 :    
150 :     sub load_stats
151 :     {
152 :     my($self) = @_;
153 :     my $file = "$self->{dir}/statistics.yml";
154 :    
155 :     if (-s $file)
156 :     {
157 :     $self->{statistics} = LoadFile($file);
158 :     }
159 :     else
160 :     {
161 :     $self->compute_and_save_stats();
162 :     }
163 :     }
164 :    
165 :     sub get_statistics
166 :     {
167 :     my($self) = @_;
168 :     return $self->{statistics};
169 : olson 1.1 }
170 :    
171 :     sub name
172 :     {
173 :     my($self) = @_;
174 :     my $n = &SeedUtils::file_head("$self->{dir}/NAME",1 );
175 :     chomp $n;
176 :     return $n;
177 :     }
178 :    
179 :     sub description
180 :     {
181 :     my($self) = @_;
182 : olson 1.2 my $n = &SeedUtils::file_read("$self->{dir}/DESCRIPTION");
183 : olson 1.1 return $n;
184 :     }
185 :    
186 :     sub perform_basic_analysis
187 :     {
188 :     my($self, %args) = @_;
189 :    
190 :     my %params = %kmer_defaults;
191 :     for my $k (keys %args)
192 :     {
193 :     if (defined($params{$k}))
194 :     {
195 :     $params{$k} = $args{$k};
196 :     }
197 :     }
198 :    
199 :     #
200 :     # If dataset not specified, look up the current
201 :     # default and specify it with all future calls.
202 :     #
203 :     if (!defined($params{-kmerDataset}))
204 :     {
205 :     my $ds = $self->anno->get_dataset();
206 :     $params{-kmerDataset} = $ds;
207 :     print "Using default dataset $ds\n";
208 :     }
209 :    
210 :     #
211 :     # Create the dataset directory if not already present.
212 :     #
213 :    
214 :     my $ds = $params{-kmerDataset};
215 :     my $ds_dir = "$self->{dir}/$ds";
216 :     if (!-d $ds_dir)
217 :     {
218 :     mkdir($ds_dir) or die "Cannot mkdir $ds_dir: $!";
219 :     }
220 :    
221 :     #
222 :     # Find a new analysis dir pathname. Start at zero.
223 :     #
224 :     my $analysis_num = 0;
225 :     my $analysis_dir = "$ds_dir/$analysis_num";
226 :     while (-e $analysis_dir)
227 :     {
228 :     $analysis_num++;
229 :     $analysis_dir = "$ds_dir/$analysis_num";
230 :     }
231 :    
232 :     mkdir($analysis_dir) or die "Cannot mkdir $analysis_dir: $!";
233 :    
234 :     print STDERR "Using analysis dir $analysis_dir\n";
235 :    
236 :     DumpFile("$analysis_dir/params.yml", \%params);
237 :    
238 :     my $fh;
239 :     my $sample = "$self->{dir}/sample.fa";
240 :     if (!open($fh, "<", $sample))
241 :     {
242 :     die "Cannot open sample file $sample: $!";
243 :     }
244 :     my %otu_summary;
245 :     my %fn_summary;
246 :    
247 :     my $details_fh;
248 :     open($details_fh, ">", "$analysis_dir/sample.out") or die "Cannot write $analysis_dir/sample.out: $!";
249 :    
250 :     my $otu_sum_fh;
251 :     open($otu_sum_fh, ">", "$analysis_dir/sample.otu.sum") or die "Cannot write $analysis_dir/sample.otu.sum: $!";
252 :    
253 :     my $fn_sum_fh;
254 :     open($fn_sum_fh, ">", "$analysis_dir/sample.fn.sum") or die "Cannot write $analysis_dir/sample.fn.sum: $!";
255 :    
256 : olson 1.2 my $start_time = gettimeofday;
257 : olson 1.1 my $rh = $self->anno->assign_functions_to_dna(-input => $fh, %params);
258 :     my $totF = 0;
259 :     my $totO = 0;
260 :     while (my $res = $rh->get_next())
261 :     {
262 :     my($id, $tuple) = @$res;
263 :    
264 :     my($count, $start, $stop, $func, $otu) = @$tuple;
265 :    
266 :     my $loc = "${id}_${start}_${stop}";
267 :     print $details_fh "$id\t$count\t$loc\t$func\t$otu\n";
268 :    
269 :     if (defined($otu))
270 :     {
271 :     $otu_summary{$otu}++;
272 :     $totO++;
273 :     }
274 :     if (defined($func))
275 :     {
276 :     $fn_summary{$func}++;
277 :     $totF++;
278 :     }
279 :     }
280 : olson 1.2 my $end_time = gettimeofday;
281 : olson 1.1
282 :     close($details_fh);
283 :    
284 :     for my $fn (sort { $fn_summary{$b} <=> $fn_summary{$a} } keys %fn_summary)
285 :     {
286 :     print $fn_sum_fh join("\t",
287 :     $fn_summary{$fn},
288 :     sprintf("%0.6f", $fn_summary{$fn} / $totF),
289 :     $fn), "\n";
290 :     }
291 :     close($fn_sum_fh);
292 :    
293 :     for my $otu (sort { $otu_summary{$b} <=> $otu_summary{$a} } keys %otu_summary)
294 :     {
295 :     print $otu_sum_fh join("\t",
296 :     $otu_summary{$otu},
297 :     sprintf("%0.6f", $otu_summary{$otu} / $totO),
298 :     $otu), "\n";
299 :     }
300 :     close($otu_sum_fh);
301 :    
302 : olson 1.2 #
303 :     # Summary data.
304 :     #
305 :     my $sum = {
306 :     hits_with_function => $totF,
307 :     hits_with_otu => $totO,
308 :     distinct_functions => scalar(keys %fn_summary),
309 :     distinct_otus => scalar(keys %otu_summary),
310 :     elapsed_time => ($end_time - $start_time),
311 :     };
312 :     DumpFile("$analysis_dir/summary.yml", $sum);
313 :    
314 :     return ($ds, $analysis_num);
315 :     }
316 :    
317 :     sub get_analysis_list
318 :     {
319 :     my($self) = @_;
320 :    
321 :     my $out = {};
322 :     if (opendir(my $dh1, $self->dir))
323 :     {
324 :     while (my $f1 = readdir($dh1))
325 :     {
326 :     my $p1 = "$self->{dir}/$f1";
327 :     if (opendir(my $dh2, $p1))
328 :     {
329 :     while (my $f2 = readdir($dh2))
330 :     {
331 :     my $p2 = "$p1/$f2";
332 :     if ($f2 =~ /^\d+$/ && -d $p2)
333 :     {
334 :     push(@{$out->{$f1}}, $f2);
335 :     }
336 :     }
337 :     closedir($dh2);
338 :     }
339 :     }
340 :     closedir($dh1);
341 :     }
342 :     return $out;
343 :     }
344 :    
345 :     sub get_analysis
346 :     {
347 :     my($self, $which, $n) = @_;
348 :     my $dir = "$self->{dir}/$which/$n";
349 :     if (-d $dir)
350 :     {
351 :     return SampleAnalysis->new($self, $which, $n, $dir);
352 :     }
353 :     else
354 :     {
355 :     return undef;
356 :     }
357 :     }
358 :    
359 :     package SampleAnalysis;
360 :     use strict;
361 :     use YAML::Any qw(Dump Load DumpFile LoadFile);
362 : olson 1.4 use POSIX;
363 : olson 1.2
364 :     use base 'Class::Accessor';
365 :    
366 :     __PACKAGE__->mk_accessors(qw(sample dataset index dir));
367 :    
368 :     sub new
369 :     {
370 :     my($class, $sample, $dataset, $index, $dir) = @_;
371 :    
372 : olson 1.4 my $params_file = "$dir/params.yml";
373 :     my $summary_file = "$dir/summary.yml";
374 :    
375 :     my $params = eval { LoadFile($params_file); };
376 :     my $summary = eval { LoadFile($summary_file); };
377 : olson 1.2
378 :     my $self = {
379 :     sample => $sample,
380 :     dataset => $dataset,
381 :     index => $index,
382 :     dir => $dir,
383 :     params => ($params || {}),
384 :     summary => ($summary || {}),
385 : olson 1.4 params_file => $params_file,
386 :     summary_file => $summary_file,
387 : olson 1.2 };
388 :     return bless $self, $class;
389 :     }
390 :    
391 :     sub get_parameters
392 :     {
393 :     my($self) = @_;
394 :     return $self->{params};
395 : olson 1.1 }
396 : olson 1.2
397 : olson 1.4 sub save_parameters
398 :     {
399 :     my($self) = @_;
400 :     my $now = strftime("%Y-%m-%d-%H-%M-%S", localtime);
401 :     my $bak = "$self->{params_file}.bak.$now";
402 :     if (!rename($self->{params_file}, $bak))
403 :     {
404 :     warn "Could not rename $self->{params_file} to $bak: $!";
405 :     }
406 :     DumpFile($self->{params_file}, $self->{params});
407 :     }
408 :    
409 : olson 1.2 sub get_summary
410 :     {
411 :     my($self) = @_;
412 :     return $self->{summary};
413 :     }
414 :    
415 :     sub get_function_file
416 :     {
417 :     my($self) = @_;
418 :     return "$self->{dir}/sample.fn.sum";
419 :     }
420 :    
421 :     sub get_otu_file
422 :     {
423 :     my($self) = @_;
424 :     return "$self->{dir}/sample.otu.sum";
425 :     }
426 :    
427 :     sub get_all_hits_file
428 :     {
429 :     my($self) = @_;
430 :     return "$self->{dir}/sample.out";
431 :    
432 :     my $grid = Wx::Grid->new($self, -1);
433 :     }
434 :    
435 :     sub get_function_vector
436 :     {
437 :     my($self, $basis) = @_;
438 :    
439 :     return $basis->create_vector($self->get_function_file);
440 :     }
441 :    
442 :     sub get_otu_vector
443 :     {
444 :     my($self, $basis) = @_;
445 :    
446 :     return $basis->create_vector($self->get_otu_file);
447 :     }
448 :    
449 : olson 1.4 #
450 :     # Recreate the summarization of data based on the
451 :     # current function exclusion list.
452 :     #
453 :     sub rerun_analysis
454 :     {
455 :     my($self) = @_;
456 :    
457 :     my $analysis_dir = $self->dir;
458 :    
459 :     my $raw_fh;
460 :     if (!open($raw_fh, "<", "$analysis_dir/sample.out"))
461 :     {
462 :     warn "Cannot open $analysis_dir/sample.out: $!";
463 :     return;
464 :     }
465 :    
466 :     my $otu_sum_fh;
467 :     open($otu_sum_fh, ">", "$analysis_dir/sample.otu.sum") or die "Cannot write $analysis_dir/sample.otu.sum: $!";
468 :    
469 :     my $fn_sum_fh;
470 :     open($fn_sum_fh, ">", "$analysis_dir/sample.fn.sum") or die "Cannot write $analysis_dir/sample.fn.sum: $!";
471 :    
472 :     my $totF = 0;
473 :     my $totO = 0;
474 :     my $excluded = 0;
475 :    
476 :     my %otu_summary;
477 :     my %fn_summary;
478 :    
479 :     my $exclusions = $self->{params}->{excluded_functions};
480 :     while (<$raw_fh>)
481 :     {
482 :     chomp;
483 :     my($id, $count, $loc, $func, $otu) = split(/\t/);
484 :    
485 :     if ($exclusions->{$func})
486 :     {
487 :     $excluded++;
488 :     next;
489 :     }
490 :    
491 :     if ($otu ne '')
492 :     {
493 :     $otu_summary{$otu}++;
494 :     $totO++;
495 :     }
496 :     if ($func ne '')
497 :     {
498 :     $fn_summary{$func}++;
499 :     $totF++;
500 :     }
501 :     }
502 :     close($raw_fh);
503 :    
504 :     for my $fn (sort { $fn_summary{$b} <=> $fn_summary{$a} } keys %fn_summary)
505 :     {
506 :     print $fn_sum_fh join("\t",
507 :     $fn_summary{$fn},
508 :     sprintf("%0.6f", $fn_summary{$fn} / $totF),
509 :     $fn), "\n";
510 :     }
511 :     close($fn_sum_fh);
512 :    
513 :     for my $otu (sort { $otu_summary{$b} <=> $otu_summary{$a} } keys %otu_summary)
514 :     {
515 :     print $otu_sum_fh join("\t",
516 :     $otu_summary{$otu},
517 :     sprintf("%0.6f", $otu_summary{$otu} / $totO),
518 :     $otu), "\n";
519 :     }
520 :     close($otu_sum_fh);
521 :    
522 :     #
523 :     # Summary data.
524 :     #
525 :     my $sum = {};
526 :     %$sum = %{$self->{summary}};
527 :    
528 :     my $sum2 = {
529 :     hits_with_function => $totF,
530 :     hits_with_otu => $totO,
531 :     distinct_functions => scalar(keys %fn_summary),
532 :     distinct_otus => scalar(keys %otu_summary),
533 :     excluded_hits => $excluded,
534 :     };
535 :     $sum->{$_} = $sum2->{$_} for keys %$sum2;
536 :    
537 :     $self->{summary} = $sum;
538 :    
539 :     DumpFile("$analysis_dir/summary.yml", $sum);
540 :     }
541 :    
542 : olson 1.2 package VectorBasis;
543 :     use strict;
544 :     use PDL;
545 :    
546 :     =head1 DESCRIPTION
547 :    
548 :     A vector basis is used to define the indexes in the function or OTU vectors
549 :     to be used for each function or OTU string. The definition of these is
550 :     kept on the SEED servers and is retrieved via the get_vector_basis_sets
551 :     call. The resulting data file contains tab-separated pairs
552 :     (index, name). The indexes here are 1-based, so we convert them to be
553 :     0-based.
554 :    
555 :     =cut
556 :    
557 :     sub new
558 :     {
559 :     my($class, $file) = @_;
560 :    
561 :     my $fh;
562 :     if (!open($fh, "<", $file))
563 :     {
564 :     die "VectorBasis::new: cannot open $file: $!";
565 :     }
566 :    
567 :     my $by_idx = [];
568 :     my $by_name = {};
569 :     my $n = 0;
570 :     while (<$fh>)
571 :     {
572 :     chomp;
573 :     my($idx, $str) = split(/\t/);
574 :     $idx--;
575 :     $by_idx->[$idx] = $str;
576 :     $by_name->{$str} = $idx;
577 : olson 1.4 $n = $idx if $idx > $n;
578 : olson 1.2 }
579 : olson 1.4 $n++;
580 : olson 1.2
581 :     my $self = {
582 :     file => $file,
583 :     by_index => $by_idx,
584 :     by_name => $by_name,
585 :     length => $n,
586 :     };
587 :    
588 :     return bless $self, $class;
589 :     }
590 :    
591 :     =head2 create_vector
592 : olson 1.1
593 : olson 1.2 Create a PDL vector from the given file. It is to be of the format of the
594 :     function and OTU summary files - tab separated triples (count, fraction, name).
595 :    
596 :     =cut
597 :    
598 :     sub create_vector
599 :     {
600 :     my($self, $file) = @_;
601 :     my $fh;
602 :     if (!open($fh, "<", $file))
603 :     {
604 :     die "VectorBasis::create_vector: cannot open $file: $!";
605 :     }
606 :    
607 :     my $vec = zeroes($self->{length});
608 :    
609 :     while (<$fh>)
610 :     {
611 :     chomp;
612 :     my($count, $frac, $str) = split(/\t/);
613 :     my $idx = $self->{by_name}->{$str};
614 :     if (defined($idx))
615 :     {
616 :     $vec->index($idx) += $count;
617 :     }
618 :     }
619 :     return $vec;
620 :     }
621 :    
622 : olson 1.1
623 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3