#!/usr/bin/perl -w

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

$TEST = 0;
$a_script_per_PF = 1;
if($TEST){$test_number = 20;}
if($SGE_CLUSTER){ # nb of CPUs left available on the SGE cluster
	$CPUs_left_available = 3;
	$CPUs_left_thresh = .25; # proportion of available CPUs (left to other users)
}


##################################################################################
# WARNING: this script REQUIRES YN00 & codeml (from PAML) to be installed on     #
# your computer. One file has to sit in $curdir: Pandit6.2.txt (the DB file)    #
# for a given PF.                                                                #
#                                                                                #
# Oct. 2002 - Bioinformatics Research Center - North Carolina State University   #
#                                                                                #
# Comments to sarisbro at uottawa dot ca - Thanks.                                #
##################################################################################




#################################################################################
use POSIX;
# use Bio::DB::GenBank;
# use Bio::SeqIO;
# my $gb = new Bio::DB::GenBank;

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

# $start_at_PF =
# $end_at_PF =

########################### MAIN ###############################################
# generates file system from DB file ($infile)
Create_TreeFS();

# Initiates things
if($SGE_CLUSTER) {
	print "\nGenerating scripts for SGE cluster...";
	GenerateScriptsForMosix();
 	print "\nSpawning scripts on SGE cluster...\n";
	spawn_scripts_on_SGE();
}else{
	mkdir("$curdir/res",0777) || warn "\ncannot create results dir: $!";
	open(YN00FILE,">$curdir/res/YN00PF.dat") || warn "\ncould not open YN00PF file... $!\n";
	print YN00FILE "PF_name   seq  seq      S       N        t   kappa   omega     dN    SEdN    dS    SEdS\n";
	open(LRTFILE,">$curdir/res/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\n";
	open(PARAMFILE,">$curdir/res/paramSITES.dat") || warn "\ncould not open params file... $!\n";
	print PARAMFILE "PF_name\tM0-w\tM1-p0\tM1-p1\tM1-w0\tM1-w1\tM2-p0\tM2-p1\tM2-p2\tM2-w0\tM2-w1\tM2-w2\tM3-p0\tM3-p1\tM3-p2\tM3-w0\tM3-w1\tM3-w2\n";
	open(PAIRm0FILE,">$curdir/res/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/paramPairM1SITES.dat") || warn "\ncould not open params pairs file... $!\n";
	print PAIRm1FILE "PF_name\tSP_pair\tlnL-M1\tt-M1\tS-M1\tN-M1\tomega-M1\tdN-M1\tdS-M1\n";

	# dN/dS ratios estimation procedures
	Set_running_PF();
	Yang_Nielsen2000();
	Site_Models();
	Branch_Models;
	# Branch_Site_Models()
}

# dN/dS ratios estimation procedures
# do NOT uncomment the next lines if commented out!
# Yang_Nielsen2000();
# Site_Models();
# Branch_Models;
# Branch_Site_Models()

# The next lines are just to make sure jobs stop on SGE when they're done;
# make use of check_4_core()!
if($SGE_CLUSTER && $a_script_per_PF) {
	system("touch core.fake");
}

# summarize the outfiles
summarize();

close YN00FILE;
close DATAFILE;
close LRTFILE;
close PARAMFILE;
close PAIRm0FILE;
close PAIRm1FILE;

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

#################################################################################
#################################################################################
#################################################################################
#################################################################################
sub spawn_scripts_on_SGE{
	if($a_script_per_PF){
		$nb_of_scripts = $total_pf_number;
		$size_of_lumps = 1;
	}else{
		$nb_of_scripts = floor( $total_pf_number / $size_of_lumps);
	}

	mkdir("$curdir/outdump",0777) || warn "\nCould not create outdump dir: $!";
	chdir "$curdir/outdump" || warn "can't cd to $curdir/outdump: $!";
	for($i=0;$i<$nb_of_scripts;$i++){
		$start = 1 + $i * $size_of_lumps;
		$end = ($i + 1) * $size_of_lumps;
		LOOP_BACK:$spawn_new_script = Control_cpu_state();
		if($spawn_new_script){
			chmod(0700,"$curdir/scripts/PF$start-$end.pl");
			system("qsub -cwd -S /usr/bin/perl -j y $curdir/scripts/PF$start-$end.pl");
			sleep(30);
		}
		else{
			sleep(300);
			check_4_core();
			goto LOOP_BACK;
		}
	}
}
#################################################################################
sub check_4_core{
	system("find $curdir -type f -name \"core*\" > $curdir/kore_dump.out");
	@COREFILE = get_file_data("$curdir/kore_dump.out");
	foreach $line (@COREFILE) {
		($corepath,$idnum) = split(/\/core/,$line);
		($dum,$pfid) = split(/PF/,$corepath);
		# get script id on SGE & kill corresponding SGE job
		system("qstat -u stephane > $curdir/qstat.out");
		@PFIDFILE = get_file_data("$curdir/qstat.out");
		foreach $pfid_line (@PFIDFILE) {
			($job_id,$dumdum) = split(/     0 PF/,$pfid_line);
			($start_end,$dumdumdum) = split(/.pl stephane/,$dumdum);
			($mystart,$myend) = split(/-/,$start_end);
			if(($pfid>=$mystart)&&($pfid<=$myend)){
				# print "\nJob $job_id deleted for core at $pfid";
				system("qdel $job_id");
			}
		}
		print "$line";
		system("rm $line"); # erase core file
		unlink("$curdir/qstat.out");
	}
	unlink("$curdir/kore_dump.out");
}
#################################################################################
sub Control_cpu_state{
	$nb_of_free_cpu = 0;
	$nb_of_running_cpus = 0;
	$current_cpu_load = 0;
	system("qstat -f > $curdir/qstat.out");
	@QSTATFILE = get_file_data("$curdir/qstat.out");
	foreach $line (@QSTATFILE) {
		if(($line =~ /^node/) || ($line =~ /^Node/)){
			($left_of_BIP,$right_of_BIP) = split(/BIP/,$line);
			($left_of_slash,$right_of_slash) = split(/\//,$right_of_BIP);
			$left_of_slash =~ s/\s*//g;
			$nb_of_free_cpu += (2 - $left_of_slash);
			$current_cpu_load += $left_of_slash;
			# check how many nodes are not running properly (!!!)
			($l_of_glinux,$r_of_glinux) = split(/glinux/,$right_of_slash);
			$r_of_glinux =~ s/\s*//;
			if(!$r_of_glinux){
				$nb_of_running_cpus += 2;
			}
		}
	}
	# print "=>$nb_of_free_cpu ($nb_of_running_cpus) Load = $current_cpu_load";
	unlink("$curdir/qstat.out");

	# this ensures that $CPUs_left_available CPUs are left available on the cluster
	if($nb_of_running_cpus>3){
		if(($nb_of_running_cpus - $current_cpu_load) > $CPUs_left_available){
			$can_spawn_new_script = 1;
		}else{
			$can_spawn_new_script = 0;
		}
	}else{
		# if(($nb_of_running_cpus - $current_cpu_load) > $CPUs_left_available){
		if( (1.-($current_cpu_load / $nb_of_running_cpus)) > $CPUs_left_thresh){
			$can_spawn_new_script = 1;
		}else{
			$can_spawn_new_script = 0;
		}
	}
	return $can_spawn_new_script;
}
#################################################################################
sub kill_SGE_job{
}
#################################################################################
sub Yang_Nielsen2000{

	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");

		# generates control file for YN00
		# $icode = get_icode($pfid);
		$icode = ' 0 1 2 3 4 5 6 7 8 9 10';
		open(CTLFILE,">$curdir/data/PF$count/yn00.ctl") || warn "\ncannot create YN00 ctl file for PF$count: $!";
		print CTLFILE "      seqfile = $curdir/data/PF$count/dna-aln.nuc\n      outfile = $curdir/data/PF$count/yn.out\n      verbose = 0\n        icode = $icode  \n
    weighting = 1  * weighting pathways between codons (0/1)?
   commonf3x4 = 1  * use one set of codon freqs for all pairs (0/1)?\n* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,\n* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt.,\n* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt.,\n* 10: blepharisma nu. These codes correspond to transl_table 1 to 11 of GENEBANK.";
		close CTLFILE;

		# runs YN00 and parse result files
		if($SGE_CLUSTER){
			system("$curdir/yn00 $curdir/data/PF$count/yn00.ctl");
		}else{
			system("yn00 $curdir/data/PF$count/yn00.ctl");
		}
		rename("$curdir/data/PF$count/rst","$curdir/data/PF$count/YN00.rst") || warn "\nCannot rename rst file for YN00:$!";
		unlink("$curdir/data/PF$count/rub","$curdir/data/PF$count/rst1");
		Parse_YN00_OUT("$curdir/data/PF$count/yn.out","$pfid");
		PrintTimeStamp("TimeYN00","YN00 analysis for",$pfid);
	}
}
#################################################################################
sub Site_Models{

	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");

		# generates control files for codeml & run analyses
		for($id=0;$id<4;$id++) {
			# site models - global
			# $icode = get_icode($pfid);
			$icode = ' 0 1 2 3 4 5 6 7 8 9 10';
			open(CTLFILE,">$curdir/data/PF$count/codeml.M$id.ctl") || warn "\ncannot create ctl M$id file for $pfid: $!";
			print CTLFILE "      seqfile = $curdir/data/PF$count/dna-aln.nuc \n     treefile = dnatree.txt \n      outfile = $curdir/data/PF$count/PF$pfid-M$id.out \n        noisy = 0 \n      verbose = 0 \n      runmode = 0 \n      seqtype = 1 \n    CodonFreq = 2
        model = 0  * models for codons:
                   * 0:one, 1:b, 2:2 or more dN/dS ratios for branches
      NSsites = $id  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
                   * 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
                   * 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                   * 13:3normal>0
    fix_kappa = 0 \n        kappa = 2 \n    fix_omega = 0 \n        omega = 2 \n    fix_alpha = 1 \n        alpha = 0. \n       Malpha = 0 \n        ncatG = 3 \n        clock = 0 \n        getSE = 0 \n RateAncestor = 0 \n   Small_Diff = .5e-6 \n     cleandata = 1\n     icode = $icode";
			close CTLFILE;
			if($SGE_CLUSTER){
				system("$curdir/codeml $curdir/data/PF$count/codeml.M$id.ctl");
			}else{
				system("codeml $curdir/data/PF$count/codeml.M$id.ctl");
			}
			rename("$curdir/data/PF$count/rst","$curdir/data/PF$count/M$id.rst") || warn "\nCannot rename rst file for M$id:$!";
			rename("$curdir/data/PF$count/rst1","$curdir/data/PF$count/M$id.rst1") || warn "\nCannot rename rst1 file for M$id:$!";
			rename("$curdir/data/PF$count/2NG.dN","$curdir/data/PF$count/M$id.2NG.dN") || warn "\nCannot rename 2NG.dN file for M$id:$!";
			rename("$curdir/data/PF$count/2NG.dS","$curdir/data/PF$count/M$id.2NG.dS") || warn "\nCannot rename 2NG.dS file for M$id:$!";
			rename("$curdir/data/PF$count/2NG.t","$curdir/data/PF$count/M$id.2NG.t") || warn "\nCannot rename 2NG.t file for M$id:$!";
			unlink("$curdir/data/PF$count/rub");

			# site models - pairwise for M0 and M1
			if(($id==0) || ($id==1)){
				open(CTLFILE,">$curdir/data/PF$count/codeml.pM$id.ctl") || warn "\ncannot create ctl pM$id file for $pfid: $!";
				if($id){ # M1
					print CTLFILE "      seqfile = $curdir/data/PF$count/dna-aln.nuc \n     treefile = dnatree.txt \n      outfile = $curdir/data/PF$count/PF$pfid-pM$id.out \n        noisy = 0 \n      verbose = 0 \n      runmode = -2 \n      seqtype = 1 \n    CodonFreq = 2
        model = 0  * models for codons:
                   * 0:one, 1:b, 2:2 or more dN/dS ratios for branches
      NSsites = 0  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
                   * 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
                   * 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                   * 13:3normal>0
    fix_kappa = 0 \n        kappa = 2 \n    fix_omega = 1 \n        omega = 1 \n    fix_alpha = 1 \n        alpha = 0. \n       Malpha = 0 \n        ncatG = 3 \n        clock = 0 \n        getSE = 0 \n RateAncestor = 0 \n   Small_Diff = .5e-6 \n     cleandata = 1\n     icode = $icode";
				}else{ # M0
					print CTLFILE "      seqfile = $curdir/data/PF$count/dna-aln.nuc \n     treefile = dnatree.txt \n      outfile = $curdir/data/PF$count/PF$pfid-pM$id.out \n        noisy = 0 \n      verbose = 0 \n      runmode = -2 \n      seqtype = 1 \n    CodonFreq = 2
        model = 0  * models for codons:
                   * 0:one, 1:b, 2:2 or more dN/dS ratios for branches
      NSsites = 0  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
                   * 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
                   * 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                   * 13:3normal>0
    fix_kappa = 0 \n        kappa = 2 \n    fix_omega = 0 \n        omega = 2 \n    fix_alpha = 1 \n        alpha = 0. \n       Malpha = 0 \n        ncatG = 3 \n        clock = 0 \n        getSE = 0 \n RateAncestor = 0 \n   Small_Diff = .5e-6 \n     cleandata = 1\n     icode = $icode";
				}
				close CTLFILE;
				if($SGE_CLUSTER){
					system("$curdir/codeml $curdir/data/PF$count/codeml.pM$id.ctl");
				}else{
					system("codeml $curdir/data/PF$count/codeml.pM$id.ctl");
				}
				unlink("$curdir/data/PF$count/rub");
			}
		}

		# print result summary files
		print LRTFILE "$pfid\t";
		print PARAMFILE "$pfid\t";
		for($id=0;$id<4;$id++) {
			# Global analyses
			($lnL,$np) = get_Lnl_from_outfile("$curdir/data/PF$count/PF$pfid-M$id.out");
			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");
			}
			print PARAMFILE "$parms\t";
			# 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);
			}
		}
		print PARAMFILE "\n";
		print LRTFILE "\n";
		PrintTimeStamp("TimeSiteModels","Site models for",$pfid);
	}
}
#################################################################################
sub Branch_Models{

	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");

		# Branch-specific model (1 rate per branch)
		# AS this is quite intensive, only done for small data sets
		$do_it = KeepDataSet("$curdir/data/PF$count/dna-aln.nuc");
		if($do_it){
			# $icode = get_icode($pfid);
			$icode = ' 0 1 2 3 4 5 6 7 8 9 10';
			open(CTLFILE,">$curdir/data/PF$count/codeml.Branch.ctl") || warn "\ncannot create Branch model ctl file for $pfid: $!";
			print CTLFILE "      seqfile = dna-aln.nuc \n     treefile = dnatree.txt \n      outfile = PF$pfid-Branch.out \n        noisy = 0 \n      verbose = 0 \n      runmode = 0 \n      seqtype = 1 \n    CodonFreq = 2
        model = 1  * models for codons:
                   * 0:one, 1:b, 2:2 or more dN/dS ratios for branches
      NSsites = 0  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
                   * 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
                   * 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                   * 13:3normal>0
    fix_kappa = 0 \n        kappa = 2 \n    fix_omega = 0 \n        omega = 2 \n    fix_alpha = 1 \n        alpha = 0. \n       Malpha = 0 \n        ncatG = 3 \n        clock = 0 \n        getSE = 0 \n RateAncestor = 0 \n   Small_Diff = .5e-6 \n     cleandata = 1\n   icode = $icode";
			close CTLFILE;
			if($SGE_CLUSTER){
				system("$curdir/codeml $curdir/data/PF$count/codeml.Branch.ctl");
			}else{
				system("codeml $curdir/data/PF$count/codeml.Branch.ctl");
			}
			($lnL,$np) = get_Lnl_from_outfile("$curdir/data/PF$count/PF$pfid-Branch.out");
			print LRTFILE "$lnL\t$np\n";
			rename("$curdir/data/PF$count/rst","$curdir/data/PF$count/Branch.rst") || warn "\nCannot rename rst file for Branch model:$!";
			rename("$curdir/data/PF$count/rst1","$curdir/data/PF$count/Branch.rst1") || warn "\nCannot rename rst1 file for Branch model:$!";
			rename("$curdir/data/PF$count/2NG.dN","$curdir/data/PF$count/Branch.2NG.dN") || warn "\nCannot rename 2NG.dN file for Branch:$!";
			rename("$curdir/data/PF$count/2NG.dS","$curdir/data/PF$count/Branch.2NG.dS") || warn "\nCannot rename 2NG.dS file for Branch:$!";
			rename("$curdir/data/PF$count/2NG.t","$curdir/data/PF$count/Branch.2NG.t") || warn "\nCannot rename 2NG.t file for Branch:$!";
			unlink("$curdir/data/PF$count/rub");
		}
	}
	$pfid = Get_PF_id("$curdir/data/PF$count/TimeStamp.txt");
	PrintTimeStamp("TimeBranchModels","Branch models for",$pfid);

}
#################################################################################
sub Branch_Site_Models{

	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 Create_TreeFS {
	$pfid = '';
	$in_PFrecord1 = 0;
	$in_PFrecord2 = 0;
	open(DATAFILE,"$curdir/$infile") || warn "\ncould not open $infile... $!\n";
	mkdir("$curdir/data",0777) || warn "\ncannot create data dir: $!";
	print "\n\n";

	while(<DATAFILE>){
		chomp;
		($beginning,$rest) = split(/  /);
		if($beginning =~ /FAM/) {      # protein family name
			$count++;
			$total_pf_number = $count;
			$rest =~ s/\s//g;
			$pfid = $rest;
			if($count%2==0) {
				print " \rCounting families:   $count";
			}
			mkdir("$curdir/data/PF$count",0777);
			chdir "$curdir/data/PF$count" || warn "can't cd to $curdir/data/PF$count: $!";
			open(PFNAME,">$curdir/data/PF$count/$pfid.txt") || warn "\ncannot create $pfid: $!";
			close PFNAME;

		}
		elsif($beginning =~ /DES/) {      # protein family description
			open(PFNAME,">$curdir/data/PF$count/$pfid.txt");
			print PFNAME "$rest";
			close PFNAME;
		}
		elsif($beginning =~ /DNO/) { # extracts number of species
			$nb_species = $rest;
		}
		elsif($beginning =~ /DLN/) { # extracts DNA sequence length
			$seq_len = $rest;
		}
		elsif($beginning =~ /DPH/) { # tree with branch lengths
			open(TREEFILE,">$curdir/data/PF$count/dna-brlens.tre") || warn "\ncannot create dnatree brlens for $pfid: $!";
			print TREEFILE "$rest";
			close TREEFILE;
		}
		elsif($beginning =~ /DTP/) { # extracts DNA tree w/o branch lengths
			$dnatree = $rest;
			# create tree file
			open(TREEFILE,">$curdir/data/PF$count/dnatree.txt") || warn "\ncannot create dnatree for $pfid: $!";
			print TREEFILE "$nb_species  1\n\n$dnatree";
			close TREEFILE;
		}
		elsif($beginning =~ /SPT/) { # extracts accession numbers
			($name,$acc)= split(/:/,$rest);
			$acc =~ s/\s//g;
			if(!$in_PFrecord1) {
				open(SPNAME,">>$curdir/data/PF$count/spnames.txt") || warn "\ncannot create name for $pfid: $!";
				open(ACCFILE,">>$curdir/data/PF$count/acc.txt") || warn "\ncannot create acc for $pfid: $!";
				open(GBFILE,">$curdir/data/PF$count/data.GB") || warn "\ncannot create GB file for $pfid: $!";
				print SPNAME "$pfid: \n$name";
				print ACCFILE "$pfid: \n$acc";
				$in_PFrecord1 = 1;
				# $seqout = new Bio::SeqIO(-fh => \*GBFILE, -format => 'genbank');
				# $seqout = new Bio::SeqIO(-fh => \*STDOUT, -format => 'genbank');
				@acc_numb = qw( );
				@acc_numb = (@acc_numb, "$acc");
			}
			else {
				print SPNAME ",$name";
				print ACCFILE " $acc";
				# print "\r$pfid -> $acc";
				@acc_numb = (@acc_numb,  "$acc");
			}
			# $seq = $gb->get_Seq_by_acc('P28222');
			# $seqout->write_seq($seq);
		}
		elsif($beginning =~ /NAM/) { # extracts sp name of DNA seq file
			# this terminates recording sp. names & acc nb.
			$in_PFrecord1 = 0;
			close SPNAME;
			close ACCFILE;
			close GBFILE;
			# Extract_Taxonomy("$curdir/data/PF$count");
			# starts working on sp names for nuc file
			$dna_seq = '';
			if(!$in_PFrecord2) {
				open(DNAFILE,">>$curdir/data/PF$count/dna-aln.nuc") || warn "\ncannot create DNA data file for $pfid: $!";
				print DNAFILE "$nb_species  $seq_len\n\n";
				$in_PFrecord2 = 1;
			}
			$dna_sp_name = $rest;
		}
		elsif($beginning =~ /DSQ/) { # extracts corresponding DNA seq file
				$dna_seq = $rest;
		}
		elsif($beginning =~ /TRN/) { # this is to avoid entries with no DNA seq
			if($dna_seq) {
				print DNAFILE "$dna_sp_name     $dna_seq\n";
				$dna_seq = '';
			}
		}
		elsif($beginning =~ /\W/) { # marks the end of the entry for a given PF
			# this terminates recording data file for the PF just read
			$in_PFrecord2 = 0;
			close DNAFILE;

			# end of this PF record treatment
			$pfid =~ s/PF//;
			PrintTimeStamp("TimeStamp","Parsing flat file for",$pfid);
		}
		elsif($TEST && $count>$test_number) {
			$total_pf_number--;
			print "\nTotal = $total_pf_number families in $infile.\n";
			unlink <*.txt>;
			# rmdir("$curdir/data/PF$count") || warn "\nCannot rmdir last dir (PF$count):$!";
			print "\a\nTerminates reading flatfile here for test!\n";
			last;
		}
	}
}
#################################################################################
sub Summarize_PairwiseML{ # site models under M0 and M1
	$pfid = shift @_;
	$paiwiseMLoufile = shift @_;
	$id = shift @_;
	$start_reading_here = 0;

	@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;
			$filehandle = "PAIRm".$id."FILE";
			print $filehandle "$pfid  $pair  $lnL  $estimates";
		}
	}
}
#################################################################################
sub get_icode{
	$pfname = shift @_;
	$this_ic = "0"; # default value: is that reasonnable????
	$pfname =~ s/0*//;
	open(ICODE,"$curdir/icode.txt") || die "\ncould not open ICODE file... $!\n";

	while(<ICODE>){
		chomp;
		($pfcode,$ic) = split(/ icode = /);
		$pfcode =~ s/PF//;
		$pfcode =~ s/0*//;
		if($pfname == $pfcode){
			$this_ic = $ic;
		}
	}
	close ICODE;
	return $this_ic;
}
#################################################################################
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 @_;

	@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;
		}
	}
	return $id;
}
#################################################################################
sub Parse_YN00_OUT{
	$infile = shift @_;
	$pfid = shift @_;
	$start_next_line = 0;

	@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 "+-"
			print YN00FILE "$pfid    ".$line if $line;
		}
	}
}
#################################################################################
sub get_dNdS_rate{ # under M0
	$infile = shift @_;
	$after_dNdS_line = 0;

	@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;
		}
	}
	return $omega;
}
#################################################################################
sub get_params_from_outfile { # under M1 to M3
	$infile = shift @_;
	$after_dNdS_line = 0;
	$params = '';

	@INFILE = get_file_data($infile);
	foreach $line (@INFILE) {
		if($after_dNdS_line){
			if($line =~ /^p:/){
				$line =~ s/^p://;
				$line =~ s/\n/ /;
				$params .= $line;
			}
			elsif($line =~ /^w:/){
				$line =~ s/^w://;
				$line =~ s/\n/ /;
				$params .= $line;
			}
		}
		elsif($line =~ /^dN\/dS/) {
			$after_dNdS_line = 1;
		}
	}
	return $params;
}
#################################################################################
sub get_Lnl_from_outfile {
	$infile = shift @_;

	@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);
		}
	}
	return ($lnl,$np);
}
#################################################################################
sub KeepDataSet{
	$data_file = shift @_;
	
	$max_ns = 20;
	$max_nb_codons = 150;

	open(SEQFILE,"$data_file") || warn "\ncould not open nuc file... $!\n";
	chomp($line = <SEQFILE>);
	($ns,$ls)=split(/\s+/, $line);
	if($ns<$max_ns || $ls<$max_nb_codons*3) {
		$okay = 1;
	}else{
		$okay = 0;
	}
	close NUCSEQ;
	return $okay;
}
#################################################################################
sub Extract_Taxonomy { # ...from GB file
	$path2taxo_file = shift @_;
	$infile = "$path2taxo_file/data.GB";
	open(TAXOFILE,">$path2taxo_file/taxonomy.txt") || warn "\ncannot create taxo file for $pfid: $!";
	@GB_FILE = get_file_data($infile);
	foreach $line (@GB_FILE) {
		if($line =~ /^DEFINITION/) {
			$line =~ s/^DEFINITION\s*//;
			$line =~ s/ /_/;
			$line =~ s/\s//g;
			$definition = $line;
		}
		if($line =~ /^ACCESSION/) {
			$line =~ s/^ACCESSION\s*//;
			$line =~ s/ /_/;
			$line =~ s/\s//g;
			$accession = $line;
			$inside_record = 1;
		}
		elsif(($line =~ /^            /)&&($inside_record)) {
			$line =~ s/^            \s*//;
			$line =~ s/\s//g;
			$bino_sp_name .= $line;
		}
		elsif($line =~ /^REFERENCE/) {
			$inside_record = 0;
			$bino_sp_name .= ")";
		}
		elsif($line =~ /^  ORGANISM/) {
			$line =~ s/^  ORGANISM\s*//;
			$line =~ s/ /_/;
			$line =~ s/\s//g;
			$bino_sp_name = $line." (";
		}
		elsif($beginning =~ /\W/) { # end of a given record in the GB file: "//"
			print TAXOFILE "$accession: bino_sp_name {$definition}";
		}
	} # foreach
	close TAXOFILE;
}
#################################################################################
sub Check_ICODE{
	for($i=1;$i<$total_pf_number;$i++){
		chdir "$curdir/data/PF$i" || warn "can't cd to $curdir/data/PF$i: $!";
		$pfid = Get_PF_id("$curdir/data/PF$i/spnames.txt");
		$icode = get_icode($pfid);
		print "\nPF$i -> PF$pfid => $icode";
	}
}
#################################################################################
sub summarize {
}
#################################################################################
# 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;
}
#################################################################################
# CUT HERE ######################################################################
sub GenerateScriptsForMosix{

	# General to any plateform
	mkdir("$curdir/res",0777) || warn "\ncannot create results dir: $!";
	open(YN00FILE,">$curdir/res/YN00PF.dat") || warn "\ncould not open YN00PF file... $!\n";
	print YN00FILE "PF_name   seq  seq      S       N        t   kappa   omega     dN    SEdN    dS    SEdS\n";
	open(LRTFILE,">$curdir/res/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\n";
	open(PARAMFILE,">$curdir/res/paramSITES.dat") || warn "\ncould not open params file... $!\n";
	print PARAMFILE "PF_name\tM0-w\tM1-p0\tM1-p1\tM1-w0\tM1-w1\tM2-p0\tM2-p1\tM2-p2\tM2-w0\tM2-w1\tM2-w2\tM3-p0\tM3-p1\tM3-p2\tM3-w0\tM3-w1\tM3-w2\n";
	open(PAIRm0FILE,">$curdir/res/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/paramPairM1SITES.dat") || warn "\ncould not open params pairs file... $!\n";
	print PAIRm1FILE "PF_name\tSP_pair\tlnL-M1\tt-M1\tS-M1\tN-M1\tomega-M1\tdN-M1\tdS-M1\n";

	# Specific to SGE plateform
	mkdir("$curdir/scripts",0777) || warn "\ncannot create scripts dir: $!";
	@INFILE = get_file_data("$curdir/pundit-v2.pl");
	if($a_script_per_PF){
		$nb_of_scripts = $total_pf_number;
		$size_of_lumps = 1;
	}else{
		$nb_of_scripts = floor( $total_pf_number / $size_of_lumps);
	}
	for($i=0;$i<$nb_of_scripts;$i++){
		$start = 1 + $i * $size_of_lumps;
		$end = ($i + 1) * $size_of_lumps;
		# create perl scripts
		open(REPLICATEE,">$curdir/scripts/PF$start-$end.pl") || warn "\ncould not open PF$start-$end.pl... $!\n";
		foreach $line(@INFILE){
			if($line =~ /CUT HERE/){
				print "\n  --> Script PF$start-$end.pl printed to file!";
				chmod(0700,"$curdir/scripts/PF$start-$end.pl");
				close REPLICATEE;
				last;
			}
			elsif($line =~ /start_at_PF =/){# set start
				$line = "\$start_at_PF = ";
				$line .= " $start;\n";
			}
			elsif($line =~ /end_at_PF =/){# set end
				$line = "\$end_at_PF = ";
				$line .= " $end;\n";
			}
			elsif($line =~ /^Create_TreeFS/){# comment out Create_TreeFS function in the main
				$line = "# Create_TreeFS();\n";
			}
			elsif($line =~ /^	GenerateScriptsForMosix/){# comment out GenerateScriptsForMosix function in the main
				$line = "# GenerateScriptsForMosix();\n";
			}
			elsif($line =~ /^	spawn_scripts_on_SGE/){# comment out spawn_scripts_on_SGE function in the main
				$line = "# spawn_scripts_on_SGE();\n";
			}
			elsif($line =~ /^\# Yang_Nielsen2000/){# uncomment out YN00 analyses!
				$line = "Yang_Nielsen2000();\n";
			}
			elsif($line =~ /^\# Site_Models/){# uncomment out site model analyses!
				$line = "Site_Models();\n";
			}
			elsif($line =~ /^\# Branch_Models/){# uncomment out branch model analyses!
				$line = "Branch_Models();\n";
			}
			print REPLICATEE $line;
		}
	}
}
#################################################################################
sub Set_running_PF{
	$start_at_PF = 1;
	$end_at_PF = $total_pf_number;
}
#################################################################################
