#!/usr/bin/perl -w

# MACHINES CONFIGURATION; ONLY ONE OF THESE TWO -- OR NONE
$SGE_CLUSTER = 0;
$ABACUS_CPU = 0;
$ABACUS_CPU2 = 1;


##################################################################################
# Summarizes the outputs of pundit-vx.pl    (x=1..3)                             #
# Oct. 2002 - Bioinformatics Research Center - North Carolina State University   #
#                                                                                #
# Comments to sarisbro at uottawa dot ca - Thanks.                               #
##################################################################################




#################################################################################
use POSIX;

if($SGE_CLUSTER){
	$curdir = "/home/stephane/work/Pandit62";
}elsif($ABACUS_CPU){
	$curdir = "/home/stephane/work/Pandit62/FlatFiles";
}elsif($ABACUS_CPU2){
	$curdir = "/home/stephane/work/Pandit62/Pandit62";
}else{
	$curdir = "D:/data/Pandit62";
}
$infile = "Pandit6.2.txt";
$size_of_lumps = 5;
$total_pf_number = 0;
$min_ns=2;  $min_nb_codons=20;  $min_ns_nc=150;

########################### MAIN ###############################################

mkdir("$curdir/res_sum",0777) || warn "\ncannot create res_sum dir: $!";
open(YN00FILE,">$curdir/res_sum/YN00PF.dat") || warn "\ncould not create YN00PF file... $!\n";
print YN00FILE "PF_name   seq  seq      S       N        t   kappa   omega     dN    SEdN    dS    SEdS   PF_real_name  NS   LS\n";
open(YN00FILEname,">$curdir/res_sum/YN00PFname.dat") || warn "\ncould not create YN00PFname file... $!\n";
print YN00FILEname "PFid\tPFns\tPFls\tPF_name\n";
open(YN00FILEname50,">$curdir/res_sum/YN00PFname50.dat") || warn "\ncould not create YN00PFname50 file... $!\n";
print YN00FILEname50 "PFid\tPFns\tPFls\tPF_name\n";
open(LRTFILE,">$curdir/res_sum/LRT.dat") || warn "\ncould not open LRT file... $!\n";
print LRTFILE "PF_name\tLnl-M0\tnp-M0\tLnl-M1\tnp-M1\tLnl-M2\tnp-M2\tLnl-M3\tnp-M3\tLnl-Branch\tnp-Branch\tdeltaM0_B\tdf\tp-value(XLS:M0.B)\tdeltaM0_M3\tdf\tp-value(XLS:M0.M3)\tdeltaM2_M1\tdf\tp-value(XLS:M2.M1)\tPFname\n";
open(PARAMFILE,">$curdir/res_sum/paramSITES.dat") || warn "\ncould not open params file... $!\n";
print PARAMFILE "PF_id\tM0-w\tM1-p0\tM1-p1\tM1-w0\tM1-w1\tM2-p0\tM2-p1\tM2-p2\tM2-w0\tM2-w1\tM2-w2\tM2-PSsites\tM3-p0\tM3-p1\tM3-p2\tM3-w0\tM3-w1\tM3-w2\tM3-PSsites\tPF_name\n";
open(PAIRm0FILE,">$curdir/res_sum/paramPairM0SITES.dat") || warn "\ncould not open params pairs file... $!\n";
print PAIRm0FILE "PF_name\tSP_pair\tlnL-M0\tt-M0\tS-M0\tN-M0\tomega-M0\tdN-M0\tdS-M0\n";
open(PAIRm1FILE,">$curdir/res_sum/paramPairM1SITES.dat") || warn "\ncould not create params pairs file... $!\n";
print PAIRm1FILE "PF_name\tSP_pair\tns\tls\tlnL-M1\tPFname\tt-M1\tS-M1\tN-M1\tomega-M1\tdN-M1\tdS-M1\n";
open(BrFILE,">$curdir/res_sum/Branch.dat") || warn "\ncould not creaet Branch file... $!\n";
print BrFILE "PF_name\tdelta\tnp\tp-value\tbranch\t t\t S\t N\t dN/dS\t dN\t dS\t S*dS\t N*dN\n";
open(BrFILE_parsed,">$curdir/res_sum/parsed_Branch.dat") || warn "\ncould not create Branch file... $!\n";
print BrFILE_parsed "PF_name\tdelta\tnp\tp-value\tbranch\t t\t S\t N\t dN/dS\t dN\t dS\t S*dS\t N*dN\n";
open(BrFILE_name,">$curdir/res_sum/Branch_name.dat") || warn "\ncould not create Branch name file... $!\n";
print BrFILE_name "PFid\tPFns\tPFls\tPFname\tdelta\tdf\n";
open(BrFILE_name50,">$curdir/res_sum/Branch_name50.dat") || warn "\ncould not create Branch name 50 file... $!\n";
print BrFILE_name50 "PFid\tPFns\tPFls\tPFname\tdelta\tdf\n";
open(PFname_LIST,">$curdir/res_sum/pf_names.dat") || warn "\ncould not create PF name list... $!\n";
open(LOGFILE,">$curdir/res_sum/log.txt") || warn "\ncould not create log file... $!\n";

