#!/usr/bin/perl

#================================================================
# 0.	contact wei if there is any warning/error output
#		normal output will like this
#
#		total number of z file analyzing: 2
#		now mapping the expression control for YPR149W
#		now mapping the expression control for YKL129C
#				...
#
# 1.	find Roman.pm(from cpan, or wei)
# 		map file used by qtlcartographer 
#		saccharomyces_cerevisiae.gff 
#		hash.repeat (built hash table between probe index and ORF system names, from wei)
#		genotypes_forRelease_12_5_03.dat (from brem,Rachel Brem rbrem@fhcrc.org)
# 		
# 2.	check the user adjust needed section, to see whether you like the settings
# 3.	output file list
#		trait.out			majority of the file is for debug purpose. the end of the file is worth looking	
#		exact_match.list	trait with cis-acting QTL
#		numbers.list		lists for each trait
#		qtl.len				lists for each QTL, including QTL physical boundaries
#		forEachORF.list		lists of eORF for each pORF
# 4.	modification is needed to run it in DOS\windows
#=================================================================

#===================================================================
#           user adjust needed
#===================================================================

my $threshold=16; #what is your threshold to claim a QTL?

my $stepsize =2;   #what is your stepsize in IM scanning?

my $CI=1.5/0.217; #1.5 LOD support interval

my $currentMethod='MIM'; 
# it will be used to differentiate output files by IM or CIM or MIM
# IM CIM 's output is same, differs with those from MIM. 
# This code can handle the difference as long as $currentMethod is specified correctly

my $currentFolder='/Volumes/pooleraid/users/wzou/jessica'; 

#path of working directory, do not attach / at the end
#the following files are supposed to be in this working directory
#as well as the perl code itself

my $map_file_name='gene.map'; # file name of .map file required by qtlcartographer
$map_file_name='q.map' if($currentMethod eq 'MIM');
my $probe_id_to_gene_name='hash.repeat'; #file name records the correspondence between probe id and gene names
my $marker_physical_position='genotypes_forRelease_1_20_05.dat';

# marker physical position is specified in $marker_physical_position
# marker genetic position is specified in $map_file_name
# thus, we need assign each marker an id to connect physical and genetic position info
# when they are loaded from different files
# 
# marker id in MIM: M+lineNumber
#	for example, the 1st marker in genotypes_forRelease_1_20_05.dat is 11265_i_at_x01, coded as M27
#	the 2nd marker is coded as M28
#
# marker id in CIM and IM: m_lineNumber-27
#	for example, the 1st marker in genotypes_forRelease_1_20_05.dat is 11265_i_at_x01, coded as m_0
# 	the 2nd marker is coded as m_1	
# check sub snpfile() for more detail

my $zFolder='/Volumes/pooleraid/users/wzou/jessica'; 
#path for IM, CIM output files, do not attach / at the end
#output file should contain number to index traits
#output file is either with file extension 'z' or 'mim'

my $outputFolder='>'.'/Volumes/pooleraid/users/wzou/jessica'; 
# keep the begining > when you updating the path, do not attach / at the end




#====================================================================


use Roman;
use strict;
use warnings;



#===default parameters===============================================
my $begin=0.0100; # this number is default to IM and CIM
$begin=0.000100 if($currentMethod eq 'MIM');

my $minimalPercentageAsHit=0.5; 
#arbiturily chosen, to determine whether to claim a gene is "in" a physical interval or not
#0.5 means the code will declare the qtl hits the gene if half of the gene in the qtl.
#when the gene is only partly inside the physical boundaries
my $vicinity=10000; #the distance uplimit to claim cis-acting QTL, used by brem 02

my %output_option=();

$output_option{pORFlist}=($currentFolder); 	
# give a output path name if want gene lists in each QTL
# do not attach > at the begining, do not attach / at the end

