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


# ------------------------
# 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'
	     ) 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 @a_fastq = @ARGV;
  my %hh_res;
  my %h_len;
  my %h_fq;
  my %h_fqCuta;
  foreach my $fq ( @a_fastq ) {
    open (FQ, $fq) or die "Can't open file $fq : $!";
    my $fqname = basename($fq);
    $fqname =~ s/\..*fastq.*$//;
    while (my $line=<FQ>) {
      my $ac   = $line;
      my $seq  = <FQ>;
      my $ac2  = <FQ>;
      my $qual = <FQ>;
      chomp $seq;
      $hh_res{$fqname}{length($seq)}++;
      $h_len{length($seq)}++;
      $h_fqCuta{$fqname}++;
    }
    close FQ;

    # Catch number of sequence in the fastqc output
    #$h_fq{$fqname} = `grep "Total Sequences" 1_fastqc/*$fqname\_*/fastqc_data.txt`;
    $h_fq{$fqname} = `grep "Total Sequences" $fqname\_*/fastqc_data.txt`;
    $h_fq{$fqname} =~ s/Total Sequences\s+//;
    $h_fq{$fqname} =~ s/\s+//;
  }


  my @a_len = sort( {$a<=>$b} keys(%h_len) );
  print "##Length distribution\n#fastq\t".join("\t",@a_len)."\n";
  for my $k_fq ( keys %hh_res ) {
    print "$k_fq";
    for my $k_len  ( sort keys %{ $hh_res{$k_fq} } ) {
      print "\t$hh_res{$k_fq}{$k_len}";
    }
    print "\n";
  }
  print "ALL";
  for my $k_len  ( sort keys %h_len ) {
    print "\t$h_len{$k_len}";
  }
  print "\n";



  print "##Sequence number  by \"fastq\"\n";
  print "#Step";
  for my $k_fq ( sort keys %h_fqCuta ) {
    print "\t$k_fq";
  }
  print "\n";
  print "Initial";
  for my $k_nb ( sort keys %h_fq ) {
    print "\t$h_fq{$k_nb}";
  }
  print "\n";
  print "Cutadapt";
  for my $k_nb ( sort keys %h_fqCuta ) {
    print "\t$h_fqCuta{$k_nb}";
  }
  print "\n";
}

=pod

=head1 NAME

        stat1.pl - Statistics from fastq file(s) and FastQC output

=head1 DESCRIPTION

        Extract from  fastqc  output  dir  (1_fastqc/)  sequence  number  for
        each related fastq file provided in input.
        Compute from fastq file(s) (cutadapt output) the length distribution.

        Warning: each sequence of fastq file must be on 4 lines

=head1 SYNOPSIS

        stat1.pl FASTQ1 [FASTQ2 ...] > stat1.txt

=head1 OPTIONS

        -h,--help
	
=cut
