#!/usr/bin/perl -w

#
# tasuke_sort_by_distance.pl
# http://tasuke.dna.affrc.go.jp/
#
# Author: Ryutaro Itoh
#
# Copyright (c) 2013 National Institute of Agrobiological Sciences   
#
# This software is released under the MIT License.
#
# http://opensource.org/licenses/mit-license.php
#

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

my $DB      = '';
my $HOST    = 'localhost';
my $USER    = '';
my $PASS    = '';
my $RDBMS = '';
my $TABLESPACE = '';
my $OUTFILE = "order.conf";
my $MODE = "";
my $FILE = "TASUKE_sequences";
my $HOME = "";
my $CRTDIR = Cwd::getcwd();
my $OUTDIR = $CRTDIR."/".basename($FILE)."_tmp";
my $BIN = "fasta.php";
my @ALLFILES = ("cmd_dnadist","cmd_tree","header","infile","infile_dnadist","order.conf","outfile","outfile_dnadist","outtree","phylip_wo_header", $FILE);

################################################################
# Checking parameters and database connection
################################################################
GetOptions( 
            'db|database=s' => \$DB,
            'h|host=s'      => \$HOST,
            'u|user=s'      => \$USER,
            'b|backend=s'   => \$RDBMS,
            'o|outfile=s'   => \$OUTFILE,
            's|space=s'     => \$TABLESPACE,
            'm|mode=s'     => \$MODE,
            'p|password=s'  => \$PASS
);

if ( !$DB || !$USER || !$OUTFILE) {
    _error_method();
    exit;
}

my $connect;
if($RDBMS eq "oracle"){
    if(!$TABLESPACE){
        _error_method();
        exit;
    }
    $connect = DBI->connect('DBI:Oracle:host='.$HOST.';SERVICE_NAME='.$DB.';',$USER,$PASS)
    or die "Can't connect to database: $DBI::errstr\n";
    if(!_check_exists_tablespace($TABLESPACE)){
        print "Tablespace: '$TABLESPACE' is not found.";
        _exit(1);
    }
    $BIN = "/oracle/fasta_oracle.php";
}else{
    $connect= DBI->connect('DBI:mysql:database='.$DB.';host='.$HOST,$USER,$PASS)
    or die "Can't connect to database: $DBI::errstr\n";
}

if($MODE eq "phylip"){
    if(!defined $ENV{"TASUKE_HOME"}){
        print "Set the directory of your tasuke instance to environment variable 'TASUKE_HOME'. \n";
        exit;
    }
    if(system("type dnadist >/dev/null 2>&1") != 0){
        print "Set the path of PHYLIP to environment variable. \n";
        exit;
    }
    if(-e "$CRTDIR/$OUTFILE"){
        print "$CRTDIR/$OUTFILE already exists.\n";
        exit;
    }
    if(!-e $OUTDIR){
        mkdir("$OUTDIR") or die "$!: $OUTDIR";
        chdir $OUTDIR or die "$!: $OUTDIR";
    } else{
        print "$OUTDIR already exists.\n";
        exit;
    }
    my $accessions_ref = _get_accessions();
    _create_fasta($accessions_ref);
    _conv_compdel_seq();
    _prep_tree();
    _make_tree();
    _make_orderFile();
    move("order.conf","$CRTDIR/$OUTFILE");
    _clean_dir("./",\@ALLFILES,$CRTDIR);
    rmdir($OUTDIR);
    exit;
}else{
    print "\nGetting distances between each accession and reference genome.\n";
    my $accessions_ref = _get_accessions();
    _create_order_file($accessions_ref);
    print "Finished!\n";
    _exit();
}

###########################################################
# Functions
###########################################################

### Creating order.conf based on the accession distance ###
sub _create_order_file{
    my $accessions_ref = shift(@_);
    my @accessions = @$accessions_ref;
    my %distarray;
    
    open OUT, '>', $OUTFILE or die "file open error: $!";
    foreach my $l (@accessions) {
        $distarray{$l}= _get_distance($l);
    }
    foreach my $key (sort {$distarray{$a} <=> $distarray{$b} || $a cmp $b} keys %distarray){
        print OUT _get_accession_id($key)."\n";
    }
    close(OUT);
}

### getting accessions  #########################################
sub _get_accessions{
    my $val="";
    my @id;
    my $query = "SELECT accession_db FROM sys_accession";
    my $query_handle = $connect->prepare($query);
    
    $query_handle->execute();
    $query_handle->bind_columns(undef, \$val);
    while($query_handle->fetch()) {
        push @id, $val;
    }
    return \@id;
}

### getting accession id #########################################
sub _get_accession_id{
    my $val="";
    my @id;
    my $str = shift(@_);
    my $query = "SELECT accession_id FROM sys_accession where accession_db='$str'";
    my $query_handle = $connect->prepare($query);
    
    $query_handle->execute();
    $query_handle->bind_columns(undef, \$val);
    $query_handle->fetch();
    return $val;
}