close(LRTFILE);
close(PARAMFILE);
close(PAIRm0FILE);
close(PAIRm1FILE);


opendir(PANDIT, "$curdir/data")  || warn "\nCan't open data directory: $!";
@tmp = sort readdir(PANDIT);
closedir(PANDIT);
foreach(<@tmp>) { if(/^PF/) { push(@folders, $_); } }   @tmp=();
$n_folders=@folders;

print ("\n$n_folders data sets to be summarized:\n @folders\n");


# summarizes raw outfiles
SumYang_Nielsen2000();
SumSite_Models();
SumBranch_Models();
# SumBranch_Site_Models()

# parse summary files for possible errors
parse_YN00();

close YN00FILE;
close YN00FILEname;
close YN00FILEname50;
close DATAFILE;
close BrFILE;
close BrFILE_parsed;
close BrFILE_name;
close BrFILE_name50;
close PFname_LIST;
close LOGFILE;

print "\n\n\aAll done!\n";

#################################################################################
#################################################################################
#################################################################################
sub parse_YN00{
	open(RESFILE,"$curdir/res_sum/YN00PF.dat") || die "\ncannot open YN00 res file: $!";
	open(OUTFILE,">$curdir/res_sum/parsed_YN00PF.dat") || die "\ncannot create YN00 outfile: $!";
	# print OUTFILE "PFname\tPFns\tPFls\tt\tkappa\tomega\tdN\tdS\n";
	$nb_lines = 0;
	$PSS_thisPF = 0;
	$new_PF_100 = 0;
	$new_PF_50 = 0;
	$this_PF = 0;

	while(<RESFILE>){
		chomp;
		$_ =~ s/  */:/g;
		($PFname,$seq1,$seq2,$S,$N,$t,$kappa,$omega,$dN,$SEdN,$dS,$SEdS,$PFreal_name,$PFns,$PFls) = split(/:/);
		if($PFname != $this_PF){
			$this_PF = $PFname;
			$PSS_thisPF = 0;
			$new_PF_100 = $new_PF_50 = 1;

		}
		if($omega<99.0 && $t<39.0 && $dN<10 && $dS<10 && ($PFname && $PFns && $PFls && $t>0 && $kappa>0 && $omega>0 && $dN>0 && $dS>0)){
			# print OUTFILE "$PFname\t$PFns\t$PFls\t$t\t$kappa\t$omega\t$dN\t$dS\n";
			print OUTFILE "$PFname  $PFns  $PFls  $t  $kappa  $omega  $dN  $dS\n";
			if($omega>1){
				$PSS_thisPF ++;
				if($new_PF_100){
					# print YN00FILEname "$PFname\t$PFns\t$PFls\t$PFreal_name\n";
					print YN00FILEname "$PFname  $PFns  $PFls  $PFreal_name\n";
					$new_PF_100 = 0;
				}
				if(($PSS_thisPF > 0.5*($PFns*($PFns-1)/2))&&($new_PF_50)){
					# print YN00FILEname50 "$PFname\t$PFns\t$PFls\t$PFreal_name\n";
					print YN00FILEname50 "$PFname  $PFns  $PFls  $PFreal_name\n";
					$new_PF_50 = 0;
				}
			}
			$nb_lines++;
		}
	}
	print "\a\nYN00 parsed! ($nb_lines lines in fine)";

	close RESFILE;
	close OUTFILE;
}
#################################################################################
sub KeepDataSet{
	$data_file = shift @_;
	open(SEQFILE,"$data_file") || warn "\ncould not open YN00PF file... $!\n";

	chomp($line = <SEQFILE>);
	($ns,$ls)=split(/\s+/, $line);
	if($ns<$min_ns || $ls<$min_nb_codons*3 || $ns*$ls < $min_ns_nc*3) {
		$okay = 0;
	}else{
		$okay = 1;
	}
	close NUCSEQ;
	return $okay;
}
#################################################################################
sub Get_PF_ns{
	$infile = shift @_;

	open(MYDNATREE,"$infile") || warn "\ncould not open $infile ... $!\n";
	chomp($ns = <MYDNATREE>);
	$ns =~ s/\s//g;
	close(MYDNATREE);
	return $ns;
}
#################################################################################
sub Get_PF_ls{
	$infile = shift @_;

	open(MYDNATREE,"$infile") || warn "\ncould not open $infile ... $!\n";
	chomp($line = <MYDNATREE>);
	($ns,$ls)=split(/\s+/, $line);
	close(MYDNATREE);
	return $ls;
}
#################################################################################
sub Get_PF_name{
	$infile = shift @_;

	open(MYSPNAME,"$infile") || warn "\ncould not open $infile ... $!\n";
	chomp($real_name = <MYSPNAME>);
	$real_name =~ s/\s/_/g;
	$real_name =~ s/\W*//g;
	$real_name =~ s/\n*//g;
	close MYSPNAME;
	return $real_name;
}
#################################################################################
sub SumYang_Nielsen2000{

	print "\n";
	for($count=1; $count<$n_folders; $count++) {
	# for($count=$start_at_PF;$count<=$end_at_PF;$count++){
		print "\rNow doing SumYang_Nielsen2000: PF$count";
		chdir "$curdir/data/PF$count" || warn "can't cd to $curdir/data/PF$count: $!";
		if(-e "TimeYN00.txt"){
			$pfid = Get_PF_id("$curdir/data/PF$count/spnames.txt");
			$PFname = Get_PF_name("$curdir/data/PF$count/PF$pfid.txt");
			$ns = Get_PF_ns("$curdir/data/PF$count/dnatree.txt");
			$seq_len = Get_PF_ls("$curdir/data/PF$count/dna-aln.nuc");
			print PFname_LIST "$pfid\t$PFname\n";

			# runs YN00 and parse result files
			$keep_it = KeepDataSet("$curdir/data/PF$count/dna-aln.nuc");
			if($keep_it){
				Parse_YN00_OUT("$curdir/data/PF$count/yn.out","$pfid","$PFname","$ns","$seq_len");
			}else{
				print LOGFILE "PF$pfid (in dir PF$count)too short - discard entries!\n";
			}
		}
	}
}
#################################################################################
sub SumSite_Models{

	print "\n";
	for($count=1; $count<$n_folders; $count++) {
	#for($count=1; $count<floor($n_folders/2); $count++) {
	#for($count=floor($n_folders/2); $count<$n_folders; $count++) {
	# for($count=$start_at_PF;$count<=$end_at_PF;$count++){
		open(LRTFILE,">>$curdir/res_sum/LRT.dat") || warn "\ncould not open LRT file... $!\n";
		open(PARAMFILE,">>$curdir/res_sum/paramSITES.dat") || warn "\ncould not open params file... $!\n";
		open(PAIRm0FILE,">>$curdir/res_sum/paramPairM0SITES.dat") || warn "\ncould not open params pairs file... $!\n";
		open(PAIRm1FILE,">>$curdir/res_sum/paramPairM1SITES.dat") || warn "\ncould not create params pairs file... $!\n";
		print "\rNow doing SumSite & SumBranch Models: PF$count";
		chdir "$curdir/data/PF$count" || warn "can't cd to $curdir/data/PF$count: $!";
		$pfid = Get_PF_id("$curdir/data/PF$count/spnames.txt");
		if((-e "TimeSiteModels.txt") && (-s "PF$pfid-M0.out")){

			# print result summary files
			print LRTFILE "$pfid\t";
			print PARAMFILE "$pfid\t";
			$lnLBr = 0.0; $npBr = 0.0;
			$lnLM0 = 0.0; $npM0 = 0.0;
			$lnLM1 = 0.0; $npM1 = 0.0;
			$lnLM2 = 0.0; $npM2 = 0.0;
			$lnLM3 = 0.0; $npM3 = 0.0;
			for($id=0;$id<4;$id++) {
				# Global analyses
				$keep_it = KeepDataSet("$curdir/data/PF$count/dna-aln.nuc");
				if($keep_it){
					($lnL,$np) = get_Lnl_from_outfile("$curdir/data/PF$count/PF$pfid-M$id.out");
					if($id == 1){
						$lnLM1 = $lnL; $npM1 = $np;
					}elsif($id == 2){
						$lnLM2 = $lnL; $npM2 = $np;
					}elsif($id == 3){
						$lnLM3 = $lnL; $npM3 = $np;
					}
					print LRTFILE "$lnL\t$np\t";
					if($id){
						$parms = get_params_from_outfile("$curdir/data/PF$count/PF$pfid-M$id.out");
					}else{
						$parms = get_dNdS_rate("$curdir/data/PF$count/PF$pfid-M$id.out");
						$lnLM0 = $lnL; $npM0 = $np;
					}
					$parms =~ s/  *//g;
					$parms =~ s/\t/ /g;
					if($id == 2 || $id == 3){
						$pss = getPosSelSites("$curdir/data/PF$count/PF$pfid-M$id.out");
						$pss =~ s/\n*//g;
						$parms .= ' '.$pss.' ';
					}
					print PARAMFILE "$parms";
					# site models - pairwise for M0 and M1
					if(($id==0) || ($id==1)){
						Summarize_PairwiseML($pfid,"$curdir/data/PF$count/PF$pfid-pM$id.out",$id,$count);
					}
				}
			}
			# Branch models and LRT M0 vs. Branch model
			if(-e "TimeBranchModels.txt"){
				if($keep_it){
					($lnL,$np) = get_Lnl_from_outfile("$curdir/data/PF$count/PF$pfid-Branch.out");
					$lnLBr = $lnL; $npBr = $np;
					print LRTFILE "$lnL\t$np";
				}
			}else{
				print LRTFILE "0\t0";
			}
			$PFname = Get_PF_name("$curdir/data/PF$count/PF$pfid.txt");
			print PARAMFILE "\t$PFname\n";
			if($lnLBr && $lnLM0){
				# Branch ("free-ratio") model vs. M0 ("one-ratio"): a test for variable selective pressure among lineage
				$delta = 2.0*($lnLBr - $lnLM0);
				$df = $npBr - $npM0;
				print LRTFILE "\t$delta\t$df\t";
			}else{
				print LRTFILE "\t0\t0\t";
			}
			if($lnLM3 && $lnLM0){
				# M3 vs. M0: a test for variable selective pressure among sites
				$delta = 2.0*($lnLM3 - $lnLM0);
				$df = $npM3 - $npM0;
				print LRTFILE "\t$delta\t$df\t";
			}else{
				print LRTFILE "\t0\t0\t";
			}
			if($lnLM2 && $lnLM1){
				# M2 vs. M1: a test for positive selection
				$delta = 2.0*($lnLM2 - $lnLM1);
				$df = $npM2 - $npM1;
				print LRTFILE "\t$delta\t$df\t";
			}else{
				print LRTFILE "\t0\t0\t";
			}
			$PFname = Get_PF_name("$curdir/data/PF$count/PF$pfid.txt");
			print LRTFILE "\t$PFname\n";
		}
		close(LRTFILE);
		close(PARAMFILE);
		close(PAIRm0FILE);
		close(PAIRm1FILE);
	}
}
#################################################################################
sub SumBranch_Models{

	print "\n";
	for($count=1; $count<$n_folders; $count++) {
		print "\rNow doing SumBranch_Models: PF$count";
		chdir "$curdir/data/PF$count" || warn "can't cd to $curdir/data/PF$count: $!";
		$after_dNdS_line = 0;
		$pfid = Get_PF_id("$curdir/data/PF$count/spnames.txt");
		if(-e "PF$pfid-Branch.out"){

			$PFname = Get_PF_name("$curdir/data/PF$count/PF$pfid.txt");
			$PFns = Get_PF_ns("$curdir/data/PF$count/dnatree.txt");
			$ls = Get_PF_ls("$curdir/data/PF$count/dna-aln.nuc");
			$new_PF_100 = $new_PF_50 = 1;
			$PSS_thisPF = 0;

			# LRT for branch model vs. M0 (const. rate across sites and lineages)
			$lnLBr = 0.0; $npBr = 0.0;
			($lnLBr,$npBr) = get_Lnl_from_outfile("$curdir/data/PF$count/PF$pfid-Branch.out");
			$lnLM0 = 0.0; $npM0 = 0.0;
			($lnLM0,$npM0) = get_Lnl_from_outfile("$curdir/data/PF$count/PF$pfid-M0.out");
			$delta = 2.0*($lnLBr - $lnLM0);
			$df = $npBr - $npM0;

			# gets rates for branches
			if(-e "$curdir/data/PF$count/PF$pfid-Branch.out"){
				@INFILE = get_file_data("$curdir/data/PF$count/PF$pfid-Branch.out");
				foreach $line (@INFILE) {
					if($after_dNdS_line){
						if($line =~ /^  /){
							$line =~ s/  */:/g;
							($trash,$br,$t,$S,$N,$omega,$dN,$dS,$SdS,$NdN) = split(/:/,$line);
							print BrFILE "$pfid\t$delta\t$df\t\t$br\t$t\t$S\t$N\t$omega\t$dN\t$dS\t$SdS\t$NdN";
							if(($t>0)&&($dS>0.01)&&($delta>=0)){
								print BrFILE_parsed "$pfid\t$delta\t$df\t\t$br\t$t\t$S\t$N\t$omega\t$dN\t$dS\t$SdS\t$NdN";
								if($omega>1 && ($t<10 && $dS<10 && $dN<10)){
									$PSS_thisPF ++;
									if($new_PF_100){
										print BrFILE_name "PF$pfid\t$PFns\t$ls\t$PFname\t$delta\t$df\n";
										$new_PF_100 = 0;
									}
									if(($PSS_thisPF > 0.5*($PFns*($PFns-1)/2))&&($new_PF_50)){
										print BrFILE_name50 "PF$pfid\t$PFns\t$ls\t$PFname\t$delta\t$df\n";
										$new_PF_50 = 0;
									}
								}
							}
						}
					}
					elsif($line =~ /^ branch  /) {
						$after_dNdS_line = 1;
					}
				}
			}
		}
	}
}
#################################################################################
sub Branch_Site_Models{

	for($count=1; $count<$n_folders; $count++) {
	# for($count=$start_at_PF;$count<=$end_at_PF;$count++){
		chdir "$curdir/data/PF$count" || warn "can't cd to $curdir/data/PF$count: $!";
		$pfid = Get_PF_id("$curdir/data/PF$count/spnames.txt");
	}
	PrintTimeStamp("TimeBranchSiteModels","Branch-site models for",$pfid);

}
#################################################################################
sub Summarize_PairwiseML_old{ # site models under M0 and M1
	$pfid = shift @_;
	$paiwiseMLoufile = shift @_;
	$id = shift @_;
	$count = shift @_;
	$start_reading_here = 0;

	if(-e $paiwiseMLoufile){
		@INFILE = get_file_data("$paiwiseMLoufile");
		foreach $line (@INFILE){
			if($line =~ /^pairwise comparison/){
				$start_reading_here = 1;
			}elsif($start_reading_here && $line =~ /\(/){
				$line =~ s/\s*//g;
				$line =~ s/\.\.\./../;
				$pair = $line;
			}elsif($start_reading_here && $line =~ /lnL/){
				$line =~ s/^lnL =//;
				$line =~ s/\s*//g;
				$lnL = $line;
			}elsif($start_reading_here && $line =~ /t=/){
				$line =~ s/t=//;
				$line =~ s/S=//;
				$line =~ s/N=//;
				$line =~ s/dN\/dS=//;
				$line =~ s/dN=//;
				$line =~ s/dS=//;
				$estimates = $line;
				# $estimates =~ s/\n//g;
				$filehandle = "PAIRm".$id."FILE";
				if($id){
					$PFname = Get_PF_name("$curdir/data/PF$count/PF$pfid.txt");
					$ns = Get_PF_ns("$curdir/data/PF$count/dnatree.txt");
					$ls = Get_PF_ls("$curdir/data/PF$count/dna-aln.nuc");
					print $filehandle "$pfid  $pair  $ns  $ls  $lnL  $PFname  $estimates";
				}else{
					print $filehandle "$pfid  $pair  $lnL  $estimates";
				}
			}
		}
	}else{
		print LOGFILE "$paiwiseMLoufile does not exist!\n";
	}
}
#################################################################################
sub Summarize_PairwiseML{ # site models under M0 and M1
	$pfid = shift @_;
	$paiwiseMLoufile = shift @_;
	$id = shift @_;
	$count = shift @_;
	$start_reading_here = 0;

	if(-e $paiwiseMLoufile){
		open(INFILE,"$paiwiseMLoufile") || warn "\ncould not open $paiwiseMLoufile... $!\n";
		while(<INFILE>){
			chomp;
			$line = $_;
			if($line =~ /^pairwise comparison/){
				$start_reading_here = 1;
			}elsif($start_reading_here && $line =~ /\(/){
				$line =~ s/\s*//g;
				$line =~ s/\.\.\./../;
				$pair = $line;
			}elsif($start_reading_here && $line =~ /lnL/){
				$line =~ s/^lnL =//;
				$line =~ s/\s*//g;
				$lnL = $line;
			}elsif($start_reading_here && $line =~ /t=/){
				$line =~ s/t=//;
				$line =~ s/S=//;
				$line =~ s/N=//;
				$line =~ s/dN\/dS=//;
				$line =~ s/dN=//;
				$line =~ s/dS=//;
				$estimates = $line;
				# $estimates =~ s/\n//g;
				$filehandle = "PAIRm".$id."FILE";
				if($id){
					$PFname = Get_PF_name("$curdir/data/PF$count/PF$pfid.txt");
					$ns = Get_PF_ns("$curdir/data/PF$count/dnatree.txt");
					$ls = Get_PF_ls("$curdir/data/PF$count/dna-aln.nuc");
					print $filehandle "$pfid  $pair  $ns  $ls  $lnL  $PFname  $estimates\n";
				}else{
					print $filehandle "$pfid  $pair  $lnL  $estimates\n";
				}
			}
		}
	}else{
		print LOGFILE "$paiwiseMLoufile does not exist!\n";
	}
}
#################################################################################
sub PrintTimeStamp{
	$filename = shift @_;
	$message = shift @_;
	$pfid = shift @_;

	($sec,$min,$hour,$mday,$mon,$year,$wday) = gmtime(time);
	$year += 1900; $hour -= 5;
	open(TIMESTAMP,">$curdir/data/PF$count/$filename.txt") || warn "\ncannot create TimeStamp for $pfid: $!";
	print TIMESTAMP "$message Protein Family PF$pfid completed on\n $mon.$mday.$year ($wday-th day of the week) at $hour:$min:$sec EST(USA).";
	close TIMESTAMP;
}
#################################################################################
sub Get_PF_id{
	$SP_file = shift @_;

	if(-e $SP_file)
	{
		@INFILE = get_file_data($SP_file);
		foreach $line (@INFILE) {
			if($line =~ /^PF/){
				$line =~ s/\D//;
				$line =~ s/\s*//g;
				$line =~ s/F//;
				$line =~ s/://;
				$id = $line;
			}
		}
	}else{
		$id = 0;
	}
	return $id;
}
#################################################################################
sub Parse_YN00_OUT{
	$infile = shift @_;
	$pfid = shift @_;
	$PFname = shift @_;
	$ns = shift @_;
	$ls = shift @_;
	$start_next_line = 0;

	if(-e $infile){
		@INFILE = get_file_data($infile);
		foreach $line (@INFILE) {
			if($line =~ /^seq. seq./){
				$start_next_line = 1;
			}
			elsif($start_next_line && $line =~ /\+\-/){
				$line =~ s/\+\-/ /; # removes the 1st "+-" of each result line
				$line =~ s/\+\-/ /; # removes the 2nd "+-"
				if($line =~ /99.0000/){
					print LOGFILE "99.000 at PF$pfid in YN00\n";
				}else{
					$line =~ s/\n//g;
					print YN00FILE "$pfid    ".$line."    $PFname    $ns    $ls\n" if $line;
				}
			}
		}
	}else{
		print LOGFILE "$infile does not exist!\n";
	}

}
#################################################################################
sub get_dNdS_rate{ # under M0
	$infile = shift @_;
	$after_dNdS_line = 0;

	if(-e $infile){
		@INFILE = get_file_data($infile);
		foreach $line (@INFILE) {
			if($after_dNdS_line){
				if($line =~ /^  /){
					$line =~ s/  */:/g;
					($trash,$br,$t,$S,$N,$omega,$dN,$dS,$SdS,$NdN) = split(/:/,$line);
				}
			}
			elsif($line =~ /^ branch/) {
				$after_dNdS_line = 1;
			}
		}
	}else{
		print LOGFILE "$infile does not exist!\n";
		$omega = '';
	}
	return $omega;
}
#################################################################################
sub get_params_from_outfile { # under M1 to M3
	$infile = shift @_;
	$after_dNdS_line = 0;
	$params = '';

	if(-e $infile){
		@INFILE = get_file_data($infile);
		foreach $line (@INFILE) {
			if($after_dNdS_line){
				if($line =~ /^p:/){
					$line =~ s/^p://;
					$line =~ s/\n/ /;
					$line =~ s/  */\t/g;
					$params .= $line;
				}
				elsif($line =~ /^w:/){
					$line =~ s/^w://;
					$line =~ s/\n/ /;
					$line =~ s/  */\t/g;
					$params .= $line;
				}
			}
			elsif($line =~ /^dN\/dS/) {
				$after_dNdS_line = 1;
			}
		}
	}else{
		print LOGFILE "$infile does not exist!\n";
		$params = '';
	}
	return $params;
}
#################################################################################
sub getPosSelSites{ # under M2 and M3 only! and sites w/ post proba > .95 only.
	$infile = shift @_;
	$after_PSS_line = 0;
	$first_pss_entry = 1;
	$pss = '';

	if(-e $infile){
		@INFILE = get_file_data($infile);
		foreach $line (@INFILE) {
			if($after_PSS_line){
				if(!($line =~ /Note/)){
					if($line =~ /\*\*/ || $line =~ /\*/ ){# sites with post proba > .95
						if($line =~ / 0./){
							($left,$right) = split(/ 0\./,$line);
						}else{
							($left,$right) = split(/ 1\./,$line);
						}
						$left =~ s/  *//g;
						$right =~ s/\d//g; # to stick back the right number of stars (significance level)
						if($first_pss_entry){
							$pss = $left.$right;
							$first_pss_entry = 0;
						}else{
							$pss .= ','.$left.$right;
						}
					}elsif(!$pss){# case of sites with post proba < .95 when only one site is selected in the whole EB analysis
						$pss .= '-';
					}
				}
			}
			elsif($line =~ /^Positively selected sites/){
				$after_PSS_line = 1;
			}
		}
		if(!$after_PSS_line){# implies some PSS have not been detected
			$pss = '-';
		}
	}
	return $pss;
}
#################################################################################
sub get_Lnl_from_outfile {
	$infile = shift @_;

	if(-e $infile){
		@INFILE = get_file_data($infile);
		foreach $line (@INFILE) {
			if($line =~ /^lnL/) {
				($begin,$nt,$np,$lnlsd) = split(/:/,$line);
				$np =~ s/\D//;
				$np =~ s/\)//;
				($lnl,$sd) = split(/\+/,$lnlsd);
			}
		}
	}else{
		print LOGFILE "$infile does not exist!\n";
		$lnl = ''; $np = '';
	}
	return ($lnl,$np);
}
#################################################################################
# The following routine is adapted from:
# Beginning Perl for Bioinformatics
# by James Tisdall
# published by O'Reilly & Associates
# (c) 2001 James Tisdall
#################################################################################
sub get_file_data {
    my($filename) = @_;
    # Initialize variables
    my @filedata = (  );
    unless( open(GET_FILE_DATA, $filename) ) {
        print STDERR "Cannot open file \"$filename\"\n\n";
        exit;
    }
    @filedata = <GET_FILE_DATA>;
    close GET_FILE_DATA;
    return @filedata;
}
#################################################################################
