Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

extended RAxML's capability to read in pre-computed binary model para…

…meter files as produced

by the -f e option for the following options:

1. standard tree search
2. per-site log likelihood calculation for a bunch of trees (-f g option)
3. log-likelihood evaluation of a bunch of trees (-f n option)

also added an example script for the read analyses called pavlatoras.pl
  • Loading branch information...
commit f8cdf42e14efd0ba62b0b103e86ca119904ddbff 1 parent 9a1aa34
@stamatak authored
View
26 axml.c
@@ -7282,8 +7282,17 @@ static void computePerSiteLLs(tree *tr, analdef *adef, char *bootStrapFileName)
treeReadLen(treeFile, tr, FALSE, FALSE, FALSE, adef, TRUE);
assert(tr->ntips == tr->mxtips);
- if(i == 0)
- modOpt(tr, adef, TRUE, adef->likelihoodEpsilon, TRUE);
+ if(i == 0)
+ {
+ if(adef->useBinaryModelFile)
+ {
+ readBinaryModel(tr);
+ evaluateGenericInitrav(tr, tr->start);
+ treeEvaluate(tr, 2);
+ }
+ else
+ modOpt(tr, adef, TRUE, adef->likelihoodEpsilon, TRUE);
+ }
else
treeEvaluate(tr, 2);
@@ -8033,7 +8042,12 @@ void readBinaryModel(tree *tr)
model;
FILE
- *f = fopen(binaryModelParamsInputFileName, "r");
+ *f;
+
+
+ printBothOpen("\nRAxML is reading a binary model file and not optimizing model params\n");
+
+ f = fopen(binaryModelParamsInputFileName, "r");
/* cdta */
@@ -8064,7 +8078,7 @@ void readBinaryModel(tree *tr)
}
#ifdef _USE_PTHREADS
- /* TODO need to add stuff if we want this to work for CAT models and SP implementations as well */
+ /* TODO need to add stuff if we want this to work for CAT models as well */
masterBarrier(THREAD_RESET_MODEL, tr);
#endif
@@ -8427,9 +8441,7 @@ int main (int argc, char *argv[])
printLog(tr, adef, TRUE);
printResult(tr, adef, TRUE);
}
-
-
-
+
break;
case ANCESTRAL_STATES:
initModel(tr, rdta, cdta, adef);
View
9 multiple.c
@@ -1460,7 +1460,14 @@ void doInference(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta)
{
restoreTL(rl, tr, best);
onlyInitrav(tr, tr->start);
- modOpt(tr, adef, FALSE, adef->likelihoodEpsilon, FALSE);
+ if(!adef->useBinaryModelFile)
+ modOpt(tr, adef, FALSE, adef->likelihoodEpsilon, FALSE);
+ else
+ {
+ readBinaryModel(tr);
+ evaluateGenericInitrav(tr, tr->start);
+ treeEvaluate(tr, 2);
+ }
bestLH = tr->likelihood;
tr->likelihoods[best] = tr->likelihood;
saveTL(rl, tr, best);
View
2  optimizeModel.c
@@ -2408,7 +2408,7 @@ void modOpt(tree *tr, analdef *adef, boolean resetModel, double likelihoodEpsilo
int *unlinked = (int *)malloc(sizeof(int) * tr->NumberOfModels);
-
+ assert(!adef->useBinaryModelFile);
modelEpsilon = 0.0001;
View
44 pavlatoras.pl
@@ -0,0 +1,44 @@
+#!/usr/bin/perl
+use strict;
+use warnings;
+
+# those are the ref alignment and ref tree
+
+my $referenceData = "44.phy";
+my $referenceTree = "44.tree";
+
+# this is the subalignment, i.e., placements on a branch
+
+my $subAlignment = "10.phy";
+
+# estimate model parames on reference tree to obtain a binary model parameter file
+
+system("./raxmlHPC-SSE3 -f e -m GTRGAMMA -t ".$referenceTree." -s ".$referenceData." -n MODEL\n");
+
+# do a ML search on the subalignment using the binary model data file, that is, without optimizing
+# ML model params at all
+
+system("./raxmlHPC-SSE3 -m GTRGAMMA -R RAxML_binaryModelParameters.MODEL -p 12345 -s ".$subAlignment." -n ML_TREE\n");
+
+my $fileNames = "cat RAxML_bestTree.ML_TREE ";
+my $range = 1000000;
+
+# now generate 10 or more random trees on the subalignment of reads
+
+for(my $i = 0; $i < 10; $i++)
+{
+ my $r = rand($range);
+
+ system("./raxmlHPC-SSE3 -y -d -m GTRCAT -p ".$r." -s ".$subAlignment." -n R".$i."\n");
+
+ $fileNames = $fileNames." RAxML_randomTree.R".$i." ";
+}
+
+# concatenate all tree files, i.e., the ML tree and the 10 random trees
+
+system($fileNames. " > allTrees");
+
+# compute per site log likelihoods on all those trees, the per-site log likes can then be used with CONSEL for a significance
+# test
+
+system("./raxmlHPC-SSE3 -m GTRGAMMA -R RAxML_binaryModelParameters.MODEL -s ".$subAlignment." -z allTrees -f g -n allLikes\n");
View
25 searchAlgo.c
@@ -1193,7 +1193,16 @@ void computeBIGRAPID (tree *tr, analdef *adef, boolean estimateModel)
Thorough = 0;
if(estimateModel)
- modOpt(tr, adef, FALSE, 10.0, FALSE);
+ {
+ if(adef->useBinaryModelFile)
+ {
+ readBinaryModel(tr);
+ evaluateGenericInitrav(tr, tr->start);
+ treeEvaluate(tr, 2);
+ }
+ else
+ modOpt(tr, adef, FALSE, 10.0, FALSE);
+ }
else
treeEvaluate(tr, 2);
@@ -1208,7 +1217,12 @@ void computeBIGRAPID (tree *tr, analdef *adef, boolean estimateModel)
bestTrav = adef->bestTrav = adef->initial;
if(estimateModel)
- modOpt(tr, adef, FALSE, 5.0, FALSE);
+ {
+ if(adef->useBinaryModelFile)
+ treeEvaluate(tr, 2);
+ else
+ modOpt(tr, adef, FALSE, 5.0, FALSE);
+ }
else
treeEvaluate(tr, 1);
@@ -1303,7 +1317,12 @@ void computeBIGRAPID (tree *tr, analdef *adef, boolean estimateModel)
recallBestTree(bestT, 1, tr);
if(estimateModel)
- modOpt(tr, adef, FALSE, 1.0, FALSE);
+ {
+ if(adef->useBinaryModelFile)
+ treeEvaluate(tr, 2);
+ else
+ modOpt(tr, adef, FALSE, 1.0, FALSE);
+ }
else
treeEvaluate(tr, 1.0);
Please sign in to comment.
Something went wrong with that request. Please try again.