package xmlInputForEqtlViewer;
use strict;
use warnings;
require Exporter;
use vars qw(@ISA @EXPORT) ;
@ISA = qw(Exporter);
@EXPORT = qw(output);

=input data structure

*****************************************************************************

genome is stored in a hash of hash, eg, '%genome'
1st level key is a gene/genomic element name, 1st level value is a reference to the 2nd level hash
2nd level keys are	'chr'=>the gene's chromosome location
					'pos'=>the gene's position within the chromosome

*****************************************************************************

mapping results are stored in a hash of arrays, eg, '%mappingResults'
1st level key is the trait gene name, 1st level value is a reference to the 2nd level array
2nd level array has the same length as the number QTL of the trait

if QTL region is estimated using lod support interval, each element of the 2nd level array
is a reference to an other array of strings. 
Each string is a name of a gene which is located with the qtl region

if QTL region is estimated using an arbituray physical region around a gene/marker, each element
of the 2nd level array is a string. The string is the gene/marker name, which is the center of the QTL

If you use a genetic marker as the center of a QTL, please incorporate that marker into %genome

*****************************************************************************

=cut



sub output{
	my ($genome,$mappingResults,$setting,$vicinity,$outFileName)=@_;
	
	my $gList=sort_genome_to_chr($genome);

	open(OUT,'>'.$outFileName) or die;
	print OUT '<?xml version="1.0" ?>'."\n".'<!DOCTYPE QTLViewer SYSTEM "http://statgen.ncsu.edu/~wzou/svg/QTLViewer.dtd">'."\n".'<QTLViewer>'."\n";
	print OUT "\t",'<setting>',"\n";
	foreach my $s (sort keys %{$setting}) {
		my $v=$setting->{$s};
		if ($s eq 'urlPrefix') {
			$v=~s/(\W)/sprintf("%%%x", ord($1))/eg;
		}
		print OUT "\t\t",'<',$s,'>',$v,'</',$s,'>',"\n";
	}
	print OUT "\t",'</setting>',"\n";
	
	my %unknown=();
	
	foreach my $trait (sort keys %{$mappingResults}) {
		if (! exists $genome->{$trait}) {
			print "cannot find the genomic position of $trait...all its eQTL are skipped\n";
			next;
		}
		foreach my $qtlRegion (@{$mappingResults->{$trait}}) {
			print OUT join("\n",gen_qtl_info($trait,$qtlRegion,$genome,$gList,$vicinity,\%unknown));
		}
	}

	foreach my $u (keys %unknown) {
		print "cannot find the genomic position of gene $u, which is included in $unknown{$u} eQTL...skipped\n";
	}
			

	my %boundary=();

	foreach my $g (keys %{$genome}) {
		print OUT "\t<marker name=",'"',$g,'" chr="chromosome ',$genome->{$g}->{'chr'};
		print OUT '" position="',$genome->{$g}->{'pos'},'"/>',"\n";
		if (exists $boundary{$genome->{$g}->{'chr'}}) {
			$boundary{$genome->{$g}->{'chr'}}{min}=$genome->{$g}->{'pos'}
			if ($boundary{$genome->{$g}->{'chr'}}{min}>$genome->{$g}->{'pos'});
			$boundary{$genome->{$g}->{'chr'}}{max}=$genome->{$g}->{'pos'}
			if ($boundary{$genome->{$g}->{'chr'}}{max}<$genome->{$g}->{'pos'});
		}
		else{
			$boundary{$genome->{$g}->{'chr'}}{min}=$genome->{$g}->{'pos'};
			$boundary{$genome->{$g}->{'chr'}}{max}=$genome->{$g}->{'pos'};
		}
	}

	
	foreach my $c (sort_mixed(keys %{$gList})) {
		print OUT "\t<chr>chromosome $c</chr>\n";
	}
	
	print OUT "</QTLViewer>\n";
}




sub gen_qtl_info{
	my ($traitGeneName,$qtl,$genome,$gList,$vicinity,$unknown)=@_;
	my @r=();
	push(@r,"\t<QTL>");
	push(@r,"\t\t<ofTrait>$traitGeneName</ofTrait>");
	die unless(exists $genome->{$traitGeneName});

	if ($vicinity) {
		die "if specifying the vicinity, each qtl region is represented by the cental marker."
			if(ref($qtl)); 
		die "cannot find $qtl 's genomic information" unless exists $genome->{$qtl};
		foreach my $g (@{get_genes($genome->{$qtl},$gList,$vicinity)}) {
			push(@r,"\t\t<onMap>$g</onMap>");
		}
	}
	else{
		die "if qtl regions are estimated using lod support interval, each qtl should be an array reference to the gene list in the qtl"
			unless(ref($qtl) eq "ARRAY");	

		foreach my $g (@{$qtl}) {
			unless (exists $genome->{$g}){
				$unknown->{$g}++;
				next;
			}
			push(@r,"\t\t<onMap>$g</onMap>");
		}
	}

	push(@r,"\t</QTL>\n");
	return @r;
}

sub get_genes{

	my ($marker,$gList,$vicinity)=@_;
	my $c=$marker->{'chr'} or die;
	my $pos=$marker->{'pos'} or die;

	die unless(exists $gList->{$c});
	my @r=();
	foreach my $g (@{$gList->{$c}}) {
		if( abs($pos - $g->{'pos'}) <$vicinity ){
			push(@r,$g->{name});
		}
	}
	return \@r;
}


sub sort_genome_to_chr{
	my ($genome)=@_;
	my %gList=();

	foreach my $g (@{sort_genome($genome)}) {
		die "function sort_genome_to_chr requires $g 's physical position information" unless(exists $genome->{$g}->{'pos'});
		die "function sort_genome_to_chr requires $g 's chromosome information" unless (exists $genome->{$g}->{'chr'});
		my %r=('name'=>$g,'pos'=>$genome->{$g}->{'pos'});
		push(@{$gList{$genome->{$g}->{'chr'}}},\%r);
	}
	return \%gList;
}

sub sort_genome{
	my $genome=shift;
	my %location=();
	my @gList=();

	foreach my $g (keys %{$genome}) {
		die "function sort_genome requires $g 's physical position information" unless(exists $genome->{$g}->{'pos'});
		die "function sort_genome requires $g 's chromosome information" unless (exists $genome->{$g}->{'chr'});
		my $position=int($genome->{$g}->{'pos'});
		while(exists $location{$genome->{$g}->{'chr'}}{$position})  {
			$position++;
		}
		$genome->{$g}->{'pos'}=$position;
		$location{$genome->{$g}->{'chr'}}{$position}=$g;
	}

	foreach my $c (sort_mixed(keys %location)) {
		foreach my $pos (sort {$a <=> $b} (keys %{$location{$c}}))  {
			push(@gList,$location{$c}{$pos});
		}
	}

	print "\nfinish sorting genes in genome\n","=" x 60,"\n\n";
	
	return [@gList];
}

sub sort_mixed{
	
	my (@num,@char)=();
	foreach my $e (@_) {
		if ($e=~/\d/) {
			push(@num,$e);
		}
		else{
			push(@char,$e);
		}
	}

	my @r1= sort {$a <=> $b} @num;
	my @r2= sort @char;
	return (@r1,@r2);
}
1