$output_option{eORFlist}=(); 	
# give a output path name if want gene lists mapped onto each pORF
# do not attach > at the begining, do not attach / at the end

#=================================================================

my(%nuc_orf,%com_ponent_to_plex,@geneticMap,%markerHash,@factor,@target)=();
my(@markerHash,@complexPair,@transPair,%targethash,%factorhash)=();
my(%num_to_name,%name_to_num,%phy_pos,@featureList,@featureBegin,@featureEnd,%featurePos)=();
my($num_diff_chr,$num_same_chr,$num_approxi,$num_exact,$num_trans_QTL,$num_tf_confirmed,
	$num_interaction,$num_no_qtl,$num_no_orf,$num_unknown_location)=(0,0,0,0,0,0,0,0,0,0,0,0);

my(%orf_traitList,%alias)=();


gff();# deal with saccharomyces_cerevisiae.gff

read_qtlcart_map(); #marker<->number,gives $markerHash[$chr][$number]=$markername
read_hash_repeat();#name_to_num,num_to_name, ORF are numbered in cDNA array
snpfile();# deal with genotypes_forRelease_1_20_05.dat, get the position of markers


openFiles();

#if ($output_option{eORFlist}) {
#	open(TF_CSH,'>..\output\tf.csh') or die;
#	print TF_CSH 'perl analyze.pl gene_association.sgd 7272 ';
#}


my %trait;


opendir THISDIR, "$zFolder" or die "serious dainbramage: $!";
my @allfiles = grep /.(mim|z)$/, readdir THISDIR;
closedir THISDIR;

print "total number of z file analyzing: ",scalar @allfiles,"\n";
#die "number of z files is incorrect" if(@allfiles != 6195);

foreach my $input (@allfiles) {
	my $inTheRegion=0;
	my $startChr=0;
	my $maxLR=0;
	my @curr_marker_list=();
	%trait=();
	
	$input=~/(\d+)/ or die "I assume all output file names are gene+digit\n";
	$trait{'index'}=$1 or die;
	
unless ( open(ZMAP, $zFolder."/$input")  ){die "cannot open $zFolder/$input";}

line52: while (<ZMAP>) {
	if (/Analyzed trait/) {
			
			my @space=split(/\s+/);
			$space[4]=~ /(\d+)/ or die;
			
			if (exists $num_to_name{$1}){
				$trait{name}=$num_to_name{$1} or die;	
			}
			else{
				print 'gene'.$trait{index}.'.z is not yeast gene, overlooked',"\n";
				print NOTE 'gene'.$trait{index}.'.z is not yeast gene, overlooked',"\n";
				last line52;
			}

			init();

			next line52;
	}
	

	if ( /^[\s\d]/) {
		my($space,$chr,$markerIndex,$pos,$lr)=();
		if (/:/) {next line52;}

		if ($currentMethod eq 'MIM'){
			my $ic;
			($space,$chr,$markerIndex,$pos,$ic,$lr)=split/\s+/;
		}
		else{	
			if (/^\s/) {	($space,$chr,$markerIndex,$pos,$lr)=split/\s+/;		}

			if (/^\d/) {	($chr,$markerIndex,$pos,$lr)=split/\s+/;		}
		}

		$pos=$pos*100 unless($currentMethod eq 'MIM');
		
		if ($lr>$threshold-$CI) {
			
			if ($inTheRegion==0){
				$inTheRegion=1;
				$startChr=$chr;
				@curr_marker_list=();
				$maxLR=0;
			}
			else{	
				if ($startChr != $chr) {
					if ($maxLR>$threshold) {
#						$numQtl=$numQtl+1;
#						$currNumQtl++;
						map_marker_to_orf($maxLR,@curr_marker_list);
					}
					$startChr=$chr;
					@curr_marker_list=();
					$maxLR=0;
				}
			}

			push(@curr_marker_list,"$chr,$markerIndex,$pos,$lr");
			if ($lr>$maxLR) {
				$maxLR=$lr;
			}		
		}
		else{
			if ($inTheRegion==1) {
				if ($maxLR>$threshold) {
#					$numQtl=$numQtl+1;
#					$currNumQtl++;
					map_marker_to_orf($maxLR,@curr_marker_list);
				}
				$inTheRegion=0;
			}
		}	
	}
}

if ($inTheRegion==1) {
	if ($maxLR>$threshold) {
#		$numQtl=$numQtl+1;
#		$currNumQtl++;
		map_marker_to_orf($maxLR,@curr_marker_list);
	}
}
close(ZMAP);
OUTPUT_horizontal() if (exists $trait{numQTL});
}

