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

Annotation of /FigKernelPackages/SampleDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.2
2 : olson 1.1 # Utilities for dealing with samples, function summaries,
3 :     # OTU summaries, and the vectors that represent them.
4 :     #
5 :    
6 :     #
7 :     # A sample is stored in a directory.
8 :     #
9 :     # The raw DNA for the sample is in sample.fa
10 :     # Metadata about the sample is stored in a set of flat files:
11 :     #
12 :     # NAME
13 :     # DESCRIPTION
14 :     #
15 :     # Each analysis is stored in a directory
16 :     #
17 :     # dataset-name/idx
18 :     #
19 :     # where dataset-name is the name of the Kmer dataset used for the analysis,
20 :     # and idx is a numeric index. Each analysis with a different set of
21 :     # parameters is stored in its own directory.
22 :     #
23 :     # The processing parameters are stored in a YAML file params.yml.
24 :     # The parameters file contains a serialized hash with the following keys:
25 :     # kmer: kmer size
26 :     # max_gap: max-gap parameter
27 :     # min_size: min-size parameter
28 :     # dataset: Kmer dataset identifier
29 :     #
30 :    
31 :     package SampleDir;
32 :    
33 :     use strict;
34 :     use SeedUtils;
35 : olson 1.2 use gjoseqlib;
36 : olson 1.1 use ANNOserver;
37 :     use Data::Dumper;
38 :     use YAML::Any qw(Dump Load DumpFile LoadFile);
39 : olson 1.2 use Time::HiRes 'gettimeofday';
40 : olson 1.1
41 :     use base qw(Class::Accessor);
42 :    
43 :     __PACKAGE__->mk_accessors(qw(dir sequence_type anno));
44 :    
45 :     our %kmer_defaults = (-kmer => 8,
46 :     -minHits => 3,
47 :     -maxGap => 600);
48 :    
49 :     sub create
50 :     {
51 : olson 1.2 my($class, $dir, $sample_file, $name, $desc) = @_;
52 : olson 1.1
53 :     if (-d $dir)
54 :     {
55 :     die "SampleDir::create: directorty $dir already exists\n";
56 :     }
57 :    
58 :     if (! -f $sample_file)
59 :     {
60 :     die "SampleDir::create: sample file $sample_file does not exist\n";
61 :     }
62 :    
63 :     mkdir($dir) or die "SampleDir::create: mkdir $dir failed: $!";
64 :     my $seq_type = SeedUtils::validate_fasta_file($sample_file, "$dir/sample.fa");
65 :    
66 : olson 1.2 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 : olson 1.1 my $obj = $class->new($dir);
79 :     $obj->sequence_type($seq_type);
80 :    
81 :     return $obj;
82 :     }
83 :    
84 :     sub new
85 :     {
86 :     my($class, $dir, %args ) = @_;
87 :    
88 :     my $self = {
89 :     dir => $dir,
90 :     anno => ANNOserver->new(exists($args{-url}) ? (-url => $args{-url}) : ()),
91 :     args => \%args,
92 :     };
93 :    
94 : olson 1.2 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 : olson 1.1 }
165 :    
166 :     sub name
167 :     {
168 :     my($self) = @_;
169 :     my $n = &SeedUtils::file_head("$self->{dir}/NAME",1 );
170 :     chomp $n;
171 :     return $n;
172 :     }
173 :    
174 :     sub description
175 :     {
176 :     my($self) = @_;
177 : olson 1.2 my $n = &SeedUtils::file_read("$self->{dir}/DESCRIPTION");
178 : olson 1.1 return $n;
179 :     }
180 :    
181 :     sub perform_basic_analysis
182 :     {
183 :     my($self, %args) = @_;
184 :    
185 :     my %params = %kmer_defaults;
186 :     for my $k (keys %args)
187 :     {
188 :     if (defined($params{$k}))
189 :     {
190 :     $params{$k} = $args{$k};
191 :     }
192 :     }
193 :    
194 :     #
195 :     # If dataset not specified, look up the current
196 :     # default and specify it with all future calls.
197 :     #
198 :     if (!defined($params{-kmerDataset}))
199 :     {
200 :     my $ds = $self->anno->get_dataset();
201 :     $params{-kmerDataset} = $ds;
202 :     print "Using default dataset $ds\n";
203 :     }
204 :    
205 :     #
206 :     # Create the dataset directory if not already present.
207 :     #
208 :    
209 :     my $ds = $params{-kmerDataset};
210 :     my $ds_dir = "$self->{dir}/$ds";
211 :     if (!-d $ds_dir)
212 :     {
213 :     mkdir($ds_dir) or die "Cannot mkdir $ds_dir: $!";
214 :     }
215 :    
216 :     #
217 :     # Find a new analysis dir pathname. Start at zero.
218 :     #
219 :     my $analysis_num = 0;
220 :     my $analysis_dir = "$ds_dir/$analysis_num";
221 :     while (-e $analysis_dir)
222 :     {
223 :     $analysis_num++;
224 :     $analysis_dir = "$ds_dir/$analysis_num";
225 :     }
226 :    
227 :     mkdir($analysis_dir) or die "Cannot mkdir $analysis_dir: $!";
228 :    
229 :     print STDERR "Using analysis dir $analysis_dir\n";
230 :    
231 :     DumpFile("$analysis_dir/params.yml", \%params);
232 :    
233 :     my $fh;
234 :     my $sample = "$self->{dir}/sample.fa";
235 :     if (!open($fh, "<", $sample))
236 :     {
237 :     die "Cannot open sample file $sample: $!";
238 :     }
239 :     my %otu_summary;
240 :     my %fn_summary;
241 :    
242 :     my $details_fh;
243 :     open($details_fh, ">", "$analysis_dir/sample.out") or die "Cannot write $analysis_dir/sample.out: $!";
244 :    
245 :     my $otu_sum_fh;
246 :     open($otu_sum_fh, ">", "$analysis_dir/sample.otu.sum") or die "Cannot write $analysis_dir/sample.otu.sum: $!";
247 :    
248 :     my $fn_sum_fh;
249 :     open($fn_sum_fh, ">", "$analysis_dir/sample.fn.sum") or die "Cannot write $analysis_dir/sample.fn.sum: $!";
250 :    
251 : olson 1.2 my $start_time = gettimeofday;
252 : olson 1.1 my $rh = $self->anno->assign_functions_to_dna(-input => $fh, %params);
253 :     my $totF = 0;
254 :     my $totO = 0;
255 :     while (my $res = $rh->get_next())
256 :     {
257 :     my($id, $tuple) = @$res;
258 :    
259 :     my($count, $start, $stop, $func, $otu) = @$tuple;
260 :    
261 :     my $loc = "${id}_${start}_${stop}";
262 :     print $details_fh "$id\t$count\t$loc\t$func\t$otu\n";
263 :    
264 :     if (defined($otu))
265 :     {
266 :     $otu_summary{$otu}++;
267 :     $totO++;
268 :     }
269 :     if (defined($func))
270 :     {
271 :     $fn_summary{$func}++;
272 :     $totF++;
273 :     }
274 :     }
275 : olson 1.2 my $end_time = gettimeofday;
276 : olson 1.1
277 :     close($details_fh);
278 :    
279 :     for my $fn (sort { $fn_summary{$b} <=> $fn_summary{$a} } keys %fn_summary)
280 :     {
281 :     print $fn_sum_fh join("\t",
282 :     $fn_summary{$fn},
283 :     sprintf("%0.6f", $fn_summary{$fn} / $totF),
284 :     $fn), "\n";
285 :     }
286 :     close($fn_sum_fh);
287 :    
288 :     for my $otu (sort { $otu_summary{$b} <=> $otu_summary{$a} } keys %otu_summary)
289 :     {
290 :     print $otu_sum_fh join("\t",
291 :     $otu_summary{$otu},
292 :     sprintf("%0.6f", $otu_summary{$otu} / $totO),
293 :     $otu), "\n";
294 :     }
295 :     close($otu_sum_fh);
296 :    
297 : olson 1.2 #
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 : olson 1.1 }
385 : olson 1.2
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 : olson 1.1
476 : olson 1.2 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 : olson 1.1
506 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3