use strict;
use warnings;
use Util;
use xmlInputForEqtlViewer;
chdir 'G:\back\Downloads\doc\bystrykh';

my %cumulative=();
load_cumulative();

my %synonymous=();
load_synonymous();

my %affy=();
load_affy();

my %genome=();

load_coordinate();



my %qtl=();
load_data('ng1497-S2_wei.txt',8,0,'cis');
load_data('ng1497-S3_wei.txt',9,0,'trans');

my $vic=10*1e6;

output(\%genome,\%qtl,getSetting(),$vic,'mice_test.xml');
print "\a\a\a";

sub getSetting{
	my %r=(
		'urlPrefix'=>'http://www.informatics.jax.org/javawi2/servlet/WIFetch?page=searchTool&selectedQuery=Genes+and+Markers&query=',
		'title'=>'eQTL mapping in D2xB6');
	return \%r;
}


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{'markerName'});
			}
			else{
				$qtl{$s[$title{'symbol'}]}=[($r{'markerName'})];
			}

			
		}
	}
	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_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;
}