#!/usr/bin/perl -w

#
# tasuke_bedtotsv.pl
# http://tasuke.dna.affrc.go.jp/
#
# Author: Ryutaro Itoh
#
# Copyright (c) 2013 National Institute of Agrobiological Sciences   
#

use strict;
use Getopt::Long;
use File::Basename;
use File::Copy;
use Cwd;

my $OUTFILE = "";
my $CHRLIST = "";
my $BED = "";
my $FORMAT = "BED";
my $bedgraph;
my $VALUE_CL=4;

################################################################
# Checking parameters and database connection
################################################################
GetOptions( 
            'o|outfile=s'     => \$OUTFILE,
            'c|chr=s'        => \$CHRLIST,
            'g|graph'        => \$bedgraph,
            'i|inbam=s'      => \$BED     
);

if ($OUTFILE eq "" || $CHRLIST eq "" || $BED eq "") {
  _error_method();  
  exit;
}
if($bedgraph){
    $FORMAT = "BEDGraph";
    $VALUE_CL=3;
}

open IN, '<', $BED or die "file open error: $!";
my $columns    =0;
    while (<IN>) {
        chomp;
        my @line = split(/\t/, $_);
        $columns = @line;
        last;
    }
close(IN);
if($FORMAT eq "BED" && $columns !=6){
    print "The Input file is probably not BED formatted file.\n\n";
    exit;
}elsif($FORMAT eq "BEDGraph" && $columns !=4){
    print "The Input file is probably not BEDGraph formatted file.\n\n";
    exit;
}

###########################################################
# Creating TSV files 
###########################################################
#Set chromosomes 
my $OUTDIR=$OUTFILE."_TMP";
my $loadchr_rt = _loadChr($CHRLIST);
_print_time(" Begin converting ".$FORMAT." to TSV");
_mkTmpDir($OUTDIR);
_createTSV($BED, $OUTDIR, $loadchr_rt, $OUTFILE);
_print_time(" Complete!");

###########################################################
# Loading chromosome list
###########################################################
sub _loadChr{
    my $CHRLIST = shift;
    my @CHR;
    my $max = 0;
    open IN, '<', $CHRLIST or die "file open error: $!";
    
    while (<IN>) {
        chomp;
        my @line = split(/,/, $_);
        if($line[1] !~ /\d/){
            print "CSV file is not correct.\n";
            exit;
        }else{
            push(@CHR,$line[0]);
            if($line[1]>$max){
                $max=$line[1];
            }
        }
    }    
    close(IN);
    my @chr_return = ($max, \@CHR);
    return \@chr_return;
}


##############################################
# create bam
##############################################
sub _mkTmpDir{
    my $OUTDIR = shift;
    
    if(!-e "$OUTDIR"){
        mkdir("$OUTDIR") or die "$!: $OUTDIR";
    } else{
        print "$OUTDIR already exists.\n";
        exit(0);
    }

}

