Mercurial > repos > dereeper > pgap

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Thu, 24 Jun 2021 13:51:52 +0000
parents
children 70b7a5270968
comparison
equal deleted inserted replaced
-1:000000000000 0:83e62a1aeeeb
1 #!/usr/bin/perl
2 use strict;
3 use warnings;
4 use Getopt::Long;
5
6 ### programs from BLAST
7 my $formatdb="formatdb";
8 my $blastall="blastall";
9
10 ### programs from mcl
11 my $mcl="mcl";
12
13 ### programs from mafft
14 my $mafft="mafft";
15
16 ### programs from PHYLIP
17 my $seqboot="seqboot";
18 my $neighbor="neighbor";
19 my $consense="consense";
20 my $dnaml="dnaml";
21 my $dnadist="dnadist";
22 my $dnapars="dnapars";
23
24 my $count_tree=0;
25
26 my $sampleSize=8000; # when calculate the pan-genome size, we will sample $sampleSize combinations
27 # if the total combination number is larger than $sampleSize for specific genomes
28 # Surely, the number of $sampleSize is, the larger, the better.
29 # However, the larger the $sampleSize is, the more time would be consumed.
30 # we suggest the range: 5000 ~ 20,000
31
32 #####################################################################
33 # DOn't modify the following code, unless you know their functions
34 #####################################################################
35
36 my %opt=qw();
37 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!");
38
39 my @usage=qq(
40 ====== Pan-Genome Analysis Pipeline (PGAP) ======
41 Version 1.2.1
42
43 Usage: perl PGAP.pl [Options]
44
45 Options:
46 --strains String Input strains nicknames, and join them with '+', for example: A+B+C
47 --input String Input data directory
48 --output String Result output directory
49
50 --cluster Run homologous gene clustering
51 --pangenome Run pan-genome analysis
52 --variation Run homologous clusters variation analysis
53 --evolution Run evolution analysis
54 --function Run Function analysis
55
56 --method String GF for GeneFamily method, and MP for MultiParanoid method
57 for GF: fast, but not very accurate
58 evalue, score, indentity, coverage are employed
59 for MP: slow, but more accurate
60 score, coverage, local, global are employed
61 --thread Int Number of processors to use in blastall. [default:1]
62 --score Int Minimum score in blastall. [default:40]
63 --evalue Decimal Maximal E-value in blastall. [default:1e-10]
64 --coverage Decimal Minimum alignment coverage for two homologous proteins. [default:0.5]
65 --local Decimal Minimum local alignment overlap in MP method. [default:0.25]
66 --global Decimal Minimum global alignment overlap in MP method. [default:0.5]
67 --identity Decimal Minimum alignment indentity for two homologous proteins. [default:0.5]
68 --bootstrap Int Bootstrap times for phylogenetics tree. [default:1]
69
70 --h or help Display this message
71 );
72
73
74 ############# specified variable #############
75 my $inputDIR;
76 my $outputDIR;
77 my $run_cluster;
78 my $run_pangenome;
79 my $run_variation;
80 my $run_evolution;
81 my $run_function;
82 my $method="";
83 my $thread;
84 my $score;
85 my $identity;
86 my $evalue;
87 my $coverage;
88 my $global;
89 my $local;
90 my $bootstrap;
91
92
93 my %pep;
94 my %nuc;
95 my $spnum;
96 my @clusters;
97 my $Cluster;
98 my @SpecieCombination;
99 my @spID;
100 my %genenum;
101 my %aaAln;
102 my %ntAln;
103 my %cog;
104 my %description;
105 #my %aa4tree; ### AA sequence for Phylogenetic Tree
106 my %nt4tree; ### nucleotide sequence for Phylogenetic Tree
107 my @SNPPosition; ### SNP position
108 my $dieMessage="You did not run PGAP.pl in the program directory\n";
109 my $section;
110
111 ######### common temporary variable #############
112 my $i;
113 my $j;
114 my $line;
115 my %tmpHash;
116 my @tmp;
117 my $tmp;
118 my $key;
119 my @row;
120 my $inparacount;
121 my $ClusterID;
122 my $orth;
123 my @content;
124 my $clusterName;
125 my @xdata;
126 my @ydata;
127 my @fit;
128
129 my $fit_A;
130 my $fit_A_interval;
131 my $fit_B;
132 my $fit_C;
133 my $fit_C_interval;
134 my $fit_Rsquare;
135
136
137
138 #### check option
139
140 my $opt_error=0;
141
142 if ((scalar(keys %opt) ==0) or (exists($opt{"help"})))
143 {
144 print join("\n",@usage)."\n";
145 exit;
146 }
147
148
149 ###################### public info
150 ### strains name
151 my @species;
152 if (exists($opt{"strains"}))
153 {
154 @species=split(/\+/,$opt{"strains"});
155 $spnum=scalar(@species);
156 }else
157 {
158 print "Please assign strains nick name!\n";
159 exit;
160 }
161
162 ### input data directory
163
164 if (exists($opt{"input"}))
165 {
166 $inputDIR=$opt{"input"};
167 if ($inputDIR!~/\/$/)
168 {
169 $inputDIR=$inputDIR."/";
170 }
171 }else
172 {
173 print "Please assign input data directory!\n\n";
174 exit;
175 }
176 ### output data directory
177
178 if (exists($opt{"output"}))
179 {
180 $outputDIR=$opt{"output"};
181 if ($outputDIR!~/\/$/)
182 {
183 $outputDIR=$outputDIR."/";
184 }
185 }else
186 {
187 print "Please assign result output directory!\n\n";
188 exit;
189 }
190
191 ###################### section info
192
193 if (exists($opt{"cluster"}))
194 {
195 $run_cluster=1;
196 }else
197 {
198 $run_cluster=0;
199 }
200
201 if (exists($opt{"pangenome"}))
202 {
203 $run_pangenome=1;
204 }else
205 {
206 $run_pangenome=0;
207 }
208
209 if (exists($opt{"variation"}))
210 {
211 $run_variation=1;
212 }else
213 {
214 $run_variation=0;
215 }
216
217 if (exists($opt{"evolution"}))
218 {
219 $run_evolution=1;
220 }else
221 {
222 $run_evolution=0;
223 }
224
225 if (exists($opt{"function"}))
226 {
227 $run_function=1;
228 }else
229 {
230 $run_function=0;
231 }
232
233 if ($run_cluster)
234 {
235 ### method
236 if (exists($opt{"method"}))
237 {
238 $method=uc($opt{"method"});
239 if ($method!~/^GF$/ and $method!~/^MP$/)
240 {
241 print "Unknown method: ".$opt{"method"}."\n";
242 exit;
243 }
244 }else
245 {
246 print "Please assign the cluster method!\n\n";
247 exit;
248 }
249
250 ##thread
251 if (exists($opt{"thread"}))
252 {
253 $thread=$opt{"thread"};
254 if ($thread==0)
255 {
256 print "please assign an applicable thread value.\n";
257 exit;
258 }
259 }else
260 {
261 $thread=1;
262 }
263
264 ##score
265 if (exists($opt{"score"}))
266 {
267 $score=$opt{"score"};
268 if ($score<=0)
269 {
270 print "please assign an applicable score value.\n";
271 exit;
272 }
273 }else
274 {
275 $score=40;
276 }
277
278 if ($method eq "GF")
279 {
280 ###identity
281 if (exists($opt{"identity"}))
282 {
283 $identity=$opt{"identity"};
284 if ($identity>1 or $identity<=0)
285 {
286 print "identity should be 0 ~ 1 \n";
287 exit;
288 }
289 }else
290 {
291 $identity=0.5;
292 }
293
294 ###evalue
295 if (exists($opt{"evalue"}))
296 {
297 $evalue=$opt{"evalue"};
298 }else
299 {
300 $evalue=1e-10;
301 }
302
303 ###coverage
304 if (exists($opt{"coverage"}))
305 {
306 $coverage=$opt{"coverage"};
307 if ($coverage>1 or $coverage<=0)
308 {
309 print "coverage should be 0 ~ 1 \n";
310 exit;
311 }
312 }else
313 {
314 $coverage=0.5;
315 }
316 }
317
318
319 if ($method eq "MP")
320 {
321 ###global
322 if (exists($opt{"global"}))
323 {
324 $global=$opt{"global"};
325 if ($global>1 or $global<=0)
326 {
327 print "global coverage should be 0 ~ 1 \n";
328 exit;
329 }
330 }else
331 {
332 $global=0.5;
333 }
334 ###local
335 if (exists($opt{"local"}))
336 {
337 $local=$opt{"local"};
338 if ($local<=0)
339 {
340 print "local coverage should be 0 ~ [global coverage value] \n";
341 exit;
342 }
343 if ($local>$global)
344 {
345 print "local coverage should be less than global coverage!\n";
346 exit;
347 }
348 }else
349 {
350 $local=0.25;
351 }
352 }
353 }
354
355 if ($run_evolution)
356 {
357 if (exists($opt{"bootstrap"}))
358 {
359 $bootstrap=$opt{"bootstrap"};
360 if ($bootstrap<=0)
361 {
362 print "please assign an applicable bootstrap value.\n";
363 }
364 }else
365 {
366 $bootstrap=1;
367 }
368 }
369
370 print "Program begin at ".localtime()."\n";
371 print "The following are the parameters for current process:\n";
372 print "Strains: ".join(",",@species)."\n";
373 print "Input directory: $inputDIR\n";
374 print "Output directory: $outputDIR\n";
375 if ($run_cluster)
376 {
377 print "Cluster analysis: yes\n";
378 print " Method: $method\n";
379 print " Thread: $thread\n";
380 if ($method eq "GF")
381 {
382 print " E-value: $evalue\n";
383 print " Identity: $identity\n";
384 print " Coverage: $coverage\n";
385 print " Score: $score\n";
386 }
387 if ($method eq "MP")
388 {
389 print " Local: $local\n";
390 print " Global: $global\n";
391 print " Score: $score\n";
392 }
393 }else
394 {
395 print "Cluster analysis: no\n";
396 }
397 if ($run_pangenome)
398 {
399 print "Pan-genome analysis: yes\n";
400 }else
401 {
402 print "Pan-genome analysis: no\n";
403 }
404 if ($run_variation)
405 {
406 print "Variation analysis: yes\n";
407 }else
408 {
409 print "Variation analysis: no\n";
410 }
411 if ($run_evolution)
412 {
413 print "Evolution analysis: yes\n";
414 print " Bootstrap: $bootstrap\n";
415 }else
416 {
417 print "Evolution analysis: no\n";
418 }
419 if ($run_function)
420 {
421 print "Function analysis: yes\n";
422 }else
423 {
424 print "Function analysis: no\n";
425 }
426
427 $section=$run_cluster.$run_pangenome.$run_variation.$run_evolution.$run_function;
428
429 ###############################################
430 # section 0) check input file and program
431 ###############################################
432 if (!(-e $outputDIR))
433 {
434 system("mkdir $outputDIR");
435 }
436 system("chmod +rw $outputDIR");
437
438 if (!(-w $outputDIR))
439 {
440 print "There is no WRITE permission in $outputDIR\n";
441 exit;
442 }
443 @tmp=qw();
444 &CheckInputFile(\@species,$inputDIR,$section,$method,\@tmp);
445 &CheckExtraProgram($section,$method,\@tmp);
446
447 if (scalar(@tmp)>0)
448 {
449 open(R,">".$outputDIR."0.error.message");
450 print R join("",@tmp)."\n";
451 close(R);
452 print "error!\nlog are saved in ${outputDIR}0.error.message\n";
453 exit;
454 }
455
456
457 ############################################
458 # section 1) cluster analysis
459 ############################################
460
461 if ($run_cluster)
462 {
463 print "\n\n############################################\n";
464 print "# section 1) cluster analysis\n";
465 print "############################################\n\n\n";
466
467 #### cluster gene and return result to the array @clusters
468
469 if ($method eq "MP")
470 {
471 print "Begin cluster gene with MP method ...\n";
472 &MP();
473 }else
474 {
475 print "Begin cluster gene with GF method ...\n";
476 &GF();
477 }
478
479 #### output normal cluster format
480
481 &FormatClusterOutPut(\@species,"${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
482
483 #### Retrieve cluster
484
485 &RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
486
487 ##### gene distribution in each strains
488 %tmpHash=();
489 &GeneDistribution(\@clusters,\%tmpHash);
490
491 open(R,">${outputDIR}1.Gene_Distribution_By_Conservation.txt");
492 print R "SharedBy_Strains\t".join("\t",@species)."\n";
493
494 for ($i=$spnum;$i>0;$i--)
495 {
496 print R $i;
497 for ($j=0;$j<$spnum;$j++)
498 {
499 if (exists($tmpHash{$i."|".$j}))
500 {
501 print R "\t".$tmpHash{$i."|".$j};
502 }else
503 {
504 print R "\t0";
505 }
506 }
507 print R "\n";
508 }
509
510 close(R);
511
512
513 }else
514 {
515 print "Homologous gene clustering is skipped!\n";
516 }
517
518
519 if ($run_pangenome)
520 {
521 print "\n\n############################################\n";
522 print "# section 2) Pan-genome analysis\n";
523 print "############################################\n\n\n";
524
525 #### Retrieve cluster
526 &RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
527 chomp(@clusters);
528
529 #### convert file into 0-1 matrix
530 for ($line=0;$line<@clusters;$line++)
531 {
532 @row=split(/\t/,$clusters[$line]);
533 splice(@row,0,1);
534 for ($i=0;$i<@row;$i++)
535 {
536 if ($row[$i] eq "-")
537 {
538 $row[$i]=0;
539 }else
540 {
541 $row[$i]=1;
542 }
543 }
544 $clusters[$line]=join("\t",@row);
545 }
546
547 #### fetch gene number of each strains
548 for ($i=0;$i<$spnum;$i++)
549 {
550 open(F,"$inputDIR$species[$i].pep");
551 @tmp=<F>;
552 close(F);
553 @tmp=grep(/^>/,@tmp);
554 $genenum{$species[$i]}=scalar(@tmp);
555 }
556
557 #### pan genome size and core genome size
558 print "Deducing pan genome size and core genome size for each composition...\n\n";
559
560 open(PAN,">${outputDIR}2.PanGenome.Data.txt");
561 print PAN "ClusterConservation\tTotalGeneNumber\tPanGenome\tCoreGenome\n";
562
563 for ($i=1;$i<=scalar(@species);$i++)
564 {
565 #@SpecieCombination=&Combination(\@species,$i);
566 #@SpecieCombination=&Combination($spnum,$i);
567 if (&ChkCombinationValue($spnum,$i) !=0) ### transfer the array reference to the subroutine
568 {
569 &Combination($spnum,$i,\@SpecieCombination); ## if the combination number is less than sampleSize, then fecth all, else sample
570 }else
571 {
572 &SampleCombination($spnum,$i,\@SpecieCombination);
573 }
574
575 foreach $key (@SpecieCombination)
576 {
577 ##### count total gene number in current combination
578 $tmp=0;
579 @spID=split(/\t/,$key); #### speices id in current combination
580 foreach (@spID)
581 {
582 $tmp=$tmp+$genenum{$species[$_]};
583 }
584 ##### scan pangenome and coregenome
585 @tmp=split(/\t/,&PanGenomeNumber(\@spID));
586 print PAN "$i\t$tmp\t".join("\t",@tmp)."\n";
587
588 }
589 }
590
591 close(PAN);
592
593 #### data fit
594
595 #### for model A
596
597 if ($spnum<3)
598 {
599 print "There are $spnum strains. For pan-genome function fitting, at least 3 strains data are required.\n";
600 }else
601 {
602 open(R,">${outputDIR}2.PanGenome.Profile.txt");
603 ##### genome number & pan-genome size
604 @xdata=qw();
605 @ydata=qw();
606 &ReadData2Array("${outputDIR}2.PanGenome.Data.txt",\@xdata,0,\@ydata,2);
607 &SumData(\@xdata,\@ydata,"mean");
608 ($fit_Rsquare, $fit_A, $fit_A_interval, $fit_B, $fit_C, $fit_C_interval)=&fit_model_A(\@xdata,\@ydata);
609 print R "The relation bewteen genome number and pan-genome size\n\n";
610 print R "Function model: y=A*x**B +C \n";
611 print R "\ty denotes pan-genome size, x denotes genome number, and A, B, C are fitting parameters.\n\n";
612 print R "Fitting result:\n";
613 print R "\ty = $fit_A *x**$fit_B + $fit_C\n";
614 print R "\tR-square = $fit_Rsquare\n";
615 print R "\tA 95% confidence interval: ($fit_A - $fit_A_interval , $fit_A + $fit_A_interval)\n";
616 print R "\tC 95% confidence interval: ($fit_C - $fit_C_interval , $fit_C + $fit_C_interval)\n\n\n\n\n";
617
618 ##### total gene number & pan-genome size
619 #@xdata=qw();
620 #@ydata=qw();
621 #&ReadData2Array("${outputDIR}2.PanGenome.Data.txt",\@xdata,1,\@ydata,2);
622 #&SumDataByMedian(\@xdata,\@ydata);
623 #($fit_Rsquare, $fit_A, $fit_B, $fit_C)=&fit_model_A(\@xdata,\@ydata);
624 #print R "The relation bewteen total gene number and pan-genome size\n\n";
625 #print R "$fit_Rsquare, $fit_A, $fit_B, $fit_C\n";
626 #print R "\ty = $fit_A *x**$fit_B + $fit_C R-square = $fit_Rsquare\n";
627 #print R "\tx: total gene number\n";
628 #print R "\ty: pan-genome size\n\n\n\n\n";
629
630 ##### genome number & core genome
631 @xdata=qw();
632 @ydata=qw();
633 &ReadData2Array("${outputDIR}2.PanGenome.Data.txt",\@xdata,0,\@ydata,3);
634 &SumData(\@xdata,\@ydata,"mean");
635 ($fit_Rsquare, $fit_A, $fit_A_interval, $fit_B, $fit_C, $fit_C_interval)=&fit_model_B(\@xdata,\@ydata);
636 print R "The relation bewteen genome number and core genome size\n\n";
637 print R "Function model: y=A*exp(B*x) +C \n";
638 print R "\ty denotes pan-genome size, x denotes genome number, and A, B, C are fitting parameters.\n\n";
639 print R "Fitting result:\n";
640 print R "\ty = $fit_A *exp($fit_B * x) + $fit_C R-square = $fit_Rsquare\n";
641 print R "\tR-square = $fit_Rsquare\n";
642 print R "\tA 95% confidence interval: ($fit_A - $fit_A_interval , $fit_A + $fit_A_interval)\n";
643 print R "\tC 95% confidence interval: ($fit_C - $fit_C_interval , $fit_C + $fit_C_interval)\n\n\n\n\n";
644 close(R);
645 }
646
647 }
648
649
650 ############################################
651 # section 3) CDS variation analysis
652 ############################################
653
654 if ($run_variation)
655 {
656 print "\n\n############################################\n";
657 print "# section 3) CDS variation analysis\n";
658 print "############################################\n\n\n";
659
660 #### Retrieve cluster
661 &RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
662 chomp(@clusters);
663
664 ## protein
665 system("rm -rf *.pep");
666 &PrepareFasta(\@species,$inputDIR,".pep"); ###prepare pep file
667 system("cat *.pep > All.faa && rm -rf *.pep && mv All.faa All.pep");
668 &ReadSequenceInToHash("All.pep",\%pep);
669 ## nucleic
670 system("rm -rf *.nuc");
671 &PrepareFasta(\@species,$inputDIR,".nuc"); ###prepare nuc file
672 system("cat *.nuc > All.ffn && rm -rf *.nuc && mv All.ffn All.nuc");
673 &ReadSequenceInToHash("All.nuc",\%nuc);
674
675 ## scanning SNP
676 %nt4tree=();
677 for ($i=0;$i<$spnum;$i++)
678 {
679 $nt4tree{"S".$i}="";
680 }
681
682 open(VAR,">${outputDIR}3.CDS.variation.txt");
683 print VAR "ClusterID\tStrains_Number\tGene_Number\tPosition\taaType\tntType\tntProfile\tVariation type\n";
684
685 open(VA,">${outputDIR}3.CDS.variation.analysis.txt");
686 print VA "ClusterID\tInDel Base\tNonsynonymous mutation\tSynonymous mutation\n";
687
688 for ($line=0;$line<@clusters;$line++)
689 {
690 @row=split(/\t|\,/,$clusters[$line]);
691 $ClusterID=$row[0];
692 splice(@row,0,1);
693 @row=grep(/^S/,@row);
694 if (scalar(@row) >=2)
695 {
696 open(PEP,">$ClusterID.pep");
697 open(NUC,">$ClusterID.nuc");
698 foreach $key (@row)
699 {
700 print PEP ">$key\n$pep{$key}\n";
701 print NUC ">$key\n$nuc{$key}\n";
702 }
703 close(PEP);
704 close(NUC);
705 system("$mafft --quiet $ClusterID.pep > $ClusterID.pal");
706 #system("perl ./pal2nal.pl $ClusterID.pal $ClusterID.nuc -output fasta > $ClusterID.nal");
707 $tmp=&pal2nal("$ClusterID.pal","$ClusterID.nuc","$ClusterID.nal");
708 if ($tmp == 0)
709 {
710 system("rm -rf $ClusterID.*");
711 next;
712 }
713
714 @tmp=&DetectSNP();
715 if (scalar(@tmp)>0)
716 {
717 print VA $ClusterID."\t".&VarAnalysis(\@tmp)."\n";
718 print VAR join("",@tmp);
719 ### core orthologs
720 @row=split(/\t/,$clusters[$line]);
721 splice(@row,0,1);
722 if ((&CountGeneInCluster(join("\t",@row)) ==$spnum) and (&CountSpeicesInCluster(join("\t",@row)) == $spnum) )
723 {
724 $count_tree++;
725 %tmpHash=();
726 foreach (@row)
727 {
728 $tmpHash{$_}="";
729 }
730 &RemoveHeadGap("$ClusterID.nal",\%tmpHash);
731 &ExtractSNP4tree(\%tmpHash,\%nt4tree);
732 }
733 }
734 system("rm -rf $ClusterID.*");
735 }
736 }
737 close(VAR);
738 close(VA);
739
740 open(R,">${outputDIR}3.CDS.variation.for.evolution.txt");
741 foreach $key (keys %nt4tree)
742 {
743 $_=$key;
744 s/s//gi;
745 print R ">$species[$_]\n$nt4tree{$key}\n";
746 }
747 close(R);
748 print $count_tree."\n\n";
749
750 ###
751 system("rm All.nuc All.pep");
752 }else
753 {
754 print "CDS variation is skipped.\n";
755 }
756
757 ############################################
758 # section 4) CDS variation analysis
759 ############################################
760
761 if ($run_evolution)
762 {
763 #### Retrieve cluster
764 &RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
765 chomp(@clusters);
766
767
768
769 ##################
770 ##
771 ## Distance based
772 ##
773 ##################
774
775 ### caculate the distance between each two strains by cluster
776 &ClusterProfil4Specie(\%tmpHash); # caculate Clusters profile for each specie
777 &DistanceMatrix($spnum,\%tmpHash); # caculate distance matrix acoording to Clusters profile
778 ###output distance
779 open(DIST,">${outputDIR}4.Species_Distance_Clusters_Based.txt");
780 ###header
781 printf DIST "%5d", $spnum;
782 print DIST "\n";
783 foreach $i (0..($spnum-1))
784 {
785 $key="sp".$i."sp";
786 printf DIST "%-10s",$key;
787 foreach $j (0..($spnum-1))
788 {
789 printf DIST " %8f",$tmpHash{$i."-".$j};
790 }
791 print DIST "\n";
792 }
793 close(DIST);
794
795 %tmpHash=();
796
797 ### based on pan genome (distance)
798 print "\nDraw pangenome based phylogenetic tree ...\n\n";
799
800 &PanBasedTree("${outputDIR}4.Species_Distance_Clusters_Based.txt","${outputDIR}4.PanBased");
801
802 ##################
803 ##
804 ## SNP based
805 ##
806 ##################
807 %tmpHash=();
808 if (!(-e "${outputDIR}3.CDS.variation.for.evolution.txt"))
809 {
810 print "Variation in core orthologs cluster is not found from ${outputDIR}3.CDS.variation.for.evolution.txt.\n";
811 print "Maybe you have skipped CDS variation analysis.\n";
812 }else
813 {
814 &ReadSequenceInToHash("${outputDIR}3.CDS.variation.for.evolution.txt",\%tmpHash);
815 open(R,">mlst.aln");
816 for ($i=0;$i<@species;$i++)
817 {
818 print R ">sp${i}sp\n".$tmpHash{$species[$i]}."\n";
819 }
820 close(R);
821 &fasta2phylip("mlst.aln","mlst.phylip");
822 system("rm -rf mlst.aln");
823
824 print "\nDraw SNP based phylogenetic tree ...\n\n";
825 &SNPBasedTree("mlst.phylip","${outputDIR}4.SNPBased");
826 system("rm -rf mlst.phylip")
827 }
828
829 ########
830 # replace speices name
831 ########
832
833 opendir(DIR,"${outputDIR}");
834 @tmp=readdir(DIR);
835 closedir(DIR);
836 @tmp=grep(/^4/,@tmp);
837 foreach $tmp (@tmp)
838 {
839 &ReplaceName(\@species,"${outputDIR}$tmp");
840 }
841 }else
842 {
843 print "Evolution analysis is skipped.\n";
844 }
845
846
847 ############################################
848 # section 4) Function analysis
849 ############################################
850
851 if ($run_function)
852 {
853 #### Retrieve cluster
854 &RetrieveClusterFromFile("${outputDIR}1.Orthologs_Cluster.txt",\@clusters);
855 chomp(@clusters);
856
857 #### prepare annotation file
858 &PrepareTable(\@species,$inputDIR,".function"); ###prepare location file
859 &ReadAnnotation(\@species,\%cog,\%description);
860
861
862 #### assign function
863 open(R,">${outputDIR}5.Orthologs_Cluster_Function.txt");
864 print R "ClusterID\tConservation_Level\tCOG\tDescription\n";
865 for ($i=0;$i<@clusters;$i++)
866 {
867 @row=split(/\t/,$clusters[$i]);
868 $ClusterID=$row[0];
869 splice(@row,0,1);
870 print R $ClusterID."\t".&CountSpeicesInCluster(join("\t",@row))."\t".&getCOG(\@row,\%cog)."\t".&getDescription(\@row,\%description)."\n";
871 }
872 close(R);
873
874 #### COG distribution
875
876 ###Whole Clusters COG Distribution
877 &outputCOGStatistic("${outputDIR}5.Orthologs_Whole_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",$spnum,1));
878
879 ###Core Clusters COG Distribution
880 &outputCOGStatistic("${outputDIR}5.Orthologs_Core_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",$spnum,$spnum));
881
882 ###Dispensable Clusters COG Distribution
883 &outputCOGStatistic("${outputDIR}5.Orthologs_Dispensable_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",($spnum-1),2));
884
885 ###strains specifc Clusters COG Distribution
886 &outputCOGStatistic("${outputDIR}5.Orthologs_specifc_Cluster_COG_Distribution.txt",&scanCOG("${outputDIR}5.Orthologs_Cluster_Function.txt",1,1));
887
888 system("rm -rf *.function");
889
890 }else
891 {
892 print "Function analysis is skipped.\n";
893 }
894
895 sub outputCOGStatistic()
896 {
897 (my $file,my $subcogcount)=@_;
898 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 -");
899 my @cogdesc=("INFORMATION STORAGE AND PROCESSING","CELLULAR PROCESSES AND SIGNALING","METABOLISM","POORLY CHARACTERIZED");
900 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 -);
901 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");
902 my %subcogdesc;
903 my $key;
904 my @cog;
905 my $i;
906 my $cognum;
907
908 for ($i=0;$i<@subcogcat;$i++)
909 {
910 $subcogdesc{$subcogcat[$i]}=$subcogdesc[$i];
911 }
912
913 open(R,">$file");
914 for ($i=0;$i<@cogcat;$i++)
915 {
916 $cognum=0;
917 foreach $key (split(" ",$cogcat[$i]))
918 {
919 $cognum=$cognum+$$subcogcount{$key};
920 }
921 print R $cogdesc[$i]." ( ".$cognum." )\n";
922 foreach $key (split(" ",$cogcat[$i]))
923 {
924 printf R "%-6d %s\n",$$subcogcount{$key},$subcogdesc{$key};
925 }
926 print R "\n";
927
928 }
929 close(R);
930
931 }
932
933 sub scanCOG()
934 {
935 (my $file,my $max_orth,my $min_orth)=@_;
936 my @row;
937 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 -);
938 my %subcogcount;
939 my $cog;
940 my $key;
941
942 foreach $key (@subcogcat)
943 {
944 $subcogcount{$key}=0;
945 }
946
947 @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);
948
949 open(F,"$file");
950 $_=<F>;
951 while (<F>)
952 {
953 @row=split(/\t/,$_);
954 if ($row[1]>=$min_orth and $row[1]<=$max_orth)
955 {
956 if ($row[2] eq "-")
957 {
958 $subcogcount{"-"}++;
959 }else
960 {
961 $_=uc($row[2]);
962 s/COG//gi;
963 $cog=$_;
964 foreach $key (@subcogcat)
965 {
966 if ($cog=~/$key/)
967 {
968 $subcogcount{$key}++;
969 }
970 }
971 }
972 }
973 }
974 close(F);
975
976 return \%subcogcount;
977 }
978
979
980
981
982 sub getCOG()
983 {
984 (my $data,my $coghash)=@_;
985 my $cog="";
986 my @cog;
987 my $key;
988 my %hash;
989 my @gene=split(/\t|\,/,join("\t",@$data));
990 @gene=grep(/^S/,@gene);
991
992 foreach $key (@gene)
993 {
994 if (($$coghash{$key} ne "-") and ($$coghash{$key} ne ""))
995 {
996 $cog=$cog.",".$$coghash{$key};
997 }
998 }
999 @cog=split(/,/,$cog);
1000 foreach $cog (@cog)
1001 {
1002 if ($cog ne "")
1003 {
1004 $hash{$cog}=1;
1005 }
1006 }
1007
1008 $cog=join(",",(keys %hash));
1009 if ($cog eq "")
1010 {
1011 $cog="-";
1012 }
1013 return $cog;
1014 }
1015
1016 sub getDescription()
1017 {
1018 (my $data,my $deschash)=@_;
1019 my $desc="";
1020 my $key;
1021 my @gene=split(/\t|\,/,join("\t",@$data));
1022 @gene=grep(/^S/,@gene);
1023
1024 foreach $key (@gene)
1025 {
1026 if ( ($$deschash{$key} ne "") and ($$deschash{$key} ne "-") and ($$deschash{$key}!~/hypothetical/))
1027 {
1028 $desc=$$deschash{$key};
1029 }
1030 }
1031
1032 if ($desc eq "")
1033 {
1034 $desc="hypothetical protein";
1035 }
1036
1037 return $desc;
1038 }
1039
1040
1041
1042 sub ReadAnnotation()
1043 {
1044 (my $species,my $cog,my $description)=@_;
1045 my $i;
1046 my @row;
1047
1048 for ($i=0;$i<@$species;$i++)
1049 {
1050 open(F,"$$species[$i].function");
1051 while (<F>)
1052 {
1053 chomp($_);
1054 @row=split(/\t/,$_);
1055 if (scalar(@row)>=2)
1056 {
1057 $$cog{$row[0]}=$row[1];
1058 }else
1059 {
1060 $$cog{$row[0]}="-";
1061 }
1062
1063 if (scalar(@row)>=3)
1064 {
1065 $$description{$row[0]}=$row[2];
1066 }else
1067 {
1068 $$description{$row[0]}="hypothetical protein";
1069 }
1070 }
1071 close(F);
1072 }
1073 }
1074
1075
1076 sub SNPBasedTree()
1077 {
1078 (my $infile,my $outfileprefix)=@_;
1079 my $tmpin=$infile;
1080 my $tmpout;
1081
1082
1083 #### boootstrap
1084 print "\n#### seqboot ...\n\n";
1085 open(R,">seqboot.cmd");
1086 #print R "$tmpin\n";
1087 print R "R\n";
1088 print R "$bootstrap\n";
1089 print R "Y\n";
1090 print R "1\n";
1091 close(R);
1092
1093 system("cp $tmpin infile");
1094 system("$seqboot < seqboot.cmd");
1095 system("mv outfile 100dnaseq");
1096 system("rm -rf infile");
1097 system("rm seqboot.cmd"); # 100dnasseq
1098
1099 #### dnaml
1100 print "\n#### dnaml ...\n\n";
1101 open(R,">dnaml.cmd");
1102 #print R "100dnaseq\n";
1103 print R "T\n";
1104 print R "25\n";
1105 if ($bootstrap>1)
1106 {
1107 print R "M\n";
1108 print R "D\n";
1109 #print R "100\n";
1110 print R "$bootstrap\n";
1111 print R "1\n"; # Random number seed (must be odd)?
1112 print R "5\n"; # Number of times to jumble?
1113 }
1114 print R "Y\n";
1115 close(R);
1116
1117 system("cp 100dnaseq infile");
1118 system("$dnaml < dnaml.cmd");
1119 system("rm -rf outfile");
1120 system("rm -rf infile");
1121 system("mv outtree 100dnaseqtree"); # 100dnaseq, 100dnaseqtree
1122
1123 #### consense
1124 print "\n#### dnaml consense ...\n\n";
1125
1126 open(R,">consense.cmd");
1127 #print R "100dnaseqtree\n";
1128 print R "Y\n";
1129 close(R);
1130
1131 system("cp 100dnaseqtree intree");
1132 system("$consense < consense.cmd");
1133 system("mv outfile ${outfileprefix}.ML.outfile");
1134 system("mv outtree ${outfileprefix}.ML.tree");
1135 system("rm -rf infile");
1136 system("rm -rf 100dnaseqtree"); # 100dnaseq
1137
1138 #### dnadist
1139 print "\n#### dnadist ...\n\n";
1140 open(R,">dnadist.cmd");
1141 #print R "100dnaseq\n";
1142 print R "T\n";
1143 print R "25\n";
1144 if ($bootstrap>1)
1145 {
1146 print R "M\n";
1147 print R "D\n";
1148 #print R "100\n";
1149 print R "$bootstrap\n";
1150 }
1151 print R "Y\n";
1152 close(R);
1153
1154 system("cp 100dnaseq infile");
1155 system("$dnadist < dnadist.cmd");
1156 system("rm -rf 100dnaseq");
1157 system("rm -rf infile");
1158 system("mv outfile 100dnadist"); # 100dnadist
1159
1160 #### Neighbor-joining tree
1161 print "\n#### Neighbor-joining ...\n\n";
1162 open(R,">NJ.cmd");
1163 if ($bootstrap>1)
1164 {
1165 #print R "100dnadist\n";
1166 print R "M\n";
1167 #print R "100\n";
1168 print R "$bootstrap\n";
1169 print R "1\n";
1170 }
1171 print R "Y\n";
1172 close(R);
1173
1174 system("cp 100dnadist infile");
1175 system("$neighbor < NJ.cmd");
1176 system("mv outtree 100dnadistNJtree");
1177 system("rm outfile");
1178 system("rm -rf infile");
1179 system("rm -rf NJ.cmd"); # 100dnadist,100dnadistNJtree
1180
1181 #### NJ-consense
1182 print "\n#### NJ-consense ...\n\n";
1183 open(R,">NJ-consense.cmd");
1184 #print R "100dnadistNJtree\n";
1185 print R "Y\n";
1186 close(R);
1187
1188 system("cp 100dnadistNJtree intree");
1189 system("$consense < NJ-consense.cmd");
1190 system("mv outfile ${outfileprefix}.Neighbor-joining.outfile");
1191 system("mv outtree ${outfileprefix}.Neighbor-joining.tree");
1192 system("rm -rf NJ-consense.cmd");
1193 system("rm -rf intree");
1194 system("rm -rf 100dnadistNJtree");
1195
1196
1197 #### UPGMA tree
1198 print "\n#### UPGMA ...\n\n";
1199 open(R,">UPGMA.cmd");
1200 #print R "100dnadist\n";
1201 print R "N\n";
1202 if ($bootstrap>1)
1203 {
1204 print R "M\n";
1205 #print R "100\n";
1206 print R "$bootstrap\n";
1207 print R "1\n";
1208 }
1209 print R "Y\n";
1210 close(R);
1211
1212 system("cp 100dnadist infile");
1213 system("$neighbor < UPGMA.cmd");
1214 system("mv outtree 100dnadistUPGMAtree");
1215 system("rm -rf outfile");
1216 system("rm -rf infile");
1217 system("rm -rf UPGMA.cmd");
1218
1219 #### UPGMA-consense
1220 print "\n#### UPGMA-consense ...\n\n";
1221 open(R,">UPGMA-consense.cmd");
1222 #print R "100dnadistUPGMAtree\n";
1223 print R "Y\n";
1224 close(R);
1225
1226 system("cp 100dnadistUPGMAtree intree");
1227 system("$consense < UPGMA-consense.cmd");
1228 system("mv outfile ${outfileprefix}.UPGMA.outfile");
1229 system("mv outtree ${outfileprefix}.UPGMA.tree");
1230 system("rm -rf UPGMA-consense.cmd");
1231 system("rm -rf 100dnadistUPGMAtree");
1232 system("rm -rf intree");
1233
1234 ###CLEAN TMP FILE
1235
1236 system("rm -rf *.cmd");
1237 system("rm -rf 100dnadist");
1238 }
1239
1240 sub PanBasedTree()
1241 {
1242 (my $infile,my $outfileprefix)=@_;
1243 my $tmpin;
1244 my $tmpout;
1245
1246 $tmpin=$infile;
1247 #### Neighbor-joining tree
1248
1249 open(R,">NJ.cmd");
1250 #print R "$tmpin\n";
1251 print R "Y\n";
1252 close(R);
1253
1254 system("cp $tmpin infile");
1255 system("$neighbor < NJ.cmd");
1256 system("mv outfile ${outfileprefix}.Neighbor-joining.outfile");
1257 system("mv outtree ${outfileprefix}.Neighbor-joining.tree");
1258 system("rm -rf NJ.cmd");
1259 system("rm -rf infile");
1260
1261 #### UPGMA tree
1262
1263 open(R,">UPGMA.cmd");
1264 #print R "$tmpin\n";
1265 print R "N\n";
1266 print R "Y\n";
1267 close(R);
1268
1269 system("cp $tmpin infile");
1270 system("$neighbor < UPGMA.cmd");
1271 system("mv outfile ${outfileprefix}.UPGMA.outfile");
1272 system("mv outtree ${outfileprefix}.UPGMA.tree");
1273 system("rm -rf UPGMA.cmd");
1274 system("rm -rf infile");
1275
1276
1277 ###CLEAN TMP FILE
1278 system("rm -rf *.cmd");
1279 }
1280
1281 sub DistanceMatrix()
1282 {
1283 (my $spnum,my $hash)=@_;
1284 my $i;
1285 my $j;
1286 my $k;
1287 my $dist;
1288 my $ref;
1289 my $query;
1290 foreach $i (0..($spnum-1))
1291 {
1292 foreach $j ($i..($spnum-1))
1293 {
1294 $ref=$$hash{$i};
1295 $query=$$hash{$j};
1296 $dist=0;
1297 for ($k=0;$k<length($ref);$k++)
1298 {
1299 if (substr($ref,$k,1) ne substr($query,$k,1))
1300 {
1301 $dist++;
1302 }
1303 }
1304 $$hash{$i."-".$j}=$dist;
1305 $$hash{$j."-".$i}=$dist;
1306 }
1307 }
1308 }
1309
1310
1311 sub ClusterProfil4Specie
1312 {
1313 (my $hash)=@_;
1314 my @row;
1315 my $i;
1316
1317 foreach (0..($spnum-1)) #initialization Hash
1318 {
1319 $$hash{$_}="";
1320 }
1321
1322 foreach (@clusters)
1323 {
1324 @row=split(/\t/,$_);
1325 splice(@row,0,1);
1326 if (&CountSpeicesInCluster(join("\t",@row))>1)
1327 {
1328 for ($i=0;$i<@row;$i++)
1329 {
1330 if ($row[$i] eq "-")
1331 {
1332 $$hash{$i}=$$hash{$i}."0";
1333 }else
1334 {
1335 $$hash{$i}=$$hash{$i}."1";
1336 }
1337 }
1338 }
1339 }
1340 }
1341
1342 # &ExtractSNP4tree(\%tmpHash,\%nt4tree);
1343
1344 sub ExtractSNP4tree()
1345 {
1346 (my $hash,my $nt4treeRef)=@_;
1347 my $key;
1348 my @row;
1349 my $i;
1350 my $len;
1351 my @tribases;
1352 foreach $key (keys %$hash)
1353 {
1354 $$hash{substr($key,0,index($key,"G"))}=$$hash{$key};
1355 delete($$hash{$key});
1356 }
1357
1358 for ($i=0;$i<$spnum;$i++)
1359 {
1360 $nt4tree{"S".$i}=$nt4tree{"S".$i}.$$hash{"S".$i};
1361 }
1362 }
1363
1364
1365 =pod
1366 sub ExtractSNP4tree()
1367 {
1368 (my $hash,my $nt4treeRef)=@_;
1369 my $key;
1370 my @row;
1371 my $i;
1372 my $len;
1373 my @tribases;
1374 foreach $key (keys %$hash)
1375 {
1376 $$hash{substr($key,0,index($key,"G"))}=$$hash{$key};
1377 delete($$hash{$key});
1378 }
1379 @_=(keys %$hash);
1380 $len=length($_[0]);
1381 for ($j=0;3*$j<$len;$j++)
1382 {
1383 ##### scanning each codon
1384 for ($i=0;$i<$spnum;$i++)
1385 {
1386 $tribases[$i]=substr($$hash{"S".$i},3*$j,3);
1387 }
1388 ##### checking each codon
1389 if (&IsTheSame(@tribases) ==0)
1390 {
1391 for ($i=0;$i<@tribases;$i++)
1392 {
1393 $nt4tree{"S".$i}=$nt4tree{"S".$i}.$tribases[$i];
1394 }
1395 }
1396 }
1397 }
1398 =cut
1399
1400
1401 sub pal2nal()
1402 {
1403 (my $pal,my $nuc, my $nal)=@_;
1404 my %aaAln=();
1405 my %ffn=();
1406 my %ntAln=();
1407 my %nt;
1408 my $dna;
1409 my $nt;
1410 my $key;
1411 my $flag=1;
1412 my $i=0;
1413 my $j;
1414
1415 ### read protein aligment result
1416 &ReadAlignmentToHash("$pal",\%aaAln);
1417 ### read nt sequences
1418 &ReadSequenceInToHash("$nuc",\%ffn);
1419 foreach $key (keys %ffn)
1420 {
1421 $dna=$ffn{$key};
1422 #if (int(length($nt{$key})/3)*3 ne length($nt{$key}))
1423 if (int(length($dna)/3)*3 ne length($dna))
1424 {
1425 $flag=0;
1426 print "The length of nucleotide sequence is not 3 integer times.\n";
1427 last;
1428 }else
1429 {
1430 for ($i=0;$i<(length($dna)/3);$i++)
1431 {
1432 $nt{$key."|".$i}=substr($dna,$i*3,3);
1433 }
1434 }
1435 }
1436
1437 if ($flag==0)
1438 {
1439 return 0;
1440 }else
1441 {
1442 foreach $key (keys %aaAln) ### replace aa with corresponding nt
1443 {
1444 $nt="";
1445 $i=0;
1446 for ($j=0;$j<length($aaAln{$key});$j++)
1447 {
1448 if (substr($aaAln{$key},$j,1) eq "-")
1449 {
1450 $nt=$nt."---";
1451 }else
1452 {
1453 $nt=$nt.$nt{$key."|".$i};
1454 $i++;
1455 }
1456 }
1457 $ntAln{$key}=$nt;
1458 }
1459
1460 ### output
1461 open(R,">$nal");
1462 foreach (keys %ntAln)
1463 {
1464 print R ">$_\n".$ntAln{$_}."\n";
1465 }
1466 close(R);
1467
1468 return 1;
1469 }
1470 }
1471
1472
1473
1474
1475 sub DetectSNP()
1476 {
1477 my %faa;
1478 my %ffn;
1479 my @row;
1480 my $count_gene;
1481 my $count_sp;
1482 my @genelist;
1483 my $i;
1484 my $j;
1485 my $pepalnlen;
1486 my @cdsvar=qw();
1487 my $cdi=0;
1488 my @tribases;
1489 my @bases;
1490 my @aa;
1491
1492
1493 ### fetch gene list
1494 open(F,"$ClusterID.pep");
1495 @genelist=<F>;
1496 close(F);
1497 @genelist=grep(/^>/,@genelist);
1498 chomp(@genelist);
1499 $_=join("\t",@genelist);
1500 s/>//g;
1501 @genelist=split(/\t/,$_);
1502
1503 ### count gene number and species number
1504 @row=split(/\t/,$clusters[$ClusterID-1]);
1505 splice(@row,0,1);
1506 $count_sp=&CountSpeicesInCluster(join("\t",@row));
1507 $count_gene=&CountGeneInCluster(join("\t",@row));
1508
1509 ### read alignment sequences
1510 &ReadAlignmentToHash("$ClusterID.pal",\%faa);
1511 &ReadAlignmentToHash("$ClusterID.nal",\%ffn);
1512
1513 @_=(keys %faa);
1514 $pepalnlen=length($faa{$_[0]});
1515 ### scan SNP
1516 for ($i=1;$i<=$pepalnlen;$i++)
1517 {
1518 @tmp=qw();
1519 @tribases=qw();
1520 for ($j=0;$j<@genelist;$j++) ### fetch triplet codon
1521 {
1522 $tribases[$j]=substr($ffn{$genelist[$j]},3*($i-1),3);
1523 }
1524 if (&IsTheSame(@tribases) ==0) ### if triplet codon is not consistent
1525 {
1526 @aa=qw();
1527 for ($j=0;$j<@genelist;$j++)
1528 {
1529 $aa[$j]=substr($faa{$genelist[$j]},($i-1),1);
1530 }
1531 if (&IsTheSame(@aa) ==0) ### aa is not consistent
1532 {
1533 if (join("",@aa) =~/-/)
1534 {
1535 $cdsvar[$cdi++]=$ClusterID."\t".$count_sp."\t".$count_gene."\t".$i."\t".&CharType(\@aa)."\t-\t-\tInDel\n";
1536 }else
1537 {
1538 #### base 1
1539 for ($j=0;$j<@genelist;$j++)
1540 {
1541 $bases[$j]=substr($ffn{$genelist[$j]},3*($i-1),1);
1542 }
1543 if (&IsTheSame(@bases) ==0)
1544 {
1545 $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";
1546 }
1547 #### base 2
1548 for ($j=0;$j<@genelist;$j++)
1549 {
1550 $bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+1,1);
1551 }
1552 if (&IsTheSame(@bases) ==0)
1553 {
1554 $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";
1555 }
1556 #### base 3
1557 for ($j=0;$j<@genelist;$j++)
1558 {
1559 $bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+2,1);
1560 }
1561 if (&IsTheSame(@bases) ==0)
1562 {
1563 $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";
1564 }
1565 }
1566 }else
1567 {
1568 #### base 1
1569 for ($j=0;$j<@genelist;$j++)
1570 {
1571 $bases[$j]=substr($ffn{$genelist[$j]},3*($i-1),1);
1572 }
1573 if (&IsTheSame(@bases) ==0)
1574 {
1575 $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";
1576 }
1577 #### base 2
1578 for ($j=0;$j<@genelist;$j++)
1579 {
1580 $bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+1,1);
1581 }
1582 if (&IsTheSame(@bases) ==0)
1583 {
1584 $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";
1585 }
1586 #### base 3
1587 for ($j=0;$j<@genelist;$j++)
1588 {
1589 $bases[$j]=substr($ffn{$genelist[$j]},3*($i-1)+2,1);
1590 }
1591 if (&IsTheSame(@bases) ==0)
1592 {
1593 $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";
1594 }
1595 }
1596 }
1597 }
1598 return @cdsvar;
1599 }
1600
1601
1602 sub VarAnalysis()
1603 {
1604 (my $data)=@_;
1605 my @data=@$data;
1606 my $indel=0;
1607 my $syn=0;
1608 my $nonsyn=0;
1609 my @tmp;
1610 $indel=scalar(grep(/InDel$/,@data));
1611 $nonsyn=scalar(grep(/Nonsynonymous mutation$/,@data));;
1612 $syn=scalar(grep(/Synonymous mutation$/,@data));
1613 return "$indel\t$nonsyn\t$syn";
1614 }
1615
1616
1617
1618 sub CharType()
1619 {
1620 (my $str)=@_;
1621 my %hash;
1622 my @data=@$str;
1623 foreach (@data)
1624 {
1625 $hash{$_}=1;
1626 }
1627 return join(",",(keys %hash));
1628 }
1629
1630 sub IsTheSame()
1631 {
1632 (my @data)=@_;
1633 my %hash;
1634 foreach (@data)
1635 {
1636 $hash{$_}=1;
1637 }
1638 if (scalar(keys %hash) ==1)
1639 {
1640 return 1;
1641 }else
1642 {
1643 return 0;
1644 }
1645 }
1646
1647
1648
1649 sub FormatClusterOutPut()
1650 {
1651 (my $speices,my $file,my $cluster)=@_;
1652 my @row;
1653 my $gid=1;
1654 my $key;
1655 my %hash;
1656 my $gene;
1657 my @tmp;
1658 my $i;
1659 my $j;
1660 open(R,">$file");
1661 print R "ClutserID\t".join("\t",@$speices)."\n";
1662 foreach $key (@$cluster)
1663 {
1664 @row=split(/\t/,$key);
1665 for ($i=0;$i<@row;$i++)
1666 {
1667 if ($row[$i] ne "-")
1668 {
1669 @tmp=split(/,/,$row[$i]);
1670 for ($j=0;$j<@tmp;$j++)
1671 {
1672 $_=$tmp[$j];
1673 s/^S[0-9]+G//;
1674 $tmp[$j]=$_;
1675 }
1676 $row[$i]=join(",",@tmp);
1677 }
1678 }
1679 print R $gid."\t".join("\t",@row)."\n";
1680 $gid++;
1681 }
1682 close(R);
1683 }
1684 sub RetrieveClusterFromFile()
1685 {
1686 (my $file,my $clusters)=@_;
1687 my @content;
1688 my @row;
1689 my $spid;
1690 my $line=0;
1691 my $i=0;
1692 my $j;
1693 my @tmp;
1694 open(F,$file) or die "Could open $file\n";
1695 @content=<F>;
1696 close(F);
1697 splice(@content,0,1);
1698 chomp(@content);
1699 foreach (@content)
1700 {
1701 @row=split(/\t/,$_);
1702 $$clusters[$line]=$row[0];
1703 splice(@row,0,1);
1704 for ($i=0;$i<@row;$i++)
1705 {
1706 if ($row[$i] ne "-")
1707 {
1708 @tmp=split(/,/,$row[$i]);
1709 for ($j=0;$j<@tmp;$j++)
1710 {
1711 $tmp[$j]="S${i}G".$tmp[$j];
1712 }
1713 $row[$i]=join(",",@tmp);
1714 }
1715 }
1716 $$clusters[$line]=$$clusters[$line]."\t".join("\t",@row)."\n";
1717 $line++;
1718 }
1719 }
1720
1721
1722
1723
1724 sub GeneDistribution()
1725 {
1726 (my $clusters,my $hash)=@_;
1727 my @row;
1728 my $spid;
1729 my $orth;
1730 my $key;
1731 foreach (@$clusters)
1732 {
1733 @row=split(/\t/,$_);
1734 splice(@row,0,1);
1735 $orth=&CountSpeicesInCluster(join("\t",@row));
1736 @row=split(/\t|\,/,join("\t",@row));
1737 foreach $key (@row)
1738 {
1739 if ($key ne "-")
1740 {
1741 $spid=substr($key,1,(index($key,'G')-1)); ###extract strains id
1742 if (exists($$hash{$orth."|".$spid}))
1743 {
1744 $$hash{$orth."|".$spid}++;
1745 }else
1746 {
1747 $$hash{$orth."|".$spid}=1;
1748 }
1749 }
1750 }
1751 }
1752 }
1753
1754 sub CountSpeicesInCluster()
1755 {
1756 (my $str)=@_;
1757 chomp($str);
1758 my @list=split(/\t/,$str);
1759 my $key;
1760 my $count=0;
1761
1762 foreach $key (@list)
1763 {
1764 if ($key ne "-")
1765 {
1766 $count++;
1767 }
1768 }
1769 return $count;
1770 }
1771
1772 sub CountGeneInCluster()
1773 {
1774 (my $str)=@_;
1775 chomp();
1776 my @list=split(/\t|\,/,$str);
1777 my $key;
1778 my $count=0;
1779 foreach $key (@list)
1780 {
1781 if ($key ne "-")
1782 {
1783 $count++;
1784 }
1785 }
1786 return $count;
1787 }
1788
1789
1790
1791
1792 sub GF()
1793 {
1794 &PrepareFasta(\@species,$inputDIR,".pep"); ###prepare pep file
1795 system("cat ".join(".pep ",@species).".pep > All.pep");
1796 system("grep '>' All.pep > genelist");
1797 system("$formatdb -p T -i All.pep");
1798 system("$blastall -p blastp -i All.pep -d All.pep -M BLOSUM45 -m9 -e $evalue -o All.blastp -a $thread");
1799 system("perl ./Blast_Filter.pl All.blastp All.pep $coverage $identity $score | $mcl - --abc -I 2.0 -o All.cluster");
1800 &FormatCluster("All.cluster","genelist",$spnum,\@clusters);
1801 #system("rm -rf *.pep* All.blastp All.cluster genelist");
1802 }
1803
1804 sub MP()
1805 {
1806 # (my $species,my $inputDIR,my $thread,my $evalue,my $score,my $coverage,my $identity)=@_;
1807 my $i;
1808 my $j;
1809 &PrepareFasta(\@species,$inputDIR,".pep"); ###prepare pep file
1810 system("cat ".join(".pep ",@species).".pep > All.pep");
1811 system("grep '>' All.pep > genelist");
1812 system("rm -rf All.pep");
1813 for ($i=0;$i<$spnum;$i++)
1814 {
1815 for ($j=$i+1;$j<$spnum;$j++)
1816 {
1817 system("perl ./inparanoid.pl $blastall $thread $formatdb $score $global $local $species[$i].pep $species[$j].pep");
1818 }
1819 }
1820 system("perl ./multiparanoid.pl -species ".join(".pep+",@species).".pep -unique 1");
1821 ###convert the MP result to table list based on gene
1822 &MP_Result_to_Table("MP.Cluster","All.cluster");
1823 &FormatCluster("All.cluster","genelist",$spnum,\@clusters);
1824 system("rm -rf sqltable.* *.pep* MP.Cluster genelist");
1825 }
1826
1827 sub fasta2phylip()
1828 {
1829 (my $input,my $output)=@_;
1830 use Bio::AlignIO;
1831 my $inputfilename = "10.aln";
1832 my $in= Bio::AlignIO->new(-file => $input ,
1833 -format => 'fasta');
1834 my $out = Bio::AlignIO->new(-file => ">$output" ,
1835 -format => 'phylip');
1836 while ( my $aln = $in->next_aln() )
1837 {
1838 $out->write_aln($aln);
1839 }
1840 }
1841
1842
1843 sub RemoveHeadGap()
1844 {
1845 (my $nal,my $hash)=@_;
1846 my %aln;
1847 my $key;
1848 my $gaplength=0;
1849 my $len1;
1850 my $len2;
1851 &ReadSequenceInToHash("$nal",\%aln);
1852 foreach $key (keys %aln)
1853 {
1854 $len1=length($aln{$key});
1855 $_=$aln{$key};
1856 s/^-+//;
1857 $len2=length($_);
1858 if (($len1-$len2)>$gaplength)
1859 {
1860 $gaplength=$len1-$len2;
1861 }
1862 }
1863 foreach $key (keys %aln)
1864 {
1865 $$hash{$key}=$$hash{$key}.substr($aln{$key},$gaplength,(length($aln{$key})-$gaplength));
1866 }
1867 }
1868
1869 sub PrepareFasta()
1870 {
1871 (my $species,my $inputDIR,my $extention)=@_;
1872 my $sp;
1873 my $file;
1874 my $i;
1875 my %hash;
1876 my $key;
1877 for ($i=0;$i<@$species;$i++)
1878 {
1879 $file=$inputDIR.$$species[$i].$extention;
1880 %hash=();
1881 &ReadSequenceInToHash($file,\%hash);
1882 open(R,">$$species[$i]${extention}") or die "Could write into $file\n";
1883 foreach $key (keys %hash)
1884 {
1885 print R ">S${i}G$key\n";
1886 print R $hash{$key}."\n";
1887 }
1888 close(R);
1889 }
1890 }
1891
1892 sub PrepareTable()
1893 {
1894 (my $species,my $inputDIR,my $extention)=@_;
1895 my @content;
1896 my $i;
1897 my @row;
1898 my $file;
1899 for ($i=0;$i<@$species;$i++)
1900 {
1901 $file=$inputDIR.$$species[$i].$extention;
1902 open(F,$file) or die "Could open $file\n";
1903 @content=<F>;
1904 close(F);
1905 chomp(@content);
1906 open(R,">$$species[$i]${extention}") or die "Could write into $file\n";
1907 foreach (@content)
1908 {
1909 @row=split(/\t/,$_);
1910 $row[0]="S${i}G$row[0]";
1911 if ($extention eq ".location")
1912 {
1913 $row[0]=$row[0]."\t".$row[0];
1914 }
1915 print R join("\t",@row)."\n";
1916 }
1917 close(R);
1918 }
1919 }
1920
1921 sub CheckExtraProgram
1922 {
1923 #(my $section, my $method, my $tmparray)=@_;
1924 my @error;
1925 my $ei=0;
1926
1927 #####cluster gene
1928 if (substr($section,0,1) eq "1")
1929 {
1930 ###MP: blastall formatdb
1931 ###GF: blastall formatdb mcl
1932
1933 if (!(-e $formatdb))
1934 {
1935 $error[$ei++]="formatdb is not found at $formatdb\n";
1936 }
1937
1938 if (!(-X $formatdb))
1939 {
1940 $error[$ei++]="there is not premission to execute $formatdb\n";
1941 }
1942
1943 if (!(-e $blastall))
1944 {
1945 $error[$ei++]="blastall is not found at $blastall\n";
1946 }
1947
1948 if (!(-X $blastall))
1949 {
1950 $error[$ei++]="there is not premission to execute $blastall\n";
1951 }
1952
1953 if ($method eq "GF")
1954 {
1955 if (!(-e $mcl))
1956 {
1957 $error[$ei++]="mcl is not found at $mcl\n";
1958 }
1959 if (!(-X $mcl))
1960 {
1961 $error[$ei++]="there is not premission to execute $mcl\n";
1962 }
1963 }
1964 }
1965
1966 #####CDS variation
1967 if (substr($section,2,1) eq "1")
1968 {
1969 if (!(-e $mafft))
1970 {
1971 $error[$ei++]="mafft is not found at $mafft\n";
1972 }
1973 if (!(-X $mafft))
1974 {
1975 $error[$ei++]="there is not premission to execute $mafft\n";
1976 }
1977 }
1978
1979 #####CDS variation
1980 if (substr($section,3,1) eq "1")
1981 {
1982 if (!(-e $mafft))
1983 {
1984 $error[$ei++]="mafft is not found at $mafft\n";
1985 }
1986 if (!(-X $mafft))
1987 {
1988 $error[$ei++]="there is not premission to execute $mafft\n";
1989 }
1990 }
1991 #####Evolution analysis
1992 if (substr($section,3,1) eq "1")
1993 {
1994 if (-e $seqboot)
1995 {
1996 $error[$ei++]="there is not premission to execute $seqboot\n" if(!(-X $seqboot));
1997 }else
1998 {
1999 $error[$ei++]="seqboot is not found at $seqboot\n";
2000 }
2001 if (-e $dnaml)
2002 {
2003 $error[$ei++]="there is not premission to execute $dnaml\n" if(!(-X $dnaml));
2004 }else
2005 {
2006 $error[$ei++]="dnaml is not found at $dnaml\n";
2007 }
2008 if (-e $dnadist)
2009 {
2010 $error[$ei++]="there is not premission to execute $dnadist\n" if(!(-X $dnadist));
2011 }else
2012 {
2013 $error[$ei++]="dnadist is not found at $dnadist\n";
2014 }
2015 if (-e $neighbor)
2016 {
2017 $error[$ei++]="there is not premission to execute $neighbor\n" if(!(-X $neighbor));
2018 }else
2019 {
2020 $error[$ei++]="neighbor is not found at $neighbor\n";
2021 }
2022 if (-e $consense)
2023 {
2024 $error[$ei++]="there is not premission to execute $consense\n" if(!(-X $consense));
2025 }else
2026 {
2027 $error[$ei++]="consense is not found at $consense\n";
2028 }
2029 if (-e $dnapars)
2030 {
2031 $error[$ei++]="there is not premission to execute $dnapars\n" if(!(-X $dnapars));
2032 }else
2033 {
2034 $error[$ei++]="dnapars is not found at $dnapars\n";
2035 }
2036 }
2037 #@$tmparray=(@$tmparray,@error);
2038 @tmp=(@tmp,@error);
2039 }
2040
2041
2042 sub CheckInputFile()
2043 {
2044 (my $species,my $inputDIR,my $section,my $method,my $tmparray)=@_;
2045 ####cluster
2046 if (substr($section,0,1) eq "1")
2047 {
2048 if ($method eq "MM")
2049 {
2050 @$tmparray=(@$tmparray,&chk2SEQ($species,$inputDIR)); ### check pep and nuc
2051 @$tmparray=(@$tmparray,&chktab($species,$inputDIR,".location"));### chk pep nuc location
2052 }else
2053 {
2054 @$tmparray=(@$tmparray,&chk1SEQ($species,$inputDIR));
2055 }
2056 }
2057 ###CDS variation
2058 if (substr($section,2,1) eq "1")
2059 {
2060 @$tmparray=(@$tmparray,&chk2SEQ($species,$inputDIR));
2061 }
2062 ###function analysis
2063 if (substr($section,4,1) eq "1")
2064 {
2065 @$tmparray=(@$tmparray,&chktab($species,$inputDIR,".function"));
2066 }
2067 }
2068
2069
2070 sub chk1SEQ()
2071 {
2072 (my $species,my $inputDIR)=@_;
2073 my @error;
2074 my $ei=0;
2075 my $sp;
2076 my $pepfile;
2077 my %pep;
2078 foreach $sp (@$species)
2079 {
2080 %pep=();
2081 $pepfile=$inputDIR.$sp.".pep";
2082 &ReadSequenceInToHash($pepfile,\%pep);
2083 if (scalar(keys %pep)<2)
2084 {
2085 $error[$ei++]="format error in $pepfile\n";
2086 }
2087 }
2088 return @error;
2089 }
2090
2091 sub chk2SEQ()
2092 {
2093 (my $species,my $inputDIR)=@_;
2094 my $sp;
2095 my %pep;
2096 my %nuc;
2097 my $pepfile;
2098 my $nucfile;
2099 my $key;
2100 my @error;
2101 my $ei=0;
2102 foreach $sp (@$species)
2103 {
2104 $pepfile=$inputDIR.$sp.".pep";
2105 $nucfile=$inputDIR.$sp.".nuc";
2106 %pep=();
2107 %nuc=();
2108 &ReadSequenceInToHash("$pepfile",\%pep);
2109 &ReadSequenceInToHash("$nucfile",\%nuc);
2110 if (scalar(keys %pep) ne scalar(keys %nuc))
2111 {
2112 $error[$ei++]="Sequences number is not consistent in the following two file:\n\t$pepfile\n\t$nucfile\n";
2113 }else
2114 {
2115 foreach $key (keys %pep)
2116 {
2117 if (exists($nuc{$key}))
2118 {
2119 if (length($nuc{$key}) ne ((length($pep{$key})+1)*3))
2120 {
2121 $error[$ei++]="the length of $key in $nucfile is not consistent with its corresponding protein length\n";
2122 }
2123 }else
2124 {
2125 $error[$ei++]="$key lost in $nucfile\n";
2126 }
2127 }
2128
2129 foreach $key (keys %nuc)
2130 {
2131 if (!exists($pep{$key}))
2132 {
2133 $error[$ei++]="1048 $key lost in $pepfile\n";
2134 }
2135 }
2136 }
2137 }
2138 return @error;
2139 }
2140
2141
2142 sub chktab()
2143 {
2144 (my $species,my $inputDIR,my $extention)=@_;
2145 my %pep;
2146 my @row;
2147 my $key;
2148 my %tab;
2149 my @error;
2150 my $ei=0;
2151 my $sp;
2152 my $tabfile;
2153 my $pepfile;
2154 foreach $sp (@$species)
2155 {
2156 %tab=();
2157 %pep=();
2158 $tabfile=$inputDIR.$sp.$extention;
2159 open(F,"$tabfile");
2160 while (<F>)
2161 {
2162 chomp();
2163 @row=split(/\t/,$_);
2164 if (scalar(@row)<3)
2165 {
2166 $error[$ei++]="format error in $tabfile\n";
2167 }else
2168 {
2169 $tab{$row[0]}=$row[1];
2170 }
2171 }
2172 close(F);
2173 $pepfile=$inputDIR.$sp.".pep";
2174 &ReadSequenceInToHash($pepfile,\%pep);
2175 foreach $key (keys %pep)
2176 {
2177 if (!exists($tab{$key}))
2178 {
2179 $error[$ei++]="sequence $key lost infomation in $tabfile\n";
2180 }
2181 }
2182 }
2183 return @error;
2184 }
2185
2186
2187
2188 sub ReadSequenceInToHash()
2189 {
2190 use Bio::SeqIO;
2191 (my $file,my $hash)=@_;
2192 my $seq;
2193 my $in=Bio::SeqIO->new(-file=>"$file",-format=>"fasta");
2194 while ($seq=$in->next_seq())
2195 {
2196 #$$hash{$id."|".$seq->id}=$seq->seq();
2197 $$hash{$seq->id}=$seq->seq();
2198 }
2199 }
2200
2201 sub ReadAlignmentToHash()
2202 {
2203 (my $file,my $hash)=@_;
2204 my $name="";
2205 my $seq="";
2206 my @content;
2207 my $line;
2208 open(F,"$file");
2209 @content=<F>;
2210 close(F);
2211 chomp(@content);
2212 for ($line=0;$line<@content;$line++)
2213 {
2214 if ($content[$line]=~/^>/)
2215 {
2216 if ($line>0)
2217 {
2218 $$hash{$name}=$seq;
2219 $name="";
2220 }
2221
2222 $_=$content[$line];
2223 s/^>//;
2224 $name=$_;
2225 $seq="";
2226 }else
2227 {
2228 if ($name ne "")
2229 {
2230 $seq=$seq.$content[$line];
2231 }
2232 }
2233 }
2234 $$hash{$name}=$seq;
2235 }
2236
2237
2238
2239 sub Combination()
2240 {
2241 (my $m,my $n,my $comRef)=@_;
2242 my $str="";
2243 my %hash;
2244 my $fpos;
2245 my $num0;
2246 my $rest;
2247 my $tmp;
2248 my $i;
2249 my $j;
2250 my $key;
2251 #my $m=scalar(@$array);
2252 my @combination;
2253
2254 for ($i=1;$i<=$n;$i++)
2255 {
2256 $str="1".$str;
2257 }
2258
2259 for ($i=1;$i<=($m-$n);$i++)
2260 {
2261 $str=$str."0";
2262 }
2263
2264 $hash{$str}=1;
2265 while ($str=~/10/)
2266 {
2267 $fpos=index($str,"10");
2268 $_=$str;
2269 s/10/01/;
2270 $str=$_;
2271 $tmp=substr($str,0,$fpos);
2272 $_=$tmp;
2273 s/0//g;
2274 $rest=$_;
2275 $num0=$fpos-length($_);
2276 for ($i=1;$i<=$num0;$i++)
2277 {
2278 $rest="$rest"."0";
2279 }
2280 $str="$rest".substr($str,$fpos,$m-$fpos);
2281 $hash{$str}=1;
2282 }
2283 $j=0;
2284 foreach $key (keys %hash)
2285 {
2286 $combination[$j]="";
2287 for ($i=0;$i<$m;$i++)
2288 {
2289 if (substr($key,$i,1) eq "1")
2290 {
2291 if ($combination[$j] ne "")
2292 {
2293 #$combination[$j]=$combination[$j]."\t".$$array[$i];
2294 $combination[$j]=$combination[$j]."\t".$i;
2295 }else
2296 {
2297 #$combination[$j]=$$array[$i]; ### For return species ID
2298 $combination[$j]=$i;
2299 }
2300 }
2301 }
2302 $j++;
2303 }
2304 @$comRef=@combination; ### update the data through the physic address
2305 }
2306
2307 sub ChkCombinationValue()
2308 {
2309 (my $m,my $n)=@_;
2310 my %hash;
2311 my %vhash;
2312 my $value=0;
2313 my $key;
2314 my @row;
2315 my @sdA;
2316 my @sdB;
2317
2318 ### initialization
2319 $hash{$m."-".$n}=1;
2320
2321 ### split combination
2322 while (scalar(keys %hash)>0 and $value<=$sampleSize)
2323 {
2324 foreach $key (keys %hash)
2325 {
2326 if ($value > $sampleSize) ### threshold
2327 {
2328 last;
2329 }
2330 if (!exists($hash{$key}))
2331 {
2332 next;
2333 }
2334 @row=split(/-/,$key);
2335 #print $row[0]."|".$row[1]."\n";
2336 if ($row[0] eq $row[1])
2337 {
2338 $value=$value+$hash{$key};
2339 }else
2340 {
2341 ##split
2342 $sdA[0]=$row[0]-1;
2343 $sdA[1]=$row[1];
2344 $sdB[0]=$row[0]-1;
2345 $sdB[1]=$row[1]-1;
2346 ##storing A
2347 if (($sdA[0] eq $sdA[1]) or $sdA[1] ==0)
2348 {
2349 $value=$value+$hash{$key};
2350 }else
2351 {
2352 if (exists($hash{$sdA[0]."-".$sdA[1]}))
2353 {
2354 $hash{$sdA[0]."-".$sdA[1]}=$hash{$sdA[0]."-".$sdA[1]}+$hash{$key};
2355 }else
2356 {
2357 $hash{$sdA[0]."-".$sdA[1]}=$hash{$key};
2358 }
2359 }
2360
2361 ##storing B
2362 if (($sdB[0] eq $sdB[1]) or $sdB[1]==0)
2363 {
2364 $value=$value+$hash{$key};
2365 }else
2366 {
2367 if (exists($hash{$sdB[0]."-".$sdB[1]}))
2368 {
2369 $hash{$sdB[0]."-".$sdB[1]}=$hash{$sdB[0]."-".$sdB[1]}+$hash{$key};
2370 }else
2371 {
2372 $hash{$sdB[0]."-".$sdB[1]}=$hash{$key};
2373 }
2374 }
2375 }
2376 #delete original combination
2377 delete($hash{$key});
2378 }
2379 }
2380
2381 if ($value>$sampleSize)
2382 {
2383 return 0;
2384 }else
2385 {
2386 return $value;
2387 }
2388 }
2389
2390
2391 sub SampleCombination()
2392 {
2393 (my $m,my $n,my $comRef)=@_;
2394 my %hash;
2395 my $sampleTimes=0;
2396 my @randNum;
2397 my @sortID;
2398 my $i;
2399 my $j;
2400 my $tmp;
2401 while ( scalar(keys %hash)<$sampleSize and $sampleTimes<($sampleSize*2))
2402 {
2403 for ($i=0;$i<$m;$i++) # generate random data
2404 {
2405 $randNum[$i]=int(100000 * rand(100));
2406 $sortID[$i]=$i;
2407 }
2408
2409 for ($i=0;$i<$m;$i++) # sorting random data
2410 {
2411 for ($j=0;$j<$m;$j++)
2412 {
2413 if ($randNum[$sortID[$i]]<$randNum[$sortID[$j]])
2414 {
2415 $tmp=$sortID[$i];
2416 $sortID[$i]=$sortID[$j];
2417 $sortID[$j]=$tmp;
2418 }
2419 }
2420 }
2421
2422 #storing data
2423 $tmp=join("\t",sort {$a<=>$b} (splice(@sortID,0,$n)));
2424 $hash{$tmp}=1;
2425 $sampleTimes++;
2426 }
2427 @$comRef=keys %hash;
2428 }
2429
2430
2431 sub PanGenomeNumber()
2432 {
2433 (my $spID)=@_;
2434 my $pan=0;
2435 my $core=0;
2436 my $count; #### counter;
2437 my @row;
2438
2439 foreach (@clusters)
2440 {
2441 $count=0;
2442 @row=split(/\t/,$_);
2443
2444 foreach (@$spID)
2445 {
2446 $count=$count+$row[$_];
2447 }
2448
2449 if ($count>0)
2450 {
2451 $pan++;
2452 if ($count == scalar(@$spID))
2453 {
2454 $core++;
2455 }
2456 }
2457 }
2458 return $pan."\t".$core;
2459 }
2460
2461 sub fit_model_A()
2462 {
2463 ### model y = A * x**B + C
2464 (my $xdata,my $ydata)=@_;
2465 my $i;
2466 my $b;
2467 my $max_B=0;
2468 my $max_R=0;
2469 my $max_A=0;
2470 my $max_A_interval;
2471 my $max_C=0;
2472 my $max_C_interval;
2473 my $R=1e-100;
2474 my $start;
2475 my $end;
2476 my $step;
2477 my @xValues;
2478 my @yValues;
2479
2480 $start=1;
2481 $step=0.001;
2482 $b=$start;
2483 $max_R=0;
2484 $R=1e-100;
2485
2486 use Statistics::LineFit;
2487 use Statistics::Distributions;
2488
2489 while ($max_R<=$R)
2490 {
2491 if (($b < 0.02) and ($b >-0.02))
2492 {
2493 $b=-0.02;
2494 }
2495
2496 for ($i=0;$i<@$xdata;$i++)
2497 {
2498 $xValues[$i]=$$xdata[$i]**$b;
2499 }
2500 @yValues=@$ydata;
2501 my $lineFit = Statistics::LineFit->new();
2502 $lineFit->setData (\@xValues, \@yValues) or die "Invalid data";
2503 (my $intercept, my $slope) = $lineFit->coefficients();
2504 my $rSquared = $lineFit->rSquared();
2505 my $meanSquaredError = $lineFit->meanSqError();
2506 my $durbinWatson = $lineFit->durbinWatson();
2507 my $sigma = $lineFit->sigma();
2508 (my $tStatIntercept, my $tStatSlope) = $lineFit->tStatistics();
2509 (my $varianceIntercept,my $varianceSlope) = $lineFit->varianceOfEstimates();
2510
2511 $max_R=$R;
2512 $R=$rSquared;
2513 if ($max_R<=$R)
2514 {
2515 $max_R=$R;
2516 ($max_C,$max_A)=$lineFit->coefficients();
2517 $max_A_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceSlope);
2518 $max_C_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceIntercept);
2519 }
2520 $b=$b-$step;
2521 }
2522 $max_B=$b;
2523 return ($max_R,$max_A,$max_A_interval,$max_B,$max_C,$max_C_interval);
2524
2525 }
2526
2527 sub fit_model_B()
2528 {
2529 ### model y = A * exp(x*B) + C
2530 (my $xdata,my $ydata)=@_;
2531 my $i;
2532 my $b;
2533 my $max_B=0;
2534 my $max_R=0;
2535 my $max_A=0;
2536 my $max_A_interval;
2537 my $max_C=0;
2538 my $max_C_interval;
2539 my $R=1e-100;
2540 my $start;
2541 my $end;
2542 my $step;
2543 my @xValues;
2544 my @yValues;
2545
2546 $start=0;
2547 $step=0.001;
2548 $b=$start;
2549 $max_R=0;
2550 $R=1e-100;
2551
2552 use Statistics::LineFit;
2553 use Statistics::Distributions;
2554
2555 while ($max_R<=$R)
2556 {
2557 if (($b < 0.02) and ($b >-0.02))
2558 {
2559 $b=-0.02;
2560 }
2561
2562 for ($i=0;$i<@$xdata;$i++)
2563 {
2564 $xValues[$i]=exp($$xdata[$i]*$b);
2565 }
2566 @yValues=@$ydata;
2567 my $lineFit = Statistics::LineFit->new();
2568 $lineFit->setData (\@xValues, \@yValues) or die "Invalid data";
2569 (my $intercept, my $slope) = $lineFit->coefficients();
2570 my $rSquared = $lineFit->rSquared();
2571 my $meanSquaredError = $lineFit->meanSqError();
2572 my $durbinWatson = $lineFit->durbinWatson();
2573 my $sigma = $lineFit->sigma();
2574 (my $tStatIntercept, my $tStatSlope) = $lineFit->tStatistics();
2575 (my $varianceIntercept,my $varianceSlope) = $lineFit->varianceOfEstimates();
2576
2577 $max_R=$R;
2578 $R=$rSquared;
2579 if ($max_R<=$R)
2580 {
2581 $max_R=$R;
2582 ($max_C,$max_A)=$lineFit->coefficients();
2583 $max_A_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceSlope);
2584 $max_C_interval=Statistics::Distributions::tdistr (($spnum-2),.025)*sqrt($varianceIntercept);
2585 }
2586 $b=$b-$step;
2587 }
2588 $max_B=$b;
2589 return ($max_R,$max_A,$max_A_interval,$max_B,$max_C,$max_C_interval);
2590 }
2591
2592
2593 sub ReadData2Array()
2594 {
2595 (my $file, my $array1,my $col1,my $array2,my $col2)=@_;
2596 my $i=0;
2597 open(F,$file);
2598 $_=<F>;
2599 while (<F>)
2600 {
2601 chomp();
2602 @_=split(/\t/,$_);
2603 $$array1[$i]=$_[$col1];
2604 $$array2[$i]=$_[$col2];
2605 $i++;
2606 }
2607 close(F);
2608 }
2609
2610 sub SumData()
2611 {
2612 (my $xdata,my $ydata,my $SumMethod)=@_;
2613 my %hash;
2614 my $i;
2615 my $key;
2616 my $max=0;
2617 for ($i=0;$i<@$xdata;$i++)
2618 {
2619 if (exists($hash{$$xdata[$i]}))
2620 {
2621 $hash{$$xdata[$i]}=$hash{$$xdata[$i]}." ".$$ydata[$i];
2622 }else
2623 {
2624 $hash{$$xdata[$i]}=$$ydata[$i];
2625 if ($$xdata[$i]>$max)
2626 {
2627 $max=$$xdata[$i];
2628 }
2629 }
2630 }
2631 @$xdata=qw();
2632 @$ydata=qw();
2633 $i=0;
2634 foreach $i (1..$max)
2635 {
2636 $$xdata[$i-1]=$i;
2637 if ($SumMethod eq "median")
2638 {
2639 $$ydata[$i-1]=&median($hash{$i});
2640 }else
2641 {
2642 $$ydata[$i-1]=&mean($hash{$i});
2643 }
2644 }
2645 #print join(",",@$xdata)."\n";
2646 #print join(",",@$ydata)."\n";
2647 }
2648
2649 sub median()
2650 {
2651 (my $data)=@_;
2652 my @data=split(/ /,$data);
2653 my $arraylen=scalar(@data);
2654 @data=sort{$a<=>$b} @data;
2655 if (int($arraylen/2)*2 == $arraylen)
2656 {
2657 return ($data[$arraylen/2]+$data[$arraylen/2-1])/2;
2658 }else
2659 {
2660 return $data[int($arraylen/2)];
2661 }
2662 }
2663
2664 sub mean()
2665 {
2666 (my $data)=@_;
2667 my @data=split(/ /,$data);
2668 my $sum=0;
2669 foreach (@data)
2670 {
2671 $sum=$sum+$_;
2672 }
2673 return int(($sum/scalar(@data))*1000)/1000;
2674 }
2675
2676 sub ReplaceName()
2677 {
2678 (my $sp,my $file)=@_;
2679 my @content;
2680 my $line;
2681 my $i;
2682 my $target;
2683 open(F,$file);
2684 @content=<F>;
2685 close(F);
2686 for ($line=0;$line<@content;$line++)
2687 {
2688 for ($i=0;$i<@$sp;$i++)
2689 {
2690 $_=$content[$line];
2691 $target="sp".$i."sp";
2692 s/$target/$$sp[$i]/;
2693 $content[$line]=$_;
2694 }
2695 }
2696 open(R,">$file");
2697 print R @content;
2698 close(R);
2699 }
2700
2701 sub MP_Result_to_Table()
2702 {
2703 (my $MPresult, my $outputfile)=@_;
2704 my %hash;
2705 my $maxid=0;
2706 my $i;
2707 my @row;
2708
2709 open(F,"$MPresult");
2710 $_=<F>;
2711 while (<F>)
2712 {
2713 @row=split(/\t/,$_);
2714 if (exists($hash{$row[0]}))
2715 {
2716 $hash{$row[0]}=$hash{$row[0]}."\t".$row[2];
2717 }else
2718 {
2719 $hash{$row[0]}=$row[2];
2720 if ($row[0]>$maxid)
2721 {
2722 $maxid=$row[0];
2723 }
2724 }
2725 }
2726 close(F);
2727
2728 open(R,">$outputfile");
2729 foreach $i (1..$maxid)
2730 {
2731 print R $hash{$i}."\n";
2732 }
2733 close(R);
2734 }
2735
2736
2737
2738 sub FormatCluster()
2739 {
2740 (my $infile,my $genelist,my $spnum,my $cluster)=@_;
2741 my %hash;
2742 my %gene;
2743 my $key;
2744 my @row;
2745 my $sp;
2746 my $line;
2747 my $i=0;
2748 my $j=0;
2749 my @content;
2750
2751 ### record gene in clusters
2752 open(F,"$infile");
2753 @content=<F>;
2754 close(F);
2755 chomp(@content);
2756 for ($line=0;$line<@content;$line++)
2757 {
2758 @row=split(/\t/,$content[$line]);
2759 foreach $key (@row)
2760 {
2761 $gene{$key}=1;
2762 }
2763 }
2764 ###retrieves gene which is not in clutsers
2765
2766 open(F,"$genelist");
2767 while ($key=<F>)
2768 {
2769 if ($key=~/^>/)
2770 {
2771 chomp($key);
2772 $_=$key;
2773 s/^>//;
2774 $key=$_;
2775 if (!exists($gene{$key}))
2776 {
2777 $content[$line]=$key;
2778 $line++;
2779 }
2780 }
2781 }
2782 close(F);
2783
2784 #### initialization @cluster
2785 @$cluster=qw();
2786 $j=0;
2787
2788 foreach $line (@content)
2789 {
2790 if ($line ne "")
2791 {
2792 %hash=();
2793 @row=split(/\t/,$line);
2794 foreach $key (@row)
2795 {
2796 $sp=substr($key,0,index($key,"G"));
2797 $gene{$key}=1;
2798 if (exists($hash{$sp}))
2799 {
2800 $hash{$sp}=$hash{$sp}.",".$key;
2801 }else
2802 {
2803 $hash{$sp}=$key;
2804 }
2805 }
2806
2807 $i=0;
2808 @row=qw();
2809
2810 foreach $i (0..($spnum-1))
2811 {
2812 if (exists($hash{"S$i"}))
2813 {
2814 $row[$i]=$hash{"S$i"};
2815 }else
2816 {
2817 $row[$i]="-";
2818 }
2819 }
2820 $$cluster[$j++]=join("\t",@row);
2821 }
2822 }
2823 }
2824
2825