#!/usr/bin/perl
use strict;
use warnings;
use Pod::Usage;
use Getopt::Long;

# main
main : {
  # options hash table
  my %h_options;

  Getopt::Long::Configure( "no_ignorecase" );
  Getopt::Long::Configure( "bundling" );

  # retrieve all options
  GetOptions (\%h_options,
	      'h|help',
	      'i=i',
	      'a=i',
	      'o|output=s'
	     ) or  pod2usage( -verbose => 1 , noperldoc => 1) ;
	
  # help
  pod2usage( -verbose => 1 , noperldoc => 1) if (exists $h_options{'h'});

  # test ARGV parameters
  pod2usage( -verbose => 2 , noperldoc => 1) if (@ARGV != 1);

  my %hh_x0nm;    # hash 1 = X0; hash 2 = NM
  my %h_x0;       # Number of each X0
  my %h_nm;       # Number of each NM
  my $nb_all = 0;
  my $max_nm = 0;
  my %h_flag;     # Number of each flag
  open(FILE, "$ARGV[0]") || die "Error enable to open $ARGV[0]";
  while (my $line = <FILE>) {
    chomp $line;
    next if($line =~ /^@/);
    $nb_all++;

    my ($x0) = $line =~ /X0:i:(\d+)/;
    my ($nm) = $line =~ /NM:i:(\d+)/;
    if(defined($x0) && defined($nm)) {
      $x0 < 11 ? $hh_x0nm{$x0}{$nm}++ : $hh_x0nm{"11"}{$nm}++;
    }

    if(defined($x0)) {
      $x0 < 11 ? $h_x0{$x0}++ : $h_x0{"11"}++;
    }

    if(defined($nm)) {
      $h_nm{$nm}++;
      $max_nm = $nm if($max_nm < $nm);
    }

    # Count by flag
    my @a_line = split('\t', $line);
    $h_flag{$a_line[1]}++;
  }
  close FILE;

  print "##Global count (by flag)\n";
  print "#Flag\tFlags for the $nb_all sequences\n";
  for my $flag ( sort{$a<=>$b} keys %h_flag ) {
    print "$flag\t$h_flag{$flag}\n";
  }

  print "##NM (edit distance) count for each X0 (number of best hits)\n";
  print "#X0\tALL";
  for my $k_x0 ( sort{$a<=>$b} keys %hh_x0nm ) {
    $k_x0 == 11 ? print "\tX0>=$k_x0" : print "\tX0=$k_x0";
  }
  print "\n";
  print "ALL\t$nb_all";
  for my $k_x0 ( sort{$a<=>$b} keys %h_x0 ) {
    print "\t$h_x0{$k_x0}";
  }
  print "\n";
  my $nm = 0;
  while($nm <= $max_nm) {
    print "NM=$nm\t";
    if(defined $h_nm{$nm}) { print "$h_nm{$nm}" }
    else                   { print "0" }
    for my $k_x0 ( sort{$a<=>$b} keys %hh_x0nm ) {
      if(defined $hh_x0nm{$k_x0}{$nm}) {
	print "\t$hh_x0nm{$k_x0}{$nm}";
      }
      else {
	print "\t0";
      }
    }
    $nm++;
    print "\n";
  }

}


=pod

=head1 NAME

        bwasam2stat.pl - Statistics from SAM file.

=head1 DESCRIPTION

        Global count by flag.
        NM (edit distance) count for each X0 (number of best hits)

        Warning: Input file must be a SAM file created by BWA

=head1 SYNOPSIS

        bwasam2stat.pl FILE.SAM > stat.txt
        OR
        samtools view FILE.BAM | bwasam2stat.pl - > stat.txt

=head1 OPTIONS

        -h,--help

=cut

