#!/usr/bin/env perl

use strict;
use Getopt::Long;

my $assembly = ''; my $sga_file = ''; my $span = 150; my $step=10; my $cutoff = 1;
# Is it just me, or do we have an option called '--span' that goes into variable
# $step, and another option (previously --fragmentSize, now --length) that
# goes into variable $span?
# NOTE: Getopt::Long automatically abbreviates long forms up to uniqueness,
# which is rather cool :-) but forces us to change --fragmentSize to --length
# because we want `-l`, not `-f`.
# TODO use =i for numeric values (instead of =s)
GetOptions('assembly=s' => \$assembly, 'length=s' => \$span, 'span=s' => \$step, 'count=s' => \$cutoff );

if ($assembly eq ''){
    print STDERR "Convert SGA format to WIG with Fixed Step.

Usage:
       sga2wigSmooth_FS.pl <options> <SGA file>

where options are:
       -a|--assembly   Assembly used in the SGA file
       -l|--length     Fragment length of ChIP-seq experiment. Default value 150
       -s|--span       Span used in the WIG sile, default set to 10
       -c|--count      Count cut-off, defoult set to 1

";
    exit;
}

my %t_high;
my $pos;
my $chr;
my $strand;
my $high;
my $oldpos=0;
my $newpos=0;
my $f;
my $start=99999999999999999999;
my $stop=0;
my $oldchr;
my $chk;
my $gap = 2000;
$gap = $span if $span > $gap;

my %chromosome;
# Most scripts do this in 2 steps, so I'm harmonizing.
my $DB = './';
my $chr_file = $DB.$assembly."/chr_NC_gi";
open(CHRN, "<$chr_file");
foreach my $line (<CHRN>){
    chomp $line;
    my @fields = split(/\t/, $line);
    $fields[0] = "M" if ($fields[0] eq "Mt");
    $chromosome{$fields[1]} = "chr".$fields[0];
#    print STDERR "$fields[1]  -> $chromosome{$fields[1]}\n"
}
close(CHRN);

# ucsctools: fetchChromSizes
system("fetchChromSizes $assembly > $assembly.chrom.sizes 2> /dev/null")  if ! -f "$assembly.chrom.sizes";
open(CHRL, "<$assembly.chrom.sizes");
my %l_chr;
foreach my $line (<CHRL>){
    chomp $line;
    my @fields = split(/\t/, $line);
    $l_chr{$fields[0]} = $fields[1];
}
close(CHRL);

#print STDERR "###    Average fragment length set to $span    ###\n\n";

#open(SGA, "<$sga_file");
#foreach my $line (<>){
while (<>){
    my $line = $_;
    chomp $line;
    ($chr, $f, $pos, $strand, $high) = split(/\s+/, $line);
    $high = $cutoff if ($high > $cutoff);
    $pos = $pos - $span + 1 if ($strand eq "-");
    $span = $span + $pos - 1 if $pos < 1;
    $pos = 1 if $pos < 1;
    $newpos = $pos;
    if ($chr eq $oldchr){
        # if the actual position is more distant than the GAP, print
        # the previous region
        if ($pos - $oldpos > $gap){
            print "fixedStep  chrom=$chromosome{$chr}  start=$start  step=$step\n";
            for (my $i = $start; $i < $stop; $i++) {
                my $module = ($i-$start) % $step;
                if ($module == 0){
                    if ($i < $l_chr{$chromosome{$chr}}){
                        print "$t_high{$i}\n"  if ($t_high{$i} ne '');
                        print "0\n"            if ($t_high{$i} eq '');
                    }
                }
            }
            undef %t_high;
            $start = sprintf('%d', $pos);
        }
    }
    else{
        # when chr change, print previous region
        print "fixedStep  chrom=$chromosome{$oldchr}  start=$start  step=$step\n" if $start != 99999999999999999999;
        for (my $i = $start; $i < $stop; $i++) {
            my $module = ($i-$start) % $step;
            if ($module == 0){
                if ($i < $l_chr{$chromosome{$oldchr}}){
                    print "$t_high{$i}\n" if ($t_high{$i} ne "");
                    print "0\n" if ($t_high{$i} eq "");
                }
            }
        }
        undef %t_high;
        $start = sprintf('%d', $pos);
    }

    # fill the region with the hight considering the span:
    until ($newpos - $pos > $span-1){
        $t_high{$newpos} += $high;
        $newpos++;
        $stop = $newpos;
        $stop = $l_chr{$chromosome{$chr}} if ($newpos >= $l_chr{$chromosome{$chr}});
    }

    $start = $pos if ($pos < $start);
    $oldpos = $pos;
    $oldpos += $span;
    $oldchr = $chr;
#    $oldhight = $high;
}

# print the last region beefore closing:
print "fixedStep  chrom=$chromosome{$chr}  start=$start  step=$step\n";
for (my $i = $start; $i < $stop; $i++) {
    my $module = ($i-$start) % $step;
    if ($module == 0){
        if ($i < $l_chr{$chromosome{$chr}}){
            print "$t_high{$i}\n" if ($t_high{$i} ne "");
            print "0\n" if ($t_high{$i} eq "");
        }
        else{
            print "$i > $l_chr{$chromosome{$chr}} $chr\n";
        }
    }
    else{
        print "module = $module\n";
    }
}

