#!/usr/bin/env perl
#
# reformats partitioning output files (SGA) into bed format
# usage: partit2bed.pl [-T track_name -d desc --chr chr_nb -b chr_start -e -chr_end] <ChIP-part.output>

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; # To avoid the "dangling /" problem - see l. 36

my %opt;
my @options = ("help", "h", "desc=s", "d=s", "track=s", "color=s",
        "T=s", "t=s", "chr=s", "b=i", "e=i", "db=s", "i=s");

my $desc = "ChIP-Seq";
my $track_name = "Test-BED";
my $color= "200,100,0";
my $chr_nb = "chr0";
my $chr_start = -1;
my $chr_end = -1;
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 '') {   # See gffsga.pl, ll. 34-44
    $DB = $opt{'i'};
}
if ($opt{'db'} ne '') {
    $DB = $opt{'db'};
}

my $nstorage_path = File::Spec->catfile($DB, "chro_idx.nstorage");
my $chroac = retrieve($nstorage_path);

if ($opt{'d'} ne '') {
    $desc = $opt{'d'};
}
if ($opt{'desc'} ne '') {
    $desc = $opt{'desc'};
}
if ($opt{'t'} ne '') {
  warn "WARNING: option -t is deprecated, please use -T.\n";
    $track_name = $opt{'t'};
}
if ($opt{'T'} ne '') {
    $track_name = $opt{'T'};
}
if ($opt{'track'} ne '') {
    $track_name = $opt{'track'};
}
if ($opt{'color'} ne '') {
    $color = $opt{'color'};
}
if ($opt{'chr'} ne '') {
    $chr_nb = "chr".$opt{'chr'};
}
if ($opt{'b'} ne '') {
    $chr_start = $opt{'b'};
}
if ($opt{'e'} ne '') {
    $chr_end = $opt{'e'};
}

my $file = $ARGV[0];

my $start_reg = -1;
my $end_reg = -1;

my $chr_start_flag = 0;
my $chr_end_flag = 0;
my $chr_nb_flag =  0;
my $chr_name_flag =  1;
my $assembly = "";

my @pos = ();

# Read First line of Partition File
# Print out BED File Header
open(PART, "<$file") || die "can't read $file : $!";
my $lin=<PART>;
my @field=split(/\t/,$lin);

if ( $chr_start != -1) {
    $start_reg = $chr_start;
    $chr_start_flag = 1;
} else {
    $start_reg = $field[2];
}
if ( $chr_end != -1) {
    $end_reg = $chr_end;
    $chr_end_flag = 1;
} else {
    $end_reg = $field[2] + 1000000;
}
if ($chr_nb ne "chr0") {
    $chr_nb_flag = 1;
}
if (exists($chroac->{$field[0]})) {
    $chr_name_flag = 0;
}
#if ($chr_nb_flag) {
#    print "browser position $chr_nb:$start_reg-$end_reg\n";
#} elsif (!$chr_name_flag) {
#    print "browser position $chroac->{$field[0]}:$start_reg-$end_reg\n";
#} else {
#    print "browser position $field[0]:$start_reg-$end_reg\n";
#}
#print "browser full refGene\n";
chomp $track_name;
chomp $desc;

print "track name=$track_name description=\"$desc\" visibility=1 color=$color\n";
#if ($chr_name_flag) {
#    print "track name=$track_name description=\"ChIP-Seq-$experiment/$exp_data Partitioning\" visibility=1 color=200,100,0\n";
#} else {
#    print "track name=$track_name description=\"",$chroac->{'assembly'}->{$field[0]},"-ChIP-Seq-$experiment/$exp_data Partitioning\" visibility=1 color=200,100,0\n";
#}

close(PART);

