Gene Trees: RAxML

A faster alternative to running MrBayes on each gene and BUCKy is to run instead RAxML on each gene (with bootstrap for later use), followed by counting the number of gene trees that have each quartet (which will be done within PhyloNetworks).

The advantage of MrBayes+BUCKy is that uncertainty in gene trees is reduced and integrated out to estimate quartet concordance factors. But it's slower.

To run RAxML on each gene, we can use the script raxml.pl in the script/ folder. To do so, still from the data folder baseline.gamma0.3_n300/, we first uncompress the tarball that has the alignments in nexus format. The script won't do this for us unfortunately. Also, to keep (thousands of) files organized, we uncompress these files in a folder nexus/ within our inut/ folder:

mkdir input/nexus
tar -xzf input/1_seqgen.tar.gz -C input/nexus

We can then run the script, which is not parallelized by the way and cannot send jobs to a cluster (unlike the scripts from TICR). RAxML sends a lot of output to the screen, 300 times if we have 300 genes, so we redirect this overwhelming output (and any possible error message) to a file mylog like this:

../scripts/raxml.pl --seqdir=input/nexus --raxmldir=raxml --astraldir=astral > mylog 2>&1 &

Let's check that everything ran smoothly and that we have our desired output. We have a new raxml/ directory with a bunch of things, including one file containing the best tree from each gene:

$ ls raxml/
besttrees.tgz	besttrees.tre	bootstrap	contrees.tgz	raxml.pl.log

$ head -n 4 raxml/besttrees.tre
((3:0.03722741374604016107,((1:0.00200161809596385333,2:0.00201350363534626934):0.02244234032100871773,4:0.01019471115416328323):0.00892507159858969140):0.01108038091781409287,5:0.05741592338682129787,6:0.06219449652861407801):0.0;
(((2:0.01227567249720099701,1:0.01840652224788371163):0.01469961352294601314,(3:0.00599379238893725259,4:0.00409652386305139546):0.02927896686177304530):0.00889065450733303973,5:0.03875872627548612032,6:0.07170201496079345316):0.0;
(((4:0.01028714366807356555,3:0.00183212712076471769):0.01214201212353699204,(1:0.00403084968085533086,2:0.00605548655541476161):0.01076366509971924193):0.01687768622690864057,5:0.06258702709089380978,6:0.08608073372761951281):0.0;
((1:0.00390601028837478034,2:0.01433216163343447279):0.02778992481020635397,((4:0.01423422190143536215,3:0.00409243847219765196):0.03425648595931906487,5:0.03106362878596600346):0.01154330379556680469,6:0.04409154556286316168):0.0;

We also have a new folder raxml/boostrap/ containing one bootstrap tree file per gene:

$ ls raxml/bootstrap/ | head
RAxML_bootstrap.1_seqgen1
RAxML_bootstrap.1_seqgen10
RAxML_bootstrap.1_seqgen100
RAxML_bootstrap.1_seqgen101
RAxML_bootstrap.1_seqgen102
RAxML_bootstrap.1_seqgen103
RAxML_bootstrap.1_seqgen104
RAxML_bootstrap.1_seqgen105
RAxML_bootstrap.1_seqgen106
RAxML_bootstrap.1_seqgen107

and we have a new folder astral/ that has the results from ASTRAL (which the script runs after RAxML by default):

$ ls astral/
BSlistfiles		astral.screenlog	astral.tre

$ head -3 astral/astral.tre 
(3,(4,((5,6)1:0.5764721782579237,(2,1)1:2.483246898369636)1:1.154550029033082)); 
(3,(4,((5,6)1:0.6311555045319776,(2,1)1:2.68060633252813)1:1.1234191104379085)); 
(3,(4,((5,6)1:0.6172017897581119,(2,1)1:2.68060633252813)1:1.1157854855828375)); 

$ tail -n 5 astral/astral.tre 
(3,(4,((5,6)1:0.6405673387143238,(2,1)1:2.610402073854883)1:1.1082096917743796)); 
(3,(4,((5,6)1:0.5545414747639535,(2,1)1:2.68060633252813)1:1.0932281381587627)); 
(3,(4,((5,6)1:0.5898652355943617,(2,1)1:2.7561138850362754)1:1.1006908593603526)); 
(3,(4,((5,6)100.0,(1,2)100.0)100.0)); 
(3,(4,((5,6)100.0:0.6548532959618005,(2,1)100.0:2.6448882499260518)100.0:1.1388634328653822)); 

In this astral/ folder, we also have the list of all bootstrap RAxML files, which was needed for ASTRAL, and which will be reused in PhyloNetworks to get bootstrap species networks:

$ head astral/BSlistfiles 
raxml/bootstrap/RAxML_bootstrap.1_seqgen1
raxml/bootstrap/RAxML_bootstrap.1_seqgen10
raxml/bootstrap/RAxML_bootstrap.1_seqgen100
raxml/bootstrap/RAxML_bootstrap.1_seqgen101
raxml/bootstrap/RAxML_bootstrap.1_seqgen102
raxml/bootstrap/RAxML_bootstrap.1_seqgen103
raxml/bootstrap/RAxML_bootstrap.1_seqgen104
raxml/bootstrap/RAxML_bootstrap.1_seqgen105
raxml/bootstrap/RAxML_bootstrap.1_seqgen106
raxml/bootstrap/RAxML_bootstrap.1_seqgen107

Now that we got all our results with no error, we can remove the uninteresting output from RAxML in mylog:

rm mylog

Next: from quartet CFs or gene trees to phylogenetic networks with PhyloNetworks.