OUTPUT_vertical();
closeFiles();



########################################################
#OUTPUT started here
#################################################

sub OUTPUT_vertical{

	open(ORF,$outputFolder."/$currentMethod".'_forEachORF.list') or die;
	print ORF "gene\tisTarget\tisFactor\tnumTarget\ttotalTrait\tcorrectTarget\tcorrectFactor\n";
	
	foreach my $pORF  (keys(%orf_traitList)) {
		if ($output_option{eORFlist} ) {
	
			open(forTF,">$output_option{eORFlist}/$pORF") or die;
			print forTF join("\n",keys(%{$orf_traitList{$pORF}})) ;
			close(forTF);
		}

		
		my ($correctTarget	,$correctFactor, $numTarget)=(0,0,0);
		
	if(@factor>0){			
		foreach my $eORF (keys (%{$orf_traitList{$pORF}})) {
			if ($orf_traitList{$pORF}{$eORF} and $orf_traitList{$pORF}{$eORF}=~/TARGET/) {
				$correctTarget++;
			}
			if ($orf_traitList{$pORF}{$eORF} and $orf_traitList{$pORF}{$eORF}=~/FACTOR/) {
				$correctFactor++;
			}
		}
		


		if (get_target($pORF)) {
			my @s=split(/\s+/,get_target($pORF));
			$numTarget=@s;
		}
	}
		print ORF $pORF."\t".((exists($targethash{$pORF}))?1:0)."\t".((exists($factorhash{$pORF}))?1:0)."\t";
		print ORF $numTarget."\t";
		print ORF scalar(keys(%{$orf_traitList{$pORF}}))."\t".$correctTarget."\t".$correctFactor."\n";
	}

	close(ORF);
}


sub OUTPUT_horizontal{
	
	my ($transIsCis,$targetIsCis)=(0,0);
	my $numORF=0;


	for (my $k=0;$k< $trait{numQTL};$k++) {
		$numORF=$numORF+scalar(@{$trait{qtl}->[$k]});
	
	if(@factor>0){	
		if ($trait{qtlID}->[$k]=~/TRANS/){
			$num_trans_QTL++
		}
			
		if ($trait{qtlID}->[$k]=~/TRANS/ and $trait{qtlID}=~/CIS/) {
			$transIsCis++;
		}

		if ($trait{qtlID}->[$k]=~/TARGET/ and $trait{qtlID}=~/CIS/) {
			$targetIsCis++;
		}

	}
	}

	print OUT "$trait{numQTL} Qtls are detected, containing $numORF ORFs \n";
	
	$num_no_orf++ if ($numORF==0) 	;
	$num_no_qtl++ if ($trait{numQTL}==0) ;

	if ($trait{numQTL}  && !$numORF) {
		print NOTE "$trait{name} QTL detected without ORF in it\n";
	}
	
	if (!$trait{numQTL}  && $numORF) {
		die "it is impossble";
	}
	

	print NUMBERS "$trait{name}\t$trait{numQTL}\t".(@{$trait{position}}>0)."\t$numORF\t";
#                    total num of trans factor       num of orf in trans qtl      num of orf in the qtl is in nuclear               			
	print NUMBERS $trait{numTarget}."\t".$trait{numFactor}."\t";
	print NUMBERS $trait{mapped_to_cis}."\t".$trait{mapped_to_trans}."\t".$trait{mapped_to_target}."\t";
	print NUMBERS $transIsCis."\t".$targetIsCis."\n";

	if (@{$trait{position}}>0) {
		if ($trait{closest}<10000000-0.1) {
			print OUT "the closest position to the boundary of $trait{name} is $trait{closest} \n";
					
			if ($trait{closest} == 0) {
				$num_exact++;
				print exact_match "$trait{name} \n";
				if ($trait{mapped_to_cis}==0) {
					die "inconsistency about cis mapping\n";
				}
			}
					
			$num_approxi++ if ($trait{closest}<=$vicinity && $trait{closest}>0  ) ;
			$num_same_chr++ if ($trait{closest}>$vicinity) 		;
				
		}
#############################################
#notice that when $numORF==0, it is still possible that closest <10000000
#closest is the distance between QTL boundary and unknow orf postion, it might be no ORF in the QTL
# on the other hand, $num_diff_chr might count less than necessary
		else{
			if ($trait{numQTL}>0){	
				$num_diff_chr++;
				print OUT "no QTL on the cis chromosome, but there are some ORF on other chromosome \n";
				print OUT "\n";
			}
		}
	}
	else{
		print OUT "$trait{name} position unknonw, because not in GFF \n \n";
		$num_unknown_location++;				
	}
}
###################################################
#end of OUTPUT()
##############################################


