#!/usr/bin/perl

use strict;
use warnings;
use DBI;
use Getopt::Long;
use File::Basename;
use Cwd;
use File::Temp qw/tempfile tempdir/;
use File::Path qw/mkpath/;
use Data::Dumper;
use FindBin qw($Bin);
use lib "$Bin";
use CommonLibTASUKE;

my $DB      = '';
my $HOST    = 'localhost';
my $USER    = '';
my $PASS    = '';
#my $TASUKE_WWW_DIR = '';
my $REMOVE = 0;

# Checking parameters and database connection
GetOptions( 'db|database=s' => \$DB,
            'h|host=s'      => \$HOST,
            'u|user=s'      => \$USER,
            'r|remove'    => \$REMOVE,
#	    't|tasuke=s'    => \$TASUKE_WWW_DIR,
);

#if(!$TASUKE_WWW_DIR || !-d $TASUKE_WWW_DIR){
#	die "The specified tasuke_www dir is invalid.\n";
#}

system("stty -echo");
print "Enter db password: ";
chomp($PASS = <STDIN>);
print "\n";
system("stty echo");

my $connect = DBI->connect('DBI:mysql:database='.$DB.';host='.$HOST,$USER,$PASS)
  or die "Can't connect to database: $DBI::errstr\n";
#&delete_phenotype_tables();
_initTable();
_checkUpdate();
my $orderFlg = _checkOrder();
my $dbs = "";
my $accs = "";
my @DBary;
my @ACCary;
if ($orderFlg) {
	($dbs,$accs) = _collectDBs();
	@DBary = split(/,/,$dbs);
	@ACCary = split(/,/,$accs);
}

my $commonMysql = new CommonMysql();
my $MODE_IS_MULTIVCF = $commonMysql->modeIsPlusMultiVcf($connect);

my $phenotype_hash_ref = &make_phenotype_hash_from_csv();
my $phenotype_id_hash_ref = &set_phenotype_to_db($phenotype_hash_ref);
my $snp_site_hash_hash_ref = &make_snp_site_hash_hash($phenotype_hash_ref);
&set_phenotype_prob_table_and_make_plot_images($phenotype_hash_ref, $phenotype_id_hash_ref, $snp_site_hash_hash_ref);
print "clearing maps...\n";
undef $phenotype_hash_ref;
undef $phenotype_id_hash_ref;
undef $snp_site_hash_hash_ref;
print "done.\n";
&set_phenotype_gene_table();
#&set_phenotype_paper_table($phenotype_id_hash_ref);
#&set_phenotype_paper_table();
print "All done.\n";

sub _initTable {
	if ($REMOVE) {
        	&exec_sql("DROP TABLE IF EXISTS sys_phenotype;");
        	&exec_sql("DROP TABLE IF EXISTS sys_phenotype_expression;");
        	&exec_sql("DROP TABLE IF EXISTS sys_phenotype_gene;");
        	&exec_sql("DROP TABLE IF EXISTS sys_phenotype_prob;");
        	&exec_sql("DROP TABLE IF EXISTS id2eff;");
        print "Done Remove Data!!\n";
        exit;
	} elsif ( ! defined _check_exists_table("sys_phenotype") ) {
                &exec_sql("CREATE TABLE `sys_phenotype` (
                        `phenotype_id` varchar(40) NOT NULL,
                        `phenotype` varchar(40) NOT NULL,
                        PRIMARY KEY (`phenotype_id`)
                        )");
                &exec_sql("CREATE TABLE `sys_phenotype_expression` (
                        `accession_id` varchar(40) NOT NULL,
                        `phenotype_id` varchar(40) NOT NULL,
                        `expression` varchar(100) NOT NULL,
                        PRIMARY KEY (`accession_id`,`phenotype_id`)
                        )");
                &exec_sql("CREATE TABLE `sys_phenotype_prob` (
                        `id` bigint(20) unsigned NOT NULL AUTO_INCREMENT,
                        `phenotype_id` varchar(40) NOT NULL DEFAULT '',
                        `chr` varchar(50) NOT NULL DEFAULT '',
                        `pos` int(11) NOT NULL DEFAULT '0',
                        `prob` float DEFAULT NULL,
                        PRIMARY KEY (`id`),
                        UNIQUE KEY `id` (`id`),
                        KEY `phenotype_id` (`phenotype_id`,`chr`,`pos`),
                        KEY `pos_idx` (`pos`),
                        KEY `chr_pos_idx` (`chr`,`pos`)
                        )");
                &exec_sql("CREATE TABLE `sys_phenotype_gene` (
                        `id` bigint(20) unsigned NOT NULL,
                        `phenotype_prob_id` int(10) NOT NULL,
                        `phenotype_id` varchar(40) DEFAULT NULL,
                        `chr` varchar(50) DEFAULT NULL,
                        `pos` int(11) DEFAULT NULL,
                        `prob` float DEFAULT NULL,
                        `Effect` text DEFAULT NULL,
                        `note` text DEFAULT NULL,
                        UNIQUE KEY `id` (`id`),
                        KEY `phenotype_id` (`phenotype_id`),
                        KEY `prob` (`prob`)
                        )");
                &exec_sql("CREATE TABLE `id2eff` (
                        `id` bigint(20) NOT NULL,
                        `eff_grp` smallint(6) NOT NULL,
                        PRIMARY KEY (`id`,`eff_grp`)
                        )");
=pod
        	&exec_sql("UNLOCK TABLES;");
		print "Done Remove Data!!\n";
		exit;
	        &exec_sql("DROP TABLE IF EXISTS sys_phenotype_paper");
	        &exec_sql("CREATE TABLE `sys_phenotype_paper` (
                `accession_id` varchar(40) NOT NULL DEFAULT '',
                `phenotype_id` varchar(40) NOT NULL,
                `PMID` int(10) NOT NULL,
                `title` text NOT NULL,
                `author` text NOT NULL,
                `abstract` text,
                PRIMARY KEY (`accession_id`,`phenotype_id`,`PMID`)
                )");
=cut
	}
}

