#! /usr/bin/perl

use POSIX;
use strict;
use File::Basename;
use Getopt::Std;

my %opts=( k=>undef, n=>undef, r=>undef, s=>100, f=>".", u=>undef, q=>undef, c=>25);
getopts('u:k:r:g:f:n:s:c:',\%opts);

if (@ARGV <1)  {
    print "usage: sim_imrdenom.pl [-k -r -g -f -n] <ref>", "\n";
    exit;
}

my $key_id=$opts{k};   #
my $ofolder=$opts{f};  # set the upper folder. final folder will be $ofolder/$key_id
my $readlen=$opts{s};
my $cov=$opts{c};

my $setref=$ARGV[0];

if ( ! $setref ) {
	print "please set the refernce file\n";
}

my @paired;
my @single;
my @paired_uneq;


if (!defined $key_id) {
	my $BATCHW="sim_";
	my ($BATCHI)= split(/\./, basename($setref));
	$key_id=$BATCHW.$BATCHI;
}


if ( ! -e "$setref.fai" ) {
	die "Please index the reference first\n"
}

open my $fh, "<$setref.fai" or die "cannot open the reference index file, please index the reference first\n";
my $chrlen;
while (<$fh>)  {
   	my $gp=$_;
	my @s=split;
	$chrlen += $s[1];
}
#print STDERR "reference genome size is $chrlen\n";

my $nread= ceil($chrlen* $cov / 2 / $readlen);

my $outread1="./outread_1.fa";
my $outread2="./outread_2.fa";


my $cmd="wgsim -N $nread -1 $readlen -2 $readlen -r 0 -R 0 $setref $outread1 $outread2";
print $cmd, "\n";
system($cmd);


my $tfile=$ofolder."/$key_id.t";
open my $tfh, ">$tfile" or die "cannot write to file $tfile\n";
print $tfh "outputfolder $ofolder/$key_id\n";
print $tfh "reference $setref\n";
print $tfh "iterations 1\nthread 4\nmaxreads 2000000\n";
print $tfh "loaddata $outread1 $outread2\n";