#########################################################################################
sub init{

	$trait{closest}=10000000;		
	$trait{numQTL}=0;
	$trait{mapped_to_cis}=0;
	$trait{mapped_to_trans}=0;
	$trait{mapped_to_target}=0;
	$trait{numFactor}	=0;
	$trait{numTarget}=0;

	if(exists($featurePos{$trait{name}})){
		$trait{position}=[split(/\s+/,$featurePos{$trait{name}})];
		if (@{$trait{position}}<3) {die;
		}
	}
	else{
		$trait{position}=[()];
	}

if(@factor>0){		
	if (get_factor($trait{name})) {
		$trait{factor}=get_factor($trait{name}) ;
		my @s=split(/\s+/,$trait{factor});
		$trait{numFactor}=@s;
	}

	if (get_target($trait{name})) {
		$trait{target}=get_target($trait{name});
		my @s=split(/\s+/,$trait{target});
		$trait{numTarget}=@s;
	}
		
}
	print "\n"."now mapping the expression control for $trait{name}";
	print OUT "\n"."now mapping the expression control for $trait{name} \n";
	print OUT "which is supposed to be contronlled by $trait{factor}\n" if($trait{factor});

}







#########################################################################################
sub map_marker_to_orf{
	my $peak=shift;
	my $currChr;
	
	my ($inTheRegion,@qtl_start,@qtl_stop)=();
	
	foreach  (@_) {
		my ($chr,$markerIndex,$pos,$lr)=split(/,/);
		if(!$currChr){
			$currChr=$chr;
		}
		else{
			if ($currChr != $chr) {die;	}
		}
	

		if ($lr>$peak-$CI) {
		
			if (!$inTheRegion) {
				my @section=getPosition($chr,$markerIndex,$pos);
				push(@qtl_start,$section[0]);
				if(@qtl_start>@qtl_stop){
					push(@qtl_stop,$section[0]);
				}
				else{
					die;
				}
				$inTheRegion=1;
			}
			else{
				if (@qtl_start==@qtl_stop){
					my @section=getPosition($chr,$markerIndex,$pos);
					$qtl_stop[-1]=$section[0];
				}
				else{
					die;
				}
			}
		}
		else{
			if ($inTheRegion) {
				$inTheRegion='';
#				my @section=getPosition($chr,$markerIndex,$pos);
#				push(@qtl_stop,$section[1]);
			}
		}

	}


	
	if ($inTheRegion) {
		my ($chr,$markerIndex,$pos,$lr)=split(/,/,$_[-1]);
		my @section=getPosition($chr,$markerIndex,$pos);
		push(@qtl_stop,$section[1]);
	}

	if (@qtl_start != @qtl_stop) {die;
	}

	my $qtl_length=$qtl_stop[-1]-$qtl_start[0];
	
	my $info_for_mim=$currChr."\t".$qtl_start[0]."\t".$qtl_stop[-1]."\t".$qtl_length."\t".$peak;

#get physical postion, start to obtain orf within the physical boundaries

	my $at_cis_now=0;
	my $at_trans_now=0;
	my $at_target_now=0;

	my %orf_at_this_qtl=();
	my $category;
	


	die unless ($qtl_length>0);

	for (my $k=0;$k<@qtl_start ;$k++) {
#		discontinous QTL is not allowed
#		$qtl_length=$qtl_length+$qtl_stop[$k]-$qtl_start[$k];

		
		foreach  (getFeature($qtl_start[$k],$qtl_stop[$k],$currChr)) {

		if(@factor>0){
			if (isFactorOfTrait($trait{name},$_)) {
				$at_trans_now=1;
				$trait{mapped_to_trans}++;
				$category.='TRANS';

				push(@transPair,"$trait{name} $_");
				print OUT "known regulation between $trait{name} and $_ is confirmed \n";
				$num_tf_confirmed++;
			}

			if (isTargetOfTrait($trait{name},$_)) {
				$at_target_now=1;
				$trait{mapped_to_target}++;
				$category.='TARGET';
			}
		}
		
			if ($_ eq $trait{name}) {
				$at_cis_now=1;
				$trait{mapped_to_cis}++;
				$category.='CIS';
			}
	
			$orf_at_this_qtl{$_}='';
		}
	}

	foreach  (keys %orf_at_this_qtl) {
		$orf_traitList{$_}{$trait{name}}=$category;
	}


	push(@{$trait{qtl}},[keys %orf_at_this_qtl]); # an empty list is allowed to be pushed into

	$trait{numQTL}++;

	if ($category) {
		push (@{$trait{qtlID}},$category);
	}
	else{
		push (@{$trait{qtlID}},$trait{numQTL});
	}
		
	print LEN "$trait{index}\t$trait{name}\t$trait{qtlID}->[-1]\t$at_cis_now\t$at_trans_now\t$at_target_now\t",
		$info_for_mim,"\n";
	
	if(scalar (keys (%orf_at_this_qtl))==0){
		print NOTE "trait $trait{name} 's qtl $trait{qtlID}->[-1] contains no ORF\n";
		return;
	}
		
	if ($output_option{pORFlist}) {
		output_list( $trait{qtl}->[-1], "$output_option{pORFlist}/$trait{name}".'_'.$trait{qtlID}->[-1] );
	}

}