sub _checkUpdate {
	##check whether data is installed##
	my $valupd = $connect->selectrow_array("select * from sys_phenotype");
	if (defined $valupd) {
		##check if data is updated new one.##
		print "Caution: The phenotype data you already installed will update.\n";
		print "Do you want to update phenotype value?[y/N]\n#";
		while(defined(my $line = <STDIN>)){
                	$line =~ s/[\r\n]+\z//;
                        if($line eq "y" || $line eq "Y"){
				return 1;
                        }elsif($line eq "n" || $line eq "N" || $line eq ""){
      				print "Installation was stopped.\n";
                             	exit;
                       	}else{
                        	print "Caution: The phenotype data you already installed will update.\n";
                                print "Do you want to update phenotype value?[y/N]\n#";
                        }
                }
        }
}
sub _checkOrder {
	print "Do you use selected accession(require order.conf using TASUKE)?[y/N]\n#";
	while(defined(my $line = <STDIN>)){
		$line =~ s/[\r\n]+\z//;
		if($line eq "y" || $line eq "Y"){
			my $yFlg;
			if (-e "order.conf") {
				if ( -s "order.conf" > 0) {
                			$yFlg = 1;
				}
			}
			if (!$yFlg) {
	        	        print "Caution: 'order.conf' is not found.\n";
		                print "Use all accessions.\n";
				$yFlg = 0;
			}
			return $yFlg;
                }else {
	                print "Caution: 'order.conf' is not found.\n";
         	        print "Use all accessions.\n";	
                	return 0;
                }
        }
}

sub _collectDBs {
	if (-e "order.conf") {
        	my $dbFlg = 1;
        	my $accFlg = 1;
        	open my $IN, "<order.conf";
        	while (<$IN>) {
                	my $line = $_;
                	$line =~ s/[\r\n]//g;
                	my $query = "SELECT accession_db from sys_accession where accession_id = '$line'";
                	my ($sth, $rc) = &exec_sql($query);
                	while(my @row = $sth->fetchrow){
                       		if ($dbFlg) {
                        	      	$dbs = "$row[0]";
                        	       	$dbFlg = 0;
                        	} else {
                        	       	$dbs .= ",$row[0]";
                        	}
			}
			if ($accFlg) {
                                $accs = "$line";
                                $accFlg = 0;
                        } else {
                                $accs .= ",$line";
                        }
                }
        	return $dbs,$accs;
	}

}

sub make_phenotype_hash_from_csv(){
    my $file_path;

	print "> Where is phenotype data csv file?\n";
	#print "(File format: Breed Name,Phenotype,[1|0])\n# ";
	print "(File format: Breed Name,Phenotype,Phenotype Value)\n# ";
	while(defined(my $line = <STDIN>)){
    		$line =~ s/[\r\n]+\z//;
		$line =~ s/^\s*(.*?)\s*$/$1/;
	     	$file_path = $line;
     		if($file_path ne "" && -e $file_path) {
			if (&check_csv($file_path, 3)) {
         			print "wrong format.\n";
			}else{
         			last;
			}
     		}elsif($file_path ne "" && !(-e $file_path)) {
			print "File does not exists.\n";
		}
     		print "> Where is phenotype data csv file?\n# ";
	}
	
	my %phenotype_exists_hash = ();
	my %phenotype_hash = ();
	open my $fh, '<', $file_path or die "file open error: $!";
	while(<$fh>) {
		$_ =~ s/[\r\n]+\z//;
		my @line = split(/,/, $_);
		$phenotype_exists_hash{$line[0] . "," . $line[1]} = $line[2];
		$phenotype_hash{$line[1]} = 1;
	}
	close $fh;
	
	my $sql = "SELECT accession_id FROM sys_accession";
	my ($sth, $rc) = &exec_sql($sql);
	while(my @row = $sth->fetchrow){
		for my $phenotype (keys %phenotype_hash){
			my $key = $row[0] . "," . $phenotype;
			if(!exists $phenotype_exists_hash{$key}){
				$phenotype_exists_hash{$key} = 0;
			}
		}
	}
	
	return \%phenotype_exists_hash;
}

sub check_csv(){
	my $file_path = shift;
	my $columns_size = shift;
	my $err_flg = 0;
	open my $fh, '<', $file_path or die "file open error: $!";
	while (<$fh>) {
		$_ =~ s/[\r\n]+\z//;
		my @line = split(/,/, $_, -1);
		$err_flg = 1 if @line != $columns_size; 
	}
	close $fh;
	return $err_flg;
}

sub check_tsv(){
	my $file_path = shift;
	my $columns_size = shift;
	my $err_flg = 0;
	open my $fh, '<', $file_path or die "file open error: $!";
	while (<$fh>) {
		$_ =~ s/[\r\n]+\z//;
		my @line = split(/\t/, $_, -1);
		$err_flg = 1 if @line != $columns_size; 
	}
	close $fh;
	return $err_flg;
}

