#!/usr/bin/perl
############################################################
# file: index.cgi, in tools/MeltMap
# desc: cgi-script for PubGene MeltMap viewer
# auth: Tor-Kristian Jenssen, PubGene AS
# comm:
#===========================================================
#
# Copyright 2003 PubGene AS, all rights reserved
#
#===========================================================
# who | dd/mm/yyyy | what
#-----------------------------------------------------------
# tkj | 03/03/2003 | created
############################################################

use lib '../../lib';

use strict;
use Carp;
use Benchmark;

require PubGene;
require MeltMap::Configuration;
require MeltMap::Database;
require MeltMap::Plot;
require MeltMap::MySQL;

# file globals
my $DEBUG=0;

# start program
my $pg = PubGene::Web::Page->new("Melting Map Exon Viewer",
																 'melting_map');
unless (defined $pg){
	exit(PubGene::Web::Page::FATAL_ERROR("cannot create Page Object"));
}
#$DEBUG && $pg->list_params();
my $q = $pg->{'cgi'};
if ($q->param('gethelp') eq 'yes'){
	exit(PubGene::Web::Page::Help::layout($pg));
}

$pg->add("With this tool you may view melting maps for individual exons.<br>");
$pg->add("The melting maps are generated with temperature step of .1 degrees Celcius.<br>");

$pg->add("NOTE: exon files and database data ready for all exons from both the Ensembl and RefSeq exon prediction sets.<br>");

$pg->add("NOTE: Exons only for June 28 2002 freeze.<br>");

#if ($q->param('ListExons')){
#	$pg->add("List Exons selected");
#}
$pg->add(chrmenu($pg));

if ($q->param('GetMeltMap')){
	$pg->add(results($pg));
}
$pg->layout($DEBUG);
exit;