sub output_list{
	my ($list,$fileName)=@_;
	open(OUTPUTLIST,">$fileName") || die;
	print OUTPUTLIST join("\n",@{$list});
	close(OUTPUTLIST);
}	



#############################################################################
			
#######################################################################
sub getFeature{
	my @thisQtl=();
	my($posB,$posE,$currC,$k)="";
	($posB,$posE,$currC)=@_;
	
	for ($k=0;$k<@{$featureList[$currC]};$k++) {
		if ( hit($featureBegin[$currC]->[$k],$featureEnd[$currC]->[$k],$posB,$posE)==1){
			push (@thisQtl,$featureList[$currC]->[$k]);
		}
	}


	if ($trait{position}->[0] && $trait{position}->[0] eq $currC) {
		
		if (hit($trait{position}->[1],$trait{position}->[2],$posB,$posE)==1) {
			$trait{closest}=0;
#at cis now will de determined after getting pORF list in a qtl region
#			$at_cis_now=1;
#			print "intersect with $r1 $r2 \n";
		}
		else{
			$trait{closest}=min($trait{closest},abs($posB-$trait{position}->[1]), abs($posE-$trait{position}->[2]), 
				abs($posB-$trait{position}->[2]),abs($posE-$trait{position}->[1]));		
		}
	}


	return @thisQtl;
}