sub set_phenotype_to_db(){
	my $phenotype_exists_hash_ref = shift;
	
	my %phenotype_id_hash = ();
	my $id_count = 0;
	my $sql = "SELECT count(*) FROM sys_phenotype";
	my ($sth, $rc) = &exec_sql($sql);
	while(my @row = $sth->fetchrow){
		$id_count = $row[0]
	}
	for my $key (keys(%$phenotype_exists_hash_ref)){
		my ($accession_id, $phenotype) = split(/,/, $key);
		if(!exists $phenotype_id_hash{$phenotype}){
			my $sql = $connect->selectrow_array("SELECT * FROM sys_phenotype where phenotype = '$phenotype'");
                        if(! defined $sql) {
				++$id_count;
				my $id_str = "phenotype";
				if($id_count < 10){
					$id_str .= "0$id_count";
				}else{
					$id_str .= $id_count;
				}
				$phenotype_id_hash{$phenotype} = $id_str;
                        	&exec_sql("LOCK TABLES sys_phenotype WRITE;");
                        	my @sql_params = ($phenotype_id_hash{$phenotype}, $phenotype);
                        	&exec_sql("INSERT INTO sys_phenotype (phenotype_id, phenotype) VALUES(?, ?)", @sql_params);
                        	&exec_sql("UNLOCK TABLES;");
			} else {	
				my $sql = "SELECT phenotype_id FROM sys_phenotype where phenotype = ?";
                        	my ($sth, $rc) = &exec_sql($sql, $phenotype);
                        	if ($rc eq '0E0') {
                               		next;
                        	}else{
                               		my @row = $sth->fetchrow;
					$phenotype_id_hash{$phenotype} = $row[0];
				}
			}
		}
	}
	for my $key (keys(%$phenotype_exists_hash_ref)){
		my ($accession_id, $phenotype) = split(/,/, $key);
		my @sql_params = ($accession_id, $phenotype_id_hash{$phenotype}, $$phenotype_exists_hash_ref{$key});
		my $valexp = $connect->selectrow_array("select * from sys_phenotype_expression where accession_id='$accession_id' and phenotype_id='$phenotype_id_hash{$phenotype}'");
		if(! defined $valexp) {
			&exec_sql("LOCK TABLES sys_phenotype_expression WRITE;");
			&exec_sql("INSERT INTO sys_phenotype_expression (accession_id, phenotype_id, expression) VALUES(?, ?, ?)", @sql_params);
        		&exec_sql("UNLOCK TABLES;");
		} else {
			&exec_sql("LOCK TABLES sys_phenotype_expression WRITE;");
            &exec_sql("UPDATE sys_phenotype_expression SET expression='$$phenotype_exists_hash_ref{$key}' where accession_id='$accession_id' and phenotype_id='$phenotype_id_hash{$phenotype}'");
            &exec_sql("UNLOCK TABLES;");
		}
	}
	
	return \%phenotype_id_hash;
}

sub exec_sql() {
	my ($sql, @data) = @_;
	my $rc = 0;
	
	my $sth = $connect->prepare($sql);
	$rc = $sth->execute(@data);
	die "db error: " . $sth->errstr unless $rc;

	return ($sth, $rc);
}

sub make_snp_site_hash_hash(){
	my $phenotype_hash_ref = shift;
	
	my %accession_id_hash = ();
	for my $key (keys(%$phenotype_hash_ref)){
		my ($accession_id, $phenotype) = split(/,/, $key);
		$accession_id_hash{$accession_id} = 1;
	}
		
	my %accession_db_hash = ();
	for my $accession_id (keys(%accession_id_hash)){
		my $sql = "SELECT accession_db FROM sys_accession WHERE accession_id = ?";
		my ($sth, $rc) = &exec_sql($sql, $accession_id);
		if ($rc eq '0E0') {
			next;	
		}else{
			my @row = $sth->fetchrow;
			$accession_db_hash{$accession_id} = $row[0];
		}
	}
	
	print "making a snp sites map...\n";
	my %snp_site_hash_hash = ();
	my @accs;
	if ($dbs eq "") {
		@accs = keys %accession_db_hash;
	} else {
		@accs = @ACCary;	# Deep copy
	}
	foreach my $accession_id (@accs) {
		my $sql = "SELECT chr,pos FROM ". $accession_db_hash{$accession_id}. "_variant WHERE variant = 's' AND ". getGtIsNotZeroSql();
		my ($sth, $rc) = &exec_sql($sql);
		if ($rc eq '0E0') {
			next;
		}else{
			my %snp_site_hash = ();
			while (my @row = $sth->fetchrow) {
				$snp_site_hash{$row[0] . ":" . $row[1]} = 1;
			}
			$snp_site_hash_hash{$accession_id} = \%snp_site_hash;
		}
	}
	print "done.\n";
	
	return \%snp_site_hash_hash;	
}

sub set_phenotype_prob_table_and_make_plot_images(){
	my ($phenotype_exists_hash_ref, $phenotype_id_hash_ref, $snp_site_hash_hash_ref) = @_;
	my %chr_size_hash = ();
	my @chr_list = ();
	my $sql = "SELECT id,size FROM sys_chromosome";
	my ($sth, $rc) = &exec_sql($sql);
	if ($rc eq '0E0') {
		next;	
	}else{
		my $chr_count = 1;
		while (my @row = $sth->fetchrow) {
			$chr_size_hash{$chr_count} = $row[1];
			push @chr_list, $row[0];
			$chr_count++;
		}
	}
	&set_prob_table_and_make_plots_from_qqman_file($phenotype_id_hash_ref, \@chr_list, \%chr_size_hash);
}

