#!/usr/bin/env perl

# converts gff format to sga format
# usage: ./gff2sga.pl [<-f feature> <-s species> <-C centered> <-u unoriented> <-x extended SGA>] <file.gff>

use strict;
use Getopt::Long qw(:config no_ignore_case); # silences warnings due to -c being (wrongly) equated with -C
use Storable qw (retrieve nstore);         # package to store persistently variables in files [http://search.cpan.org/author/AMS/Storable-2.07/Storable.pm]

#my %opts;
#getopt('asf:', \%opts);  # -a, -s, & -f take arg. Values in %opts, Hash keys will be the switch names

my %opt;
my @options = ( "help", "h", "species=s", "a=s", "s=s", "feature=s", "f=s", "c", "C", "u", "x", "db=s", "i=s");

my $file = "";
my $species = "";
my $feature = "";
my $centered = 0;
my $unoriented = 0;
my $extended = 0;

my $ext_sc = 0;
my $score = 0;

my $DB = './';

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

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

# NOTE: IIRC the idiomatic way to test the existence of a key in a
# hash is:
#          if (exists $hash{$key}) ...
# but since there are a dozen tests below that all use the "ne ''" form (which,
# BTW, is rather strange since it equates the empty string to nil/None/nothing
# and would likely NOT work in most other languages (only checked Python) - how
# then can you store an empty string in a Perl hash?), and since I'm here to
# harmonize options and not to generally modify the code, I'll follow the
# style already in use.
#
if ($opt{'i'} ne '') {   # should be if (exists $opt{'i'}, etc.
    $DB = $opt{'i'};
}
if ($opt{'db'} ne '') {
    $DB = $opt{'db'};
}
open FH, $DB."chro_idx.nstorage" or die "Wrong Chrom Id Storable file $DB.\"chro_idx.nstorage\": $!";
my $chr2SV = retrieve($DB."chro_idx.nstorage");
close FH;

if ($opt{'f'} ne '') {
    $feature = $opt{'f'};
}
if ($opt{'feature'} ne '') {
    $feature = $opt{'feature'};
}
# If more than one of -a, --species and -s are supplied, the last one counts.
if ($opt{'s'} ne '') {
    warn("WARNING: option -s is deprecated, please use -a.\n");
    $species = $opt{'s'};
}
if ($opt{'a'} ne '') {
    $species = $opt{'a'};
}
if ($opt{'species'} ne '') {
    $species = $opt{'species'};
}
if ($opt{'C'} ne '') {
    $centered = 1;
} elsif ($opt{'c'} ne '') {
    warn "WARNING: option -c is deprecated, please use -C.";
    $centered = 1;
}
if ($opt{'u'} ne '') {
    $unoriented = 1;
}
if ($opt{'x'} ne '') {
    $extended = 1;
}

$file = $ARGV[0];

#print "GFF file : $file\n";
#print "Options: feature : $feature,  Species $species, Centered: $centered, Unoriented : $unoriented\n";

my $wrong_seqs = 0;

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

if ($centered) {
    if ($unoriented) {
        if ($extended) {
            print_peak_sga_ext($GFF3);
        } else {
            print_peak_sga($GFF3);
        }
    } else {
        if ($extended) {
            print_peak_sga_oriented_ext($GFF3);
        } else {
            print_peak_sga_oriented($GFF3);
        }
    }
} else {
    if ($extended) {
        print_sga_ext($GFF3);
    } else {
        print_sga($GFF3);
    }
}
close ($GFF3);