sub getPosition{
	my($relative_genetic_postion,$gL,$pL)=();
	my($dist_from_currM,$dist_from_nextM,$r1,$r2)=();
	my ($currC,$currM,$currP)=@_;
	my $nextM=$currM+1;

	my $phy_zmap_begin=$phy_pos{$markerHash[$currC]->[$currM]} or die "chr $currC marker $currM";
	my $phy_zmap_end=$phy_pos{$markerHash[$currC]->[$nextM]} or die;
#gL for interval length in genetic map, pL for interval length in physical map

	$gL=$geneticMap[$currC]->[$nextM]-$geneticMap[$currC]->[$currM];
	$pL=$phy_zmap_end-$phy_zmap_begin;

	if ($gL<0 || $pL<0) {die;
	}
	
	$dist_from_currM=($currP-$geneticMap[$currC]->[$currM])/$gL;
	return ($phy_zmap_begin+$dist_from_currM*$pL);

=the following code is obsolete
	if ( ($currP-$geneticMap[$currC]->[$currM])>$stepsize/2.) {
		$dist_from_currM=$currP-$stepsize/2-($geneticMap[$currC]->[$currM]);
	}
	else{
		$dist_from_currM=($currP-$geneticMap[$currC]->[$currM])/2.;
	}
	
	if ( ($geneticMap[$currC]->[$nextM]-$currP)>$stepsize/2.) {
		$dist_from_nextM=$geneticMap[$currC]->[$nextM]-$currP-$stepsize/2.;
	}
	else{
		$dist_from_nextM=($geneticMap[$currC]->[$nextM]-$currP)/2.;
	}

	if ($dist_from_currM<-0.01){
		die $geneticMap[$currC]->[$currM]." ".$currP." ".$geneticMap[$currC]->[$nextM]; 
	}

	if($dist_from_nextM<-0.01) {
		die $geneticMap[$currC]->[$currM]." ".$currP." ".$geneticMap[$currC]->[$nextM]; 
	}

	$dist_from_currM=$dist_from_currM/$gL;
	$dist_from_nextM=$dist_from_nextM/$gL;
	
	$r1=$phy_zmap_begin+$dist_from_currM*$pL;
	$r2=$phy_zmap_end-$dist_from_nextM*$pL;
	
	return ($r1,$r2);
=cut
}


sub hit{ # this subroutine will judge whether two intervals have sufficient amount in common
	my($b1,$e1,$b2,$e2,$len1,$len2)="";
	
	($b1,$e1,$b2,$e2)=@_;
	
	if ( ($e1-$e2)*($b1-$b2)<0) { return 1;}
	
	$len1=$e1-$b1;$len2=$e2-$b2;

	if ($len1<0 || $len2<0) {die;
	}

	if ($e1>$e2 && ($e2-$b1)>min($len1,$len2)*$minimalPercentageAsHit ) {
		return 1;
	}

	if ($b1<$b2 && ($e1-$b2)>min($len1,$len2)*$minimalPercentageAsHit ) {
		return 1;
	}
	
	return 0;

}


################################################################################
#the following subroutines are to be called at the begining of the program to load informations
#########################################################################