sub set_prob_table_and_make_plots_from_qqman_file(){
	my ($phenotype_id_hash_ref, $chr_list_ref, $chr_size_hash_ref) = @_;
	
	my $qqman_file;
	my $phenotype;
	print "> Where is qqman output file?\n# ";
	while(defined(my $line = <STDIN>)){
    		$line =~ s/[\r\n]+\z//;
		$line =~ s/^\s*(.*?)\s*$/$1/;
    		$qqman_file = $line;
    		if($qqman_file ne "" && -e $qqman_file) {
    			my $sys_phenotype_prob_insert_count = 0;
     		my $sys_phenotype_prob_insert_value = "";
        		print "> What is a phenotype of the file?\n# ";
        		while(defined($line = <STDIN>)){
        			$line =~ s/[\r\n]+\z//;
				$line =~ s/^\s*(.*?)\s*$/$1/;
     			$phenotype = $line;
     			if(!exists $$phenotype_id_hash_ref{$line}){
					print "Invalid phenotype: $line\n";
					print "> What is a phenotype of the file?\n# ";
					next;
				}
     			&exec_sql("LOCK TABLES sys_phenotype_prob WRITE;");	
     			open my $qqman_fh, '<', $qqman_file or die "file open error: $!";
     			my ($tempfh, $tempfile) = tempfile;
     			<$qqman_fh>;
     			my $snp_count = 0;
     			my $chr_count = 0;
     			my $prev_chr = "";  			
     			my ($sth, $rc) = &exec_sql("select chr,pos from sys_phenotype_prob where phenotype_id='$$phenotype_id_hash_ref{$phenotype}'");
     			my %chrpos_map = ();
     			if($rc ne '0E0'){
     				while(my @row = $sth->fetchrow){
     					$chrpos_map{"$row[0]:$row[1]"} = 1;
     				}
     			}
     			while (<$qqman_fh>) {
					$_ =~ s/[\r\n]+\z//;
					my @line_elements = split(/\t/, $_, -1);
					my $minus_log10_p_value = -log($line_elements[2]) / log(10);
					if(!exists $chrpos_map{"$line_elements[0]:$line_elements[1]"}) {
						if($sys_phenotype_prob_insert_value ne ""){
							$sys_phenotype_prob_insert_value .= ",";
						}
						$sys_phenotype_prob_insert_value .= "('$$phenotype_id_hash_ref{$phenotype}', '$line_elements[0]', $line_elements[1], $minus_log10_p_value)";
						$sys_phenotype_prob_insert_count++;
						if($sys_phenotype_prob_insert_count == 1000){
							&exec_sql("INSERT INTO sys_phenotype_prob (phenotype_id, chr, pos, prob) VALUES $sys_phenotype_prob_insert_value");
							$sys_phenotype_prob_insert_count = 0;
							$sys_phenotype_prob_insert_value = "";
						}
						$chrpos_map{"$line_elements[0]:$line_elements[1]"} = 1;
     					} else {
						&exec_sql("UPDATE sys_phenotype_prob SET prob='$minus_log10_p_value' where phenotype_id='$$phenotype_id_hash_ref{$phenotype}' and chr='$line_elements[0]' and pos='$line_elements[1]'");
						&exec_sql("LOCK TABLES sys_phenotype_gene WRITE;");
                    				&exec_sql("UPDATE sys_phenotype_gene SET prob='$minus_log10_p_value' where phenotype_id='$$phenotype_id_hash_ref{$phenotype}' and chr='$line_elements[0]' and pos='$line_elements[1]'");
						&exec_sql("UNLOCK TABLES;");
					}	
     				$snp_count++;
     				if($prev_chr ne $line_elements[0]){
     					$prev_chr = $line_elements[0];
     					$chr_count++;
     				}
     				print $tempfh "snp$snp_count\t$chr_count\t$line_elements[1]\t$minus_log10_p_value\n";
     			}
     			close $tempfh;
     			close $qqman_fh;
     			&exec_sql("UNLOCK TABLES;");
         		last;
         	}
         	if($sys_phenotype_prob_insert_value ne ""){
         		&exec_sql("INSERT INTO sys_phenotype_prob (phenotype_id, chr, pos, prob) VALUES $sys_phenotype_prob_insert_value");
         	}
		#&make_manhattan_plot_images($qqman_file, $phenotype, $chr_list_ref, $chr_size_hash_ref);
         	print "> Completed.\n> Read the next file? [y/N]\n# ";
         	my $is_last = 0;
         	while(defined($line = <STDIN>)){
         		$line =~ s/[\r\n]+\z//;
				$line =~ s/^\s*(.*?)\s*$/$1/;
				if($line eq "y" || $line eq "Y"){
					print "> Where is qqman output file?\n# ";
					last;
				}elsif($line eq "n" || $line eq "N" || $line eq ""){
					$is_last = 1;
					last;
				}else{
					print "Read the next file? [y/N]\n# ";
				}
         	}
         	if($is_last == 1){
         		last;
         	}
     	}elsif($qqman_file eq "" || !(-e $qqman_file)) {
			print "File does not exists.\n";
			print "> Where is qqman output file?\n# ";
		}
	}
}
=head
sub make_manhattan_plot_images(){
        my $p_values_file = shift;
        my $phenotype = shift;
        my $chr_list_ref = shift;
        my $chr_size_hash_ref = shift;

        my $plot_file = "plot_$phenotype.png";
        $plot_file =~ s/\s+/_/g;
        my $out_dir = "$TASUKE_WWW_DIR/images/plot";
        if(!-d $out_dir){
                mkpath($out_dir);
        }
        my $r_file = &make_manhattan_plot_r_file($p_values_file, $phenotype, "$out_dir/$plot_file", $chr_list_ref, $chr_size_hash_ref);
        system("Rscript --vanilla --slave $r_file");
}
=cut
sub make_manhattan_plot_r_file(){
        my $p_values_file = shift;
        my $phenotype = shift;
        my $plot_file = shift;
        my $chr_list_ref = shift;
        my $chr_size_hash_ref = shift;
        
	my %chr_pos_hash = ();
        my $sql = "SELECT id,size FROM sys_chromosome";
	my ($sth, $rc) = &exec_sql($sql);
	my $chr_name_vector_str = "chrName <- c(";
	my $chr_names_vector_str = "names(Pos) <- c(";
	while(my @row = $sth->fetchrow){
		$chr_pos_hash{$row[0]} = $row[1];
		$chr_name_vector_str .= "\"$row[0]\",";
		$chr_names_vector_str .= "\"$row[0]\",";
	}
	my $all_chr_num = keys %chr_pos_hash;
	chop($chr_name_vector_str);
	chop($chr_names_vector_str);
	$chr_name_vector_str .= ")";
	$chr_names_vector_str .= ")";
	print "$chr_name_vector_str\n";

	my $chr_size_vector_str = "chrSize <- c(";
	my $pos_vector_str = "";
	for my $key (keys %chr_pos_hash) {
		$chr_size_vector_str .= $chr_pos_hash{$key} . ",";
		$pos_vector_str .= "Pos[\"$key\"] <- $chr_pos_hash{$key}\;";
	}
	chop($chr_size_vector_str);
	$chr_size_vector_str .= ")";
	chop($pos_vector_str);
	print "$chr_size_vector_str\n";
	print "$pos_vector_str\n";

        my $script = <<"EOF";
                rcMpPvals <- read.table('$p_values_file', sep="\t",header=T)
		colnames(rcMpPvals) <- c('chr', 'pos', 'logp')
		chrNum=$all_chr_num
                $chr_size_vector_str
		$chr_name_vector_str
		rcMpPvals<-rcMpPvals[order(rcMpPvals\$chr),]

		chrBeginPos <- vector(length=chrNum)
		chrBeginPos[1] <- 1
		chrEndPos <- vector(length=chrNum)

		Pos <- vector("list",length=$all_chr_num)
		$chr_names_vector_str
		$pos_vector_str
		for(i in 1:chrNum) {
			if (!is.na(rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i+1])),"pos"][1])) {
		 		 rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i+1])),"pos"]<-rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i+1])),"pos"]+max(rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i])),"pos"])
			}
		}
		for(i in 1:chrNum){
			rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i])), "Tick"] <- (min(rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i])),"pos"]) + max(rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i])),"pos"]))/2
		}
		rcMpPvals\$log10p <- - log10(rcMpPvals\$logp)
		for (i in 1:chrNum) {
		  	rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i])), "Color_Discovery"] <- i %% 2
		  	rcMpPvals[rcMpPvals\$chr==as.character(names(Pos[i])), "Color_Replication"] <- (i+1) %% 2
		}

		bpMidVec <- vector(length=chrNum)
		for (i in 1:chrNum){
			bpMidVec[i] <- ((chrEndPos[i] - chrBeginPos[i] + 1)/2) + chrBeginPos[i] - 1
		}
		library(ggplot2)
		library(grid)
		p <- ggplot(rcMpPvals,aes(pos,log10p))
		p <- p + geom_point(size=1,aes(colour = as.factor(Color_Discovery)))
		p <- p + scale_color_manual(values=rep(c('#87cefa', '#6495ed'), 6))
		p <- p + theme_bw(base_size=15)
		p <- p + scale_x_continuous(labels=as.character(names(Pos[1:chrNum])), breaks=unique(rcMpPvals\$Tick), expand=c(0,0))
		p <- p + ggtitle('$phenotype') + xlab('chromosome') + ylab('-log10(P)')
		p <- p + theme(legend.position='none', plot.title=element_text(hjust=0.5), text=element_text(size=15))
		#, panel.spacing=unit(c(0, 0, 0, 0), "mm"), plot.margin=unit(c(0, 0, 0, 0), "mm"), axis.title.x = element_blank())
		ggsave(file='$plot_file', dpi = 100, width= 8, height = 4)

