Merge pull request LANL-Bioinformatics#9 from mshakya/master

minor bug fixes

Author: mshakya

.gitignore

@@ -26,4 +26,5 @@ thirdParty/evolver.out
 thirdParty/errors.log
 messages.log
 *.log
-docs/_build/*
+docs/_build/*
+.vscode/

README.md

@@ -77,11 +77,9 @@ Then test phame to see if the installation was complete using a short test run:
 
 ## CITATION
 
-From raw reads to trees: Whole genome SNP phylogenetics across the tree of life.
+Shakya, M., Ahmed. S.A, Davenport K.W., Flynn M.C., Lo. C-C, Chain P.S.G. Standardized phylogenetic and molecular evolutionary analysis applied to species across the microbial tree of life. Sci. Rep. 10, 1723 (2020).
 
-Sanaa Afroz Ahmed, Chien-Chi Lo, Po-E Li, Karen W Davenport, Patrick S.G. Chain
-
-bioRxiv doi: http://dx.doi.org/10.1101/032250
+doi:  https://doi.org/10.1038/s41598-020-58356-1
 
 --------------------------------------------------------------
 ## Contact

lib/PhaME.pm

@@ -593,6 +593,7 @@ sub buildTree {
     my $thread = shift;
     my $tree   = shift;
     my $name   = shift;
+    my $bsignal = shift; # signal showing if bootstrap should be done or not.
     my $bootstrap = shift;
     my $error  = shift;
     my $log    = shift;
@@ -616,11 +617,25 @@ sub buildTree {
             # = "raxmlHPC-PTHREADS -p 10 -T $thread -m GTRGAMMAI -s $outdir/$name\_snp_alignment.fna -w $outdir -n $name 2>>$error >> $log\n\n";
         print OUT $raxml;
         if ( system($raxml) ) { die "Error running $raxml.\n"; }
+        if ( $bsignal == 1 ) {
+            open( OUT, ">>$log" );
+            my $bootTrees
+                = "raxmlHPC-PTHREADS -p 10 -T $thread -m GTRGAMMAI -b 10000 -t $outdir/RAxML_bestTree.$name -s $outdir/$name\_all_alignment.fna -w $outdir -N $bootstrap -n $name\_b -k 2>>$error >> $log\n\n";
+            print OUT $bootTrees;
+            if ( system($bootTrees) ) { die "Error running $bootTrees.\n"; }
+            my $bestTree
+            = "raxmlHPC-PTHREADS -p 10 -T $thread -f b -m GTRGAMMAI -t $outdir/RAxML_bestTree.$name -s $outdir/$name\_all_alignment.fna -z $outdir/RAxML_bootstrap.$name\_b -w $outdir -n $name\_best 2>>$error >> $log\n\n";
+            print OUT $bestTree;
+            if ( system($bestTree) ) { die "Error running $bestTree.\n"; }
+
+        return "Bootstrap complete";
+        close OUT;
+    }
     }
 
     if ( $tree == 3 || $tree == 4 ) {
         print OUT "Reconstructing phylogeny using IQ-tree after finding the best model\n";
-        if ($bootstrap > 0){
+        if ($bsignal == 1){
             print OUT "Also bootstraping IQ-Trees trees\n";
             print OUT "\n";
             my $iqtree
@@ -667,11 +682,11 @@ sub bootstrap {
     open( OUT, ">>$log" );
     if ( $tree == 2 || $tree == 4 ) {
         my $bootTrees
-            = "raxmlHPC-PTHREADS -p 10 -T $thread -m GTRGAMMAI -b 10000 -t $outdir/RAxML_bestTree.$name -s $outdir/$name\_snp_alignment.fna -w $outdir -N $bootstrap -n $name\_b -k 2>>$error >> $log\n\n";
+            = "raxmlHPC-PTHREADS -p 10 -T $thread -m GTRGAMMAI -b 10000 -t $outdir/RAxML_bestTree.$name -s $outdir/$name\_all_alignment.fna -w $outdir -N $bootstrap -n $name\_b -k 2>>$error >> $log\n\n";
         print OUT $bootTrees;
         if ( system($bootTrees) ) { die "Error running $bootTrees.\n"; }
         my $bestTree
-            = "raxmlHPC-PTHREADS -p 10 -T $thread -f b -m GTRGAMMAI -t $outdir/RAxML_bestTree.$name -s $outdir/$name\_snp_alignment.fna -z $outdir/RAxML_bootstrap.$name\_b -w $outdir -n $name\_best 2>>$error >> $log\n\n";
+            = "raxmlHPC-PTHREADS -p 10 -T $thread -f b -m GTRGAMMAI -t $outdir/RAxML_bestTree.$name -s $outdir/$name\_all_alignment.fna -z $outdir/RAxML_bootstrap.$name\_b -w $outdir -n $name\_best 2>>$error >> $log\n\n";
         print OUT $bestTree;
         if ( system($bestTree) ) { die "Error running $bestTree.\n"; }
 
@@ -739,14 +754,15 @@ sub translateGenes {
     my $bindir  = shift;
     my $thread  = shift;
     my $program = shift;
+    my $gencode = shift;
     my $error   = shift;
     my $log     = shift;
     my $genedir = $dir . '/PSgenes';
 
     open( OUT, ">>$log" );
     print OUT "\n";
     my $translate
-        = "parallel_run.pl -d $genedir -t $thread -m $program 2>>$error >> $log\n\n";
+        = "parallel_run.pl -d $genedir -t $thread -m $program -g $gencode 2>>$error >> $log\n\n";
     print OUT $translate;
     if ( system($translate) ) { die "Error running $translate.\n"; }
     close OUT;
@@ -758,14 +774,15 @@ sub alignGenes {
     my $bindir  = shift;
     my $thread  = shift;
     my $program = shift;
+    my $gencode = shift;
     my $error   = shift;
     my $log     = shift;
     my $genedir = $dir . '/PSgenes';
 
     open( OUT, ">>$log" );
     print OUT "\n";
     my $align
-        = "parallel_run.pl -d $genedir -t $thread -m $program 2>>$error >> $log\n\n";
+        = "parallel_run.pl -d $genedir -t $thread -m $program -g $gencode 2>>$error >> $log\n\n";
     print OUT $align;
     if ( system($align) ) { die "Error running $align.\n"; }
     close OUT;
@@ -777,14 +794,15 @@ sub revTransGenes {
     my $bindir  = shift;
     my $thread  = shift;
     my $program = shift;
+    my $gencode = shift;
     my $error   = shift;
     my $log     = shift;
     my $genedir = $dir . '/PSgenes';
 
     open( OUT, ">>$log" );
     print OUT "\n";    
     my $revTrans
-        = "parallel_run.pl -d $genedir -t $thread -m $program 2>>$error >> $log\n\n";
+        = "parallel_run.pl -d $genedir -t $thread -m $program -g $gencode 2>>$error >> $log\n\n";
     print OUT $revTrans;
     if ( system($revTrans) ) { die "Error running $revTrans.\n"; }
     close OUT;

src/parallel_run.pl

@@ -20,6 +20,7 @@
    'd=s'      => \$dir,
    't=i'      => \$thread,
    'm=s'      => \$program,
+   'g=i'      => \$gencode,
    'h|help'   => sub {usage()},
 );
 
@@ -74,7 +75,7 @@ sub usage
    my ($ref_file_name,$ref_file_path,$ref_file_suffix)=fileparse("$files[$i]", qr/\.[^.]*/);
    if ($program=~/oneline/){`fasta_oneline.pl $files[$i] > $ref_file_path/$ref_file_name.faa`}
    if ($program=~ /mafft/i){`unset MAFFT_BINARIES; mafft --quiet --auto $mafft_options $files[$i] > $ref_file_path/$ref_file_name.msa`;}
-   if ($program=~/pal2nal/){`pal2nal.pl $files[$i] $ref_file_path/$ref_file_name.fna -output fasta > $ref_file_path/$ref_file_name.cdn`;}
+   if ($program=~/pal2nal/){`pal2nal.pl $files[$i] $ref_file_path/$ref_file_name.fna -codontable $gencode -output fasta > $ref_file_path/$ref_file_name.cdn`;}
 #   if ($program=~/pal2nal/){`pal2nal.pl $files[$i] $ref_file_path/$ref_file_name.fna -output fasta > $ref_file_path/$ref_file_name.cdn`;}
    if ($program=~ /translate/){
     print "Translating $files[$i] to $ref_file_path/$ref_file_name.faa\n";