use strict;
use warnings;
use Util;

chdir 'G:\back\Downloads\doc\bystrykh';

my %cumulative=();
load_cumulative();

my %synonymous=();
load_synonymous();

my %affy=();
load_affy();

my %genome=();

#$genome{'Mcmd6'}{chr}=1;
#$genome{'Mcmd6'}{pos}=(130159137+130187210)/2;

#$genome{'AI840675'}{'chr'}=4;
#$genome{'AI840675'}{'pos'}=(43403079+43448188)/2;

#'Kiaa0251-hp'

load_coordinate();



my %qtl=();
load_data('ng1497-S2_wei.txt',8,0,'cis');
load_data('ng1497-S3_wei.txt',9,0,'trans');

my @gList=();
sort_genome();
my %gList=();
sort_genome_to_chr();
my $vicinity=10*1e6;
output();
print "\a\a\a";


sub output{
	my $url='http://www.informatics.jax.org/javawi2/servlet/WIFetch?page=searchTool&selectedQuery=Genes+and+Markers&query=';
	$url=~s/(\W)/sprintf("%%%x", ord($1))/eg;
	open(OUT,'>mice.xml') 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","\t\t",'<title>','eQTL mapping in D2xB6','</title>',"\n";
	print OUT "\t\t",'<urlPrefix>'.$url.'</urlPrefix>',"\n";
	print OUT "\t",'</setting>',"\n";
	
	foreach my $trait (sort keys %qtl) {
		foreach my $qtl (@{$qtl{$trait}}) {
			print OUT join("\n",gen_qtl_info($qtl));
		}
	}
			

	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 (keys %boundary)) {
		print $c,"\t",(scalar @{$gList{$c}}),"\t",($boundary{$c}{max}-$boundary{$c}{min}),"\n";
	}
	
	foreach my $c (sort_mixed(keys %gList)) {
		print OUT "\t<chr>chromosome $c</chr>\n";
	}
	
	print OUT "</QTLViewer>\n";
}




sub gen_qtl_info{
	my $qtl=shift;
	my @r=();
	push(@r,"\t<QTL>");
	push(@r,"\t\t<ofTrait>$qtl->{traitName}</ofTrait>");
	die unless(exists $genome{$qtl->{traitName}} and exists $genome{$qtl->{markerName}});
	foreach my $g (@{get_genes($genome{$qtl->{markerName}})}) {
		push(@r,"\t\t<onMap>$g</onMap>");
	}	
	push(@r,"\t</QTL>\n");
	return @r;
}

sub get_genes{

	my $marker=shift;
	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{
	foreach my $g (@gList) {
		my %r=('name'=>$g,'pos'=>$genome{$g}{'pos'});
		push(@{$gList{$genome{$g}{'chr'}}},\%r);
	}
}

sub sort_genome{

=the following two genes
are found with genes of different names but exactly same genomic boundarys
it is found their "entrez gene id" is null, there are replaced.

Olfr1354
Olfr767
=cut

	my %location=();
	
	foreach my $g (keys %genome) {
		
		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";
}

# sort_mixed sorts chromosomes, numbers are at the begining, sorted numerically
# characters like X Y are at the end, sorted alphabetically

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);
}


sub load_data{
	my ($inFile,$ncol,$keyPos,$qtl_type)=@_;
	my %title=();
	open (IN,$inFile) or die;
	while (<IN>) {
		chomp;
		my @s=split(",") or die;
		if ($.==1) {
			die "check $inFile" unless(@s == $ncol);
			for (my $i=0;$i<@s ;$i++) {
				$title{$s[$i]}=$i;
			}
		}
		else{
# handling marker
			if (!exists $genome{$s[$title{'QTL marker'}]}) {
				$genome{$s[$title{'QTL marker'}]}{'chr'}=$s[$title{'QTL chromosome'}] or die;
				$genome{$s[$title{'QTL marker'}]}{'pos'}=$s[$title{'Marker position (Megabase)'}]-$cumulative{$s[$title{'QTL chromosome'}]}*1e6;
				die "check $inFile line $." if($genome{$s[$title{'QTL marker'}]}{'pos'}<0);
			}

# change trait name if possible
			if(!$s[$title{'symbol'}]){
				if(exists $affy{$s[$title{'Probeset ID'}]}) {
					$s[$title{'symbol'}]=$affy{$s[$title{'Probeset ID'}]}{'symbol'} or die;
				}
				else{
					next;
				}
			}
			
			if (exists $synonymous{$s[$title{'symbol'}]}) {
				$s[$title{'symbol'}]=$synonymous{$s[$title{'symbol'}]} or die;
			}

# qtl informaiton
			my %r=();
			$r{traitName}=$s[$title{'symbol'}] or die;
			$r{'markerName'}=$s[$title{'QTL marker'}] or die;

			if (!exists $genome{$r{traitName}}) {
				if (exists $affy{$s[$title{'Probeset ID'}]}) {
					$genome{$r{traitName}}{'chr'}=$affy{$s[$title{'Probeset ID'}]}{'chr'} or die "check $inFile line $.";
					$genome{$r{traitName}}{'pos'}=$affy{$s[$title{'Probeset ID'}]}{'pos'} or die;
				}
				else{
					print "no information for $r{traitName} from $inFile line $.-1\n";
					next;
				}
			}

			if (exists $qtl{$s[$title{'symbol'}]}) {
				push(@{$qtl{$s[$title{'symbol'}]}},\%r);
			}
			else{
				$qtl{$s[$title{'symbol'}]}=[(\%r)];
			}

			
		}
	}
	close (IN);
}