EOF
        my ($tempfh, $tempfile) = tempfile;
        print $tempfh $script;
        close($tempfh);

        return $tempfile;
}
sub make_fisher_r_file(){
	my $fisher_counts_file = shift;
	
	my $script = <<"EOF";
z <- read.csv("$fisher_counts_file", header=FALSE)
f <- function(x) fisher.test(matrix(c(x),nrow=2))
y <- apply(z,1,f)
print(y)
EOF
	my ($tempfh, $tempfile) = tempfile;
	print $tempfh $script;
	close($tempfh);
	
	return $tempfile;
}

sub make_p_values_hash_by_r(){
	my $counts_file = shift;
	
	my @n11_n12_list = ();
	open my $fh, '<', $counts_file or die "file open error: $!";
	while(<$fh>){
		$_ =~ s/[\r\n]+\z//;
		my @line_elements = split(/,/);
		push(@n11_n12_list, $line_elements[0] . ',' . $line_elements[2]);
	}
	close($fh);
	
	my $r_script_file = &make_fisher_r_file($counts_file);
	my ($tempfh, $tempfile) = tempfile;
	system("Rscript --vanilla --slave $r_script_file > $tempfile");
	#system("/usr/local/bin/Rscript --vanilla --slave $r_script_file > $tempfile");
	my %p_values_hash;
	my $n11_n12_list_pos = 0;
	open $fh, '<', $tempfile or die "file open error: $!";
	while (<$fh>){
		$_ =~ s/[\r\n]+\z//;
		if(/^p-value = (\S+)$/){
			$p_values_hash{$n11_n12_list[$n11_n12_list_pos]} = $1;
			$n11_n12_list_pos++;	
		}else{
			next;
		}
	}
	close($fh);
	
	return \%p_values_hash;
}