sub gff{
my(@space)=();
my($chr)=();
unless (open (GFF,$currentFolder.'/saccharomyces_cerevisiae.gff')){die;}
while (<GFF>) {
	if (/FASTA/) {last;
	}

	if (/ID.*Name/) {
	
		@space=split/\t/;
		
		$space[0]=~/chr(\w+)/;
		if (isroman($1)) {
			$chr= arabic($1)
		}
		else{
			if ($1 ne "Mito") {
			print $_;
			die ;
			}
			next;
		}

		if ($chr==44) {
			die "line 204 $_";
		}

		$space[8]=~/ID=([-\.\)\(\w]+);.*Name=([-\.\)\(\w]+);/;
		if ($1 ne $2) {
			print $_."\n";
			die $1." ".$2;
		}
		
		my $id=$1;

		push(@{$featureList[$chr]},$id);
		push(@{$featureBegin[$chr]},$space[3]);
		push(@{$featureEnd[$chr]},$space[4]);
		$featurePos{$id}=join(" ",($chr,$space[3],$space[4]));

		if ($space[8]=~/Alias=([-,\.\)\(\w]+);/) {
			my $alias=$1;
			foreach my $orf (split(/,/,$alias)) {

				if ($orf=~ /Y[A-P][LR]\d{3}[WC]/) {
					if ($orf=~ /Y[A-P][LR]\d{3}[WC]\w/) {
						$orf=~ s/(Y[A-P][LR]\d{3}[WC])(\w)/$1-$2/;
					}

		#			$featurePos{$orf}=join(" ",($chr,$space[3],$space[4]));
					
				}
				$alias{$orf}=$id;
			}
		}

	}
}
close(GFF);
}

################################################################################

