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

# main
main : {
  # options hash table
  my %h_options;
  $h_options{'i'} = 2;
  $h_options{'a'} = 5;

  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);

  # output must be defined
  pod2usage( -verbose => 1 , noperldoc => 1) if (!defined $h_options{'o'});

  my $allseq = 0;     # Total number of sequences
  my $filtseq= 0;     # Number of filtered sequences
  my $n      = 0;     # Number of sequence with at least one N
  my $n_f    = 0;     # Number of filtered sequence with at least one N
  my %h_inXindiv;     # How many seq suported by X indiv. (key=X, value=seqnumber)
  my %h_inXindivf;    # How many seq suported by X indiv. after filtering (key=X, value=seqnumber)
  my %h_inXcopy;      # How many seq suported by X copy in all indiv (key=X, value=copynumber)
  my %h_inXcopyf;     # How many seq suported by X copy in all indiv after fitering (key=X, value=copynumber)
  my %h_histlen;      # How many seq by length (key=seqLen, value=seqnumber)

  my %h_id2f;         # Column index to filename
  my %hh_histlenbyind;# How many seq by indiv by length (key1=indiv, key2=seqLen, value=seqnumber)
  my %h_histlenf;     # How many seq by length after filtering (key=seqLen, value=seqnumber)

  my %h_spe;          # By fastq, number of specific sequence
  my %h_spef;         # By fastq, number of specific sequence after filtering

  print STDERR "Matrix ($ARGV[0]) scan started:\n";

  open(STAT, ">$h_options{'o'}.filter.txt") || die "Error enable to create $h_options{'o'}.filter.txt";
  open(TFA,  ">$h_options{'o'}.filter.tfa") || die "Error enable to create $h_options{'o'}.filter.tfa";
  open(MAT,  ">$h_options{'o'}.filter.csv") || die "Error enable to create $h_options{'o'}.filter.csv";

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

    #Header
    if($a_line[0]=~/^#/) {
      print MAT "#seq\tSupportedBy\tTotal";
      for(my $i=1; $i<=$#a_line; $i++) {
	$h_id2f{$i} = $a_line[$i];
	$hh_histlenbyind{$a_line[$i]} = {};
	print MAT "\t$a_line[$i]";
      }
      print MAT "\n";
      next;
    }

    $allseq++;
    $n++ if($a_line[0]=~/N/);
    my $nbindiv = 0;
    my $nbcopy  = 0;
    my $whichfastq = "";
    for (my $i=1; $i<=$#a_line; $i++) {
      if($a_line[$i]>0) {
	$nbindiv++;
	$hh_histlenbyind{$h_id2f{$i}}{length($a_line[0])}++;
	$whichfastq = $h_id2f{$i};
      }
      $nbcopy += $a_line[$i];
    }
    $h_inXindiv{$nbindiv}++;
    $h_inXcopy{$nbcopy}++;
    $h_histlen{length($a_line[0])}++;
    $h_spe{$whichfastq}++ if($nbindiv == 1);

    if($nbindiv >= $h_options{'i'} || $nbcopy >= $h_options{'a'}) {
      $h_spef{$whichfastq}++ if($nbindiv == 1);
      $h_inXindivf{$nbindiv}++;
      $h_inXcopyf{$nbcopy}++;
      $h_histlenf{length($a_line[0])}++;
      $filtseq++;
      $n_f++ if($a_line[0]=~/N/);
      $line =~ s/^[ATCGN]+\t//;
      #$line =~ s/\t/_/g;
      print TFA ">seq$filtseq#$nbindiv#$nbcopy\n$a_line[0]\n";
      print MAT "seq$filtseq#$nbindiv#$nbcopy\t$nbindiv\t$nbcopy\t$line\n";
    }

    print STDERR "$allseq...done\n" if(($allseq%50000) == 0)
  }
  close FILE;

  print STDERR "End of the analysis ($filtseq pass filter - ".($filtseq/$allseq*100)."%)\n";


  print STAT "##Global count\n";
  print STAT "#Step\t \n";
  print STAT "ALL_uniq\t$allseq\n";
  print STAT "Pass_filter\t$filtseq\n";
  print STAT "ALL_uniq_N\t$n\n";
  print STAT "Pass_filter_N\t$n_f\n";


  print STAT "##Number of sequences supported by X fastq\n";
  print STAT "#IndivNumber";
  foreach my $k (sort { $a <=> $b} keys %h_inXindiv) {
    print STAT "\t$k";
  }
  print STAT "\nBefore filtering";
  foreach my $k (sort { $a <=> $b} keys %h_inXindiv) {
    print STAT "\t$h_inXindiv{$k}";
  }
  print STAT "\nAfter filtering";
  foreach my $k (sort { $a <=> $b} keys %h_inXindivf) {
    print STAT "\t$h_inXindivf{$k}";
  }
  print STAT "\n";


  print STAT "##Number of sequences in X copy\n";
  print STAT "#CopyNumber";
  foreach my $k (sort { $a <=> $b} keys %h_inXcopy) {
    print STAT "\t$k";
    last if($k >= 25);
  }
  print STAT "\nBefore filtering";
  foreach my $k (sort { $a <=> $b} keys %h_inXcopy) {
    print STAT "\t$h_inXcopy{$k}";
    last if($k >= 25);
  }
  print STAT "\nAfter filtering";
  foreach my $k (sort { $a <=> $b} keys %h_inXcopy) {
    if(defined $h_inXcopyf{$k}) { print STAT "\t$h_inXcopyf{$k}"; }
    else                        { print STAT "\t0";}
    last if($k >= 25);
  }
  print STAT "\n";


  my @a_len = sort( {$a<=>$b} keys(%h_histlen) );
  print STAT "##Length distribution after redundancy filter\n#fastq\t".join("\t",@a_len)."\n";
  for my $k_fq ( keys %hh_histlenbyind ) {
    print STAT "$k_fq";
    for my $k_len  ( sort keys %{ $hh_histlenbyind{$k_fq} } ) {
      print STAT "\t$hh_histlenbyind{$k_fq}{$k_len}";
    }
    print STAT "\n";
  }
  print STAT "ALL";
  for my $k_len  ( sort keys %h_histlen ) {
    print STAT "\t$h_histlen{$k_len}";
  }
  print STAT "\n";
  print STAT "ALL_fil.";
  foreach my $k (sort { $a <=> $b} keys %h_histlenf) {
    print STAT "\t$h_histlenf{$k}";
  }
  print STAT "\n";


  print STAT "##Number of specific sequence of each \"fastq\"\n";
  print STAT "#";
  for my $k_fq ( keys %h_spe ) {
    print STAT "\t$k_fq";
  }
  print STAT "\nBefore filtering";
  for my $k_fq ( keys %h_spe ) {
    print STAT "\t$h_spe{$k_fq}";
  }
  print STAT "\nAfter filtering";
  for my $k_fq ( keys %h_spe ) {
    if(defined $h_spef{$k_fq}) { print STAT "\t$h_spef{$k_fq}"; }
    else                       { print STAT "\t0";}
  }
  print STAT "\n";

  close STAT;
  close TFA;
  close MAT;
}


=pod

=head1 NAME

        matrix.pl - Build a filtered formated  fasta a  file, a
                    statistics file and filtered matrix file.

=head1 DESCRIPTION

        From matrix file provided by join.pl build three files:
          - multifasta file of the filtered sequences
          - statistics file
          - filtered matrix file

        Sequence from matrix input file can be filter by:
          - at least in N different columns (-i)
          - global minimum number of copy (-a)
        Sequence conserved if appears in  -i  different columns
        OR global minimum number of copy is -a.

        Warning: Input matrix file must be sequence(<TAB>int)xN
                 see join.pl for details.

=head1 SYNOPSIS

        matrix.pl ALL.matrix -o OUTPUT_PREFIX
        matrix.pl ALL.matrix -i 1 -a 10 -o OUTPUT_PREFIX

=head1 OPTIONS

        -h,--help

        -i INT
          sequence must be at least in N different columns [2]

        -a INT
          global minimum number of copy [5]

        -o, --output STRING
          MANDATORY: output prefix for each output files
	
=cut