###########################################################
# samtools depth
###########################################################
sub _createTSV{
    my $BED = shift;
    my $OUTDIR = shift;
    my $rt = shift;
    my @rt = @$rt;
    my $max= shift @rt;
    my $chr_sc = shift @rt;
    my $OUTFILE = shift;
    my @CHR= @$chr_sc;
    my @CHR_files;
    my @CHR_header;
    my @CHR_DEL;
    my $CD = Cwd::getcwd();

    open OUT, '>', $OUTDIR."/"."position" or die "file open error: $!";
    for(my $i=1; $i< $max+1; $i++){
        print OUT $i."\n";
    }
    close(OUT);
    foreach my $x(@CHR){
        my $min=1;
        my $searched=0;
        my $end=1;
        open IN, '<', $BED or die "file open error: $!";
        open OUT, '>', $OUTDIR."/".$x or die "file open error: $!";
        READLINE: while (<IN>) {
            chomp;
            my @line = split(/\t/, $_);# chr, pos , end , val
            $line[1]++;
            if($line[0] eq $x){
                if($searched==1){#blank between 2 rows
                    for(my $b=$end+1; $b< $line[1]; $b++){
                        if($b > $max){last READLINE;}
                        print OUT "\\N\n";
                    }
                }else{# first line of each chr
                    for(my $a=1; $a< $line[1]; $a++){
                        if($a > $max){last READLINE;}
                        print OUT "\\N\n";
                    }
                }
                for(my $c=$line[1]; $c<= $line[2]; $c++){
                    if($c > $max){last READLINE;}
                    print OUT $line[$VALUE_CL]."\n";
                }
                $end=$line[2];
                $searched=1;
            }elsif(($line[0] ne $x) and $searched == 1){
                if($end <= $max){
                    for(my $i=$end+1; $i< $max+1; $i++){
                        print OUT "\\N\n";
                    }
                }
                last;
            }
            if(eof IN and $searched == 1){
                if($end <= $max){
                    for(my $i=$end+1; $i< $max+1; $i++){
                        print OUT "\\N\n";
                    }
                }
            }
        }
        if($searched==0){
            for(my $d=1; $d< $max+1; $d++){
                print OUT "\\N\n";
            }
        }
        close(IN);
        close(OUT);
    }
    
    foreach my $x(@CHR){
        my $val=0;
        my $isnull=0;
        my $line=1;
        my $searched=0;
        my $prtline=0;
        
        open IN, '<', $OUTDIR."/".$x or die "file open error: $!";
        open OUT, '>', $OUTDIR."/".$x."_K" or die "file open error: $!";
        while (<IN>) {
            chomp;
            if ($_ =~ /^[0-9]+$/) {
                $val+=$_;
            }else{
                $isnull++;
            }
            if($line % 1000 == 0 ){
                if($isnull != 1000){
                    print OUT int($val/1000)."\n";
                }else{
                    print OUT "\\N\n";
                }
                $val=0;
                $isnull=0;
                $prtline++;
            }
            if(eof IN){
                if($isnull != 1000){
                    print OUT int($val/1000)."\n";
                }else{
                    print OUT "\\N\n";
                }
                $prtline++;
                for(my $i=$prtline+1; $i< $max+1; $i++){
                    print OUT "\\N\n";
                }
            }
            $line++;
        }
        close(IN);
        close(OUT);
    }

    push(@CHR_header, "#position#reads=none");
    foreach my $x(@CHR){
        push(@CHR_files, $OUTDIR."/".$x);
        push(@CHR_files, $OUTDIR."/".$x."_K");
        push(@CHR_header, $x);
        push(@CHR_header, $x."_K");
        push(@CHR_DEL, $x);
        push(@CHR_DEL, $x."_K");
    }
    
    my $param1= join(" ", @CHR_files);
    my $param2 = $OUTDIR."/".basename($BED).".depth_tmp";
    my $param3 = $OUTDIR."/".basename($BED).".depth";
    my $param4= join("\t", @CHR_header);
    push(@CHR_DEL, "position");
    push(@CHR_DEL, "header");
    push(@CHR_DEL, basename($BED).".depth_tmp");
    
    my $cmd2 = "paste ".$OUTDIR."/position $param1 > $param2";
    $cmd2= sprintf("%s", $cmd2);
    if (system($cmd2)){
        print  "Error: merge depth file.\n";
        exit;
    }
    my $cmd3="echo \"".$param4."\" > ".$OUTDIR."/header";
    $cmd3= sprintf("%s", $cmd3);
    if (system($cmd3)){
        print  "Error: creating header file.\n ";
        exit;
    }
    my $cmd4="cat ".$OUTDIR."/header ".$param2." | head -n  ".($max+1)." > ".$param3;
    $cmd4= sprintf("%s", $cmd4);
    if (system($cmd4)){
        print  "Error: creating depth file.\n ";
        exit;
    }
    my $result = basename($BED).".depth";
    move($OUTDIR."/".$result, $CD."/".$OUTFILE);
    clean_dir($OUTDIR, \@CHR_DEL, $CD);
    rmdir($OUTDIR);
}

#### time #######################################
sub _print_time {
  my ($str) = @_;
  my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime;
  $year += 1900;
  $mon += 1;
  my $datetime = sprintf "%04d-%02d-%02d %02d:%02d:%02d", $year, $mon, $mday, $hour, $min, $sec;
  print "$str : $datetime\n";
}

### Removing temp files ###############################
sub clean_dir{
    my $work_dir = shift;
    my $chr_sc = shift;
    my $CD = shift;
    my @files= @$chr_sc;
    chdir("$work_dir") or die "$! : $work_dir";
    foreach my $file(@files){    
        if(-e $file ){
            unlink $file;
        }
    }
    chdir($CD) or die "$! : $CD";
}

#### print usage #######################################
sub _error_method {
  print "\nConvert BED or BedGraph to TSV file.\n";
  print "tasuke_bedtotsv.pl [Options] -i [bed] -c [chromosome_information] -o [output_file]  \n\n";
  print "Options : \n";
  print "-g    Loading as BEDGraph format.\n\n";
}
