#!/usr/bin/perl

use strict;
use warnings;

#
# usage:
#   make_table.pl [-table TABLE] [-listdir CHR_LIST_DIR] LISTS...
#

sub usage () {
  printf (STDERR "Usage:\n make_table.pl [-table TABLE] [-chr_list_dir CHR_LIST_DIR] [-chr_list_suffix SUFFIX] [-wordlen WORDLEN] [-list_dir LIST_DIR] LISTS|DIRS...\n");
  printf (STDERR "  TABLE - existing table (if sample is already present in table, it is not analyzed again)\n");
  printf (STDERR "  CHR_LIST_DIR - directory where per-chromosome whitelists are found\n");
  printf (STDERR "  SUFFIX - suffix of per-chromosome whitelists (filenames are contructed by adding suffix to chromosome number)\n");
  printf (STDERR "  WORDLEN - k-mer lengths (only needed if new lists are created from fastq files)\n");
  printf (STDERR "  LIST_DIR - directory where newly created sample lists will be stored\n");
  printf (STDERR "  LISTS - full pathnames of sample lists\n");
  printf (STDERR "  DIRS - per sample directories of fastq files\n");
  exit (1);
}

# Program locations
my $gt4dir = ".";
my $glm = "$gt4dir/glistmaker";
my $glq = "$gt4dir/glistquery";
my $glc = "$gt4dir/glistcompare";
# Chromosome whitelists
my @chrs = ("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "X", "Y");

my $table = "";
my $chr_list_dir = ".";
my $chr_list_suffix = "_cleaned_25.list";
my $list_dir = ".";
my $wlen = 25;
my @lists = ();
my @dirs = ();

my $debug = "";

my $i = 0;
while ($i < scalar (@ARGV)) {
  if ($ARGV[$i] eq "-table") {
    $i += 1;
    if (!defined ($ARGV[$i])) { usage (); }
    $table = $ARGV[$i];
  } elsif ($ARGV[$i] eq "-chr_list_dir") {
    $i += 1;
    if (!defined ($ARGV[$i])) { usage (); }
    $chr_list_dir = $ARGV[$i];
  } elsif ($ARGV[$i] eq "-chr_list_suffix") {
    $i += 1;
    if (!defined ($ARGV[$i])) { usage (); }
    $chr_list_suffix = $ARGV[$i];
  } elsif ($ARGV[$i] eq "-list_dir") {
    $i += 1;
    if (!defined ($ARGV[$i])) { usage (); }
    $list_dir = $ARGV[$i];
  } elsif ($ARGV[$i] eq "-wordlen") {
    $i += 1;
    if (!defined ($ARGV[$i])) { usage (); }
    $wlen = $ARGV[$i];
  } elsif ($ARGV[$i] eq "-D") {
    $debug = "-D";
  } elsif (-d $ARGV[$i]) {
    push (@dirs, $ARGV[$i]);
  } elsif (-f $ARGV[$i]) {
    push (@lists, $ARGV[$i]);
  } else {
    usage ();
  }
  $i += 1;
}

if (!@lists && !@dirs) {
  usage ();
}

printf (STDERR "table:           $table\n");
printf (STDERR "chr_list_dir:    $chr_list_dir\n");
printf (STDERR "chr_list_suffix: $chr_list_suffix\n");
printf (STDERR "list_dir:        $list_dir\n");
printf (STDERR "wordlen:         $wlen\n");

if (!(-d $list_dir)) {
  system ("mkdir $list_dir");
}

# Read existing table
my @chr_total = ();
my @chr_gc = ();
my %known;
my @known;
if ($table ne "") {
  open (my $ifs, $table) or die $!;
  # Header
  my $line = <$ifs>;
  # Totals
  $line = <$ifs>;
  chomp ($line);
  my @t = split ("\t", $line);
  @chr_total = splice (@t, 4, @t - 4);
  $line = <$ifs>;
  chomp ($line);
  @t = split ("\t", $line);
  @chr_gc = splice (@t, 4, @t - 4);
  while (my $line = <$ifs>) {
    chomp ($line);
    my @t = split ("\t", $line);
    my $id = $t[0];
    if (!defined ($known{$id})) {
      #printf (STDOUT "Registering %s\n", $id);
      $known{$id} = $line;
      push (@known, $line);
    }
  }
}

# If we have plain directories, try to create lists first
foreach (@dirs) {
  my $cmd = "ls $_/*.fastq $_/*.fastq.gz";
  my $files = qx($cmd);
  $files =~ s/\n/\ /g;
  printf (STDERR "Files: %s\n", $files);
  my @t = split ("/", $_);
  my $name = $t[-1];
  $cmd = "$glm $files -w $wlen $debug -o $list_dir/$name";
  printf (STDERR "Creating list %s\n", $name);
  if (system ($cmd)) {
    printf ("Failed: $cmd\n");
    exit (1);
  }
  my $list = "${list_dir}/${name}_25.list";
  push (@lists, $list);
}
#exit (1);
#printf (STDERR "Lists:\n");
#foreach (@lists) {
#  printf (STDERR "%s\n", $_);
#}

if ((scalar (@chr_total) == 0) || (scalar (@chr_gc == 0))) {
  foreach (@chrs) {
    my $chr = $_;
    my $list = "${chr}${chr_list_suffix}";
    my $cmd = "$glq ${chr_list_dir}/${list} -stat";
    #printf (STDOUT "%s\n", $cmd);
    my $out = qx($cmd);
    # Wordlength      32
    # NUnique 33574272
    # NTotal  39292371
    my @t = split ("\n", $out);
    my @tt = split ("\t", $t[2]);
    my $unique = $tt[1];
    @tt = split ("\t", $t[3]);
    my $total = $tt[1];
    my $gc = gc ("${chr_list_dir}/${list}");
    push (@chr_total, $total);
    push (@chr_gc, $gc);
  }
}

#Header row
printf (STDOUT "CHR\tTOTAL\tUNIQUE\tGC");
foreach (@chrs) {
  printf (STDOUT "\t%s", $_);
}
printf (STDOUT "\n");
#Totals
printf (STDOUT "TOTAL\t\t\t");
for (my $i = 0; $i < @chrs; $i++) {
  printf (STDOUT "\t%s", $chr_total[$i]);
}
printf (STDOUT "\n");
#GC
printf (STDOUT "GC\t\t\t");
for (my $i = 0; $i < @chrs; $i++) {
  printf (STDOUT "\t%.3f", $chr_gc[$i]);
}
printf (STDOUT "\n");

# Known entries
foreach (@known) {
  printf (STDOUT "%s\n", $_);
}

foreach (@lists) {
  my $list = $_;
  my @t = split ("/", $list);
  my $name = $t[-1];
  @t = split ("\\.", $name);
  $name = $t[0];
  $name = substr ($name, 0, length ($name) - 3);
  next if (defined ($known {$name}));
  printf (STDOUT "%s", $name);
  my $cmd = "$glq ${list} -stat";
  #printf (STDOUT "%s\n", $cmd);
  my $out = qx($cmd);
  # Wordlength      32
  # NUnique 33574272
  # NTotal  39292371
  @t = split ("\n", $out);
  my @tt = split ("\t", $t[1]);
  my $wlen = $tt[1];
  @tt = split ("\t", $t[2]);
  my $unique = $tt[1];
  @tt = split ("\t", $t[3]);
  my $total = $tt[1];
  my $gc = gc ("${list}");
  printf (STDOUT "\t%s", $total);
  printf (STDOUT "\t%s", $unique);
  printf (STDOUT "\t%.3f", $gc);
  
  foreach (@chrs) {
    my $chr = $_;
    my $chr_list = "${chr}${chr_list_suffix}";
    my $cmd = "$glc $list ${chr_list_dir}/${chr_list} -i --count_only";
    my $out = qx($cmd);
    #printf (STDERR "%s\n%s\n", $cmd, $out);
    # NUnique 8443982
    # NTotal  8443982
    @t = split ("\n", $out);
    my @tt = split ("\t", $t[1]);
    my $b = $tt[1];
    #printf (STDOUT "%s\t%.2f\t%s\t%s\n", $_, ($a + 0.0) / $b, $a, $b);
    printf (STDOUT "\t%s", $b);
  }
  printf (STDOUT "\n");
}

sub gc {
  my $list = $_[0];
  my $cmd = "$glq $list -gc";
  #printf (STDERR "%s\n", $cmd);
  my $out = qx($cmd);
  #printf (STDERR "%s\n", $out);
  my @t = split ("\t", $out);
  my $gc = $t[1];
  #printf (STDERR "%s\n", $gc);
  return $gc;
}