######
# results
######
sub results {
	my $pg = shift;
	my $q = $pg->{'cgi'};
    
	my $error = 0; # error code
	my $msg;
	my $t0 = Benchmark->new();

	my $results = "";  # result buffer (string)

	my $organism = 'hs';
	my $chr = $q->param('chr');
	my $exon_set = $q->param('exon_set');
	my $exon = $q->param('exon');
	my $transcript = $q->param('transcript');
	my $exon_num = $q->param('exon_num');

	if ($exon){
		($transcript, $exon_num) = split(/_/, $exon);
	}else{
		if ($transcript && $exon_num){
			$exon = $transcript . '_' . $exon_num;
		}else{
			$results .= $pg->error("exon name is empty");
			return $results;
		}
	}

	# lookup transcript in db
	my $rh_transcript_info;
	my $tmp_msg;
	($error,
	 $tmp_msg,
	 $rh_transcript_info) = MeltMap::MySQL::get_transcript($transcript,
																												 $exon_set);
	if ($error){
		$msg .= "Failed to retrieve transcript info for $transcript:<br>";
		$msg .= $tmp_msg;
		return $results .= $msg;
	}
	
	if (0){
		$results .= $tmp_msg;
		foreach (sort keys %$rh_transcript_info){
			$results .= "$_ = " . $rh_transcript_info->{$_} . "<br>";
		}
		return $results;
	}

	# lookup gene in db
	$tmp_msg = '';
	my $gene;
	($error,
	 $tmp_msg,
	 $gene) = MeltMap::MySQL::get_gene_by_transcript($transcript,
																									 $exon_set);
	if ($error){
		$msg .= "Failed to retrieve gene info for $transcript:<br>";
		$msg .= $tmp_msg;
		return $results .= $msg;
	}
	$rh_transcript_info->{'gene'} = $gene;
	if (0){
		$results .= $tmp_msg;
		$results .= "gene: >$gene<";
		return $results;
	}

	my $other_set;
	if ($exon_set eq 'ensembl'){
		$other_set = 'refseq';
	}elsif ($exon_set eq 'refseq'){
		$other_set = 'ensembl';
	}

	my $rh_other_tr_info;
	if ($gene){
		my $other_transcript;
		($error, $tmp_msg, $other_transcript) =
		 MeltMap::MySQL::get_transcript_by_gene($gene, $other_set);
		if ($error){
			$msg .= "Failed to retrieve transcript for gene $gene in $other_set:<br>";
			$msg .= $tmp_msg;
			return $results .= $msg;
		}
		if (0){
			$results .= $tmp_msg;
			$results .= "name of transcript for $gene in $other_set: >$other_transcript<";
			#return $results;
		}

		($error,
		 $tmp_msg,
		 $rh_other_tr_info) = MeltMap::MySQL::get_transcript($other_transcript,
																												 $other_set);
		if ($error){
			$msg .= "Failed to retrieve transcript info for $other_transcript in $other_set:<br>";
			$msg .= $tmp_msg;
			return $results .= $msg;
		}
		$rh_other_tr_info->{'gene'} = $gene;
	
		if (0){
			$results .= $tmp_msg;
			foreach (sort keys %$rh_transcript_info){
				$results .= "$_ = " . $rh_transcript_info->{$_} . "<br>";
			}
			return $results;
		}
	}
	my $fname = MeltMap::Configuration::get_exon_fname($exon,
																										 $chr,
																										 $exon_set,
																										 'bin',
																										 'maps',
																										 'golden');
	unless ($fname){
		$results .= $pg->error("cannot get exon melt map file name");
		return $results;
	}
	unless (-e $fname){
		$results .= $pg->warning("melt map file for exon is not ready, yet",
														 $fname,
														 "please select another exon");
		return $results;
	}
	unless (-r $fname){
		$results .= $pg->warning("melt map file for exon is not readable, yet",
														 $fname,
														 "please select another exon");
		return $results;
	}

	# ASSERT: file exists and is readable, start and end are OK
	
	$results .= "Melting map for exon $exon";
	$results .= "&nbsp;from chromosome $chr,";
	$results .= "&nbsp;gene prediction set is $exon_set<br>";
	#$results .= $fname ."<br>";

	my $start = 1;
	my $end = 10000;
	my $avg_check;
	my $avg_method;
	my $avg_wsize;
	my $num_bps = $end - $start + 1;

	# check averaging params
	my $avg_lim = $MeltMap::Database::AVG_MAX_NUM_POINTS;
	if ($num_bps > $avg_lim){
		unless ($avg_check eq 'on'){
			$results .= $pg->warning("sequence length is larger than $avg_lim",
															 "but averaging was not selected",
															 "averaging will be performed with auto size");
			$avg_check = 'on';
			$avg_wsize = 'auto';
		}
	}

	my @plot_points;
	my $window_size;

	my $max_lim = $MeltMap::Database::MAX_LENGTH;
	if ($num_bps > $max_lim){
		$results .= $pg->warning("sequence length is too large; max is $max_lim");
		return $results;
	}

	# check read length
	my $read_lim = $MeltMap::Database::MAX_READ_LENGTH;
	if ($num_bps > $read_lim){
		#$results .= $pg->warning("sequence length is larger than $read_lim",
		#												 "sequence reading will be performed per window");
		my $tmp_msg;
		if ($avg_wsize eq 'auto'){
			($window_size, $tmp_msg) = 
			 MeltMap::Database::auto_wsize($num_bps);
			$results .= "auto window size determined to $window_size<br>";
		}else{
			$window_size = $avg_wsize;
			my $num_windows = $num_bps / $window_size;
			if ($num_windows > $MeltMap::Database::AVG_MAX_NUM_POINTS){
				$results .= $pg->warning("window size too small, for sequence length",
																 "would require $num_windows points",
																 "please increase the window size");
				return $results;
			}
		}

		if ($window_size > $read_lim){
			$results .= $pg->warning("the window size is too large: max is $read_lim");
			return $results;
		}

		my $this_start_bp = $start;
		while ($this_start_bp < $end){
			my $this_end_bp = $this_start_bp + $window_size - 1;
			$this_end_bp = $end if ($this_end_bp > $end);
			#$results .= "$this_start_bp, $this_end_bp<br>";
			my @temps;
			($error,$msg) = 
			 MeltMap::Database::get_melt_map_chr_segment($fname,
																									 'bin',
																									 $this_start_bp,
																									 $this_end_bp,
																									 \@temps);
			if ($error){
				$results .= $pg->error("cannot get melt map for segment",
															 "message was: $msg");
				return $results;
			}
			unless (@temps){
				$results .= $pg->error("could not read temperatures from $fname",
															 "wanted to read bases $this_start_bp to $this_end_bp");
				return $results;
			}
			my $tmp = MeltMap::Database::calc_avg(\@temps, $avg_method);
			unless ($tmp){
				$results .= $pg->error("failed to calc average ($avg_method)");
				return $results;
			}
			push @plot_points, $tmp;
			$this_start_bp = $this_end_bp + 1;
		}
		my $n_points = scalar(@plot_points);
	}else{

		my @temps;
		($error,$msg) = 
		 MeltMap::Database::get_melt_map_chr_segment($fname,
																								 'bin',
																								 $start,
																								 $end,
																								 \@temps);
		if ($error){
			$results .= $pg->error("cannot get melt map for segment",
														 "message was: $msg");
			return $results;
		}
		unless (scalar(@temps)){
			$results .= $pg->error("could not read temperatures from $fname",
														 "wanted to read bases $start to $end");
			return $results;
		}

		my $num_tmps = scalar(@temps);
		#$results .= "read $num_tmps temperatures from $fname<br>";
		my $first = $temps[0];
		#$results .= "first temperature is $first.<br>";
		if ($temps[$#temps] == 0){
			pop @temps;
		}
	
		if ($avg_check eq 'on'){
			($error, $msg, $window_size) =
			 MeltMap::Database::average(\@temps,
																	$avg_method,
																	$avg_wsize);
			if ($error){
				$results .= $pg->error("failed to average over windows",
															 $msg);
				return $results;
			}
			$results .= "window size is $window_size<br>";
			$results .= $msg;
			$results .= "<br>";
			unless (@temps){
				$results .= $pg->error("no points after averaging<br>");
				return $results;
			}
		}
		push @plot_points, @temps;
	}
	# ASSERT: points to plot are in @plot_points, $window_size is set

	my $n_points = scalar(@plot_points);
	#$results .= "got $n_points points to plot<br>";

	my $tmp_fname = $pg->create_image_fname('svg');
	my $image_fname = $PubGene::TMP_DIR . $tmp_fname;
	my $image_url = $PubGene::Web::TMP_DIR . $tmp_fname;
	($error, $msg) = MeltMap::Plot::svg(\@plot_points,
																			$start,
																			$n_points,
																			$window_size,
																			600,
																			200,
																			$image_fname,
																			);
	if ($error){
		$results .= $pg->error("failed to create svg file",
													 $msg);
		return $results;
	}

	$results .= format_exon_info($rh_transcript_info, $exon_num);

	my $svg_html = $pg->svgimage($image_url,600,200);
	$results .= $svg_html;

	if ($rh_other_tr_info && $other_set){
		$results .= "<hr>";
		$results .= "<br>";
		$results .= "Transcipt information for this gene in $other_set:<br>";
		$results .= format_exon_info($rh_other_tr_info, $exon_num);
	}

	my $tstr = timestr(timediff(Benchmark->new(), $t0));
	$results .= "<br>time used on server: $tstr<br>";

	return $results;
}

######
# chrmenu
######
sub chrmenu {
	my $pg = shift;

	my $q = $pg->{'cgi'};

	my $error;
	my $msg;
	my $nmenu;

	my @exons;
	my %transcripts;
	my %genes;
	my $chr = $q->param('chr') || 'Y';
	my $exon_set = $q->param('exon_set') || 'ensembl';
	my $exon = $q->param('exon');
	my $transcript = $q->param('transcript');
	my $exon_num = $q->param('exon_num');

	unless ($exon){
		if ($transcript && $exon_num){
			$exon = $transcript . '_' . $exon_num;
		}
	}

	my $list_exons = ($q->param('ListExons') ? 1 : 0);
	my $list_transcripts = ($q->param('ListTranscripts') ? 1 : 0);
	my $list_genes = ($q->param('ListGenes') ? 1 : 0);

	if ($list_exons){
		($error, $msg, @exons) = 
		 MeltMap::Configuration::list_exons($chr,
																				$exon_set);
		if ($error){
			$nmenu .= "Failed to get list of exons for $chr, $exon_set<br>";
			$nmenu .= $msg;
			return $nmenu;
		}
		$nmenu .= $msg;
	}elsif ($list_transcripts){
		($error, $msg, %transcripts) = 
		 MeltMap::Configuration::list_transcripts($chr,
																							$exon_set);
		if ($error){
			$nmenu .= "Failed to get list of transcripts for $chr, $exon_set<br>";
			$nmenu .= $msg;
			return $nmenu;
		}
		$nmenu .= $msg;
	}elsif ($list_genes){
		($error, $msg, %genes) =
		 MeltMap::MySQL::list_genes($chr,
																$exon_set);
		if ($error){
			$nmenu .= "Failed to get list of genes for $chr, $exon_set<br>";
			$nmenu .= $msg;
			return $nmenu;
		}
		#$nmenu .= $msg;
	}
    
	my $organism = 'hs';
	my @chrs = MeltMap::Configuration::get_chromosomes($organism);
	my %chrs;
	@chrs{@chrs} = @chrs;
	my @exon_sets = MeltMap::Configuration::get_exon_sets($organism);
	my %exon_sets;
	@exon_sets{@exon_sets} = @exon_sets;
	$nmenu .= $q->startform('POST');
	$nmenu .= "Select chromosome:&nbsp;";
	$nmenu .= $q->popup_menu('-name' => 'chr',
													 '-values' => \@chrs,
													 '-labels' => \%chrs,
													 '-default' => '1');
	$nmenu .= "Select exon prediction set:&nbsp;";
	$nmenu .= $q->popup_menu('-name' => 'exon_set',
													 '-values' => \@exon_sets,
													 '-labels' => \%exon_sets,
													 '-default' => 'ensembl');
	
	$nmenu .= "<br>";
	$nmenu .= $q->submit('-name' => 'ListGenes',
											 '-value' => 'List Genes (named)');
	$nmenu .= $q->submit('-name' => 'ListTranscripts',
											 '-value' => 'List Transcripts');
	$nmenu .= $q->submit('-name' => 'ListExons',
											 '-value' => 'List Exons');
	$nmenu .= "<br>";
	if ($list_exons){
		$nmenu .= "Select exon from popup list and get melt map (";
		$nmenu .= scalar(@exons) . " exons on chromosome $chr, acc. to $exon_set)";
		$nmenu .= "<br>";
		$nmenu .= $q->popup_menu('-name' => 'exon',
														 '-values' => \@exons);
	}elsif ($list_transcripts){
		my @transcripts = sort keys %transcripts;
		my %labels;
		foreach (@transcripts){
			$labels{$_} .= $_ . ' - exons 1 to ' . $transcripts{$_};
		}
		$nmenu .= "Select transcript from popup list, enter exon number, and get melt map (";
		$nmenu .= scalar(@transcripts) . " transcripts on chromosome $chr, acc. to $exon_set)";
		$nmenu .= "<br>";
		$nmenu .= $q->popup_menu('-name' => 'transcript',
														 '-values' => \@transcripts,
														 '-labels' => \%labels);
		$nmenu .= "Select exon number from gene (transcript): ";
		$nmenu .= $q->textfield('-name' => 'exon_num',
														'-default' => '1',
														#'-align' => 'right',
														'-size' => 3);
	}elsif ($list_genes){
		my @genes = sort keys %genes;
		my @transcripts;
		my %labels;
		foreach (@genes){
			my $tr = $genes{$_}->[0];
			push @transcripts, $tr;
			$labels{$tr} = $_ . ' - exons 1 to ' . $genes{$_}->[1];
		}
		
		$nmenu .= "Select gene from popup list, enter exon number, and get melt map (";
		$nmenu .= scalar(@genes) . " genes on chromosome $chr, acc. to $exon_set)";
		$nmenu .= "<br>";
		$nmenu .= $q->popup_menu('-name' => 'transcript',
														 '-values' => \@transcripts,
														 '-labels' => \%labels);
		$nmenu .= "Select exon number from gene (transcript): ";
		$nmenu .= $q->textfield('-name' => 'exon_num',
														'-default' => '1',
														#'-align' => 'right',
														'-size' => 3);
	}else{
		$nmenu .= "Input exon name directly: ";
		$nmenu .= $q->textfield('-name' => 'exon',
														'-default' => $exon,
														#'-default' => 'ENST00000013913_1',
														#'-align' => 'right',
														'-size' => 25);
	}
	$nmenu .= "<br>";
	$nmenu .= $q->submit('-name' => 'GetMeltMap',
											 '-value' => 'Get MeltMap');

	$nmenu .= "<br>";
	$nmenu .= $q->endform();
	
	return $nmenu;
}

sub format_exon_info {
	my $rh_tr_info = shift;
	my $exon_num = shift;

	my $info_str;
	#foreach (sort keys %$rh_tr_info){
		#$info_str .= "$_ = " . $rh_tr_info->{$_} . "<br>";
	#}
	$info_str .= "Melting map for exon $exon_num of transcript ";
	$info_str .= $rh_tr_info->{'name'} . ".<br>";
	my $gene;
	if ($gene = $rh_tr_info->{'gene'}){
		$info_str .= "This transcript corresponds to the $gene gene.<br>";
	}else{
		$info_str .= "This transcript does not correspond to a named gene.<br>";
	}
	if ($rh_tr_info->{'strand'} eq '+'){
		$info_str .= "This transcript is on the sense strand.<br>";
	}else{
		$info_str .= "This transcript is on the anti-sense strand.<br>";
	}
	$info_str .= "This transcript starts at: " . $rh_tr_info->{'txStart'};
	$info_str .= ", and ends at " . $rh_tr_info->{'txEnd'};
	my $num = $rh_tr_info->{'txEnd'} - $rh_tr_info->{'txStart'};
	$info_str .= " ($num bases).<br>";
	$info_str .= "The coding region starts at: " . $rh_tr_info->{'cdsStart'};
	$info_str .= ", and ends at " . $rh_tr_info->{'cdsEnd'};
	$num = $rh_tr_info->{'cdsEnd'} - $rh_tr_info->{'cdsStart'};
	$info_str .= " ($num bases).<br>";
	$info_str .= "This transcript has " . $rh_tr_info->{'exonCount'} . " exons<br>";
	my @starts = split(/,/,$rh_tr_info->{'exonStarts'});
	my @ends = split(/,/,$rh_tr_info->{'exonEnds'});
	$info_str .= "This exon starts at " . $starts[$exon_num-1];
	$info_str .= ", and ends at " . $ends[$exon_num-1];
	$num = $ends[$exon_num-1] - $starts[$exon_num-1];
	$info_str .= " ($num bases).<br>";
	$info_str .= "The melting map for this exon is generated from the exon sequence with 1000 bases upstream and 500 bases downstream.<br>";

	return $info_str;
}