sub snpfile{

#_snp means hash is constructed from data read out of snp_processed.txt
unless(open (SNP,$currentFolder.'/'.$marker_physical_position)){die;}
while (<SNP>) {
	chomp;
	next if(/^\#/);
	next if(/probePairKey/);
	
	my @s=split(/\t/);

	my $id;

	if ($currentMethod eq 'MIM'){
		$id='M'.$.
	}

	else{
		$id='m_';
		$id.=$.-27;
	}

	$phy_pos{$id}=$s[-3];
}
close(SNP);
}

################################################################################
sub read_hash_repeat{
	open (HASH,$currentFolder.'/'.$probe_id_to_gene_name) or die;

	while (<HASH>){
		chomp;
		my @space=split(/ /);
		foreach (@space) {
			s/-$//;
		}

		if ($alias{$space[0]}) {
			$space[0]=$alias{$space[0]};
		}

		$name_to_num{$space[0]}=$space[1];
		for (my $k=1;$k<@space;$k++) {
			$num_to_name{$space[$k]}=$space[0];
#			$num_to_num{$space[$k]}=$space[1];
		}

		
	}	
	close(HASH);
#	print "total number of genes in hash.repeat: ".scalar(keys(%name_to_num))."\n";
}

sub read_qtlcart_map{

my $started =0;
my (@num_marker,$numChr,$numMarker)=();
unless(open(MPS,$currentFolder.'/'.$map_file_name)){die;}
while (<MPS>) {
	chomp;
	if (/^-c/) {
		my @s=split(/\s+/);
		$numChr=$s[1] or die;

		for (my $chr=1;$chr<=$numChr ;$chr++) {
			push(@{$geneticMap[$chr]},-1);
			push(@{$geneticMap[$chr]},0);
		}
	}

	if (/^-i/) {
		my @s=split(/\s+/);
		$numMarker=$s[1] or die;
	}

	if (/^-Number/) {
		s/(.*)\|//;
		die 'check the index of marker' unless (/^ /);
		my @num_marker=split(/\s+/);
		
		$geneticMap[0]=[@num_marker]; #$geneticMap[0][0]=0
		
		for (my $chr=1;$chr <= $numChr ;$chr++) {
			pop(@{$geneticMap[$chr]});
		}
		die "check total number of chr " unless(((scalar @num_marker) -1) == $numChr);

		for (my $chr=1;$chr<=$numChr;$chr++){
			die "inconsistent marker number for chr $chr\n"
			 unless(@{$geneticMap[$chr]}==(1+$num_marker[$chr]));
		}
	}

	if (/^-l\s+(\d+)\s+/) {
		die unless($numChr);
		next if ($1==0);
		s/(.*)\|//;		
		die 'check the index of marker' unless (/^ /);
		for (my $chr=1;$chr <= $numChr ;$chr++) {
			my $n=substr($_,8*($chr-1)+2,8);
			next unless($n=~/\./);
			$geneticMap[$chr][-1]=$begin if($geneticMap[$chr][-1]<0.00005 and $n>0.00005);
			push(@{$geneticMap[$chr]},$n+$geneticMap[$chr][-1]);
		}	
	}	
		
	if (/-b MarkerNames/) { 
		$started=1;
		next;
	}
	
	if ($started==0) {next;
	}
	
	if (/-e MarkerNames/) {last;
	}


	my @space=split/\s+/;

	$markerHash[$space[1]]->[$space[2]]=$space[3];
}
close (MPS);
}

sub openFiles{

#begin to handle the output file of Zqtlmap
open(OUT,$outputFolder."/$currentMethod".'_trait.out') or die;
open(exact_match,$outputFolder."/$currentMethod".'_exact_match.list') or die;
open(NUMBERS,$outputFolder."/$currentMethod".'_numbers.list');
	print NUMBERS "trait\tnumQtl\tknownPos\tnumORF\t";
	print NUMBERS "numTarget\tnumFactor\t";
	print NUMBERS "mappedToCis\tmappedToTrans\tmappedToTarget\t";
	print NUMBERS "transIsCis\ttargetIsCis\n";

open(LEN,$outputFolder."/$currentMethod".'_qtl.len') or die;
print LEN "index\ttrait\tqtlID\tat_cis_now\tat_trans_now\tat_target_now\tqtl_chr\tqtl_start\tqtl_stop\tqtl_length\tmaxLR\n";

open(NOTE,$outputFolder."/$currentMethod".'_mapping.log') or die;
}


sub closeFiles{

print OUT "number of trait for which no QTL is mapped to cis chromosome: $num_diff_chr \n";
print OUT "number of trait for which a QTL is mapped faraway but on the same chromosome: $num_same_chr \n";
print OUT "number of trait for which a QTL is mapped to a vicinity of $vicinity:  $num_approxi \n";
print OUT "number of trait for which a QTL is mapped to itself: $num_exact \n";
print OUT "number of trait with no QTL mapped: $num_no_qtl \n";
print OUT "number of trait with no ORF mapped: $num_no_orf \n";
print OUT "number of trait not in GFF and thus location unknown: $num_unknown_location \n";

close(OUT);

close(exact_match);
close(NUMBERS);

close(LEN);
close(NOTE);
}

sub min {
    my ($runningMin, @rest) = @_;
    foreach my $x (@rest) {
        $runningMin = $x if $x < $runningMin;
    }
    return $runningMin;
}

#####################
# the following subroutines are out-dated ones.
################################################################################

sub isTargetOfTrait{
	my ($eORF,$pORF)=@_;
	if ($trait{target} and "$trait{target} "=~ /$pORF /) {
			return 1;
	}
	return 0;
}

################################################################################
sub isFactorOfTrait{
	my ($eORF,$pORF)=@_;
	if ($trait{factor} and "$trait{factor} "=~ /$pORF /) {
	   return 1;
	}
	return 0;
}

################################################################################
sub get_factor{
	my $target=$_[0];
	my $k ="";
	my @list =();
	for ($k=0;$k<@target;$k++) {
			if ($target eq $target[$k]) {
			push(@list,$factor[$k]);
		}
	}
	return join(" ",@list);
}

################################################################################
sub get_target{
	my $factor=$_[0];
	my $k ="";
	my @list=();
	for ($k=0;$k<@factor;$k++) {
		if ($factor eq $factor[$k]) {
			push(@list,$target[$k]);
		}
	}
	return join(" ",@list);
}


