File Coverage

blib/lib/Bio/ViennaNGS/BamStat.pm
Criterion Covered Total %
statement 4 6 66.6
branch n/a
condition n/a
subroutine 2 2 100.0
pod n/a
total 6 8 75.0


line stmt bran cond sub pod time code
1             # -*-CPerl-*-
2             # Last changed Time-stamp: <2015-02-09 09:17:24 mtw>
3              
4             package Bio::ViennaNGS::BamStat;
5              
6 1     1   1179 use version; our $VERSION = qv('0.12_16');
  1         2  
  1         6  
7 1     1   243 use Bio::DB::Sam 1.37;
  0            
  0            
8             use Moose;
9             use Carp;
10             use Data::Dumper;
11             use namespace::autoclean;
12              
13             has 'bam' => (
14             is => 'rw',
15             isa => 'Str',
16             required => 1,
17             predicate => 'has_bam',
18             );
19              
20             has 'data' => (
21             is => 'ro',
22             isa => 'HashRef',
23             predicate => 'has_data',
24             );
25              
26             has 'control_match' => ( # provides stats how many mapped bases match the reference genome
27             is => 'rw',
28             isa => 'Bool',
29             default => '1',
30             clearer => 'clear_control_match',
31             predicate => 'has_control_match',
32             );
33              
34             has 'control_clip' => ( # provides stats how many bases are soft or hard clipped
35             is => 'rw',
36             isa => 'Bool',
37             default => '1',
38             clearer => 'clear_control_clip',
39             predicate => 'has_control_clip',
40             );
41              
42             has 'control_split' => ( # provides stats how many/often mapped reads are split
43             is => 'rw',
44             isa => 'Bool',
45             default => '1',
46             clearer => 'clear_control_split',
47             predicate => 'has_control_split',
48             );
49              
50             has 'control_qual' => ( # provides stats on quality of the match
51             is => 'rw',
52             isa => 'Bool',
53             default => '1',
54             clearer => 'clear_control_qual',
55             predicate => 'has_control_qual',
56             );
57              
58             has 'control_edit' => ( # provides stats on the edit distance between read and mapped reference
59             is => 'rw',
60             isa => 'Bool',
61             default => '1',
62             clearer => 'clear_control_edit',
63             predicate => 'has_control_edit',
64             );
65              
66             has 'control_flag' => ( # analyses the sam bit flag for qual/strands/pair_vs_single reads
67             is => 'rw',
68             isa => 'Bool',
69             default => '1',
70             clearer => 'clear_control_flag',
71             predicate => 'has_control_flag',
72             );
73              
74             has 'control_score' => ( # provides stats on per-base quality scores
75             is => 'rw',
76             isa => 'Bool',
77             default => '1',
78             clearer => 'clear_control_score',
79             predicate => 'has_control_score',
80             );
81              
82             has 'control_uniq' => ( # gives number and stats of multiplicity of readaligments
83             is => 'rw',
84             isa => 'Bool',
85             default => '1',
86             clearer => 'clear_control_uniq',
87             predicate => 'has_control_uniq',
88             );
89              
90             has 'is_segemehl' => ( # toggles to consider segemehl specific bam feature
91             is => 'rw',
92             isa => 'Bool',
93             default => '0',
94             predicate => 'has_is_segemehl',
95             );
96              
97             has 'data_out' => (
98             is => 'rw',
99             isa => 'HashRef',
100             default => sub { {} },
101             predicate => 'has_data_out',
102             );
103              
104             has 'data_match' => (
105             is => 'rw',
106             isa => 'ArrayRef',
107             default => sub { [] },
108             predicate => 'has_data_match',
109             );
110              
111             has 'data_qual' => (
112             is => 'rw',
113             isa => 'ArrayRef',
114             default => sub { [] },
115             predicate => 'has_data_qual',
116             );
117              
118             has 'data_edit' => (
119             is => 'rw',
120             isa => 'ArrayRef',
121             default => sub { [] },
122             predicate => 'has_data_edit',
123             );
124              
125             has 'data_score' => (
126             is => 'rw',
127             isa => 'ArrayRef',
128             default => sub { [] },
129             predicate => 'has_data_score',
130             );
131              
132             has 'data_clip' => (
133             is => 'rw',
134             isa => 'ArrayRef',
135             default => sub { [] },
136             predicate => 'has_data_clip',
137             );
138              
139             has 'data_split' => (
140             is => 'rw',
141             isa => 'HashRef',
142             default => sub { {} },
143             predicate => 'has_data_split',
144             );
145              
146             has 'data_flag' => (
147             is => 'rw',
148             isa => 'HashRef',
149             default => sub { {} },
150             predicate => 'has_data_flag',
151             );
152              
153             has 'data_uniq' => (
154             is => 'rw',
155             isa => 'HashRef',
156             default => sub { {} },
157             predicate => 'has_data_uniq',
158             );
159              
160             sub stat_singleBam {
161             my ($self) = @_;
162             my $this_function = (caller(0))[3];
163             confess "ERROR [$this_function] Attribute 'bam' not found $!"
164             unless ($self->has_bam);
165              
166             ## read in bam file using Bio::DB::Sam
167             #my $sam = Bio::DB::Sam->new(-bam => $self->bam);
168             my $bam = Bio::DB::Bam->open($self->bam);
169              
170             my $header = $bam->header;
171             my $target_count = $header->n_targets;
172             my @chromosomes = @{$header->target_name};
173             my $target_names = \@chromosomes;
174              
175             while (my $align = $bam->read1) {
176             my %flags = ();
177             my $qname = $align->qname;
178             my $seqid = $target_names->[$align->tid];
179             my $start = $align->pos+1;
180             my $end = $align->calend;
181             my $strand = $align->strand;
182             my $cigar = $align->cigar_str;
183             my @scores = $align->qscore;
184             my $match_qual = $align->qual;
185             my $flag = $align->flag;
186             my $attributes = $align->aux;
187             # print "\$seqid: $seqid\t \$qname: $qname\t $start: $start\t \$end: $end\t \$cigar: $cigar\t\$strand: $strand\t \$match_qual: $match_qual\t \$attributes: $attributes\n";
188             my @tags = $align->get_all_tags;
189             # print ">> tags: @tags\n";
190              
191             #############################
192             ### flag
193             if ($self->has_control_flag){
194              
195             my $multisplit_weight = 1;
196             if ($self->has_is_segemehl && $align->has_tag('XL') ) {
197             $multisplit_weight = (1/($align->aux_get('XL')));
198             }
199              
200             if( $align->get_tag_values('PAIRED')) {
201             ${$self->data_flag}{'paired'}->{'paired-end'}+=$multisplit_weight;
202             }else{
203             ${$self->data_flag}{'paired'}->{'single-end'}+=$multisplit_weight;
204             }
205              
206             if( $align->get_tag_values('REVERSED')) {
207             ${$self->data_flag}{'strand'}->{'reverse'}+=$multisplit_weight;
208             }else{
209             ${$self->data_flag}{'strand'}->{'forward'}+=$multisplit_weight;
210             }
211              
212             if( $align->get_tag_values('PAIRED') && $align->get_tag_values('M_UNMAPPED')){
213             ${$self->data_flag}{'pairs'}->{'unmapped_pair'}+=$multisplit_weight;
214             }
215             elsif( $align->get_tag_values('PAIRED') && !$align->get_tag_values('M_UNMAPPED')) {
216             ${$self->data_flag}{'pairs'}->{'mapped_pair'}+=$multisplit_weight;
217             }
218              
219             if( $align->get_tag_values('DUPLICATE') ){
220             ${$self->data_flag}{'unmapped'}->{'duplicated'}+=$multisplit_weight;
221             }
222             if( $align->get_tag_values('QC_FAILED') ){
223             ${$self->data_flag}{'unmapped'}->{'qual_failed'}+=$multisplit_weight;
224             }
225             if( $align->get_tag_values('UNMAPPED') ){
226             ${$self->data_flag}{'unmapped'}->{'unmapped'}+=$multisplit_weight;
227             }
228             }
229              
230             #############################
231             ### uniq
232             if ($self->has_control_uniq){
233              
234             my $multisplit_weight = 1;
235             if ($self->has_is_segemehl && $align->has_tag('XL') ) {
236             $multisplit_weight = ($align->aux_get('XL'));
237             }
238              
239             my $multimap_weight = 1;
240             if ( $align->has_tag('NH') ) {
241             $multimap_weight = ($align->aux_get('NH'));
242             }
243              
244             ${$self->data_uniq}{'uniq'}+=(1/$multisplit_weight) if($multimap_weight == 1);
245             ${$self->data_uniq}{'mapped'}+=1/($multimap_weight*$multisplit_weight);
246             ${$self->data_uniq}{'multiplicity'}->{$multimap_weight}+=(1/$multisplit_weight);
247              
248             }
249              
250             #############################
251             ### Quality Score read
252             if($self->has_control_qual){
253             if($self->has_is_segemehl && $match_qual == 255){
254             #carp "WARN [$this_function] no match_qual for segemehl";
255             $self->clear_control_qual;
256             }
257             elsif($match_qual){
258             push @{$self->data_qual}, $match_qual;
259             }
260             else{
261             carp "WARN [$this_function] No matchqual for read $qname";
262             carp "WARN [$this_function] Setting \$self->has_control_qual to zero";
263             $self->clear_control_qual;
264             }
265             }
266              
267             #############################
268             #### Edit distance
269             if($self->has_control_edit){
270             if( $align->has_tag('NM') ){
271             push @{$self->data_edit}, $align->aux_get('NM');
272             }
273             else{
274             carp "WARN [$this_function] No <NM> attribute for read $qname";
275             carp "WARN [$this_function] Setting \$self->has_control_edit to zero";
276             $self->clear_control_edit;
277             }
278             }
279              
280             #############################
281             ### Match bases
282             if($self->has_control_match){
283             if($cigar){
284             push @{$self->data_match},
285             sprintf("%.2f", 100*(&cigarmatch($cigar))/(&cigarlength($cigar)));
286             }
287             else{
288             carp "WARN [$this_function] No CIGAR string for read $qname";
289             carp "WARN [$this_function] Setting \$self->has_control_match to zero";
290             $self->clear_control_match;
291             }
292             }
293              
294             #############################
295             ### Clip reads
296             if($self->control_clip){
297             if($cigar){
298             push @{$self->data_clip}, sprintf("%.2f", 100*(&cigarclip($cigar)));
299             }
300             else{
301             carp "WARN [$this_function] No CIGAR string for read $qname";
302             carp "WARN [$this_function] Setting \$self->has_control_clip to zero";
303             $self->clear_control_clip;
304             }
305             }
306              
307             #############################
308             ### Split reads
309             if($self->control_split){
310             if ($self->has_is_segemehl && $align->has_tag('XL') ) {
311             my $split_counts = $align->aux_get('XL');
312             ${$self->data_split}{$split_counts}+=1/($split_counts);
313             ${$self->data_split}{'total'}+=1/($split_counts);
314             }
315             elsif ( $self->has_is_segemehl ){ }
316             else{
317             if($cigar){
318             my $split_counts = cigarsplit($cigar);
319             ${$self->data_split}{$split_counts}++;
320             ${$self->split_data}{'total'}++;
321             }
322             else{
323             carp "WARN [$this_function] No CIGAR string for read $qname";
324             carp "WARN [$this_function] Setting \$self->has_control_split to zero";
325             $self->clear_control_split;
326             }
327             }
328             }
329              
330             } # end while
331              
332              
333             #############################
334             ### Extract reference genome/chromosome sizes from BAM header
335             for (my $i=0;$i<$header->n_targets;$i++){
336             my $seqid = $header->target_name->[$i];
337             my $length = $header->target_len->[$i];
338             ${$self->data_out}{'chrom'}->{$seqid}=$length;
339             }
340              
341             #############################
342             #### uniq
343             if ($self->has_control_uniq){
344             ${$self->data_out}{'uniq'}->{'uniq_mapped_reads'}=${$self->data_uniq}{'uniq'};
345             ${$self->data_out}{'uniq'}->{'mapped_reads'}=sprintf("%d", ${$self->data_uniq}{'mapped'});
346             ${$self->data_out}{'uniq'}->{'uniq_mapped_reads_percent'}=
347             sprintf("%.2f", (100*${$self->data_uniq}{'uniq'}/${$self->data_uniq}{'mapped'}));
348             foreach my $multiplicity (sort {$a cmp $b} keys %{${$self->data_uniq}{'multiplicity'}}) {
349             ${$self->data_out}{'uniq'}->{'distribution_of_multimapper'}->{"<$multiplicity>"}=
350             (${$self->data_uniq}{'multiplicity'}->{$multiplicity}/$multiplicity);
351             }
352             }
353              
354             #############################
355             ##### quality stats
356             if($self->has_control_qual){
357             my %qual_stats=%{stats(@${$self->data_qual})};
358             if ($qual_stats{'min'} == 255 && $qual_stats{'max'} == 255){
359             carp "WARN [$this_function] No matchqual (all values set to 255)";
360             carp "WARN [$this_function] Setting \$self->has_control_qual to zero";
361             $self->clear_control_qual;
362             }
363             else{
364             ${$self->data_out}{'qual'}->{'min'} = $qual_stats{'min'};
365             ${$self->data_out}{'qual'}->{'1q'} = $qual_stats{'1q'};
366             ${$self->data_out}{'qual'}->{'mean'} = $qual_stats{'mean'};
367             ${$self->data_out}{'qual'}->{'med'} = $qual_stats{'med'};
368             ${$self->data_out}{'qual'}->{'3q'} = $qual_stats{'3q'};
369             ${$self->data_out}{'qual'}->{'max'} = $qual_stats{'max'};
370             }
371             }
372              
373             #############################
374             ### Clip data
375             if($self->has_control_clip){
376             my %clip_stats=%{stats(@{$self->data_clip})};
377              
378             ${$self->data_out}{'clip'}->{'min'} = $clip_stats{'min'};
379             ${$self->data_out}{'clip'}->{'1q'} = $clip_stats{'1q'};
380             ${$self->data_out}{'clip'}->{'mean'} = $clip_stats{'mean'};
381             ${$self->data_out}{'clip'}->{'med'} = $clip_stats{'med'};
382             ${$self->data_out}{'clip'}->{'3q'} = $clip_stats{'3q'};
383             ${$self->data_out}{'clip'}->{'max'} = $clip_stats{'max'};
384             }
385              
386             #############################
387             ##### Match data
388             if($self->has_control_match){
389             my %match_stats=%{stats(@{$self->data_match})};
390              
391             ${$self->data_out}{'match'}->{'min'} = $match_stats{'min'};
392             ${$self->data_out}{'match'}->{'1q'} = $match_stats{'1q'};
393             ${$self->data_out}{'match'}->{'mean'} = $match_stats{'mean'};
394             ${$self->data_out}{'match'}->{'med'} = $match_stats{'med'};
395             ${$self->data_out}{'match'}->{'3q'} = $match_stats{'3q'};
396             ${$self->data_out}{'match'}->{'max'} = $match_stats{'max'};
397             }
398              
399             #############################
400             ##### Edit distance
401             if($self->has_control_edit){
402             my %edit_stats=%{stats(@{$self->data_edit})};
403              
404             ${$self->data_out}{'edit'}->{'min'} = $edit_stats{'min'};
405             ${$self->data_out}{'edit'}->{'1q'} = $edit_stats{'1q'};
406             ${$self->data_out}{'edit'}->{'mean'} = $edit_stats{'mean'};
407             ${$self->data_out}{'edit'}->{'med'} = $edit_stats{'med'};
408             ${$self->data_out}{'edit'}->{'3q'} = $edit_stats{'3q'};
409             ${$self->data_out}{'edit'}->{'max'} = $edit_stats{'max'};
410             }
411              
412             #############################
413             ###### flags
414             if ($self->has_control_flag){
415              
416             ${$self->data_flag}{'pairs'}->{'mapped_pair'} = 0
417             unless ( defined(${$self->data_flag}{'pairs'}->{'mapped_pair'}) );
418             ${$self->data_flag}{'pairs'}->{'unmapped_pair'} = 0
419             unless ( defined(${$self->data_flag}{'pairs'}->{'unmapped_pair'}) );
420              
421             ${$self->data_flag}{'paired'}->{'paired-end'} = 0
422             unless ( defined(${$self->data_flag}{'paired'}->{'paired-end'}) );
423             ${$self->data_flag}{'paired'}->{'single-end'} = 0
424             unless ( defined(${$self->data_flag}{'paired'}->{'single-end'}) );
425              
426             my $total_mapped = ${$self->data_flag}{'paired'}->{'paired-end'} +
427             ${$self->data_flag}{'paired'}->{'single-end'};
428              
429             ${$self->data_out}{'aln_count'}->{'mapped_pair'} =
430             ${$self->data_flag}{'pairs'}->{'mapped_pair'};
431             ${$self->data_out}{'aln_count'}->{'unmapped_pair'} =
432             ${$self->data_flag}{'pairs'}->{'unmapped_pair'};
433             ${$self->data_out}{'aln_count'}->{'mapped_single'} =
434             ${$self->data_flag}{'paired'}->{'single-end'};
435             ${$self->data_out}{'aln_count'}->{'total'} = $total_mapped;
436              
437             ${$self->data_flag}{'strand'}->{'forward'} = 0
438             unless ( defined(${$self->data_flag}{'strand'}->{'forward'}) );
439             ${$self->data_flag}{'strand'}->{'reverse'} = 0
440             unless ( defined(${$self->data_flag}{'strand'}->{'reverse'}) );
441              
442             my $total_strands = ${$self->data_flag}{'strand'}->{'forward'} +
443             ${$self->data_flag}{'strand'}->{'reverse'};
444              
445             ${$self->data_out}{'strand'}->{'forward'} = ${$self->data_flag}{'strand'}->{'forward'};
446             ${$self->data_out}{'strand'}->{'reverse'} = ${$self->data_flag}{'strand'}->{'reverse'};
447             }
448              
449             #############################
450             ###### split
451             if ($self->has_control_split){
452             foreach my $splits (sort {$a cmp $b} keys %{$self->data_split}) {
453             my $split_counts=( defined(${$self->data_split}{$splits}) )?(${$self->data_split}{$splits}):(0);
454             ${$self->data_out}{'split'}->{'distribution_of_multisplit'}->{"$splits"}=$split_counts;
455             }
456             }
457              
458             }
459              
460             __PACKAGE__->meta->make_immutable;
461              
462              
463             no Moose;
464              
465             sub stats{
466             # usage: %h = %{stats(@a)};
467             my @vals = sort {$a <=> $b} @_;
468             my %stats = ();
469             my $median = '';
470            
471             if(@vals%2){$stats{'med'} = $vals[int(@vals/2)];} #odd median
472             else{$stats{'med'} = ($vals[int(@vals/2)-1] + $vals[int(@vals/2)])/2;} #even median
473             $stats{'med'} = sprintf("%.2f", $stats{'med'});
474             $stats{'mean'} = sprintf("%.2f", &mean(\@vals)); ## mean
475             $stats{'min'} = sprintf("%.2f", &min(\@vals)); ## min
476             $stats{'max'} = sprintf("%.2f", &max(\@vals)); ## max
477             $stats{'1q'} = sprintf("%.2f", $vals[int(@vals/4)]); ## 1.quartile
478             $stats{'3q'} = sprintf("%.2f", $vals[int((@vals*3)/4)]); ## 3.quartile
479              
480             return(\%stats);
481             }
482              
483             sub mean { # usage: $h = %{mean(\@a)};
484             my ($arrayref) = @_;
485             my $sum;
486             foreach (@$arrayref) {$sum += $_}
487             return $sum / @$arrayref;
488             }
489              
490             sub max { # usage: $h = %{max(\@a)};
491             my ($arrayref) = @_;
492             my $max = $arrayref->[0];
493             foreach (@$arrayref) {$max = $_ if $_ > $max}
494             return $max;
495             }
496              
497             sub min { # usage: $h = %{min(\@a)};
498             my ($arrayref) = @_;
499             my $min = $arrayref->[0];
500             foreach (@$arrayref) {$min = $_ if $_ < $min}
501             return $min;
502             }
503              
504             sub cigarlength { #usage: &cigarlength($cigarstring)
505             my $cigar_string = shift;
506             my $cigar_length = 0;
507             while($cigar_string=~m/(\d+)[MIX=]/g){$cigar_length+=$1}
508             return($cigar_length);
509             }
510              
511             sub cigarmatch { #usage: cigarmatch($cigarstring)
512             my $cigar_string = shift;
513             my $cigar_match = 0;
514             while($cigar_string=~m/(\d+)[M=]/g){$cigar_match+=$1}
515             return($cigar_match);
516             }
517              
518             sub cigarsplit { #usage: cigarsplit($cigarstring)
519             my $cigar_string = shift;
520             my $cigar_split = 0;
521             while($cigar_string=~m/N/g){$cigar_split+=1}
522             return($cigar_split);
523             }
524              
525             sub cigarclip { #usage: cigarclip($cigarstring)
526             my $cigar_string = shift;
527             my $cigar_length = 0;
528             my $cigar_clip = 0;
529             while($cigar_string=~m/(\d+)[SH]/g) {$cigar_clip+=$1}
530             while($cigar_string=~m/(\d+)\D/g){$cigar_length+=$1}
531             return ($cigar_length==0)?(0):($cigar_clip/$cigar_length);
532             }
533              
534              
535              
536             1;
537             __END__
538              
539              
540             =head1 NAME
541              
542             Bio::ViennaNGS::BamStat - Moose interface to BAM mapping statistics
543              
544             =head1 SYNOPSIS
545              
546             use Bio::ViennaNGS::BamStat;
547              
548             my $bamsummary = Bio::ViennaNGS::BamStatSummary->new(files => "path/to/file.bam",
549             outpath => "path/to/output.directory",
550             rlib => "path/to/R_installation",
551             is_segemehl => 0|1,
552             control_match => 0|1,
553             control_clip => 0|1,
554             control_split => 0|1,
555             control_qual => 0|1,
556             control_edit => 0|1,
557             control_flag => 0|1,
558             control_score => 0|1,
559             );
560              
561             =head1 DESCRIPTION
562              
563             This module provides a L<Moose> interface to the mapping statistics of
564             a single BAM file. It builds on L<Bio::DB::Sam>.
565              
566             =head1 DEPENDENCIES
567              
568             =over
569              
570             =item L<Bio::DB::Sam> >= 1.37
571              
572             =item L<Moose>
573              
574             =back
575              
576             =head1 SEE ALSO
577              
578             =over
579              
580             =item L<Bio::ViennaNGS>
581              
582             =item L<Bio::ViennaNGS::BamStatSummary>
583              
584             =back
585              
586             =head1 AUTHORS
587              
588             =over
589              
590             =item Fabian Amman E<lt>fabian@tbi.univie.ac.atE<gt>
591              
592             =item Michael T. Wolfinger E<lt>michael@wolfinger.euE<gt>
593              
594             =back
595              
596             =head1 COPYRIGHT AND LICENSE
597              
598             Copyright (C) 2015 by Michael T. Wolfinger
599              
600             This library is free software; you can redistribute it and/or modify
601             it under the same terms as Perl itself, either Perl version 5.10.0 or,
602             at your option, any later version of Perl 5 you may have available.
603              
604             This software is distributed in the hope that it will be useful, but
605             WITHOUT ANY WARRANTY; without even the implied warranty of
606             MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
607              
608             =cut