#!/usr/bin/env perl

# converts sga format to gff format
# usage: ./sga2gff [-l taglen -x -X] file.sga

use strict;
use Getopt::Long qw(:config no_ignore_case); # silences "duplicate specification" warnings
use Storable qw (retrieve nstore);         # package to store persistently variables in files [http://search.cpan.org/author/AMS/Storable-2.07/Storable.pm]
use File::Spec;

#my %opts;
#getopt('lf:', \%opts);
my %opt;
my @options = ("help", "h", "taglen=i", "l=i", "x", "ext", "X", "i=s", "db=s");

my $file = "";
my $taglen = 0;
my $expand = 0;
my $ext_flag = 0;
my $DB = './';
#my $DB = "/db/genome/";


if( ! GetOptions( \%opt, @options ) ) { &Usage(); }

&Usage()  if defined($opt{'help'}) || defined($opt{'h'});
&Usage()  if $#ARGV < 0;

#define options
if ($opt{'i'} ne '') {
    $DB = $opt{i};
} elsif ($opt{'db'} ne '') {
    $DB = $opt{'db'};
}
my $nstorage_path = File::Spec->catfile($DB, "chro_idx.nstorage");
open FH, $nstorage_path or die "Wrong Chrom Id Storable file \"$nstorage_path\": $!";
my $chr2SV = retrieve($nstorage_path);
close FH;

if ($opt{'l'} ne '') {
    $taglen = $opt{'l'};
}
if ($opt{'taglen'} ne '') {
    $taglen = $opt{'taglen'};
}
if ($opt{'x'} ne '') {
    $expand = 1;
}
if ($opt{'X'} ne '') {
    $ext_flag = 1;
}
if ($opt{'ext'} ne '') {
    $ext_flag = 1;
}

$file = $ARGV[0];

#print "SGA file : $file\n";
#print "Options: taglen : $taglen,  Expand : $expand\n";

# open the SGA file
open (my $SGA, "$file") || die "can't open $file : $!";

my $start;
my $end;
my $source = "ChIPSeq";

my $firstline = <$SGA>;
my @f = split(/\s+/,$firstline);
chomp $f[0];
chomp $f[4];
if ($ext_flag) {
    chomp $f[5];
    if (not defined($f[5])) {
        print STDERR "The extension option is incompatible with the format of your SGA file : it must have the 'description' field defined!/n";
        exit(1);
    }
}
if ($f[0] =~ /^chr/) {
    print "##gff-version   3\n";
    if ($f[3] eq '0') {
        $f[3] = '.';
        $start=$f[2];
        $end=$start;
    }
    if ($f[3] eq '+') {
        $start=$f[2];
        $end=$start + $taglen;
    } elsif($f[3] eq '-') {
        $start=$f[2] - $taglen;
        $end=$f[2];
    }
    if ($ext_flag) {
        if ($expand) {
            for (my $i = $f[4]; $i > 0; $i--) {
                print "$f[0]\t$f[1]\t$f[5]\t$start\t$end\t.\t$f[3]\t.\t1\n";
            }
        } else {
            print "$f[0]\t$f[1]\t$f[5]\t$start\t$end\t$f[4]\t$f[3]\t.\t1\n";
        }
        print_gff_1_ext($SGA);
    } else {
        if ($expand) {
            for (my $i = $f[4]; $i > 0; $i--) {
                print "$f[0]\t$source\t$f[1]\t$start\t$end\t.\t$f[3]\t.\t1\n";
            }
        } else {
            print "$f[0]\t$source\t$f[1]\t$start\t$end\t$f[4]\t$f[3]\t.\t1\n";
        }
        print_gff_1($SGA);
    }
} elsif ($f[0] =~ /N[CT]_\S+\.\d+/ && exists($chr2SV->{$f[0]})) {
    print "##gff-version   3\n";
    if ($f[3] eq '0') {
        $f[3] = '.';
        $start=$f[2];
        $end=$start;
    }
    if ($f[3] eq '+') {
        $start=$f[2];
        $end=$start + $taglen;
    } elsif($f[3] eq '-') {
        $start=$f[2] - $taglen;
        $end=$f[2];
    }
    if ($ext_flag) {
        if ($expand) {
            for (my $i = $f[4]; $i > 0; $i--) {
                print $chr2SV->{$f[0]},"\t$f[1]\t$f[5]\t$start\t$end\t.\t$f[3]\t.\t1\n";
            }
        } else {
            print $chr2SV->{$f[0]},"\t$f[1]\t$f[5]\t$start\t$end\t$f[4]\t$f[3]\t.\t1\n";
        }
        print_gff_2_ext($SGA);
    } else {
        if ($expand) {
            for (my $i = $f[4]; $i > 0; $i--) {
                print $chr2SV->{$f[0]},"\t$source\t$f[1]\t$start\t$end\t.\t$f[3]\t.\t1\n";
            }
        } else {
            print $chr2SV->{$f[0]},"\t$source\t$f[1]\t$start\t$end\t$f[4]\t$f[3]\t.\t1\n";
        }
        print_gff_2($SGA);
    }
} else {
    print STDERR "Unrecognized sequence version $f[0] : please, check the chromosome identifier!\n";
    exit(1);
}
close ($SGA);