sub load_synonymous{
	open (IN,'MRK_Synonym.sql.rpt') or die;
	while (<IN>) {
		chomp;
		next unless (/^\s*MGI:/);
		my $std=substr($_,53,50);
		$std=~/(.*\S)\s{5}/ or die;
		$std=$1 or die;
		my $old=substr($_,104);
		$old=~/(.*\S)\s{5}/ or die;
		$old=$1 or die;
		$synonymous{$old}=$std;
	}
	close (IN);
}
sub load_coordinate{
	my $titleLine=1;
	my %title=();

	my %throwChr=();	 # number of features thrown away for each chr

	open (IN, 'MGI_Coordinate.rpt') or die;
	while (<IN>) {
		chomp;
		my @s=split(/\t/);
		if ($titleLine == $. ) {
			for (my $i=0;$i<@s ;$i++) {
				$title{$s[$i]}=$i;
			}
			
		}
		else{
#			next unless($s[1] eq 'Gene');

			if($s[4]!~/[XY\d]/){
	#			$positionUnknown{$s[2]}='';
				$throwChr{$s[4]}++;
				next;
			}
			die "$s[2] already exists"  if(exists $genome{$s[2]});

			$genome{$s[2]}{'chr'}=$s[4];
			$genome{$s[2]}{'pos'}=($s[5]+$s[6])/2;
		}
	}
	close (IN);

	print "total number of genes from coordinate : ", scalar keys %genome,"\n";
	
	foreach my $c (keys %throwChr) {
		print "$throwChr{$c} genes in chr $c are not included\n";
	}

	print "\nfinish reading genomic information\n","=" x 60,"\n\n";
}


sub load_affy{

	
	my %title=();

	open (IN, 'MG_U74Av2_annot.txt') or die ;
	while (<IN>) {
	
		chomp;
		my @s=split(/\t/);
		
		if ($.==1) {
			for (my $i=0;$i<@s ;$i++) {
				$title{$s[$i]}=$i;
			}
			next;
		}
		else{
			my $info=$s[$title{'Alignments'}] or die;
			next if($info eq '---');
			$info=~/chr([0-9XY]+)\:(\d+)\-(\d+)\s*\(/ or next;
			$affy{$s[$title{'Probe Set ID'}]}{'symbol'}=$s[$title{'Gene Symbol'}] or die;
			if($affy{$s[$title{'Probe Set ID'}]}{'symbol'} eq '---') {
				$affy{$s[$title{'Probe Set ID'}]}{'symbol'}=$s[$title{'Probe Set ID'}];
			}
			$affy{$s[$title{'Probe Set ID'}]}{'chr'}=$1 or die;
			$affy{$s[$title{'Probe Set ID'}]}{'pos'}=($2+$3)/2 or die;
		}

	}
	close (IN);
}


sub load_EntrezGene{
	open (IN,'MGI_EntrezGene.rpt') or die;
	my %title=();
	while (<IN>) {
		chomp;
		$_.="\twei";
		my @s=split("\t");
		pop(@s);
		if ($.==1) {
			for (my $i=0;$i<@s ;$i++) {
				$title{$s[$i]}=$i
			}
		}
		else{
			die unless(@s==scalar keys %title);
			$s[$title{'cM Position'}]=~s/\s//g;
			if($s[$title{Synonyms}] and $s[$title{'cM Position'}] 
				and $s[$title{'cM Position'}]>0) {

				my @n=split(/\|/,$s[$title{Synonyms}]);
				foreach my $n (@n) {
					$genome{$n}{'chr'}=$s[$title{Chromosome}] or die "check $.";
					$genome{$n}{cm}=$s[$title{'cM Position'}] or die "check $.";
					
				}
			}
			if ($s[$title{'cM Position'}] and $s[$title{'cM Position'}]>0 and 
				(! exists $genome{$s[$title{'Marker Symbol'}]})) {
				$genome{$s[$title{'Marker Symbol'}]}{'chr'}=$s[$title{Chromosome}] or die "check $.";
				$genome{$s[$title{'Marker Symbol'}]}{cm}=$s[$title{'cM Position'}] or die "check $.";
			}
		}
	}
	close(IN);
}

sub load_cumulative{
	$cumulative{1}=0;
	my @len=();
	open (IN,'chr.TXT') or die;

	while (<IN>) {
		chomp;

		next if ($.==1);
		my @s=split("\t");
		if (scalar @len>0) {
			$cumulative{$s[0]}=$len[-1]+$cumulative{getPrevChr($s[0])};
		}
		push(@len,$s[2]);
	}
	close (IN);

	foreach my $c (sort keys %cumulative) {
		print $c,"\t",$cumulative{$c},"\n";
	}
}

sub getPrevChr{
	return 19	if ($_[0] eq 'X') ;
	return 'X' if ($_[0] eq 'Y');
	return $_[0]-1;
}