sub print_sga {
    my ($fh) = @_;
    if ($species ne '') {
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
### skip comments
            next if (/^\#/);
### skip blank lines
            next if (/^\s*$/);
###################
## syntax checks ##
###################
            my @cols = split(/\t/);
### we should have at least 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
#$wrong_seqs++;
                next;
            }
#           $cols[5]=1;
            if (($cols[5] != int($cols[5])) || ($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            unless (exists($chr2SV->{$species}->{$cols[0]}) || $cols[0] =~ /N[CT]_\S+\.\d+/){
                print STDERR "Chromosome $cols[0] without sequence accession, line $.\n";
                next;
            }
            my $id = "";
            if ($cols[0] =~ /N[CT]_\S+\.\d+/) {
                $id = $cols[0];
            } else {
                $id = $chr2SV->{$species}->{$cols[0]};
            }
            if ($feature ne '') {
                $cols[2] = $feature;
            }
            if ($cols[6] eq '+'){
                if ($ext_sc) {
                    print $id,"\t$cols[2]\t$cols[3]\t$cols[6]\t$cols[5]\t$score\n";
                } else {
                    print $id,"\t$cols[2]\t$cols[3]\t$cols[6]\t$cols[5]\n";
                }
            } elsif ($cols[6] eq '-'){
                if ($ext_sc) {
                    print $id,"\t$cols[2]\t$cols[4]\t$cols[6]\t$cols[5]\t$score\n";
                } else {
                    print $id,"\t$cols[2]\t$cols[4]\t$cols[6]\t$cols[5]\n";
                }
            } elsif ($cols[6] eq '.' || $cols[6] eq '?'){
                if ($ext_sc) {
                    print $id,"\t$cols[2]\t$cols[4]\t0\t1\t$score\n";
                } else {
                    print $id,"\t$cols[2]\t$cols[4]\t0\t1\n";
                }
            }
            undef @cols;
        }
    } else { #if (species) not defined
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
### skip comments
            next if (/^\#/);
### skip blank lines
            next if (/^\s*$/);

###################
## syntax checks ##
###################
            my @cols = split(/\t/);
### we should have 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
                next;
            }
            if (($cols[5] != int($cols[5])) || ($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            if ($feature ne '') {
                $cols[2] = $feature;
            }
            if ($cols[6] eq '+'){
                if ($ext_sc) {
                    print $cols[0],"\t$cols[2]\t$cols[3]\t$cols[6]\t$cols[5]\t$score\n";
                } else {
                    print $cols[0],"\t$cols[2]\t$cols[3]\t$cols[6]\t$cols[5]\n";
                }
            } elsif ($cols[6] eq '-'){
                if ($ext_sc) {
                    print $cols[0],"\t$cols[2]\t$cols[4]\t$cols[6]\t$cols[5]\t$score\n";
                } else {
                    print $cols[0],"\t$cols[2]\t$cols[4]\t$cols[6]\t$cols[5]\n";
                }
            } elsif ($cols[6] eq '.' || $cols[6] eq '?'){
                if ($ext_sc) {
                    print $cols[0],"\t$cols[2]\t$cols[4]\t0\t1\t$score\n";
                } else {
                    print $cols[0],"\t$cols[2]\t$cols[4]\t0\t1\n";
                }
            }
            undef @cols;
        }
    }
}

sub print_sga_ext {
    my ($fh) = @_;
    if ($species ne '') {
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
### skip comments
            next if (/^\#/);
### skip blank lines
            next if (/^\s*$/);
###################
## syntax checks ##
###################
            my @cols = split(/\t/);
### we should have at least 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
#$wrong_seqs++;
                next;
            }
#           $cols[5]=1;
            if (($cols[5] != int($cols[5])) || ($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            unless (exists($chr2SV->{$species}->{$cols[0]}) || $cols[0] =~ /N[CT]_\S+\.\d+/){
                print STDERR "Chromosome $cols[0] without sequence accession, line $.\n";
                next;
            }
            my $id = "";
            if ($cols[0] =~ /N[CT]_\S+\.\d+/) {
                $id = $cols[0];
            } else {
                $id = $chr2SV->{$species}->{$cols[0]};
            }
            if ($feature ne '') {
                $cols[1] = $feature;
            }
            if ($cols[6] eq '+'){
                print $id,"\t$cols[1]\t$cols[3]\t$cols[6]\t$cols[5]\t$cols[2] $score\n";
            } elsif ($cols[6] eq '-'){
                print $id,"\t$cols[1]\t$cols[4]\t$cols[6]\t$cols[5]\t$cols[2] $score\n";
            } elsif ($cols[6] eq '.' || $cols[6] eq '?'){
                print $id,"\t$cols[1]\t$cols[4]\t0\t1\t$cols[2] $score\n";
            }
            undef @cols;
        }
    } else { #if (species) not defined
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
### skip comments
            next if (/^\#/);
### skip blank lines
            next if (/^\s*$/);

###################
## syntax checks ##
###################
            my @cols = split(/\t/);
### we should have 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
                next;
            }
            if (($cols[5] != int($cols[5])) || ($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            if ($feature ne '') {
                $cols[1] = $feature;
            }
            if ($cols[6] eq '+'){
                print $cols[0],"\t$cols[1]\t$cols[3]\t$cols[6]\t$cols[5]\t$cols[2] $score\n";
            } elsif ($cols[6] eq '-'){
                print $cols[0],"\t$cols[1]\t$cols[4]\t$cols[6]\t$cols[5]\t$cols[2] $score\n";
            } elsif ($cols[6] eq '.' || $cols[6] eq '?'){
                print $cols[0],"\t$cols[1]\t$cols[4]\t0\t1\t$cols[2] $score\n";
            }
            undef @cols;
        }
    }
}

sub print_peak_sga {
    my ($fh) = @_;
    if ($species ne '') {
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have at least 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
#$wrong_seqs++;
                next;
            }
            if (($cols[5] != int($cols[5])) || ($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            unless (exists($chr2SV->{$species}->{$cols[0]}) || $cols[0] =~ /N[CT]_\S+\.\d+/){
                print STDERR "Chromosome $cols[0] without sequence accession, line $.\n";
                next;
            }
            my $id = "";
            if ($cols[0] =~ /N[CT]_\S+\.\d+/) {
                $id = $cols[0];
            } else {
                $id = $chr2SV->{$species}->{$cols[0]};
            }
            if ($feature ne '') {
                $cols[2] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            if ($ext_sc) {
                print $id,"\t$cols[2]\t$center\t0\t$cols[5]\t$score\n";
            } else {
                print $id,"\t$cols[2]\t$center\t0\t$cols[5]\n";
            }
            undef @cols;
        }
    } else { #if (species) not defined
        while(<$fh>){
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            if ($feature ne '') {
                $cols[2] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            if ($ext_sc) {
                print "$cols[0]\t$cols[2]\t$center\t0\t$cols[5]\t$score\n";
            } else {
                print "$cols[0]\t$cols[2]\t$center\t0\t$cols[5]\n";
            }
            undef @cols;
        }
    }
}

sub print_peak_sga_ext {
    my ($fh) = @_;
    if ($species ne '') {
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have at least 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
#$wrong_seqs++;
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            unless (exists($chr2SV->{$species}->{$cols[0]}) || $cols[0] =~ /N[CT]_\S+\.\d+/){
                print STDERR "Chromosome $cols[0] without sequence accession, line $.\n";
                next;
            }
            my $id = "";
            if ($cols[0] =~ /N[CT]_\S+\.\d+/) {
                $id = $cols[0];
            } else {
                $id = $chr2SV->{$species}->{$cols[0]};
            }
            if ($feature ne '') {
                $cols[1] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            print $id,"\t$cols[1]\t$center\t0\t$cols[5]\t$cols[2] $score\n";
            undef @cols;
        }
    } else { #if (species) not defined
        while(<$fh>){
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            if ($feature ne '') {
                $cols[1] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            print "$cols[0]\t$cols[1]\t$center\t0\t$cols[5]\t$cols[2] $score\n";
            undef @cols;
        }
    }
}

sub print_peak_sga_oriented {
    my ($fh) = @_;
    if ($species ne '') {
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have at least 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
#$wrong_seqs++;
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            unless (exists($chr2SV->{$species}->{$cols[0]}) || $cols[0] =~ /N[CT]_\S+\.\d+/){
                print STDERR "Chromosome $cols[0] without sequence accession, line $.\n";
                next;
            }
            my $id = "";
            if ($cols[0] =~ /N[CT]_\S+\.\d+/) {
                $id = $cols[0];
            } else {
                $id = $chr2SV->{$species}->{$cols[0]};
            }
            if ($feature ne '') {
                $cols[2] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            if ($ext_sc) {
                print $id,"\t$cols[2]\t$center\t$cols[6]\t$cols[5]\t$score\n";
            } else {
                print $id,"\t$cols[2]\t$center\t$cols[6]\t$cols[5]\n";
            }
            undef @cols;
        }
    } else { #if (species) not defined
        while(<$fh>){
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            if ($feature ne '') {
                $cols[2] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            if ($ext_sc) {
                print "$cols[0]\t$cols[2]\t$center\t$cols[6]\t$cols[5]\t$score\n";
            } else {
                print "$cols[0]\t$cols[2]\t$center\t$cols[6]\t$cols[5]\n";
            }
            undef @cols;
        }
    }
}

sub print_peak_sga_oriented_ext {
    my ($fh) = @_;
    if ($species ne '') {
        while(<$fh>){
# chr2    NimbleScan      1859802:RajiIPH3/SJOIPH3:BLOCK1 132843487       132843531       -0.11   +       .       seq_id=18S_CHR2:132843487-132849513;probe_id=CHR0200P132843487;count=1
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have at least 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
#$wrong_seqs++;
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            unless (exists($chr2SV->{$species}->{$cols[0]}) || $cols[0] =~ /N[CT]_\S+\.\d+/){
                print STDERR "Chromosome $cols[0] without sequence accession, line $.\n";
                next;
            }
            my $id = "";
            if ($cols[0] =~ /N[CT]_\S+\.\d+/) {
                $id = $cols[0];
            } else {
                $id = $chr2SV->{$species}->{$cols[0]};
            }
            if ($feature ne '') {
                $cols[1] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            print $id,"\t$cols[1]\t$center\t$cols[6]\t$cols[5]\t$cols[2] $score\n";
            undef @cols;
        }
    } else { #if (species) not defined
        while(<$fh>){
            chomp;
            next if (/^\#/);
            next if (/^\s*$/);

            my @cols = split(/\t/);
### we should have 8 columns.
            if (scalar(@cols) < 8) {
                print STDERR "Incorrect column count, line $.\n";
                next;
            }
            if (($cols[5] != int($cols[5]))||($cols[5] eq '.') || ($cols[5] < 0) || ($ext_sc)) {
                $score = $cols[5];
                $cols[5] = 1;
                $ext_sc = 1;
            }
            if ($feature ne '') {
                $cols[1] = $feature;
            }
            my $center = int(($cols[3]+$cols[4])/2);
            print "$cols[0]\t$cols[1]\t$center\t$cols[6]\t$cols[5]\t$cols[2] $score\n";
            undef @cols;
        }
    }
}

sub Usage {
    print STDERR <<"_USAGE_";

    gff2sga.pl [options] <GFF file>

        where options are:
        -h|--help                Show this stuff
        -i|--db <path>           Use <path> to locate Chrom Id Storable File 'chro_idx.nstorage'
        -f|--feature <feature>   Set Feature name <feature>
        -a|--species <species>   Assembly <species> (i.e hg18)
        -s                       Deprecated, please use -a.
        -C                       Generate a Centered SGA file
        -u                       Generate an Unoriented SGA file
        -x                       Generate an extended SGA file with the
                                 6th field equal to the GFF 'type' field,
                                 and the feature field equal to the GFF 'source'

_USAGE_
    exit(1);
}

1;

