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