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

# main
main : {
  # options hash table
  my %h_options;
  $h_options{'n'} = 2;
  $h_options{'0'} = 4;
  $h_options{'1'} = 50;

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

  # retrieve all options
  GetOptions (\%h_options,
	      'h|help',
	      'n|NM=i',
	      '0|X0=i',
	      '1|X1=i',
	      'I|INDEL',
	      'o|onlybest'
	     ) 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);

  #Info log
  my $all       = 0;
  my $all_0_16  = 0;
  my $dueto_n   = 0;
  my $dueto_0   = 0;
  my $dueto_1   = 0;
  my $dueto_I   = 0;
  my $dueto_XAmissing = 0;
  my $valid     = 0;

  open(FILE, "$ARGV[0]") || die "Error enable to open $ARGV[0]";
  while (my $line = <FILE>) {
    chomp $line;
    my @a_line = split('\t', $line);

    if   ($line =~ /^@/) {
      print "$line\n";
    }
    elsif($a_line[1] != 0 && $a_line[1] != 16) {
      $all++;
      $all_0_16++;
      next;
    }
    else {
      $all++;
      my ($x0) = $line =~ /X0:i:(\d+)/;
      my ($x1) = $line =~ /X1:i:(\d+)/;
      my ($nm) = $line =~ /NM:i:(\d+)/;

      my $ok = 1;
      if($nm > $h_options{'n'})                   { $dueto_n++; $ok=0; }
      if(defined $x0 && $x0 > $h_options{'0'})    { $dueto_0++; $ok=0; }
      if(defined $x1 && $x1 > $h_options{'1'})    { $dueto_1++; $ok=0; }
      if($a_line[5]!~/^\d+M$/ && $h_options{'I'}) { $dueto_I++; $ok=0; }

      # Remove from XA tag suboptimal hits
      if($h_options{'o'}) {
	my ($xa) = $line =~ /XA:Z:(.+)/;
	if(defined $xa) {
	  my @a_xa = split(';', $xa);
	  my $xafilter = "";
	  foreach my $alt ( @a_xa ) {
	    my ($altnm) = $alt =~ /,(\d+)$/;
	    $xafilter .= "$alt;" if($altnm <= $nm);
	  }
	  if($xafilter ne "") { $line =~ s/XA:Z:.+$/XA:Z:$xafilter/; }
	  else                { $line =~ s/\s+XA:Z:.+$//;               }
	}
      }

      # Valid XA tag... due to -n arg of bwa samse
      my ($xa) = $line =~ /XA:Z:(.+)/;
      if($x0!=1 && !defined($xa)) { $dueto_XAmissing++; $ok=0; }

      if($ok) {
	print "$line\n";
	$valid++;
      }
    }
  }
  close FILE;
  print STDERR "##Filter the alignment\n";
  print STDERR "#Label\t \n";
  print STDERR "ALL READS\t$all\n";
  print STDERR "Flag is not 0 or 16\t$all_0_16\n";
  print STDERR "NM>$h_options{'n'}\t$dueto_n\n";
  print STDERR "X0>$h_options{'0'}\t$dueto_0\n";
  print STDERR "X1>$h_options{'1'}\t$dueto_1\n";
  print STDERR "INDEL\t$dueto_I\n";
  print STDERR "XA missing\t$dueto_XAmissing\n";
  print STDERR "VALID READS\t$valid\n";
}


=pod

=head1 NAME

        bwasamXfilter.pl - Filter SAM file using NM X0 and X1 tags.

=head1 DESCRIPTION

        From SAM alignment file provided by BWA samse build two files:
          - filtered SAM file (STDOUT)
          - log file (STDERR)

        Sequence from SAM input file can be filter by:
          - max NM   (-n)
          - max X0   (-0)
          - max X1   (-1)
          - No INDEL (-I)

        Warning: apdapt the -n option of bwa samse...

=head1 SYNOPSIS

        samtools view -h FILE.BAM | ( bwasamXfilter.pl - > OUTPUT.SAM ) >& OUTPUT.LOG
        OR
        ( bwasamXfilter.pl - -n 3 -0 5 -I -o > OUTPUT.SAM ) >& OUTPUT.LOG

=head1 OPTIONS

        -h,--help

        -n, --NM INT
          max NM value [2].
          Sequences aligned with more than 2 for NM tag will be filtered.

        -0, --X0 INT
          max X0 value [4].
          Sequences aligned with more than 4 for X0 tag will be filtered.

        -1, --X1 INT
          max X1 value [50].
          Sequences aligned with more than 50 for X1 tag will be filtered.

        -I, --INDEL
          Filter alignment with INDEL.

        -o, --onlybest
          removed from XA tag suboptimal hit.
	
=cut