sub set_phenotype_gene_table(){
	my ($effname_ref,$effgrp_ref) = &set_phenotype_gene_snpeff();
	my %effinfos;
	foreach my $key(keys(%$effname_ref)) {
        	my $chrpos = $key;
        	my $flg = 1;
        	foreach my $key2(keys(%{$$effname_ref{$chrpos}})) {
        	    my $idname =$key2;
        	    #print $OUT "$key:$effname_ref{$key}\n";
        	    if($flg) {
        	    	$effinfos{$chrpos} = "$idname|$$effname_ref{$chrpos}{$idname}";
              		$flg=0;
            	    } else {
                	$effinfos{$chrpos} .= "---DIFFID---$idname|$$effname_ref{$chrpos}{$idname}";
            	    }
        	}
	}
	undef $effname_ref;
    	print "registering infos in the DB...\n";
	my %posflg;
	my $sys_phenotype_gene_insert_count_1 = 0;
	my $sys_phenotype_gene_insert_value_1 = "";
	my $sys_phenotype_gene_insert_count_2 = 0;
	my $sys_phenotype_gene_insert_value_2 = "";
	my $id2eff_insert_count = 0;
	my $id2eff_insert_value = "";
	my $sys_phenotype_gene_id = 0;	
	my %valphenoflg = ();
	my %sys_phenotype_gene_chrpos_flg = ();
	my ($sth,$rc) = &exec_sql("select id,phenotype_id,chr,pos from sys_phenotype_gene");
	if($rc ne '0E0'){
		while (my @rows = $sth->fetchrow_array()) {
			$valphenoflg{"$rows[1]:$rows[2]:$rows[3]"} = 1;
			$sys_phenotype_gene_chrpos_flg{"$rows[2]:$rows[3]"} = 1;
			if($sys_phenotype_gene_id < $rows[0]){
				$sys_phenotype_gene_id = $rows[0];
			}
		}
	}
	my %id2eff_id_flg = ();
	($sth,$rc) = &exec_sql("select id from id2eff");
	if($rc ne '0E0'){
		while (my @rows = $sth->fetchrow_array()) {
			$id2eff_id_flg{$rows[0]} = 1;
		}
	}
	#for my $type ("five_prime_UTR", "CDS", "three_prime_UTR", "exon"){
	#($sth, $rc) = &exec_sql("SELECT a.phenotype_prob_id, a.phenotype_id, a.chr, a.pos, a.prob FROM tmp_table2 a, sys_track b WHERE a.gene_id = b.id");
	($sth, $rc) = &exec_sql("SELECT id, phenotype_id, chr, pos, prob FROM sys_phenotype_prob");
	while (my @rows = $sth->fetchrow_array()) {
		my $posinfo = "$rows[2]:$rows[3]";
		my $valphenoinfo = "$rows[1]:$rows[2]:$rows[3]";
		if(! exists $valphenoflg{$valphenoinfo}) {	
			my @effgrps;
			my $effgrp;
			if (defined $$effgrp_ref{$posinfo}) {
				@effgrps = split(/,/,$$effgrp_ref{$posinfo});
				@effgrps = sort {$b <=> $a} @effgrps;
				$effgrp = join(",",@effgrps);
				$sys_phenotype_gene_id++;					
				if($sys_phenotype_gene_insert_value_1 ne ""){
					$sys_phenotype_gene_insert_value_1 .= ",";
				}
				$sys_phenotype_gene_insert_value_1 .= "($sys_phenotype_gene_id, '$rows[0]','$rows[1]','$rows[2]','$rows[3]','$rows[4]','$effinfos{$posinfo}','$effgrp')";
				$sys_phenotype_gene_insert_count_1++;
				if($sys_phenotype_gene_insert_count_1 == 1000){
					&exec_sql("INSERT INTO sys_phenotype_gene (id, phenotype_prob_id, phenotype_id, chr, pos, prob, Effect,note) values $sys_phenotype_gene_insert_value_1");
					$sys_phenotype_gene_insert_count_1 = 0;
					$sys_phenotype_gene_insert_value_1 = "";
				}
				my @id2eff_insert_rows = ($sys_phenotype_gene_id, $effgrp);
				my @id2eff_insert_rows_ref_array = ();
				push(@id2eff_insert_rows_ref_array, \@id2eff_insert_rows);
				if(exists $sys_phenotype_gene_chrpos_flg{"$rows[2]:$rows[3]"}){
					my ($sth2,$rc2) = &exec_sql("select id,note from sys_phenotype_gene where chr='$rows[2]' and pos='$rows[3]'");
        				while (my @rows2 = $sth2->fetchrow_array()) {
        					push(@id2eff_insert_rows_ref_array, \@rows2);
        				}
				}
				for my $id2eff_insert_rows_ref (@id2eff_insert_rows_ref_array){
					my $id2eff_id = $$id2eff_insert_rows_ref[0];
					if(!exists $id2eff_id_flg{$id2eff_id}){
						my $id2eff_effgrp = $$id2eff_insert_rows_ref[1];
	                			my @effnums = split(/,/,$id2eff_effgrp);
	                			for my $effnum (@effnums) {
	                				if($id2eff_insert_value ne ""){
								$id2eff_insert_value .= ",";
							}
							$id2eff_insert_value .= "('$sys_phenotype_gene_id','$effnum')";
							$id2eff_insert_count++;
							if($id2eff_insert_count == 1000){
								&exec_sql("insert into id2eff (id,eff_grp) values $id2eff_insert_value");
								$id2eff_insert_count = 0;
								$id2eff_insert_value = "";
							}
			                	}
        					$id2eff_id_flg{$id2eff_id} = 1;
					}
				}
			} else {
				$effgrp = "";
				$sys_phenotype_gene_id++;
				if($sys_phenotype_gene_insert_value_2 ne ""){
					$sys_phenotype_gene_insert_value_2 .= ",";
				}
				$sys_phenotype_gene_insert_value_2 .= "($sys_phenotype_gene_id, '$rows[0]','$rows[1]','$rows[2]','$rows[3]','$rows[4]')";
				$sys_phenotype_gene_insert_count_2++;
				if($sys_phenotype_gene_insert_count_2 == 1000){
					&exec_sql("INSERT INTO sys_phenotype_gene (id, phenotype_prob_id, phenotype_id, chr, pos, prob) values $sys_phenotype_gene_insert_value_2");
					$sys_phenotype_gene_insert_count_2 = 0;
					$sys_phenotype_gene_insert_value_2 = "";
				}
			}
			$valphenoflg{$valphenoinfo} = 1;
		}
	}
	if($sys_phenotype_gene_insert_value_1 ne ""){
		&exec_sql("INSERT INTO sys_phenotype_gene (id, phenotype_prob_id, phenotype_id, chr, pos, prob, Effect,note) values $sys_phenotype_gene_insert_value_1");
	}
	if($sys_phenotype_gene_insert_value_2 ne ""){
		&exec_sql("INSERT INTO sys_phenotype_gene (id, phenotype_prob_id, phenotype_id, chr, pos, prob) values $sys_phenotype_gene_insert_value_2");
	}
	if($id2eff_insert_value ne ""){
		&exec_sql("insert into id2eff (id,eff_grp) values $id2eff_insert_value");
	}
	print "done.\n";	
}
sub set_phenotype_paper_table(){
        print "> Import papers data? [y/N]\n# ";
        while(defined(my $line = <STDIN>)){
                $line =~ s/[\r\n]+\z//;
                $line =~ s/^\s*(.*?)\s*$/$1/;
                my $answer = $line;
                if($answer eq 'Y' || $answer eq 'y') {
                        &insert_paperdata($phenotype_id_hash_ref);
                        return;
                } elsif($answer eq 'N' || $answer eq 'n' || $answer eq '') {
                        return;
                }else{
                        print "> Import papers data? [y/N]\n# ";
                }
        }
}