### getting accessions  #########################################
sub _get_distance{
    my $val="";
    my @id;
    my $str = shift(@_);
    $str = $str."_variant";
    my $query = "SELECT count(seq) FROM $str";
    my $query_handle = $connect->prepare($query);
    
    $query_handle->execute();
    $query_handle->bind_columns(undef, \$val);
    $query_handle->fetch();
    return $val;
}

### getting sequences from database  #########################################
sub _create_fasta{
    my $list_ref = shift(@_);
    my @list = @$list_ref;
    my $accessions = join(",",@list);
    my $home = $ENV{"TASUKE_HOME"};
    my $len = 1;
    
    foreach my $x(_chrs()){
        $len = _max($x);
        if(system("php $home/bin/$BIN 'prm_fa_chr=$x
                                    &prm_fa_s=1
                                    &prm_fa_e=$len
                                    &prm_fa_tgt=1
                                    &prm_fa_qual=20
                                    &prm_fa_dp=2
                                    &prm_fa_list=$accessions
                                    &prm_fa_name=id
                                    &prm_fa_nameS=subvariety
                                    &prm_fa_ref=Reference
                                    &prm_fa_width=80
                                    &prm_fa_avd=0
                                    &prm_fa_p3prm=
                                    &prm_fa_Ns=0
                                    &prm_type=a
                                    &path=$home/bin/' >> ".$FILE."_$x") !=0)
        {
            print "Error: Can not execute fasta.php for getting each sequences.\n";
            exit;
        }
    }
    foreach my $y(@list){
        open OUT, '>>', $y or die "file open error: $!";
        print OUT ">".$y."\n";
        push(@ALLFILES, $y);
        foreach my $z(_chrs()){
            my $rd_flag = 0;
            open IN, '<', $FILE."_$z" or die "file open error: $!";
            push(@ALLFILES, $FILE."_$z");
            while (<IN>) {
                chomp;
                if($_ =~ /^>$y/){
                    $rd_flag = 1;
                }else{
                    if($rd_flag == 1){
                        if($_ =~ /^>/){
                            last;
                        }else{
                            print OUT $_."\n";
                        }
                    }
                }
            }
            close(IN);
        }
        close(OUT);
    }
    
    open OUTREF, '>>', "REF" or die "file open error: $!";
    print OUTREF ">"."REF"."\n";
    foreach my $z(_chrs()){
        my $rd_flag = 0;
        open INREF, '<', $FILE."_$z" or die "file open error: $!";
        push(@ALLFILES, "REF");
        while (<INREF>) {
            chomp;
            if($_ =~ /^>REF/){
                $rd_flag = 1;
            }else{
                if($rd_flag == 1){
                    if($_ =~ /^>/){
                        last;
                    }else{
                        print OUTREF $_."\n";
                    }
                }
            }
        }
        close(INREF);
    }
    close(OUTREF);
    unshift(@list, "REF");
    $accessions = join(" ",@list);
    system("cat $accessions > $FILE");
}

### remove disused sites from each seq
sub _conv_compdel_seq{
    my $in = $FILE;
    my $out = $FILE."_compdel";
    my $home = $ENV{"TASUKE_HOME"};
    
    $BIN = "dnadist.php";
    print "$in:$out:$OUTDIR:$home/bin/";
    if(system("php $home/bin/$BIN 'prm_phy_chr=none
                                &prm_phy_s=1
                                &prm_phy_e=1
                                &prm_phy_qual=20
                                &prm_phy_dp=2
                                &prm_phy_list=none
                                &prm_phy_cri=2
                                &prm_phy_gene=0
                                &prm_phy_objid=none
                                &inseq=$in
                                &outseq=$out
                                &tmpdir=$OUTDIR
                                &path=$home/bin/'") !=0)
    {
        print "Error: Can not execute dnadist.php for getting each sequences.\n";
        exit;
    }
    unlink $in;
    #rename $in, $in."_raw"  or die $!;
    rename $out, $FILE  or die $!;
}

### making tree using phylip #########################################
sub _prep_tree{
    #Converting fasta file to phylip formatted file (infile).
    my $nuc_length = 0;
    my $num_seq = 0;
    open IN, '<', $FILE or die "file open error: $!";
    open OUT, '>', "phylip_wo_header" or die "file open error: $!";
    while (<IN>) {
        chomp;
        if($_ =~ /^>/){
            $num_seq++;
            $nuc_length = 0;
            my $ident = $_;
            $ident =~ s/^>(.*)/$1/;
            if(length($ident) > 10){
                $ident = substr($ident, 0, 10);
            }else{
                for(my $i=length($ident)+1; $i<11; $i++){
                    $ident .= " ";
                }
            }
            $ident .= " ";
            print OUT "$ident\n";
        }else{
            print OUT $_."\n";
            $nuc_length += length($_);
        }
    }
    close(IN);
    close(OUT);
    
    open HD_OUT, '>', "header" or die "file open error: $!";
    print HD_OUT "$num_seq $nuc_length\n";
    close(HD_OUT);
    
    open CMD_DNA_OUT, '>', "cmd_dnadist" or die "file open error: $!";
    print CMD_DNA_OUT "i\nY\n";
    close(CMD_DNA_OUT);
    
    open CMD_TR_OUT, '>', "cmd_tree" or die "file open error: $!";
    print CMD_TR_OUT "o\n1\nY\n";
    close(CMD_TR_OUT);
    
    if (system("cat header phylip_wo_header> infile")){
        print  "Creating phylip formatted file failed.\n Exit.\n";
        exit;
    }   
}

sub _make_tree{
    if (system("dnadist < cmd_dnadist >/dev/null")){
        print  "PHYLIP/dnadist: failed\n Exit.\n";
        exit;
    }
    rename "infile", "infile_dnadist" or die $!;
    copy "outfile","infile" or die $!;
    rename "outfile", "outfile_dnadist" or die $!;
    if (system("neighbor < cmd_tree >/dev/null")){
        print  "PHYLIP/neighbor: failed\n Exit.\n";
        exit;
    }
}

sub _make_orderFile{
    my @order;
    open IN, '<', "outfile" or die "file open error: $!";
    while (<IN>) {
        if($_ =~ /^Between\s+And\s+Length/){last;}
        chomp;
        if($_ =~ /.*\+d\d+/){
            $_ =~ s/.*\+(d\d+).*/$1/;
            unshift(@order, $_);
        }
    }
    close(IN);
    open OUT, '>', "order.conf" or die "file open error: $!";
    foreach my $x(@order){
        my $str = _get_accessionid($x);
        print OUT $str."\n";
    }
    close(OUT);
}


###########################################################
# Functions
###########################################################

### Execute DB ################################################
sub _dbexe {
    my ($query) = @_;
    my $query_handle = $connect->prepare($query);
    $query_handle->execute() or warn "Excute error: $DBI::errstr\n";
}

#### Usage #######################################
sub _error_method {
    print "\nSorting registered accessions by sequence distance\n";
    print "tasuke_sort_by_distance.pl -db [database] -u [user] -p [password] -o [output] \n( For oracle database: -b oracle -s [tablespace])\n\n";
}

#### 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";
}

### Checking table ####################################
sub _check_exists_table {
    my ($t) = @_;
    my $val = $connect->selectrow_array("show tables like '$t'");
    return $val;
}

### Checking tablespace ###############################
sub _check_exists_tablespace {
    my ($t) = @_;
    $t = uc($t);
    my $val = $connect->selectrow_array("select tablespace_name from user_tablespaces where tablespace_name = '$t'");
    return $val;
}

### getting accession id #############################
sub _get_accessionid{
    my ($t) = @_;
    my $val = "";
    my $query = "SELECT accession_id FROM sys_accession where accession_db = '$t'";
    my $query_handle = $connect->prepare($query);
    
    $query_handle->execute();
    $query_handle->bind_columns(undef, \$val);
    $query_handle->fetch();
    return $val;
}

sub _chrs{
    my $val="";
    my @c;
    my $query = "SELECT id FROM sys_chromosome";
    my $query_handle = $connect->prepare($query);
    
    $query_handle->execute();
    $query_handle->bind_columns(undef, \$val);
    while($query_handle->fetch()) {
        push @c, $val;
    }
    return @c;
}

sub _max{
    my ($t) = @_;
    my $val;
    my $query;
    if($RDBMS eq "oracle"){
        $query = "SELECT seq_size FROM sys_chromosome where id='$t'";
    }else{
        $query = "SELECT size FROM sys_chromosome where id='$t'";
    }
    my $query_handle = $connect->prepare($query);
    $query_handle->execute();
    $val = $query_handle->fetchrow_hashref;
    if($RDBMS eq "oracle"){
        return $val->{uc('seq_size')};
    }else{
        return $val->{'size'};
    }
}

#### Disconnection and exit #######################################
sub _exit {
    my ($arg) = @_;
    my $code = 0;
    if(defined $arg){$code = $arg;}
    $connect->disconnect or warn "Disconnection error: $DBI::errstr\n";
    exit($code);
}

sub _clean_dir{
    my $work_dir = shift;
    my $list = shift;
    my $CD = shift;
    my @files = @$list;
    
    chdir("$work_dir") or die "$! : $work_dir";
    foreach my $file(@files){    
        if(-e $file ){
            unlink $file;
        }
    }
    chdir($CD) or die "$! : $CD";
}