sub print_gff_1 {
    my ($fh) = @_;
    my $count = 1;
    if ($expand) {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

            $count++;
            if ($ar[3] eq '0') {
                $ar[3] = '.';
                $start=$ar[2];
                $end=$start;
            }
            if ($ar[3] eq '+') {
                $start=$ar[2];
                $end=$start + $taglen;
            } elsif($ar[3] eq '-') {
                $start=$ar[2] - $taglen;
                $end=$ar[2];
            }
            for (my $i = $ar[4]; $i > 0; $i--) {
                print "$ar[0]\t$source\t$ar[1]\t$start\t$end\t.\t$ar[3]\t.\t$count\n";
            }
        }
    } else {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

#           if (exists($chr2SV->{$species}->{$ar[0]})){
    $count++;
    if ($ar[3] eq '0') {
        $ar[3] = '.';
        $start=$ar[2];
        $end=$start;
    }
    if ($ar[3] eq '+') {
        $start=$ar[2];
        $end=$start + $taglen;
    } elsif($ar[3] eq '-') {
        $start=$ar[2] - $taglen;
        $end=$ar[2];
    }
    print "$ar[0]\t$source\t$ar[1]\t$start\t$end\t$ar[4]\t$ar[3]\t.\t$count\n";
#            }
#            else {
#                my $line = $count + 1;
#                print STDERR "line: $line - Sequence version $ar[0] not in current genome assembly!\n";
#            }
        }
    }
}

sub print_gff_1_ext {
    my ($fh) = @_;
    my $count = 1;
    if ($expand) {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

            $count++;
            if ($ar[3] eq '0') {
                $ar[3] = '.';
                $start=$ar[2];
                $end=$start;
            }
            if ($ar[3] eq '+') {
                $start=$ar[2];
                $end=$start + $taglen;
            } elsif($ar[3] eq '-') {
                $start=$ar[2] - $taglen;
                $end=$ar[2];
            }
            for (my $i = $ar[4]; $i > 0; $i--) {
                print "$ar[0]\t$ar[1]\t$ar[5]\t$start\t$end\t.\t$ar[3]\t.\t$count\n";
            }
        }
    } else {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

#           if (exists($chr2SV->{$species}->{$ar[0]})){
    $count++;
    if ($ar[3] eq '0') {
        $ar[3] = '.';
        $start=$ar[2];
        $end=$start;
    }
    if ($ar[3] eq '+') {
        $start=$ar[2];
        $end=$start + $taglen;
    } elsif($ar[3] eq '-') {
        $start=$ar[2] - $taglen;
        $end=$ar[2];
    }
    print "$ar[0]\t$ar[1]\t$ar[5]\t$start\t$end\t$ar[4]\t$ar[3]\t.\t$count\n";
#            }
#            else {
#                my $line = $count + 1;
#                print STDERR "line: $line - Sequence version $ar[0] not in current genome assembly!\n";
#            }
        }
    }
}