sub insert_paperdata() {
        my $phenotype_id_hash_ref = shift;
        my $file_path;

        print "> Where is papers data tsv file?\n";
        print "(File format: Accession id[TAB]Phenotype[TAB]PMID[TAB]Title[TAB]Author[TAB]Abstract) # ";
        while(defined(my $line = <STDIN>)){
                $line =~ s/[\r\n]+\z//;
                $line =~ s/^\s*(.*?)\s*$/$1/;
                $file_path = $line;
                if($file_path ne "" && -e $file_path) {
                        if(&check_tsv($file_path, 6)) {
                                print "wrong format.\n";
                        }else{
                                last;
                        }
                }elsif($file_path ne "" && !(-e $file_path)) {
                        print "File does not exists.\n";
                } else {
                print "> Where is papers data tsv file?\n# ";
                }
        }
        my @paper_info_array_array = ();
        open my $fh, '<', $file_path or die "file open error: $!";
        while(<$fh>) {
                $_ =~ s/[\r\n]+\z//;
                my @line = split(/\t/, $_);
                if(!exists $$phenotype_id_hash_ref{$line[1]}){
                        die "Invalid phenotype: " . $line[1];
                }else{
                        $line[1] = $$phenotype_id_hash_ref{$line[1]};
                        push(@paper_info_array_array, \@line)
                }
        }
        close $fh;

        &exec_sql("LOCK TABLES sys_phenotype_paper WRITE;");
        for my $paper_info_array_ref (@paper_info_array_array){
                my $query = "insert into sys_phenotype_paper values (?, ?, ?, ?, ?, ?)";
                &exec_sql($query, @$paper_info_array_ref);
        }
        &exec_sql("UNLOCK TABLES;");
}

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

