diff PGAP-1.2.1/PGAP.pl @ 0:83e62a1aeeeb draft

Uploaded
author dereeper
date Thu, 24 Jun 2021 13:51:52 +0000 (2021-06-24)
parents
children 70b7a5270968
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/PGAP-1.2.1/PGAP.pl	Thu Jun 24 13:51:52 2021 +0000
@@ -0,0 +1,2825 @@
+#!/usr/bin/perl
+use strict;
+use warnings;
+use Getopt::Long;
+
+### programs from BLAST
+my $formatdb="formatdb";
+my $blastall="blastall";
+
+### programs from mcl
+my $mcl="mcl";
+
+### programs from mafft
+my $mafft="mafft";
+
+### programs from PHYLIP
+my $seqboot="seqboot";
+my $neighbor="neighbor";
+my $consense="consense";
+my $dnaml="dnaml";
+my $dnadist="dnadist";
+my $dnapars="dnapars";
+
+my $count_tree=0;
+
+my $sampleSize=8000; # when calculate the pan-genome size, we will sample $sampleSize combinations
+                     # if the total combination number is larger than $sampleSize for specific genomes
+                     # Surely, the number of $sampleSize is, the larger, the better. 
+                     # However, the larger the $sampleSize is, the more time would be consumed.
+                     # we suggest the range: 5000 ~ 20,000
+
+#####################################################################
+#  DOn't modify the following code, unless you know their functions
+#####################################################################
+
+my %opt=qw();
+GetOptions(\%opt,"strains:s","input:s","output:s","cluster!","pangenome!","variation!","evolution!","function!","method:s","thread:i","score:f","evalue:f","coverage:f","local:f","global:f","identity:f","bootstrap:i","help|h!");
+
+my @usage=qq(
+====== Pan-Genome Analysis Pipeline (PGAP) ======
+                   Version 1.2.1
+
+Usage:   perl PGAP.pl [Options]
+
+Options: 
+  --strains    String    Input strains nicknames, and join them with '+', for example: A+B+C
+  --input      String    Input data directory 
+  --output     String    Result output directory
+ 
+  --cluster              Run homologous gene clustering
+  --pangenome            Run pan-genome analysis
+  --variation            Run homologous clusters variation analysis
+  --evolution            Run evolution analysis
+  --function             Run Function analysis
+
+  --method     String    GF for GeneFamily method,  and MP for MultiParanoid method
+                           for GF: fast, but not very accurate
+                               evalue, score, indentity, coverage are employed
+                           for MP: slow, but more accurate
+                               score, coverage, local, global are employed
+  --thread     Int       Number of processors to use in blastall. [default:1]
+  --score      Int       Minimum score in blastall. [default:40]
+  --evalue     Decimal   Maximal E-value in blastall. [default:1e-10]
+  --coverage   Decimal   Minimum alignment coverage for two homologous proteins. [default:0.5]
+  --local      Decimal   Minimum local alignment overlap in MP method. [default:0.25]
+  --global     Decimal   Minimum global alignment overlap in MP method. [default:0.5]
+  --identity   Decimal   Minimum alignment indentity for two homologous proteins. [default:0.5]
+  --bootstrap  Int       Bootstrap times for phylogenetics tree. [default:1]
+
+  --h or help            Display this message
+);
+
+
+############# specified variable #############
+my $inputDIR;
+my $outputDIR;
+my $run_cluster;
+my $run_pangenome;
+my $run_variation;
+my $run_evolution;
+my $run_function;
+my $method="";
+my $thread;
+my $score;
+my $identity;
+my $evalue;
+my $coverage;
+my $global;
+my $local;
+my $bootstrap;
+
+
+my %pep;
+my %nuc;
+my $spnum;
+my @clusters;
+my $Cluster;
+my @SpecieCombination;
+my @spID;
+my %genenum;
+my %aaAln;
+my %ntAln;
+my %cog;
+my %description;
+#my %aa4tree;  ### AA sequence for Phylogenetic Tree
+my %nt4tree;  ### nucleotide sequence for Phylogenetic Tree
+my @SNPPosition; ### SNP position
+my $dieMessage="You did not run PGAP.pl in the program directory\n";
+my $section;
+
+######### common temporary variable #############
+my $i;
+my $j;
+my $line;
+my %tmpHash;
+my @tmp;
+my $tmp;
+my $key;
+my @row;
+my $inparacount;
+my $ClusterID;
+my $orth;
+my @content;
+my $clusterName;
+my @xdata;
+my @ydata;
+my @fit;
+
+my $fit_A;
+my $fit_A_interval;
+my $fit_B;
+my $fit_C;
+my $fit_C_interval;
+my $fit_Rsquare;
+
+
+
+#### check option
+
+my $opt_error=0;
+
+if ((scalar(keys %opt) ==0) or (exists($opt{"help"}))) 
+	{
+		print join("\n",@usage)."\n";
+		exit;
+	}
+
+
+###################### public info
+### strains name
+my @species;
+if (exists($opt{"strains"})) 
+	{
+		@species=split(/\+/,$opt{"strains"});
+		$spnum=scalar(@species);
+	}else
+		{
+			print "Please assign strains nick name!\n";
+			exit;
+		}
+
+### input data directory
+
+if (exists($opt{"input"})) 
+	{
+		$inputDIR=$opt{"input"};
+		if ($inputDIR!~/\/$/) 
+			{
+				$inputDIR=$inputDIR."/";
+			}
+	}else
+		{
+			print "Please assign input data directory!\n\n";
+			exit;
+		}
+### output data directory
+
+if (exists($opt{"output"})) 
+	{
+		$outputDIR=$opt{"output"};
+		if ($outputDIR!~/\/$/) 
+			{
+				$outputDIR=$outputDIR."/";
+			}
+	}else
+		{
+			print "Please assign result output directory!\n\n";
+			exit;
+		}
+
+###################### section info
+
+if (exists($opt{"cluster"})) 
+	{
+		$run_cluster=1;
+	}else
+		{
+			$run_cluster=0;
+		}
+
+if (exists($opt{"pangenome"})) 
+	{
+		$run_pangenome=1;
+	}else
+		{
+			$run_pangenome=0;
+		}
+
+if (exists($opt{"variation"})) 
+	{
+		$run_variation=1;
+	}else
+		{
+			$run_variation=0;
+		}
+
+if (exists($opt{"evolution"})) 
+	{
+		$run_evolution=1;
+	}else
+		{
+			$run_evolution=0;
+		}
+
+if (exists($opt{"function"})) 
+	{
+		$run_function=1;
+	}else
+		{
+			$run_function=0;
+		}
+
+if ($run_cluster) 
+	{
+		### method
+		if (exists($opt{"method"})) 
+			{
+				$method=uc($opt{"method"});
+				if ($method!~/^GF$/ and $method!~/^MP$/) 
+					{
+						print "Unknown method: ".$opt{"method"}."\n";
+						exit;
+					}
+			}else
+				{
+					print "Please assign the cluster method!\n\n";
+					exit;
+				}
+
+		##thread
+				if (exists($opt{"thread"})) 
+					{
+						$thread=$opt{"thread"};
+						if ($thread==0) 
+							{
+								print "please assign an applicable thread value.\n";
+								exit;
+							}
+					}else
+						{
+							$thread=1;
+						}
+
+		##score
+				if (exists($opt{"score"})) 
+					{
+						$score=$opt{"score"};
+						if ($score<=0) 
+							{
+								print "please assign an applicable score value.\n";
+								exit;
+							}
+					}else
+						{
+							$score=40;
+						}
+
+		if ($method eq "GF")
+			{
+		###identity
+				if (exists($opt{"identity"})) 
+					{
+						$identity=$opt{"identity"};
+						if ($identity>1 or $identity<=0) 
+							{
+								print "identity should be 0 ~ 1 \n";
+								exit;
+							}
+					}else
+						{
+							$identity=0.5;
+						}
+
+		###evalue
+				if (exists($opt{"evalue"})) 
+					{
+						$evalue=$opt{"evalue"};
+					}else
+						{
+							$evalue=1e-10;
+						}
+
+		###coverage
+				if (exists($opt{"coverage"})) 
+					{
+						$coverage=$opt{"coverage"};
+						if ($coverage>1 or $coverage<=0) 
+							{
+								print "coverage should be 0 ~ 1 \n";
+								exit;
+							}
+					}else
+						{
+							$coverage=0.5;
+						}
+			}
+
+
+		if ($method eq "MP")
+			{
+		###global
+				if (exists($opt{"global"})) 
+					{
+						$global=$opt{"global"};
+						if ($global>1 or $global<=0) 
+							{
+								print "global coverage should be 0 ~ 1 \n";
+								exit;
+							}
+					}else
+						{
+							$global=0.5;
+						}
+		###local
+				if (exists($opt{"local"})) 
+					{
+						$local=$opt{"local"};
+						if ($local<=0) 
+							{
+								print "local coverage should be 0 ~ [global coverage value] \n";
+								exit;
+							}
+						if ($local>$global) 
+							{
+								print "local coverage should be less than global coverage!\n";
+								exit;
+							}
+					}else
+						{
+							$local=0.25;
+						}
+			}
+	}
+
+if ($run_evolution) 
+	{
+		if (exists($opt{"bootstrap"})) 
+			{
+				$bootstrap=$opt{"bootstrap"};
+				if ($bootstrap<=0) 
+					{
+						print "please assign an applicable bootstrap value.\n";
+					}
+			}else
+				{
+					$bootstrap=1;
+				}
+	}
+
+print "Program begin at ".localtime()."\n";
+print "The following are the parameters for current process:\n";
+				print "Strains:                   ".join(",",@species)."\n";
+				print "Input directory:           $inputDIR\n";
+				print "Output directory:          $outputDIR\n";
+if ($run_cluster) 
+	{
+				print "Cluster analysis:          yes\n";
+				print "  Method:                  $method\n";
+				print "  Thread:                  $thread\n";
+		if ($method eq "GF") 
+			{
+				print "  E-value:                 $evalue\n";
+				print "  Identity:                $identity\n";
+				print "  Coverage:                $coverage\n";
+				print "  Score:                   $score\n";
+			}
+		if ($method eq "MP") 
+			{
+				print "  Local:                   $local\n";
+				print "  Global:                  $global\n";
+				print "  Score:                   $score\n";
+			}
+	}else
+		{
+				print "Cluster analysis:          no\n";
+		}
+if ($run_pangenome) 
+	{
+				print "Pan-genome analysis:       yes\n";
+	}else
+		{
+				print "Pan-genome analysis:       no\n";
+		}
+if ($run_variation) 
+	{
+				print "Variation analysis:        yes\n";
+	}else
+		{
+				print "Variation analysis:        no\n";
+		}
+if ($run_evolution) 
+	{
+				print "Evolution analysis:        yes\n";
+				print "  Bootstrap:               $bootstrap\n";
+	}else
+		{
+				print "Evolution analysis:        no\n";
+		}
+if ($run_function) 
+	{
+				print "Function analysis:         yes\n";
+	}else
+		{
+				print "Function analysis:         no\n";
+		}
+
+$section=$run_cluster.$run_pangenome.$run_variation.$run_evolution.$run_function;
+
+###############################################
+#	section 0) check input file and program
+###############################################
+if (!(-e $outputDIR)) 
+	{
+		system("mkdir $outputDIR");
+	}
+system("chmod +rw $outputDIR");
+
+if (!(-w $outputDIR)) 
+	{
+		print "There is no WRITE permission in $outputDIR\n";
+		exit;
+	}
+@tmp=qw();
+&CheckInputFile(\@species,$inputDIR,$section,$method,\@tmp);
+&CheckExtraProgram($section,$method,\@tmp);
+
+if (scalar(@tmp)>0) 
+	{
+		open(R,">".$outputDIR."0.error.message");
+		print R join("",@tmp)."\n";
+		close(R);
+		print "error!\nlog are saved in ${outputDIR}0.error.message\n";
+		exit;
+	}
+
+
+############################################
+#	section 1) cluster analysis
+############################################
+
+if ($run_cluster) 
+	{
+	print "\n\n############################################\n";
+	print "#	section 1) cluster analysis\n";
+	print "############################################\n\n\n";
+
+#### cluster gene and return result to the array @clusters
+
+		if ($method eq "MP") 
+			{
+				print "Begin cluster gene with MP method ...\n";
+				&MP();
+			}else
+				{
+					print "Begin cluster gene with GF method ...\n";
+					&GF();
+				}
+
+#### output normal cluster format
+
+	&FormatClusterOutPut(\@species,"${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
+
+#### Retrieve cluster
+
+	&RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
+
+##### gene distribution in each strains
+		%tmpHash=();
+		&GeneDistribution(\@clusters,\%tmpHash);
+
+		open(R,">${outputDIR}1.Gene_Distribution_By_Conservation.txt");
+		print R "SharedBy_Strains\t".join("\t",@species)."\n";
+
+		for ($i=$spnum;$i>0;$i--) 
+			{
+				print R $i;
+				for ($j=0;$j<$spnum;$j++) 
+					{
+						if (exists($tmpHash{$i."|".$j})) 
+							{
+								print R "\t".$tmpHash{$i."|".$j};
+							}else
+								{
+									print R "\t0";
+								}
+					}
+				print R "\n";
+			}
+
+		close(R);
+
+
+	}else
+		{
+			print "Homologous gene clustering is skipped!\n";
+		}
+
+
+if ($run_pangenome) 
+	{
+	print "\n\n############################################\n";
+	print "#	section 2) Pan-genome analysis\n";
+	print "############################################\n\n\n";
+
+#### Retrieve cluster
+	&RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
+	chomp(@clusters);
+
+#### convert file into 0-1 matrix
+	for ($line=0;$line<@clusters;$line++) 
+		{
+			@row=split(/\t/,$clusters[$line]);
+			splice(@row,0,1);
+			for ($i=0;$i<@row;$i++) 
+				{
+					if ($row[$i] eq "-") 
+						{
+							$row[$i]=0;
+						}else
+							{
+								$row[$i]=1;
+							}
+				}
+			$clusters[$line]=join("\t",@row);
+		}
+
+#### fetch gene number of each strains
+for ($i=0;$i<$spnum;$i++) 
+	{
+		open(F,"$inputDIR$species[$i].pep");
+		@tmp=<F>;
+		close(F);
+		@tmp=grep(/^>/,@tmp);
+		$genenum{$species[$i]}=scalar(@tmp);
+	}
+
+#### pan genome size and core genome size
+print "Deducing pan genome size and core genome size for each composition...\n\n";
+
+open(PAN,">${outputDIR}2.PanGenome.Data.txt");
+print PAN "ClusterConservation\tTotalGeneNumber\tPanGenome\tCoreGenome\n";
+
+for ($i=1;$i<=scalar(@species);$i++) 
+	{
+		#@SpecieCombination=&Combination(\@species,$i);
+		#@SpecieCombination=&Combination($spnum,$i);
+		if (&ChkCombinationValue($spnum,$i) !=0) ### transfer the array reference to the subroutine
+		{
+			&Combination($spnum,$i,\@SpecieCombination); ## if the combination number is less than sampleSize, then fecth all, else sample
+		}else
+		{
+			&SampleCombination($spnum,$i,\@SpecieCombination);
+		}
+		
+		foreach $key (@SpecieCombination) 
+			{
+		##### count total gene number in current combination
+					$tmp=0;
+					@spID=split(/\t/,$key); #### speices id  in current combination
+					foreach  (@spID) 
+						{
+							$tmp=$tmp+$genenum{$species[$_]};
+						}
+		##### scan pangenome and coregenome
+					@tmp=split(/\t/,&PanGenomeNumber(\@spID));
+					print PAN "$i\t$tmp\t".join("\t",@tmp)."\n";
+
+			}
+	}
+
+close(PAN);
+
+#### data fit
+
+#### for model A
+
+if ($spnum<3) 
+	{
+		print "There are $spnum strains. For pan-genome function fitting, at least 3 strains data are required.\n";
+	}else
+		{
+		open(R,">${outputDIR}2.PanGenome.Profile.txt");
+	##### genome number & pan-genome size
+			@xdata=qw();
+			@ydata=qw();
+			&ReadData2Array("${outputDIR}2.PanGenome.Data.txt",\@xdata,0,\@ydata,2);
+			&SumData(\@xdata,\@ydata,"mean");
+			($fit_Rsquare, $fit_A, $fit_A_interval, $fit_B, $fit_C, $fit_C_interval)=&fit_model_A(\@xdata,\@ydata);
+			print R "The relation bewteen genome number and pan-genome size\n\n";
+			print R "Function model: y=A*x**B +C \n";
+			print R "\ty denotes pan-genome size, x denotes genome number, and A, B, C are fitting parameters.\n\n";
+			print R "Fitting result:\n";
+			print R "\ty = $fit_A *x**$fit_B + $fit_C\n";
+			print R "\tR-square = $fit_Rsquare\n";
+			print R "\tA 95% confidence interval: ($fit_A - $fit_A_interval , $fit_A + $fit_A_interval)\n";
+			print R "\tC 95% confidence interval: ($fit_C - $fit_C_interval , $fit_C + $fit_C_interval)\n\n\n\n\n";
+
+	##### total gene number & pan-genome size
+			#@xdata=qw();
+			#@ydata=qw();
+			#&ReadData2Array("${outputDIR}2.PanGenome.Data.txt",\@xdata,1,\@ydata,2);
+			#&SumDataByMedian(\@xdata,\@ydata);
+			#($fit_Rsquare, $fit_A, $fit_B, $fit_C)=&fit_model_A(\@xdata,\@ydata);
+			#print R "The relation bewteen total gene number and pan-genome size\n\n";
+			#print R "$fit_Rsquare, $fit_A, $fit_B, $fit_C\n";
+			#print R "\ty = $fit_A *x**$fit_B + $fit_C   R-square = $fit_Rsquare\n";
+			#print R "\tx: total gene number\n";
+			#print R "\ty: pan-genome size\n\n\n\n\n";
+
+	##### genome number & core genome
+			@xdata=qw();
+			@ydata=qw();
+			&ReadData2Array("${outputDIR}2.PanGenome.Data.txt",\@xdata,0,\@ydata,3);
+			&SumData(\@xdata,\@ydata,"mean");
+			($fit_Rsquare, $fit_A, $fit_A_interval, $fit_B, $fit_C, $fit_C_interval)=&fit_model_B(\@xdata,\@ydata);
+			print R "The relation bewteen genome number and core genome size\n\n";
+			print R "Function model: y=A*exp(B*x) +C \n";
+			print R "\ty denotes pan-genome size, x denotes genome number, and A, B, C are fitting parameters.\n\n";
+			print R "Fitting result:\n";
+			print R "\ty = $fit_A *exp($fit_B * x) + $fit_C   R-square = $fit_Rsquare\n";
+			print R "\tR-square = $fit_Rsquare\n";
+			print R "\tA 95% confidence interval: ($fit_A - $fit_A_interval , $fit_A + $fit_A_interval)\n";
+			print R "\tC 95% confidence interval: ($fit_C - $fit_C_interval , $fit_C + $fit_C_interval)\n\n\n\n\n";
+		close(R);
+		}
+
+	}
+
+
+	############################################
+	#	section 3) CDS variation analysis
+	############################################
+
+if ($run_variation) 
+	{
+		print "\n\n############################################\n";
+		print "#	section 3) CDS variation analysis\n";
+		print "############################################\n\n\n";
+
+	#### Retrieve cluster
+		&RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
+		chomp(@clusters);
+
+	 ## protein
+		system("rm -rf *.pep");
+		&PrepareFasta(\@species,$inputDIR,".pep"); ###prepare pep file
+		system("cat *.pep > All.faa && rm -rf *.pep && mv  All.faa All.pep");
+		&ReadSequenceInToHash("All.pep",\%pep);
+	 ## nucleic
+		system("rm -rf *.nuc");
+		&PrepareFasta(\@species,$inputDIR,".nuc"); ###prepare nuc file
+		system("cat *.nuc > All.ffn && rm -rf *.nuc && mv  All.ffn All.nuc");
+		&ReadSequenceInToHash("All.nuc",\%nuc);
+
+	## scanning SNP
+		%nt4tree=();
+		for ($i=0;$i<$spnum;$i++) 
+			{
+				$nt4tree{"S".$i}="";
+			}
+
+		open(VAR,">${outputDIR}3.CDS.variation.txt");
+		print VAR "ClusterID\tStrains_Number\tGene_Number\tPosition\taaType\tntType\tntProfile\tVariation type\n";
+
+		open(VA,">${outputDIR}3.CDS.variation.analysis.txt");
+		print VA "ClusterID\tInDel Base\tNonsynonymous mutation\tSynonymous mutation\n";
+
+		for ($line=0;$line<@clusters;$line++) 
+			{
+				@row=split(/\t|\,/,$clusters[$line]);
+				$ClusterID=$row[0];
+				splice(@row,0,1);
+				@row=grep(/^S/,@row);
+				if (scalar(@row) >=2) 
+					{
+						open(PEP,">$ClusterID.pep");
+						open(NUC,">$ClusterID.nuc");
+						foreach $key (@row) 
+							{
+								print PEP ">$key\n$pep{$key}\n";
+								print NUC ">$key\n$nuc{$key}\n";
+							}
+						close(PEP);
+						close(NUC);
+						system("$mafft --quiet $ClusterID.pep > $ClusterID.pal");
+						#system("perl ./pal2nal.pl $ClusterID.pal $ClusterID.nuc -output fasta > $ClusterID.nal");
+						$tmp=&pal2nal("$ClusterID.pal","$ClusterID.nuc","$ClusterID.nal");
+						if ($tmp == 0) 
+							{
+								system("rm -rf $ClusterID.*");
+								next;
+							}
+
+						@tmp=&DetectSNP();
+						if (scalar(@tmp)>0) 
+							{
+								print VA $ClusterID."\t".&VarAnalysis(\@tmp)."\n";
+								print VAR join("",@tmp);
+						### core orthologs
+								@row=split(/\t/,$clusters[$line]);
+								splice(@row,0,1);
+								if ((&CountGeneInCluster(join("\t",@row)) ==$spnum) and (&CountSpeicesInCluster(join("\t",@row)) == $spnum) )
+									{
+										$count_tree++;
+										%tmpHash=();
+										foreach  (@row) 
+											{
+												$tmpHash{$_}="";
+											}
+										&RemoveHeadGap("$ClusterID.nal",\%tmpHash);
+										&ExtractSNP4tree(\%tmpHash,\%nt4tree);
+									}
+							}
+						system("rm -rf $ClusterID.*");
+					}
+			}
+		close(VAR);
+		close(VA);
+
+		open(R,">${outputDIR}3.CDS.variation.for.evolution.txt");
+		foreach $key (keys %nt4tree) 
+			{
+				$_=$key;
+				s/s//gi;
+				print R ">$species[$_]\n$nt4tree{$key}\n";
+			}
+		close(R);
+		print $count_tree."\n\n";
+
+###
+		system("rm All.nuc All.pep");
+	}else
+		{
+			print "CDS variation is skipped.\n";
+		}
+
+	############################################
+	#	section 4) CDS variation analysis
+	############################################
+
+if ($run_evolution) 
+	{
+#### Retrieve cluster
+		&RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
+		chomp(@clusters);
+
+
+
+##################
+##
+##  Distance based
+##
+##################
+
+### caculate the distance between each two strains by cluster
+		&ClusterProfil4Specie(\%tmpHash);   # caculate Clusters profile for each specie
+		&DistanceMatrix($spnum,\%tmpHash);   # caculate distance matrix acoording to Clusters profile
+	###output distance
+		open(DIST,">${outputDIR}4.Species_Distance_Clusters_Based.txt");
+		###header
+		printf DIST "%5d", $spnum;
+		print DIST "\n";
+		foreach $i (0..($spnum-1)) 
+			{
+				$key="sp".$i."sp";
+				printf DIST "%-10s",$key;
+				foreach $j (0..($spnum-1)) 
+					{
+						printf DIST "  %8f",$tmpHash{$i."-".$j};
+					}
+				print DIST "\n";
+			}
+		close(DIST);
+
+		%tmpHash=(); 
+
+		### based on pan genome (distance)
+		print "\nDraw pangenome based phylogenetic tree ...\n\n";
+
+		&PanBasedTree("${outputDIR}4.Species_Distance_Clusters_Based.txt","${outputDIR}4.PanBased");
+
+##################
+##
+##  SNP based
+##
+##################
+		%tmpHash=();
+		if (!(-e "${outputDIR}3.CDS.variation.for.evolution.txt")) 
+			{
+				print "Variation in core orthologs cluster is not found from ${outputDIR}3.CDS.variation.for.evolution.txt.\n";
+				print "Maybe you have skipped CDS variation analysis.\n";
+			}else
+				{
+					&ReadSequenceInToHash("${outputDIR}3.CDS.variation.for.evolution.txt",\%tmpHash);
+					open(R,">mlst.aln");
+					for ($i=0;$i<@species;$i++) 
+						{
+							print R ">sp${i}sp\n".$tmpHash{$species[$i]}."\n";
+						}
+					close(R);
+					&fasta2phylip("mlst.aln","mlst.phylip");
+					system("rm -rf mlst.aln");
+
+					print "\nDraw SNP based phylogenetic tree ...\n\n";
+					&SNPBasedTree("mlst.phylip","${outputDIR}4.SNPBased");
+					system("rm -rf mlst.phylip")
+				}
+
+########
+# replace speices name
+########
+
+			opendir(DIR,"${outputDIR}");
+			@tmp=readdir(DIR);
+			closedir(DIR);
+			@tmp=grep(/^4/,@tmp);
+			foreach $tmp (@tmp) 
+				{
+					&ReplaceName(\@species,"${outputDIR}$tmp");
+				}
+	}else
+		{
+			print "Evolution analysis is skipped.\n";
+		}
+
+
+	############################################
+	#	section 4) Function analysis
+	############################################
+
+if ($run_function) 
+	{
+#### Retrieve cluster
+		&RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
+		chomp(@clusters);
+
+#### prepare annotation file
+		&PrepareTable(\@species,$inputDIR,".function"); ###prepare location file
+		&ReadAnnotation(\@species,\%cog,\%description);
+
+
+#### assign function 
+		open(R,">${outputDIR}5.Orthologs_Cluster_Function.txt");
+		print R "ClusterID\tConservation_Level\tCOG\tDescription\n";
+		for ($i=0;$i<@clusters;$i++) 
+			{
+				@row=split(/\t/,$clusters[$i]);
+				$ClusterID=$row[0];
+				splice(@row,0,1);
+				print R $ClusterID."\t".&CountSpeicesInCluster(join("\t",@row))."\t".&getCOG(\@row,\%cog)."\t".&getDescription(\@row,\%description)."\n";
+			}
+		close(R);
+
+#### COG distribution
+
+###Whole Clusters COG Distribution
+&outputCOGStatistic("${outputDIR}5.Orthologs_Whole_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",$spnum,1));
+
+###Core Clusters COG Distribution
+&outputCOGStatistic("${outputDIR}5.Orthologs_Core_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",$spnum,$spnum));
+
+###Dispensable Clusters COG Distribution
+&outputCOGStatistic("${outputDIR}5.Orthologs_Dispensable_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",($spnum-1),2));
+
+###strains specifc Clusters COG Distribution
+&outputCOGStatistic("${outputDIR}5.Orthologs_specifc_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",1,1));
+
+system("rm -rf *.function");
+
+	}else
+		{
+			print "Function analysis is skipped.\n";
+		}
+
+sub outputCOGStatistic()
+	{
+		(my $file,my $subcogcount)=@_;
+		my @cogcat=("J A K L B","D Y V T M N Z W U O","C G E F H I P Q","R S -");
+		my @cogdesc=("INFORMATION STORAGE AND PROCESSING","CELLULAR PROCESSES AND SIGNALING","METABOLISM","POORLY CHARACTERIZED");
+		my @subcogcat=qw(J A K L B D Y V T M N Z W U O C G E F H I P Q R S -);
+		my @subcogdesc=("[J] Translation, ribosomal structure and biogenesis","[A] RNA processing and modification","[K] Transcription","[L] Replication, recombination and repair","[B] Chromatin structure and dynamics","[D] Cell cycle control, cell division, chromosome partitioning","[Y] Nuclear structure","[V] Defense mechanisms","[T] Signal transduction mechanisms","[M] Cell wall/membrane/envelope biogenesis","[N] Cell motility","[Z] Cytoskeleton","[W] Extracellular structures","[U] Intracellular trafficking, secretion, and vesicular transport","[O] Posttranslational modification, protein turnover, chaperones","[C] Energy production and conversion","[G] Carbohydrate transport and metabolism","[E] Amino acid transport and metabolism","[F] Nucleotide transport and metabolism","[H] Coenzyme transport and metabolism","[I] Lipid transport and metabolism","[P] Inorganic ion transport and metabolism","[Q] Secondary metabolites biosynthesis, transport and catabolism","[R] General function prediction only","[S] Function unknown","[-] Unclassified");
+		my %subcogdesc;
+		my $key;
+		my @cog;
+		my $i;
+		my $cognum;
+
+		for ($i=0;$i<@subcogcat;$i++) 
+			{
+				$subcogdesc{$subcogcat[$i]}=$subcogdesc[$i];
+			}
+
+		open(R,">$file");
+		for ($i=0;$i<@cogcat;$i++) 
+			{
+				$cognum=0;
+				foreach $key (split(" ",$cogcat[$i])) 
+					{
+						$cognum=$cognum+$$subcogcount{$key};
+					}
+				print R $cogdesc[$i]." ( ".$cognum." )\n";
+				foreach $key (split(" ",$cogcat[$i])) 
+					{
+						printf R "%-6d   %s\n",$$subcogcount{$key},$subcogdesc{$key};
+					}
+				print R "\n";
+				
+			}
+		close(R);
+
+	}
+
+sub scanCOG()
+	{
+		(my $file,my $max_orth,my $min_orth)=@_;
+		my @row;
+		my @subcogcat=qw(J A K L B D Y V T M N Z W U O C G E F H I P Q R S -);
+		my %subcogcount;
+		my $cog;
+		my $key;
+
+		foreach $key (@subcogcat) 
+			{
+				$subcogcount{$key}=0;
+			}
+
+		@subcogcat=qw(J A K L B D Y V T M N Z W U O C G E F H I P Q R S);
+
+		open(F,"$file");
+		$_=<F>;
+		while (<F>) 
+			{
+				@row=split(/\t/,$_);
+				if ($row[1]>=$min_orth and $row[1]<=$max_orth) 
+					{
+						if ($row[2] eq "-") 
+							{
+								$subcogcount{"-"}++;
+							}else
+								{
+									$_=uc($row[2]);
+									s/COG//gi;
+									$cog=$_;
+									foreach $key (@subcogcat) 
+										{
+											if ($cog=~/$key/) 
+												{
+													$subcogcount{$key}++;
+												}
+										}
+								}
+					}
+			}
+		close(F);
+
+		return \%subcogcount;
+	}
+
+
+
+
+sub getCOG()
+	{
+		(my $data,my $coghash)=@_;
+		my $cog="";
+		my @cog;
+		my $key;
+		my %hash;
+		my @gene=split(/\t|\,/,join("\t",@$data));
+		@gene=grep(/^S/,@gene);
+
+		foreach $key (@gene) 
+			{
+				if (($$coghash{$key} ne "-") and ($$coghash{$key} ne "")) 
+					{
+						$cog=$cog.",".$$coghash{$key};
+					}
+			}
+		@cog=split(/,/,$cog);
+		foreach $cog (@cog) 
+			{
+				if ($cog ne "") 
+					{
+						$hash{$cog}=1;
+					}
+			}
+
+		$cog=join(",",(keys %hash));
+		if ($cog eq "") 
+			{
+				$cog="-";
+			}
+		return $cog;
+	}
+
+sub getDescription()
+	{
+		(my $data,my $deschash)=@_;
+		my $desc="";
+		my $key;
+		my @gene=split(/\t|\,/,join("\t",@$data));
+		@gene=grep(/^S/,@gene);
+
+		foreach $key (@gene) 
+			{
+				if ( ($$deschash{$key} ne "") and ($$deschash{$key} ne "-") and ($$deschash{$key}!~/hypothetical/)) 
+					{
+						$desc=$$deschash{$key};
+					}
+			}
+
+		if ($desc eq "") 
+			{
+				$desc="hypothetical protein";
+			}
+
+		return $desc;
+	}
+
+
+
+sub ReadAnnotation()
+	{
+		(my $species,my $cog,my $description)=@_;
+		my $i;
+		my @row;
+
+		for ($i=0;$i<@$species;$i++) 
+			{
+				open(F,"$$species[$i].function");
+				while (<F>) 
+					{
+						chomp($_);
+						@row=split(/\t/,$_);
+						if (scalar(@row)>=2) 
+							{
+								$$cog{$row[0]}=$row[1];
+							}else
+								{
+									$$cog{$row[0]}="-";
+								}
+
+						if (scalar(@row)>=3) 
+							{
+								$$description{$row[0]}=$row[2];
+							}else
+								{
+									$$description{$row[0]}="hypothetical protein";
+								}
+					}
+				close(F);
+			}
+	}
+
+
+sub SNPBasedTree()
+	{
+		(my $infile,my $outfileprefix)=@_;
+		my $tmpin=$infile;
+		my $tmpout;
+
+
+#### boootstrap
+print "\n#### seqboot ...\n\n";
+		open(R,">seqboot.cmd");
+		#print R "$tmpin\n";
+		print R "R\n";
+		print R "$bootstrap\n";
+		print R "Y\n";
+		print R "1\n";
+		close(R);
+
+		system("cp $tmpin infile");
+		system("$seqboot < seqboot.cmd");
+		system("mv outfile 100dnaseq");
+		system("rm -rf infile");
+		system("rm seqboot.cmd");                # 100dnasseq
+
+#### dnaml
+print "\n#### dnaml ...\n\n";
+		open(R,">dnaml.cmd");
+		#print R "100dnaseq\n";
+		print R "T\n";
+		print R "25\n";
+	if ($bootstrap>1) 
+		{
+		print R "M\n";
+		print R "D\n";
+		#print R "100\n";
+		print R "$bootstrap\n";
+		print R "1\n"; # Random number seed (must be odd)?
+		print R "5\n"; # Number of times to jumble?
+		}
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnaseq infile");
+		system("$dnaml < dnaml.cmd");
+		system("rm -rf outfile");
+		system("rm -rf infile");
+		system("mv outtree 100dnaseqtree");     # 100dnaseq, 100dnaseqtree
+
+#### consense
+print "\n#### dnaml consense ...\n\n";
+
+		open(R,">consense.cmd");
+		#print R "100dnaseqtree\n";
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnaseqtree intree");
+		system("$consense < consense.cmd");
+		system("mv outfile ${outfileprefix}.ML.outfile");
+		system("mv outtree ${outfileprefix}.ML.tree");
+		system("rm -rf infile");
+		system("rm -rf 100dnaseqtree");                        # 100dnaseq
+
+#### dnadist
+print "\n#### dnadist ...\n\n";
+		open(R,">dnadist.cmd");
+		#print R "100dnaseq\n";
+		print R "T\n";
+		print R "25\n";
+	if ($bootstrap>1) 
+		{
+		print R "M\n";
+		print R "D\n";
+		#print R "100\n";
+		print R "$bootstrap\n";
+		}
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnaseq infile");
+		system("$dnadist < dnadist.cmd");
+		system("rm -rf 100dnaseq");
+		system("rm -rf infile");
+		system("mv outfile 100dnadist");                 # 100dnadist
+
+#### Neighbor-joining tree
+print "\n#### Neighbor-joining ...\n\n";
+		open(R,">NJ.cmd");
+	if ($bootstrap>1) 
+		{
+		#print R "100dnadist\n";
+		print R "M\n";
+		#print R "100\n";
+		print R "$bootstrap\n";
+		print R "1\n";
+		}
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnadist infile");
+		system("$neighbor < NJ.cmd");
+		system("mv outtree 100dnadistNJtree");
+		system("rm outfile");
+		system("rm -rf infile");
+		system("rm -rf NJ.cmd");                      # 100dnadist,100dnadistNJtree
+
+#### NJ-consense
+print "\n#### NJ-consense ...\n\n";
+		open(R,">NJ-consense.cmd");
+		#print R "100dnadistNJtree\n";
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnadistNJtree intree");
+		system("$consense < NJ-consense.cmd");
+		system("mv outfile ${outfileprefix}.Neighbor-joining.outfile");
+		system("mv outtree ${outfileprefix}.Neighbor-joining.tree");
+		system("rm -rf NJ-consense.cmd");
+		system("rm -rf intree");
+		system("rm -rf 100dnadistNJtree");
+
+
+#### UPGMA tree
+print "\n#### UPGMA ...\n\n";
+		open(R,">UPGMA.cmd");
+		#print R "100dnadist\n";
+		print R "N\n";
+	if ($bootstrap>1) 
+		{
+		print R "M\n";
+		#print R "100\n";
+		print R "$bootstrap\n";
+		print R "1\n";
+		}
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnadist infile");
+		system("$neighbor < UPGMA.cmd");
+		system("mv outtree 100dnadistUPGMAtree");
+		system("rm -rf outfile");
+		system("rm -rf infile");
+		system("rm -rf UPGMA.cmd");
+
+#### UPGMA-consense
+print "\n#### UPGMA-consense ...\n\n";
+		open(R,">UPGMA-consense.cmd");
+		#print R "100dnadistUPGMAtree\n";
+		print R "Y\n";
+		close(R);
+
+		system("cp 100dnadistUPGMAtree intree");
+		system("$consense < UPGMA-consense.cmd");
+		system("mv outfile ${outfileprefix}.UPGMA.outfile");
+		system("mv outtree ${outfileprefix}.UPGMA.tree");
+		system("rm -rf UPGMA-consense.cmd");
+		system("rm -rf 100dnadistUPGMAtree");
+		system("rm -rf intree");
+
+###CLEAN TMP FILE
+
+		system("rm -rf *.cmd");
+		system("rm -rf 100dnadist");
+	}
+
+sub PanBasedTree()
+	{
+		(my $infile,my $outfileprefix)=@_;
+		my $tmpin;
+		my $tmpout;
+
+		$tmpin=$infile;
+#### Neighbor-joining tree
+
+		open(R,">NJ.cmd");
+		#print R "$tmpin\n";
+		print R "Y\n";
+		close(R);
+
+		system("cp $tmpin infile");
+		system("$neighbor < NJ.cmd");
+		system("mv outfile ${outfileprefix}.Neighbor-joining.outfile");
+		system("mv outtree ${outfileprefix}.Neighbor-joining.tree");
+		system("rm -rf  NJ.cmd");
+		system("rm -rf infile");
+
+#### UPGMA tree
+
+		open(R,">UPGMA.cmd");
+		#print R "$tmpin\n";
+		print R "N\n";
+		print R "Y\n";
+		close(R);
+
+		system("cp $tmpin infile");
+		system("$neighbor < UPGMA.cmd");
+		system("mv outfile ${outfileprefix}.UPGMA.outfile");
+		system("mv outtree ${outfileprefix}.UPGMA.tree");
+		system("rm -rf UPGMA.cmd");
+		system("rm -rf infile");
+
+
+###CLEAN TMP FILE
+		system("rm -rf *.cmd");
+	}
+
+sub DistanceMatrix()
+	{
+		(my $spnum,my $hash)=@_;
+		my $i;
+		my $j;
+		my $k;
+		my $dist;
+		my $ref;
+		my $query;
+		foreach $i (0..($spnum-1)) 
+			{
+				foreach $j ($i..($spnum-1)) 
+					{
+						$ref=$$hash{$i};
+						$query=$$hash{$j};
+						$dist=0;
+						for ($k=0;$k<length($ref);$k++) 
+							{
+								if (substr($ref,$k,1) ne substr($query,$k,1)) 
+									{
+										$dist++;
+									}
+							}
+						$$hash{$i."-".$j}=$dist;
+						$$hash{$j."-".$i}=$dist;
+					}
+			}
+	}
+
+
+sub ClusterProfil4Specie
+	{
+		(my $hash)=@_;
+		my @row;
+		my $i;
+
+		foreach  (0..($spnum-1))   #initialization Hash
+			{
+				$$hash{$_}="";
+			}
+
+		foreach  (@clusters) 
+			{
+				@row=split(/\t/,$_);
+				splice(@row,0,1);
+				if (&CountSpeicesInCluster(join("\t",@row))>1) 
+					{
+					for ($i=0;$i<@row;$i++) 
+						{
+							if ($row[$i] eq "-") 
+								{
+									$$hash{$i}=$$hash{$i}."0";
+								}else
+									{
+										$$hash{$i}=$$hash{$i}."1";
+									}
+						}
+					}
+			}
+	}
+
+# &ExtractSNP4tree(\%tmpHash,\%nt4tree);
+
+sub ExtractSNP4tree()
+	{
+		(my $hash,my $nt4treeRef)=@_;
+		my $key;
+		my @row;
+		my $i;
+		my $len;
+		my @tribases;
+		foreach $key (keys %$hash) 
+			{
+				$$hash{substr($key,0,index($key,"G"))}=$$hash{$key};
+				delete($$hash{$key});
+			}
+
+		for ($i=0;$i<$spnum;$i++)
+		{
+			$nt4tree{"S".$i}=$nt4tree{"S".$i}.$$hash{"S".$i};
+		}
+	}
+
+
+=pod
+sub ExtractSNP4tree()
+	{
+		(my $hash,my $nt4treeRef)=@_;
+		my $key;
+		my @row;
+		my $i;
+		my $len;
+		my @tribases;
+		foreach $key (keys %$hash) 
+			{
+				$$hash{substr($key,0,index($key,"G"))}=$$hash{$key};
+				delete($$hash{$key});
+			}
+		@_=(keys %$hash);
+		$len=length($_[0]);
+		for ($j=0;3*$j<$len;$j++) 
+			{
+##### scanning each codon
+				for ($i=0;$i<$spnum;$i++) 
+					{
+						$tribases[$i]=substr($$hash{"S".$i},3*$j,3);
+					}
+##### checking each codon
+				if (&IsTheSame(@tribases) ==0) 
+					{
+					for ($i=0;$i<@tribases;$i++) 
+						{
+							$nt4tree{"S".$i}=$nt4tree{"S".$i}.$tribases[$i];
+						}
+					}
+			}
+	}
+=cut
+
+
+sub pal2nal()
+	{
+		(my $pal,my $nuc, my $nal)=@_;
+		my %aaAln=();
+		my %ffn=();
+		my %ntAln=();
+		my %nt;
+		my $dna;
+		my $nt;
+		my $key;
+		my $flag=1;
+		my $i=0;
+		my $j;
+
+### read protein aligment result
+		&ReadAlignmentToHash("$pal",\%aaAln);
+### read nt sequences
+		&ReadSequenceInToHash("$nuc",\%ffn);
+		foreach $key (keys %ffn) 
+			{
+				$dna=$ffn{$key};
+				#if (int(length($nt{$key})/3)*3 ne length($nt{$key})) 
+				if (int(length($dna)/3)*3 ne length($dna)) 
+					{
+						$flag=0;
+						print "The length of nucleotide sequence is not 3 integer times.\n";
+						last;
+					}else
+						{
+							for ($i=0;$i<(length($dna)/3);$i++) 
+								{
+									$nt{$key."|".$i}=substr($dna,$i*3,3);
+								}
+						}
+			}
+
+		if ($flag==0) 
+			{
+				return 0;
+			}else
+				{
+					foreach $key (keys %aaAln)  ### replace aa with corresponding nt
+						{
+							$nt="";
+							$i=0;
+							for ($j=0;$j<length($aaAln{$key});$j++) 
+								{
+									if (substr($aaAln{$key},$j,1) eq "-") 
+										{
+											$nt=$nt."---";
+										}else
+											{
+												$nt=$nt.$nt{$key."|".$i};
+												$i++;
+											}
+								}
+							$ntAln{$key}=$nt;
+						}
+
+					### output 
+						open(R,">$nal");
+						foreach  (keys %ntAln) 
+							{
+								print R ">$_\n".$ntAln{$_}."\n";
+							}
+						close(R);
+
+					return 1;
+				}
+	}
+
+
+
+
+sub DetectSNP()
+	{
+		my %faa;
+		my %ffn;
+		my @row;
+		my $count_gene;
+		my $count_sp;
+		my @genelist;
+		my $i;
+		my $j;
+		my $pepalnlen;
+		my @cdsvar=qw();
+		my $cdi=0;
+		my @tribases;
+		my @bases;
+		my @aa;
+
+
+### fetch gene list
+		open(F,"$ClusterID.pep");
+		@genelist=<F>;
+		close(F);
+		@genelist=grep(/^>/,@genelist);
+		chomp(@genelist);
+		$_=join("\t",@genelist);
+		s/>//g;
+		@genelist=split(/\t/,$_);
+
+### count gene number and species number
+		@row=split(/\t/,$clusters[$ClusterID-1]);
+		splice(@row,0,1);
+		$count_sp=&CountSpeicesInCluster(join("\t",@row));
+		$count_gene=&CountGeneInCluster(join("\t",@row));
+
+### read alignment sequences
+		&ReadAlignmentToHash("$ClusterID.pal",\%faa);
+		&ReadAlignmentToHash("$ClusterID.nal",\%ffn);
+
+@_=(keys %faa);
+$pepalnlen=length($faa{$_[0]});
+### scan SNP
+		for ($i=1;$i<=$pepalnlen;$i++) 
+			{
+				@tmp=qw();
+				@tribases=qw();
+				for ($j=0;$j<@genelist;$j++) ### fetch triplet codon
+					{
+						$tribases[$j]=substr($ffn{$genelist[$j]},3*($i-1),3);
+					}
+				if (&IsTheSame(@tribases) ==0) ### if triplet codon is not consistent
+					{
+						@aa=qw();
+						for ($j=0;$j<@genelist;$j++) 
+							{
+								$aa[$j]=substr($faa{$genelist[$j]},($i-1),1);
+							}
+						if (&IsTheSame(@aa) ==0) ### aa is not consistent
+							{
+								if (join("",@aa) =~/-/) 
+									{
+										$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".$i."\t".&CharType(\@aa)."\t-\t-\tInDel\n";
+									}else
+										{
+									#### base 1
+											for ($j=0;$j<@genelist;$j++) 
+												{
+													$bases[$j]=substr($ffn{$genelist[$j]},3*($i-1),1);
+												}
+											if (&IsTheSame(@bases) ==0) 
+												{
+													$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".($i+0.1)."\t".&CharType(\@aa)."\t".&CharType(\@bases)."\t".join("",@bases)."\tNonsynonymous mutation\n";
+												}
+									#### base 2
+											for ($j=0;$j<@genelist;$j++) 
+												{
+													$bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+1,1);
+												}
+											if (&IsTheSame(@bases) ==0) 
+												{
+													$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".($i+0.2)."\t".&CharType(\@aa)."\t".&CharType(\@bases)."\t".join("",@bases)."\tNonsynonymous mutation\n";
+												}
+									#### base 3
+											for ($j=0;$j<@genelist;$j++) 
+												{
+													$bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+2,1);
+												}
+											if (&IsTheSame(@bases) ==0) 
+												{
+													$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".($i+0.3)."\t".&CharType(\@aa)."\t".&CharType(\@bases)."\t".join("",@bases)."\tNonsynonymous mutation\n";
+												}
+										}
+							}else
+								{
+									#### base 1
+											for ($j=0;$j<@genelist;$j++) 
+												{
+													$bases[$j]=substr($ffn{$genelist[$j]},3*($i-1),1);
+												}
+											if (&IsTheSame(@bases) ==0) 
+												{
+													$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".($i+0.1)."\t".&CharType(\@aa)."\t".&CharType(\@bases)."\t".join("",@bases)."\tSynonymous mutation\n";
+												}
+									#### base 2
+											for ($j=0;$j<@genelist;$j++) 
+												{
+													$bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+1,1);
+												}
+											if (&IsTheSame(@bases) ==0) 
+												{
+													$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".($i+0.2)."\t".&CharType(\@aa)."\t".&CharType(\@bases)."\t".join("",@bases)."\tSynonymous mutation\n";
+												}
+									#### base 3
+											for ($j=0;$j<@genelist;$j++) 
+												{
+													$bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+2,1);
+												}
+											if (&IsTheSame(@bases) ==0) 
+												{
+													$cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".($i+0.3)."\t".&CharType(\@aa)."\t".&CharType(\@bases)."\t".join("",@bases)."\tSynonymous mutation\n";
+												}
+								}
+					}
+			}
+		return @cdsvar;
+	}
+
+
+sub VarAnalysis()
+	{
+		(my $data)=@_;
+		my @data=@$data;
+		my $indel=0;
+		my $syn=0;
+		my $nonsyn=0;
+		my @tmp;
+		$indel=scalar(grep(/InDel$/,@data));
+		$nonsyn=scalar(grep(/Nonsynonymous mutation$/,@data));;
+		$syn=scalar(grep(/Synonymous mutation$/,@data));
+		return "$indel\t$nonsyn\t$syn";
+	}
+
+
+
+sub CharType()
+	{
+		(my $str)=@_;
+		my %hash;
+		my @data=@$str;
+		foreach  (@data) 
+			{
+				$hash{$_}=1;
+			}
+		return join(",",(keys %hash));
+	}
+
+sub IsTheSame()
+	{
+		(my @data)=@_;
+		my %hash;
+		foreach  (@data) 
+			{
+				$hash{$_}=1;
+			}
+		if (scalar(keys %hash) ==1) 
+			{
+				return 1;
+			}else
+				{
+					return 0;
+				}
+	}
+
+
+
+sub FormatClusterOutPut()
+	{
+		(my $speices,my $file,my $cluster)=@_;
+		my @row;
+		my $gid=1;
+		my $key;
+		my %hash;
+		my $gene;
+		my @tmp;
+		my $i;
+		my $j;
+		open(R,">$file");
+		print R "ClutserID\t".join("\t",@$speices)."\n";
+		foreach $key (@$cluster) 
+			{
+				@row=split(/\t/,$key);
+				for ($i=0;$i<@row;$i++) 
+					{
+						if ($row[$i] ne "-") 
+							{
+								@tmp=split(/,/,$row[$i]);
+								for ($j=0;$j<@tmp;$j++) 
+									{
+										$_=$tmp[$j];
+										s/^S[0-9]+G//;
+										$tmp[$j]=$_;
+									}
+								$row[$i]=join(",",@tmp);
+							}
+					}
+				print R $gid."\t".join("\t",@row)."\n";
+				$gid++;
+			}
+		close(R);
+	}
+sub RetrieveClusterFromFile()
+	{
+		(my $file,my $clusters)=@_;
+		my @content;
+		my @row;
+		my $spid;
+		my $line=0;
+		my $i=0;
+		my $j;
+		my @tmp;
+		open(F,$file) or die "Could open $file\n";
+		@content=<F>;
+		close(F);
+		splice(@content,0,1);
+		chomp(@content);
+		foreach  (@content) 
+			{
+				@row=split(/\t/,$_);
+				$$clusters[$line]=$row[0];
+				splice(@row,0,1);
+				for ($i=0;$i<@row;$i++) 
+					{
+						if ($row[$i] ne "-") 
+							{
+								@tmp=split(/,/,$row[$i]);
+								for ($j=0;$j<@tmp;$j++) 
+									{
+										$tmp[$j]="S${i}G".$tmp[$j];
+									}
+								$row[$i]=join(",",@tmp);
+							}
+					}
+				$$clusters[$line]=$$clusters[$line]."\t".join("\t",@row)."\n";
+				$line++;
+			}
+	}
+
+
+
+
+sub GeneDistribution()
+	{
+		(my $clusters,my $hash)=@_;
+		my @row;
+		my $spid;
+		my $orth;
+		my $key;
+		foreach (@$clusters) 
+			{
+				@row=split(/\t/,$_);
+				splice(@row,0,1);
+				$orth=&CountSpeicesInCluster(join("\t",@row));
+				@row=split(/\t|\,/,join("\t",@row));
+				foreach $key (@row) 
+					{
+						if ($key ne "-") 
+							{
+							$spid=substr($key,1,(index($key,'G')-1)); ###extract strains id
+							if (exists($$hash{$orth."|".$spid})) 
+								{
+									$$hash{$orth."|".$spid}++;
+								}else
+									{
+										$$hash{$orth."|".$spid}=1;
+									}
+							}
+					}
+			}
+	}
+
+sub CountSpeicesInCluster()
+	{
+		(my $str)=@_;
+		chomp($str);
+		my @list=split(/\t/,$str);
+		my $key;
+		my $count=0;
+
+		foreach $key (@list) 
+			{
+				if ($key ne "-") 
+					{
+						$count++;
+					}
+			}
+		return $count;
+	}
+
+sub CountGeneInCluster()
+	{
+		(my $str)=@_;
+		chomp();
+		my @list=split(/\t|\,/,$str);
+		my $key;
+		my $count=0;
+		foreach $key (@list) 
+			{
+				if ($key ne "-") 
+					{
+						$count++;
+					}
+			}
+		return $count;
+	}
+
+
+
+
+sub GF()
+	{
+		&PrepareFasta(\@species,$inputDIR,".pep"); ###prepare pep file
+		system("cat ".join(".pep ",@species).".pep > All.pep");
+		system("grep '>' All.pep > genelist");
+		system("$formatdb -p T -i All.pep");
+		system("$blastall -p blastp -i All.pep -d All.pep -M BLOSUM45 -m9 -e $evalue -o All.blastp -a $thread");
+		system("perl ./Blast_Filter.pl All.blastp All.pep $coverage $identity $score | $mcl - --abc -I 2.0 -o All.cluster");
+		&FormatCluster("All.cluster","genelist",$spnum,\@clusters);
+		#system("rm -rf *.pep* All.blastp All.cluster genelist");
+	}
+
+sub MP()
+	{
+#		(my $species,my $inputDIR,my $thread,my $evalue,my $score,my $coverage,my $identity)=@_;
+		my $i;
+		my $j;
+		&PrepareFasta(\@species,$inputDIR,".pep"); ###prepare pep file
+		system("cat ".join(".pep ",@species).".pep > All.pep");
+		system("grep '>' All.pep > genelist");
+		system("rm -rf All.pep");
+		for ($i=0;$i<$spnum;$i++) 
+			{
+				for ($j=$i+1;$j<$spnum;$j++) 
+					{
+						system("perl ./inparanoid.pl $blastall $thread $formatdb $score $global $local $species[$i].pep $species[$j].pep");
+					}
+			}
+		system("perl ./multiparanoid.pl -species ".join(".pep+",@species).".pep -unique 1");
+###convert the MP result to table list based on gene
+		&MP_Result_to_Table("MP.Cluster","All.cluster");
+		&FormatCluster("All.cluster","genelist",$spnum,\@clusters);
+		system("rm -rf sqltable.* *.pep* MP.Cluster genelist");
+	}
+
+sub fasta2phylip()
+	{
+		(my $input,my $output)=@_;
+		use Bio::AlignIO;
+		my $inputfilename = "10.aln";
+		my $in= Bio::AlignIO->new(-file   => $input ,
+		                          -format => 'fasta');
+		my $out = Bio::AlignIO->new(-file   => ">$output" ,
+		                          -format => 'phylip');
+		while ( my $aln = $in->next_aln() ) 
+			{
+			$out->write_aln($aln);
+			}
+	}
+
+
+sub RemoveHeadGap()
+	{
+		(my $nal,my $hash)=@_;
+		my %aln;
+		my $key;
+		my $gaplength=0;
+		my $len1;
+		my $len2;
+		&ReadSequenceInToHash("$nal",\%aln);
+		foreach $key (keys %aln) 
+			{
+				$len1=length($aln{$key});
+				$_=$aln{$key};
+				s/^-+//;
+				$len2=length($_);
+				if (($len1-$len2)>$gaplength) 
+					{
+						$gaplength=$len1-$len2;
+					}
+			}
+		foreach $key (keys %aln) 
+			{
+				$$hash{$key}=$$hash{$key}.substr($aln{$key},$gaplength,(length($aln{$key})-$gaplength));
+			}
+	}
+
+sub PrepareFasta()
+	{
+		(my $species,my $inputDIR,my $extention)=@_;
+		my $sp;
+		my $file;
+		my $i;
+		my %hash;
+		my $key;
+		for ($i=0;$i<@$species;$i++) 
+			{
+				$file=$inputDIR.$$species[$i].$extention;
+				%hash=();
+				&ReadSequenceInToHash($file,\%hash);
+				open(R,">$$species[$i]${extention}") or die "Could write into $file\n";
+				foreach $key (keys %hash) 
+					{
+						print R ">S${i}G$key\n";
+						print R $hash{$key}."\n";
+					}
+				close(R);
+			}
+	}
+
+sub PrepareTable()
+	{
+		(my $species,my $inputDIR,my $extention)=@_;
+		my @content;
+		my $i;
+		my @row;
+		my $file;
+		for ($i=0;$i<@$species;$i++) 
+			{
+				$file=$inputDIR.$$species[$i].$extention;
+				open(F,$file) or die "Could open $file\n";
+				@content=<F>;
+				close(F);
+				chomp(@content);
+				open(R,">$$species[$i]${extention}") or die "Could write into $file\n";
+				foreach  (@content) 
+					{
+						@row=split(/\t/,$_);
+						$row[0]="S${i}G$row[0]";
+						if ($extention eq ".location") 
+							{
+								$row[0]=$row[0]."\t".$row[0];
+							}
+						print R join("\t",@row)."\n";
+					}
+				close(R);
+			}
+	}
+
+sub CheckExtraProgram
+	{
+		#(my $section, my $method, my $tmparray)=@_;
+		my @error;
+		my $ei=0;
+
+#####cluster gene 
+		if (substr($section,0,1) eq "1") 
+			{
+###MP: blastall formatdb 
+###GF: blastall formatdb mcl
+
+			if (!(-e $formatdb)) 
+				{
+					$error[$ei++]="formatdb is not found at $formatdb\n";
+				}
+
+			if (!(-X $formatdb)) 
+				{
+					$error[$ei++]="there is not premission to execute $formatdb\n";
+				}
+
+			if (!(-e $blastall)) 
+				{
+					$error[$ei++]="blastall is not found at $blastall\n";
+				}
+
+			if (!(-X $blastall)) 
+				{
+					$error[$ei++]="there is not premission to execute $blastall\n";
+				}
+
+				if ($method eq "GF") 
+					{
+						if (!(-e $mcl)) 
+							{
+								$error[$ei++]="mcl is not found at $mcl\n";
+							}
+						if (!(-X $mcl)) 
+							{
+								$error[$ei++]="there is not premission to execute $mcl\n";
+							}
+					}
+			}
+
+#####CDS variation
+		if (substr($section,2,1) eq "1") 
+			{
+				if (!(-e $mafft)) 
+					{
+						$error[$ei++]="mafft is not found at $mafft\n";
+					}
+				if (!(-X $mafft)) 
+					{
+						$error[$ei++]="there is not premission to execute $mafft\n";
+					}
+			}
+
+#####CDS variation
+		if (substr($section,3,1) eq "1") 
+			{
+				if (!(-e $mafft)) 
+					{
+						$error[$ei++]="mafft is not found at $mafft\n";
+					}
+				if (!(-X $mafft)) 
+					{
+						$error[$ei++]="there is not premission to execute $mafft\n";
+					}
+			}
+#####Evolution analysis
+		if (substr($section,3,1) eq "1") 
+			{
+				if (-e $seqboot) 
+					{
+						$error[$ei++]="there is not premission to execute $seqboot\n" if(!(-X $seqboot));
+					}else
+						{
+							$error[$ei++]="seqboot is not found at $seqboot\n";
+						}
+				if (-e $dnaml) 
+					{
+						$error[$ei++]="there is not premission to execute $dnaml\n" if(!(-X $dnaml));
+					}else
+						{
+							$error[$ei++]="dnaml is not found at $dnaml\n";
+						}
+				if (-e $dnadist) 
+					{
+						$error[$ei++]="there is not premission to execute $dnadist\n" if(!(-X $dnadist));
+					}else
+						{
+							$error[$ei++]="dnadist is not found at $dnadist\n";
+						}
+				if (-e $neighbor) 
+					{
+						$error[$ei++]="there is not premission to execute $neighbor\n" if(!(-X $neighbor));
+					}else
+						{
+							$error[$ei++]="neighbor is not found at $neighbor\n";
+						}
+				if (-e $consense) 
+					{
+						$error[$ei++]="there is not premission to execute $consense\n" if(!(-X $consense));
+					}else
+						{
+							$error[$ei++]="consense is not found at $consense\n";
+						}
+				if (-e $dnapars) 
+					{
+						$error[$ei++]="there is not premission to execute $dnapars\n" if(!(-X $dnapars));
+					}else
+						{
+							$error[$ei++]="dnapars is not found at $dnapars\n";
+						}
+			}
+		#@$tmparray=(@$tmparray,@error);
+		@tmp=(@tmp,@error);
+	}
+
+
+sub CheckInputFile()
+	{
+		(my $species,my $inputDIR,my $section,my $method,my $tmparray)=@_;
+####cluster
+		if (substr($section,0,1) eq "1") 
+			{
+				if ($method eq "MM") 
+					{
+						@$tmparray=(@$tmparray,&chk2SEQ($species,$inputDIR)); ### check pep and nuc
+						@$tmparray=(@$tmparray,&chktab($species,$inputDIR,".location"));### chk  pep nuc location
+					}else
+						{
+						@$tmparray=(@$tmparray,&chk1SEQ($species,$inputDIR));
+						}
+			}
+###CDS variation
+		if (substr($section,2,1) eq "1")
+			{
+				@$tmparray=(@$tmparray,&chk2SEQ($species,$inputDIR));
+			}
+###function analysis
+		if (substr($section,4,1) eq "1")
+			{
+				@$tmparray=(@$tmparray,&chktab($species,$inputDIR,".function"));
+			}
+	}
+
+
+sub chk1SEQ()
+	{
+		(my $species,my $inputDIR)=@_;
+		my @error;
+		my $ei=0;
+		my $sp;
+		my $pepfile;
+		my %pep;
+		foreach $sp (@$species) 
+			{
+				%pep=();
+				$pepfile=$inputDIR.$sp.".pep";
+				&ReadSequenceInToHash($pepfile,\%pep);
+				if (scalar(keys %pep)<2) 
+					{
+						$error[$ei++]="format error in $pepfile\n";
+					}
+			}
+		return @error;
+	}
+
+sub chk2SEQ()
+	{
+		(my $species,my $inputDIR)=@_;
+		my $sp;
+		my %pep;
+		my %nuc;
+		my $pepfile;
+		my $nucfile;
+		my $key;
+		my @error;
+		my $ei=0;
+		foreach $sp (@$species) 
+			{
+				$pepfile=$inputDIR.$sp.".pep";
+				$nucfile=$inputDIR.$sp.".nuc";
+				%pep=();
+				%nuc=();
+				&ReadSequenceInToHash("$pepfile",\%pep);
+				&ReadSequenceInToHash("$nucfile",\%nuc);
+				if (scalar(keys %pep) ne scalar(keys %nuc)) 
+					{
+						$error[$ei++]="Sequences number is not consistent in the following two file:\n\t$pepfile\n\t$nucfile\n";
+					}else
+						{
+							foreach $key (keys %pep) 
+								{
+									if (exists($nuc{$key})) 
+										{
+											if (length($nuc{$key}) ne ((length($pep{$key})+1)*3)) 
+												{
+													$error[$ei++]="the length of $key in $nucfile is not consistent with its corresponding protein length\n";
+												}
+										}else
+											{
+												$error[$ei++]="$key lost in $nucfile\n";
+											}
+								}
+
+							foreach $key (keys %nuc) 
+								{
+									if (!exists($pep{$key})) 
+										{
+											$error[$ei++]="1048 $key lost in $pepfile\n";
+										}
+								}
+						}
+			}
+		return @error;
+	}
+
+
+sub chktab()
+	{
+		(my $species,my $inputDIR,my $extention)=@_;
+		my %pep;
+		my @row;
+		my $key;
+		my %tab;
+		my @error;
+		my $ei=0;
+		my $sp;
+		my $tabfile;
+		my $pepfile;
+		foreach $sp (@$species) 
+			{
+				%tab=();
+				%pep=();
+				$tabfile=$inputDIR.$sp.$extention;
+				open(F,"$tabfile");
+				while (<F>) 
+					{
+						chomp();
+						@row=split(/\t/,$_);
+						if (scalar(@row)<3) 
+							{
+								$error[$ei++]="format error in $tabfile\n";
+							}else
+								{
+									$tab{$row[0]}=$row[1];
+								}
+					}
+				close(F);
+				$pepfile=$inputDIR.$sp.".pep";
+				&ReadSequenceInToHash($pepfile,\%pep);
+				foreach $key (keys %pep) 
+					{
+						if (!exists($tab{$key})) 
+							{
+								$error[$ei++]="sequence $key lost infomation in $tabfile\n";
+							}
+					}
+			}
+		return @error;
+	}
+
+
+
+sub ReadSequenceInToHash()
+	{
+		use Bio::SeqIO;
+		(my $file,my $hash)=@_;
+		my $seq;
+		my $in=Bio::SeqIO->new(-file=>"$file",-format=>"fasta");
+		while ($seq=$in->next_seq()) 
+			{
+				#$$hash{$id."|".$seq->id}=$seq->seq();
+				$$hash{$seq->id}=$seq->seq();
+			}
+	}
+
+sub ReadAlignmentToHash()
+	{
+		(my $file,my $hash)=@_;
+		my $name="";
+		my $seq="";
+		my @content;
+		my $line;
+		open(F,"$file");
+		@content=<F>;
+		close(F);
+		chomp(@content);
+		for ($line=0;$line<@content;$line++) 
+			{
+				if ($content[$line]=~/^>/) 
+					{
+						if ($line>0) 
+							{
+								$$hash{$name}=$seq;
+								$name="";
+							}
+
+						$_=$content[$line];
+						s/^>//;
+						$name=$_;
+						$seq="";
+					}else
+						{
+							if ($name ne "") 
+								{
+								$seq=$seq.$content[$line];
+								}
+						}
+			}
+		$$hash{$name}=$seq;
+	}
+
+
+
+sub Combination()
+	{
+		(my $m,my $n,my $comRef)=@_;
+		my $str="";
+		my %hash;
+		my $fpos;
+		my $num0;
+		my $rest;
+		my $tmp;
+		my $i;
+		my $j;
+		my $key;
+		#my $m=scalar(@$array);
+		my @combination;
+
+		for ($i=1;$i<=$n;$i++) 
+			{
+				$str="1".$str;
+			}
+
+		for ($i=1;$i<=($m-$n);$i++) 
+			{
+				$str=$str."0";
+			}
+
+		$hash{$str}=1;
+		while ($str=~/10/) 
+			{
+				$fpos=index($str,"10");
+				$_=$str;
+				s/10/01/;
+				$str=$_;
+				$tmp=substr($str,0,$fpos);
+				$_=$tmp;
+				s/0//g;
+				$rest=$_;
+				$num0=$fpos-length($_);
+				for ($i=1;$i<=$num0;$i++) 
+					{
+						$rest="$rest"."0";
+					}
+				$str="$rest".substr($str,$fpos,$m-$fpos);
+				$hash{$str}=1;
+			}
+		$j=0;
+		foreach $key (keys %hash) 
+			{
+				$combination[$j]="";
+				for ($i=0;$i<$m;$i++) 
+					{
+						if (substr($key,$i,1) eq "1") 
+							{
+							if ($combination[$j] ne "") 
+								{
+									#$combination[$j]=$combination[$j]."\t".$$array[$i];
+									$combination[$j]=$combination[$j]."\t".$i;
+								}else
+									{
+										#$combination[$j]=$$array[$i];  ### For return species ID
+										$combination[$j]=$i;
+									}
+							}
+					}
+				$j++;
+			}
+		@$comRef=@combination; ### update the data through the physic address
+	}
+
+sub ChkCombinationValue()
+{
+	(my $m,my $n)=@_;
+	my %hash;
+	my %vhash;
+	my $value=0;
+	my $key;
+	my @row;
+	my @sdA;
+	my @sdB;
+
+	### initialization
+	$hash{$m."-".$n}=1;
+
+	### split combination
+	while (scalar(keys %hash)>0 and $value<=$sampleSize)
+	{
+		foreach $key (keys %hash)
+		{
+			if ($value > $sampleSize) ### threshold
+			{
+				last;
+			}
+			if (!exists($hash{$key}))
+			{
+				next;
+			}
+			@row=split(/-/,$key);
+			#print $row[0]."|".$row[1]."\n";
+			if ($row[0] eq $row[1])
+			{
+				$value=$value+$hash{$key};
+			}else
+			{
+				##split
+				$sdA[0]=$row[0]-1;
+				$sdA[1]=$row[1];
+				$sdB[0]=$row[0]-1;
+				$sdB[1]=$row[1]-1;
+				##storing A
+				if (($sdA[0] eq $sdA[1]) or $sdA[1] ==0)
+				{
+					$value=$value+$hash{$key};
+				}else
+				{
+					if (exists($hash{$sdA[0]."-".$sdA[1]}))
+					{
+						$hash{$sdA[0]."-".$sdA[1]}=$hash{$sdA[0]."-".$sdA[1]}+$hash{$key};
+					}else
+					{
+						$hash{$sdA[0]."-".$sdA[1]}=$hash{$key};
+					}
+				}
+
+				##storing B
+				if (($sdB[0] eq $sdB[1]) or $sdB[1]==0)
+				{
+					$value=$value+$hash{$key};
+				}else
+				{
+					if (exists($hash{$sdB[0]."-".$sdB[1]}))
+					{
+						$hash{$sdB[0]."-".$sdB[1]}=$hash{$sdB[0]."-".$sdB[1]}+$hash{$key};
+					}else
+					{
+						$hash{$sdB[0]."-".$sdB[1]}=$hash{$key};
+					}
+				}
+			}
+			#delete original combination
+			delete($hash{$key});
+		}
+	}
+
+	if ($value>$sampleSize)
+	{
+		return 0;
+	}else
+	{
+		return $value;
+	}
+}
+
+
+sub SampleCombination()
+{
+	(my $m,my $n,my $comRef)=@_;
+	my %hash;
+	my $sampleTimes=0;
+	my @randNum;
+	my @sortID;
+	my $i;
+	my $j;
+	my $tmp;
+	while ( scalar(keys %hash)<$sampleSize and $sampleTimes<($sampleSize*2))
+	{
+		for ($i=0;$i<$m;$i++) # generate random data
+		{
+			$randNum[$i]=int(100000 * rand(100));
+			$sortID[$i]=$i;
+		}
+
+		for ($i=0;$i<$m;$i++) # sorting random data
+		{
+			for ($j=0;$j<$m;$j++)
+			{
+				if ($randNum[$sortID[$i]]<$randNum[$sortID[$j]])
+				{
+					$tmp=$sortID[$i];
+					$sortID[$i]=$sortID[$j];
+					$sortID[$j]=$tmp;
+				}
+			}
+		}
+
+		#storing data
+		$tmp=join("\t",sort {$a<=>$b} (splice(@sortID,0,$n)));
+		$hash{$tmp}=1;
+		$sampleTimes++;
+	}
+	@$comRef=keys %hash;
+}
+
+
+sub PanGenomeNumber()
+	{
+		(my $spID)=@_;
+		my $pan=0;    
+		my $core=0;
+		my $count; #### counter;
+		my @row;
+
+		foreach  (@clusters) 
+			{
+				$count=0;
+				@row=split(/\t/,$_);
+
+				foreach  (@$spID) 
+					{
+						$count=$count+$row[$_];
+					}
+				
+				if ($count>0) 
+					{
+						$pan++;
+						if ($count == scalar(@$spID)) 
+							{
+								$core++;
+							}
+					}
+			}
+		return $pan."\t".$core;
+	}
+
+sub fit_model_A()
+	{
+### model y = A * x**B + C
+		(my $xdata,my $ydata)=@_;
+		my $i;
+		my $b;
+		my $max_B=0;
+		my $max_R=0;
+		my $max_A=0;
+		my $max_A_interval;
+		my $max_C=0;
+		my $max_C_interval;
+		my $R=1e-100;
+		my $start;
+		my $end;
+		my $step;
+		my @xValues;
+		my @yValues;
+
+		$start=1;
+		$step=0.001;
+		$b=$start;
+		$max_R=0;
+		$R=1e-100;
+
+		use Statistics::LineFit;
+		use Statistics::Distributions;
+
+		while ($max_R<=$R) 
+			{
+				if (($b < 0.02) and ($b >-0.02)) 
+					{
+						$b=-0.02;
+					}
+
+				for ($i=0;$i<@$xdata;$i++) 
+					{
+						$xValues[$i]=$$xdata[$i]**$b;
+					}
+				@yValues=@$ydata;
+				my $lineFit = Statistics::LineFit->new();
+				$lineFit->setData (\@xValues, \@yValues) or die "Invalid data";
+				(my $intercept, my $slope) = $lineFit->coefficients();
+				my $rSquared = $lineFit->rSquared();
+				my $meanSquaredError = $lineFit->meanSqError();
+				my $durbinWatson = $lineFit->durbinWatson();
+				my $sigma = $lineFit->sigma();
+				(my $tStatIntercept, my $tStatSlope) = $lineFit->tStatistics();
+				(my $varianceIntercept,my $varianceSlope) = $lineFit->varianceOfEstimates();
+
+				$max_R=$R;
+				$R=$rSquared;
+				if ($max_R<=$R) 
+					{
+						$max_R=$R;
+						($max_C,$max_A)=$lineFit->coefficients();
+						$max_A_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceSlope);
+						$max_C_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceIntercept);
+					}
+				$b=$b-$step;
+			}
+		$max_B=$b;
+		return ($max_R,$max_A,$max_A_interval,$max_B,$max_C,$max_C_interval);
+
+	}
+
+sub fit_model_B()
+	{
+### model y = A * exp(x*B) + C
+		(my $xdata,my $ydata)=@_;
+		my $i;
+		my $b;
+		my $max_B=0;
+		my $max_R=0;
+		my $max_A=0;
+		my $max_A_interval;
+		my $max_C=0;
+		my $max_C_interval;
+		my $R=1e-100;
+		my $start;
+		my $end;
+		my $step;
+		my @xValues;
+		my @yValues;
+
+		$start=0;
+		$step=0.001;
+		$b=$start;
+		$max_R=0;
+		$R=1e-100;
+
+		use Statistics::LineFit;
+		use Statistics::Distributions;
+
+		while ($max_R<=$R) 
+			{
+				if (($b < 0.02) and ($b >-0.02)) 
+					{
+						$b=-0.02;
+					}
+
+				for ($i=0;$i<@$xdata;$i++) 
+					{
+						$xValues[$i]=exp($$xdata[$i]*$b);
+					}
+				@yValues=@$ydata;
+				my $lineFit = Statistics::LineFit->new();
+				$lineFit->setData (\@xValues, \@yValues) or die "Invalid data";
+				(my $intercept, my $slope) = $lineFit->coefficients();
+				my $rSquared = $lineFit->rSquared();
+				my $meanSquaredError = $lineFit->meanSqError();
+				my $durbinWatson = $lineFit->durbinWatson();
+				my $sigma = $lineFit->sigma();
+				(my $tStatIntercept, my $tStatSlope) = $lineFit->tStatistics();
+				(my $varianceIntercept,my $varianceSlope) = $lineFit->varianceOfEstimates();
+
+				$max_R=$R;
+				$R=$rSquared;
+				if ($max_R<=$R) 
+					{
+						$max_R=$R;
+						($max_C,$max_A)=$lineFit->coefficients();
+						$max_A_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceSlope);
+						$max_C_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceIntercept);
+					}
+				$b=$b-$step;
+			}
+		$max_B=$b;
+		return ($max_R,$max_A,$max_A_interval,$max_B,$max_C,$max_C_interval);
+	}
+
+
+sub ReadData2Array()
+	{
+		(my $file, my $array1,my $col1,my $array2,my $col2)=@_;
+		my $i=0;
+		open(F,$file);
+		$_=<F>;
+		while (<F>) 
+			{
+				chomp();
+				@_=split(/\t/,$_);
+				$$array1[$i]=$_[$col1];
+				$$array2[$i]=$_[$col2];
+				$i++;
+			}
+		close(F);
+	}
+
+sub SumData()
+	{
+		(my $xdata,my $ydata,my $SumMethod)=@_;
+		my %hash;
+		my $i;
+		my $key;
+		my $max=0;
+		for ($i=0;$i<@$xdata;$i++) 
+			{
+				if (exists($hash{$$xdata[$i]})) 
+					{
+						$hash{$$xdata[$i]}=$hash{$$xdata[$i]}." ".$$ydata[$i];
+					}else
+						{
+							$hash{$$xdata[$i]}=$$ydata[$i];
+							if ($$xdata[$i]>$max) 
+								{
+									$max=$$xdata[$i];
+								}
+						}
+			}
+		@$xdata=qw();
+		@$ydata=qw();
+		$i=0;
+		foreach $i (1..$max) 
+			{
+				$$xdata[$i-1]=$i;
+				if ($SumMethod eq "median") 
+					{
+					$$ydata[$i-1]=&median($hash{$i});
+					}else
+						{
+							$$ydata[$i-1]=&mean($hash{$i});
+						}
+			}
+		#print join(",",@$xdata)."\n";
+		#print join(",",@$ydata)."\n";
+	}
+
+sub median()
+	{
+		(my $data)=@_;
+		my @data=split(/ /,$data);
+		my $arraylen=scalar(@data);
+		@data=sort{$a<=>$b} @data;
+		if (int($arraylen/2)*2 == $arraylen) 
+			{
+				return ($data[$arraylen/2]+$data[$arraylen/2-1])/2;
+			}else
+				{
+					return $data[int($arraylen/2)];
+				}
+	}
+
+sub mean()
+	{
+		(my $data)=@_;
+		my @data=split(/ /,$data);
+		my $sum=0;
+		foreach  (@data) 
+			{
+				$sum=$sum+$_;
+			}
+		return int(($sum/scalar(@data))*1000)/1000;
+	}
+
+sub ReplaceName()
+	{
+		(my $sp,my $file)=@_;
+		my @content;
+		my $line;
+		my $i;
+		my $target;
+		open(F,$file);
+		@content=<F>;
+		close(F);
+		for ($line=0;$line<@content;$line++) 
+			{
+				for ($i=0;$i<@$sp;$i++) 
+					{
+						$_=$content[$line];
+						$target="sp".$i."sp";
+						s/$target/$$sp[$i]/;
+						$content[$line]=$_;
+					}
+			}
+		open(R,">$file");
+		print R @content;
+		close(R);
+	}
+
+sub MP_Result_to_Table()
+	{
+		(my $MPresult, my $outputfile)=@_;
+		my %hash;
+		my $maxid=0;
+		my $i;
+		my @row;
+
+		open(F,"$MPresult");
+		$_=<F>;
+		while (<F>) 
+			{
+				@row=split(/\t/,$_);
+				if (exists($hash{$row[0]})) 
+					{
+						$hash{$row[0]}=$hash{$row[0]}."\t".$row[2];
+					}else
+						{
+							$hash{$row[0]}=$row[2];
+							if ($row[0]>$maxid) 
+								{
+									$maxid=$row[0];
+								}
+						}
+			}
+		close(F);
+
+		open(R,">$outputfile");
+		foreach $i (1..$maxid) 
+			{
+				print R $hash{$i}."\n";
+			}
+		close(R);
+	}
+
+
+
+sub FormatCluster()
+	{
+		(my $infile,my $genelist,my $spnum,my $cluster)=@_;
+		my %hash;
+		my %gene;
+		my $key;
+		my @row;
+		my $sp;
+		my $line;
+		my $i=0;
+		my $j=0;
+		my @content;
+
+### record gene in clusters
+		open(F,"$infile");
+		@content=<F>;
+		close(F);
+		chomp(@content);
+		for ($line=0;$line<@content;$line++) 
+			{
+				@row=split(/\t/,$content[$line]);
+				foreach $key (@row) 
+					{
+						$gene{$key}=1;
+					}
+			}
+###retrieves gene which is not in clutsers
+
+		open(F,"$genelist");
+		while ($key=<F>) 
+			{
+				if ($key=~/^>/) 
+					{
+						chomp($key);
+						$_=$key;
+						s/^>//;
+						$key=$_;
+						if (!exists($gene{$key})) 
+							{
+								$content[$line]=$key;
+								$line++;
+							}
+					}
+			}
+		close(F);
+
+#### initialization @cluster
+		@$cluster=qw();
+		$j=0;
+
+		foreach $line (@content)
+			{
+			if ($line ne "") 
+				{
+					%hash=();
+					@row=split(/\t/,$line);
+					foreach $key (@row) 
+						{
+							$sp=substr($key,0,index($key,"G"));
+							$gene{$key}=1;
+							if (exists($hash{$sp})) 
+								{
+									$hash{$sp}=$hash{$sp}.",".$key;
+								}else
+									{
+										$hash{$sp}=$key;
+									}
+						}
+
+					$i=0;
+					@row=qw();
+					
+					foreach $i (0..($spnum-1)) 
+						{
+							if (exists($hash{"S$i"})) 
+								{
+									$row[$i]=$hash{"S$i"};
+								}else
+									{
+										$row[$i]="-";
+									}
+						}
+					$$cluster[$j++]=join("\t",@row);
+				}
+			}
+	}
+
+