sub print_gff_2 {
    my ($fh) = @_;
    my $count = 1;
    if ($expand) {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

            $count++;
            if ($ar[3] eq '0') {
                $ar[3] = '.';
                $start=$ar[2];
                $end=$start;
            }
            if ($ar[3] eq '+') {
                $start=$ar[2];
                $end=$start + $taglen;
            } elsif($ar[3] eq '-') {
                $start=$ar[2] - $taglen;
                $end=$ar[2];
            }
            for (my $i = $ar[4]; $i > 0; $i--) {
                print $chr2SV->{$ar[0]},"\t$source\t$ar[1]\t$start\t$end\t.\t$ar[3]\t.\t$count\n";
            }
        }
    } else {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

#           if (exists($chr2SV->{$ar[0]})){
    $count++;
    if ($ar[3] eq '0') {
        $ar[3] = '.';
        $start=$ar[2];
        $end=$start;
    }
    if ($ar[3] eq '+') {
        $start=$ar[2];
        $end=$start + $taglen;
    } elsif($ar[3] eq '-') {
        $start=$ar[2] - $taglen;
        $end=$ar[2];
    }
    print $chr2SV->{$ar[0]},"\t$source\t$ar[1]\t$start\t$end\t$ar[4]\t$ar[3]\t.\t$count\n";
#            }
#            else {
#                my $line = $count + 1;
#                print STDERR "line: $line - Sequence version $ar[0] not in current genome assembly!\n";
#            }
        }
    }
}

sub print_gff_2_ext {
    my ($fh) = @_;
    my $count = 1;
    if ($expand) {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

            $count++;
            if ($ar[3] eq '0') {
                $ar[3] = '.';
                $start=$ar[2];
                $end=$start;
            }
            if ($ar[3] eq '+') {
                $start=$ar[2];
                $end=$start + $taglen;
            } elsif($ar[3] eq '-') {
                $start=$ar[2] - $taglen;
                $end=$ar[2];
            }
            for (my $i = $ar[4]; $i > 0; $i--) {
                print $chr2SV->{$ar[0]},"\t$ar[1]\t$ar[5]\t$start\t$end\t.\t$ar[3]\t.\t$count\n";
            }
        }
    } else {
        while(<$fh>){
            my $lin=$_;
            chomp $lin;
            my @ar=split(/\s+/,$lin);

#           if (exists($chr2SV->{$ar[0]})){
    $count++;
    if ($ar[3] eq '0') {
        $ar[3] = '.';
        $start=$ar[2];
        $end=$start;
    }
    if ($ar[3] eq '+') {
        $start=$ar[2];
        $end=$start + $taglen;
    } elsif($ar[3] eq '-') {
        $start=$ar[2] - $taglen;
        $end=$ar[2];
    }
    print $chr2SV->{$ar[0]},"\t$ar[1]\t$ar[5]\t$start\t$end\t$ar[4]\t$ar[3]\t.\t$count\n";
#            }
#            else {
#                my $line = $count + 1;
#                print STDERR "line: $line - Sequence version $ar[0] not in current genome assembly!\n";
#            }
        }
    }
}

sub Usage {
    print STDERR <<"_USAGE_";

    sga2gff.pl [options] <SGA file>

        where options are:
        -h|--help                Show this stuff
        -i|--db <path>           Use <path> to locate Chrom Id Storable File 'chro_idx.nstorage'
        -l|--taglen <taglen>     Set Read length <taglen>
        -x                       Expand SGA lines into multiple GFF lines
        -X|--ext <ext_flag>      Use SGA 6th field to set GFF type field
                                 and store SGA feature into GFF source field

_USAGE_
    exit(1);
}

1;