sub set_phenotype_gene_snpeff(){
        print "making a eff infos map...\n";
        my $query = "select accession_id, accession_db, accession_name from sys_accession";
        my ($sth, $rc) = &exec_sql($query);
        my %accession_id_map = ();
        my %accession_name_map = ();
        while (my @rows = $sth->fetchrow_array()) {
                $accession_id_map{$rows[1]} = $rows[0];
                $accession_name_map{$rows[1]} = $rows[2];
        }

        $query = "select eff_name, eff_grp from sys_snpEff";
        ($sth, $rc) = &exec_sql($query);
        my %eff_grp_map = ();
        while (my @rows = $sth->fetchrow_array()) {
                $eff_grp_map{$rows[0]} = $rows[1];
        }

        my %effname = ();
        my %effgrp = ();
        my $loop_count = 0;
        my $accession_db;
	my @accs;
	if ($dbs eq "") {
		@accs = keys %accession_id_map;
	} else {
		@accs = @DBary;	# Deep copy
	}
	my $accsNum = $#accs + 1;
	foreach my $accession_db (@accs) {
		$loop_count++;
		if($loop_count % 10 == 0){
			printf "%d/%d...\n",  $loop_count, $accsNum;
		}
		if($MODE_IS_MULTIVCF) {
			$query = "select a.chr,a.pos,s.EFF from sys_variant s,". $accession_db. "_variant a,sys_phenotype_prob b ".
				"where s.chr=a.chr and s.pos_vcf=a.pos_vcf and a.chr = b.chr and a.pos=b.pos and ". getGtIsNotZeroSql();
		} else {
			$query = "select a.chr, a.pos, a.EFF from ". $accession_db. "_variant a, sys_phenotype_prob b where a.chr = b.chr and a.pos = b.pos and ". getGtIsNotZeroSql();
		}
		($sth, $rc) = &exec_sql($query);
		while (my @rows = $sth->fetchrow_array()) {
			my @effs = split(/,/, $rows[2]);
			for my $eff (@effs){
				if($eff =~ /(\S+)\((\S+)\)/){
					setEff(\%effname, \%effgrp, \%eff_grp_map, $rows[0], $rows[1], $eff);
				} else {
					setAnn(\%effname, \%effgrp, \%eff_grp_map, $rows[0], $rows[1], $eff);
				}
			}
		}
	}
	print "done.\n";
	return \%effname,\%effgrp;
}
sub setAnn {
	my ($effname, $effgrp, $eff_grp_map, $chr, $pos, $eff) = @_;

	return if($eff =~ m/[\(\)]/);

	my @rows = ($chr, $pos);
	my @eff_elements = split(/\|/, $eff);
	my $eff_ontology = $eff_elements[1];
#	my $eff_ontology = lc($eff_elements[1]);
		$eff_ontology =~ s/\s//g;

	my $eff_name = getMaxGrpAnn($eff_grp_map, $eff_ontology);
	if(!defined $eff_name) {
		printf STDERR "Undefined Ontology:%s\n", $eff_ontology;
		return;
	}
	my $eff_grp_map_key = (split(/\+/, $eff_name))[0];
	if(!exists $eff_grp_map->{$eff_grp_map_key}) {
		printf STDERR "Undefined Class:%s\n", $eff_grp_map_key;
		return;
	}

	my $geneid = $eff_elements[3];
	my $genename = $eff_elements[6];
	if (!defined $geneid) {
		$geneid = "none";
	}
	if (!defined $genename) {
		$genename = "none";
	}
	my $idname = "$geneid:$genename";
	my $re = $eff_name;
	if($eff_name =~ /\+/) {
		$re =~ s/\+/\\\+/g;
	}
        if(exists $effname->{"$rows[0]:$rows[1]"}{"$idname"}) {
                if ($effname->{"$rows[0]:$rows[1]"}{"$idname"} !~ /$re/) {
                        $effname->{"$rows[0]:$rows[1]"}{"$idname"} .= ",$eff_name";
                }
        } else {
                $effname->{"$rows[0]:$rows[1]"}{"$idname"} = "$eff_name";
        }

#        my $eff_grp_map_key = (split(/\+/, $eff_name))[0];
        if(exists $effgrp->{"$rows[0]:$rows[1]"}) {
                if ($effgrp->{"$rows[0]:$rows[1]"} !~ /$eff_grp_map->{$eff_grp_map_key}/) {
                        $effgrp->{"$rows[0]:$rows[1]"} .= ",$eff_grp_map->{$eff_grp_map_key}";
                }
        } else {
                $effgrp->{"$rows[0]:$rows[1]"} .= "$eff_grp_map->{$eff_grp_map_key}";
        }
#printf "%s\t%s\t%s\n", $eff_name, $idname, "$eff_grp_map->{$eff_grp_map_key}";
}
sub getMaxGrpAnn {
	my($eff_grp_map, $ontologyStr) = @_;
	my @ontologys = split /\&/, $ontologyStr;
	my $maxGrp = -1;
	my $maxClass = "";
	foreach my $ontology (@ontologys) {
		if($ontology && exists $eff_grp_map->{$ontology}) {
			my $grp = $eff_grp_map->{$ontology};
			if($maxGrp < $grp) {
				$maxClass = $ontology;
				$maxGrp = $grp;
			}
		} elsif(!$maxClass) {
			$maxClass = $ontology;
		}
	}
	return $maxClass;
}
sub setEff {
	my ($effname, $effgrp, $eff_grp_map, $chr, $pos, $eff) = @_;

	return if($eff !~ /(\S+)\((\S+)\)/);

	my $eff_name = $1;
	my $eff_str = $2;

	my @rows = ($chr, $pos);
#	my $idname;
	my @eff_elements = split(/\|/, $eff_str);
	my $geneid = 'none';
	my $genename = 'none';
	if(@eff_elements > 10) {
		# new EFF format
		$geneid = $eff_elements[5] if(defined $eff_elements[5]);
		$genename = $eff_elements[8] if(defined $eff_elements[8]);
	} else {
		# old EFF format
		$geneid = $eff_elements[4] if(defined $eff_elements[4]);
		$genename = $eff_elements[7] if(defined $eff_elements[7]);

	}
	my $idname = "$geneid:$genename";
=head
	if (defined $eff_elements[5] || defined $eff_elements[8]) {
		$idname = "$eff_elements[5]:$eff_elements[8]";
	} else {
		$idname = "none:none";
	}
	if (!defined $eff_elements[5] ) {
		$eff_elements[5] = "none";
	}
	if (!defined $eff_elements[8]) {
		$eff_elements[8] = "none";
	}
	my $geneid = $eff_elements[5];
	my $genename = $eff_elements[8];
	$idname = "$geneid:$genename";
=cut
	my $re = $eff_name;
	if($eff_name =~ /\+/) {
		$re =~ s/\+/\\\+/g;
	}
	if(exists $effname->{"$rows[0]:$rows[1]"}{"$idname"}) {
		if ($effname->{"$rows[0]:$rows[1]"}{"$idname"} !~ /$re/) {
			$effname->{"$rows[0]:$rows[1]"}{"$idname"} .= ",$eff_name";
		}
	} else {
		$effname->{"$rows[0]:$rows[1]"}{"$idname"} = "$eff_name";
	}

	my $eff_grp_map_key = (split(/\+/, $eff_name))[0];
	if(exists $effgrp->{"$rows[0]:$rows[1]"}) {
		if ($effgrp->{"$rows[0]:$rows[1]"} !~ /$eff_grp_map->{$eff_grp_map_key}/) {
			$effgrp->{"$rows[0]:$rows[1]"} .= ",$eff_grp_map->{$eff_grp_map_key}";
		}
	} else {
		$effgrp->{"$rows[0]:$rows[1]"} .= "$eff_grp_map->{$eff_grp_map_key}";
	}
#printf "%s\t%s\t%s\n", $eff_name, $idname, "$eff_grp_map->{$eff_grp_map_key}";	
}

sub getGtIsNotZeroSql {
	return "(GT!='0/0' AND GT!='0|0')";
}
=head
sub isGtZero {
	my($gt) = @_;

	if($gt eq '0/0' || $gt eq '0|0') {
		return 1;
	} else {
		return 0;
	} 
}
=cut