# generate BED-file
open(PART, "<$file") || die "can't read $file : $!";
if ($chr_name_flag) { # we deal with chromosome names (chrx)
    my $i = 0;
    while (<PART>){
# Partitioning algorithm reports positions of start and end of tag,
#NC_000001.9	CTCF	227595	+	1  Start Line
#NC_000001.9	CTCF	227633	-	1  End Line
        my $lin=$_;
        chomp $lin;
        my @field=split(/\t/,$lin);
        if (defined($field[0])){
            if ($field[3] eq '+' || $field[3] eq '0') {
                $pos[0] = $field[2] - 1;
            } elsif ($field[3] eq '-') {
                $i++;
                $pos[1] = $field[2];
                if ( $chr_nb_flag && !$chr_start_flag && !$chr_end_flag ) {
                    if ($field[0] eq $chr_nb) {
#print $field[0],"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                        print $field[0],"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                    }
                } elsif ($chr_nb_flag  && $chr_start_flag  && !$chr_end_flag) {
                    if ($field[0] eq $chr_nb && ($pos[0] >= $chr_start || $pos[1] >= $chr_start)) {
#print $field[0],"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                        print $field[0],"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                    }
                } elsif ($chr_nb_flag  && $chr_start_flag  && $chr_end_flag) {
                    if ($field[0] eq $chr_nb && (($pos[0] >= $chr_start || $pos[1] >= $chr_start) && ($pos[0] <= $chr_end || $pos[1] <= $chr_end))) {
#print $field[0],"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                        print $field[0],"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                    }
                } else {
#print $field[0],"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                    print $field[0],"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                }
                @pos = ();
            }
        }
    }
} else { # Chrom name Flag = 0 (NCBI IDs)
    my $i = 1;
    my $firstline = <PART>;
    chomp $firstline;
    my @field=split(/\t/,$firstline);
    if (defined($field[0])){
        $assembly = $chroac->{'assembly'}->{$field[0]};
    }
    $pos[0] = $field[2] - 1;
    my $secondline = <PART>;
    chomp $secondline;
    @field=split(/\t/,$secondline);
    $pos[1] = $field[2];
    if ( $chr_nb_flag && !$chr_start_flag && !$chr_end_flag ) {
        if ($chroac->{$field[0]} eq $chr_nb) {
            print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
        }
    } elsif ($chr_nb_flag  && $chr_start_flag  && !$chr_end_flag) {
        if ($chroac->{$field[0]} eq $chr_nb && ($pos[0] >= $chr_start || $pos[1] >= $chr_start)) {
            print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
        }
    } elsif ($chr_nb_flag  && $chr_start_flag  && $chr_end_flag) {
        if ($chroac->{$field[0]} eq $chr_nb && (($pos[0] >= $chr_start || $pos[1] >= $chr_start) && ($pos[0] <= $chr_end || $pos[1] <= $chr_end))) {
            print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
        }
    } else {
        print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
    }
    @pos = ();
    while (<PART>){
# Partitioning algorithm reports positions of start and end of tag,
#NC_000001.9	CTCF	227595	+	1  Start Line
#NC_000001.9	CTCF	227633	-	1  End Line
        my $lin=$_;
        chomp $lin;
        my @field=split(/\t/,$lin);

        if (defined($field[0])){
            if ($assembly ne $chroac->{'assembly'}->{$field[0]}) {
                if ($assembly) {
                    print STDERR "The partition output $file appears to contain entries from different species: $field[0] from ",$chroac->{'assembly'}->{$field[0]}," vs. preceeding sequences from $assembly\n";
                }
                $assembly = $chroac->{'assembly'}->{$field[0]};
            }
            if ($field[3] eq '+' || $field[3] eq '0') {
                $pos[0] = $field[2] - 1;
# check chromosome boundaries
                next if ($pos[0] > $chroac->{'length'}->{$field[0]});
            } elsif ($field[3] eq '-') {
                $i++;
                $pos[1] = $field[2];
# check chromosome boundaries
                next if ($pos[1] > $chroac->{'length'}->{$field[0]});
                if ( $chr_nb_flag && !$chr_start_flag && !$chr_end_flag ) {
                    if ($chroac->{$field[0]} eq $chr_nb) {
#print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                        print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                    }
                } elsif ($chr_nb_flag  && $chr_start_flag  && !$chr_end_flag) {
                    if ($chroac->{$field[0]} eq $chr_nb && ($pos[0] >= $chr_start || $pos[1] >= $chr_start)) {
#print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                        print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                    }
                } elsif ($chr_nb_flag  && $chr_start_flag  && $chr_end_flag) {
                    if ($chroac->{$field[0]} eq $chr_nb && (($pos[0] >= $chr_start || $pos[1] >= $chr_start) && ($pos[0] <= $chr_end || $pos[1] <= $chr_end))) {
#print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                        print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                    }
                } else {
#print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t\t\t\+\n";
                    print $chroac->{$field[0]},"\t",$pos[0],"\t",$pos[1],"\t"."R:$i"."\t$field[4]\n";
                }
                @pos = ();
            }
        }
    }
}
close (PART);

sub Usage {
    print STDERR <<"_USAGE_";

    partit2bed.pl [options] <ChIP-Part out file>

        where options are:
        -h|--help                Show this stuff
        -i|--db <path>           Use <path> to locate Chrom Id Storable File 'chro_idx.nstorage'
        -T|--track <track>       BED track name [def: $track_name]
        -t                       Deprecated, please use -T.
        -d|--desc <desc>         Description field of the BED header line [def: $desc]
        --color <RRR,GGG,BBB>    Define the main color for the annotation track [def: $color]
        --chr <str>              Chromosome number (BED declaration lines) [def: 0 - all chromosomes]
        -b <int>                 Chromosome start [def: $chr_start - entire chrom region]
        -e <int>                 Chromosome end [def: $chr_start - entire chrom region]

        The program converts the output of the ChIP-Seq partitioning program to BED format.

_USAGE_
    exit(1);
}

1;

