#!/usr/bin/perl

use strict;
use warnings;

#
# Usage
#   generate_polymorphism_lists.pl WORDLENGTH MIN_MAF REFERENCE_DIR VCF_FILES...
#
# VCF files have to be named with the same names as chromosomes (i.e. 12.fa <-> 12.vcf
# REFERENCE_DIR is directory containing chromosome FastA files (i.e. 1.fa, 2.fa)
#
# The resulting file is named:
#   polymorphic_blacklist_WORDLENGTH.list
#

if (!defined ($ARGV[3])) {
  printf (STDERR "Usage\n  generate_polymorphism_lists.pl WORDLENGTH MIN_MAF REFERENCE_DIR VCF_FILES\n");
  exit (1);
}

# Path to GenomeTester4 programs
my $gt4path = ".";
# Path for temporary files
my $tmpdir = "gpltmp";
# Print out command lines etc.
my $debug = 1;
 
my $glm = "$gt4path/glistmaker";
my $glc = "$gt4path/glistcompare";

my $wlen = $ARGV[0];
my $minmaf = $ARGV[1];
my $refdir = $ARGV[2];
my @vcf = splice (@ARGV, 3, @ARGV - 3);

system ("mkdir $tmpdir");

my @names;
my @hnames;

foreach (@vcf) {
  my @t = split ("/", $_);
  @t = split ("\\.", $t[-1]);
  my $name = $t[0];
  # Generate non-canonical haplotypes
  my $pseudoname = "$tmpdir/$name.pseudo";
  my $cmd = "./hapl_generator.pl $wlen $minmaf $_ $refdir/$name.fa > $pseudoname";
  if ($debug) {
    printf (STDERR "%s\n", $cmd);
  }
  if (system ($cmd)) {
    printf (STDERR "Cannot generate non-canonical haplotypes\n");
    exit (1);
  }
  # Generate true haplotypes
  my $kmername = "$tmpdir/$name.kmers";
  $cmd = "./kmer_creator.pl $wlen $pseudoname > $kmername";
  if ($debug) {
    printf (STDERR "%s\n", $cmd);
  }
  if (system ($cmd)) {
    printf (STDERR "Cannot generate haplotypes\n");
    exit (1);
  }
  system ("rm $pseudoname");
  push (@names, $name);
  push (@hnames, "$kmername");
}

# Generate list
my $allhtypes = join (" ", @hnames);
my $cmd = "$glm $allhtypes -w $wlen -o $tmpdir/tmpht";
if ($debug) {
  printf (STDERR "%s\n", $cmd);
}
if (system ($cmd)) {
  printf (STDERR "Cannot generate list\n");
  exit (1);
}
# Cut frequencies to 1
$cmd = "$glc $tmpdir/tmpht_$wlen.list $tmpdir/tmpht_$wlen.list -i -r 1 -o $tmpdir/tmpht2";
if ($debug) {
  printf (STDERR "%s\n", $cmd);
}
if (system ($cmd)) {
  printf (STDERR "Cannot crop list frequencies\n");
  exit (1);
}
system ("mv $tmpdir/tmpht2_${wlen}_intrsec.list plymorphic_blacklist_$wlen.list");
system ("rm $tmpdir/tmpht_$wlen.list");
foreach (@hnames) {
  my $name = $_;
  system ("rm $name");
}
system ("rmdir $tmpdir");
