HTTPS SSH
r8spermute - Runs all permutations of r8s constraints.
            (c) Simon J. Greenhill, 2009
======================================================

This is a small program to assess how robust divergence times are based on the calibrations given to the program r8s.

The program r8s by Mike Sanderson (http://loco.biosci.arizona.edu/r8s/) [1, 2] performs a number of phylogenetic dating methods on a set of phylogenetic trees. In essence, you plug in a number of calibration points from known historical evidence and the program r8s takes your estimated trees, and "smoothes" the rates of change observed on those trees, using the calibrations to convert the branches into time. 

However, it is often useful to explore the effect that certain calibration has on the various date estimates. For example, in the paper [3], I wanted to estimate the age at which a large language family originated. We had about 10 different calibration points based on historical evidence. Some of these calibrations were, however, more controversial than the others. So, I wrote this program to take a tree (or set of trees) and analyse them under all the combinations of these calibrations. This allowed us to assess the relative effects on the date estimates of the different calibrations.

If you find this program useful, please cite reference [3].

Contact me at simon@simon.net.nz if you have any problems.

Usage:

You will need two things: 1) A set of trees, 2) A r8s command block. Examples of both of these are provided in the "examples" directory. 

Your r8s block file will look something like this:

begin rates;
    blformat lengths=persite nsites=1000 ultrametric=no round=yes;
    collapse;
    
    mrca node1 Daisy Fiona;
    constrain taxon=node1 min_age=1.8 max_age=2.5; 
    
    mrca node2 Robert Tamara;
    constrain taxon=node2 min_age=1.1 max_age=1.3;
    
    unfixage taxon=Simon;
    constrain taxon=Simon min_age=0.7 max_age=1.2; 
    
    set num_restarts=5;
    set smoothing=10;
    
    divtime method=pl algorithm=tn;
    showage shownamed = yes;
    
    profile taxon=node1 parameter=age;
    profile taxon=node2 parameter=age;
    describe plot=chrono_description;
end;

... thus, we have 3 calibrations - node1, node2, and a constraint on the terminal taxon "Simon". 

We can then run r8spermute like this:

    r8spermute.py example/example.trees example/example.r8s

This will start r8s running and you should then see some output like this:

Constraints found:  3
Constraints ignored:  0
Total Constraint Combinations:  7

Running analysis 1 of 7 (0.14%)
Staging in:  /var/folders/9Z/9Z3O2T2o2RmgRU+F75TSx++++TQ/-Tmp-/tmpXn1lw9
Constraints: 
 +  constrain taxon=node1 min_age=1.8 max_age=2.5;
Logging in: example/example_100.log
r8s version 1.71 (compiled May 16 2006)
[...reading file /var/folders/9Z/9Z3O2T2o2RmgRU+F75TSx++++TQ/-Tmp-/tmpXn1lw9]
Elapsed time: 4.28s (0.07 minutes)
Average run-time so far: 4.28s (0.07 minutes)
Estimated run-time left: 25.68s (0.43 minutes)

...etc.

This tells us that we've got one constraint turned on (+  constrain taxon=node1 min_age=1.8 max_age=2.5;)' in this sub-analysis, as well as various other run-time statistics.

....and now we wait for this to finish. It could take some time..

When the analysis has finished, we use the program r8spermute_results.py to parse the results:

python r8spermute_results.py example/

This will loop over the log files and extract the relevant information. You probably want to output this information to a file, so run it like this:

python r8spermute_results.py example/ > results.txt

Now, we can look at the results, which will look something like this:

node1 node2 Simon node1 node2
0 0 1 15.22 2.84
0 1 0 54.60 1.30
0 1 1 49.45 1.30
1 0 0 2.50 4.38
1 0 1 2.50 4.18
1 1 0 2.50 1.30
1 1 1 2.50 1.30

It's probably easiest to load this into a spreadsheet. Each row is an analysis. The first 3 columns (containing 1's and 0's) are our constraints and whether they're turned on (=1) or off (=0) in the analysis. The last two columns "node1" and "node2" are our estimated ages of those nodes, under the calibrations. 

So, the first line shows that node1 is estimated to be 15.22 years old and node2 is 2.84 years old when only the calibration called "Simon" is used. 

The second line shows that when only node2 is calibrated, the estimate for node1 and node2 is 54.60 and 1.30 respectively. 

In contrast, the very last line shows that when all calibrations are turned on, then the age estimates are 2.50 and 1.30 respectively. 

References:

[1] Sanderson, M. J. 1997. A nonparametric approach to estimating divergence    
    times in the absence of rate constancy. Mol. Biol. Evol. 14:1218-1231.

[2] Sanderson, M. J. 2002. Estimating absolute rates of molecular evolution and 
    divergence times: a penalized likelihood approach. Mol. Biol. Evol.     
    19:101-109.

[3] Gray, R.D., Drummond, A.J., & Greenhill, S.J. (2009) Language Phylogenies
    Reveal Expansion Pulses and Pauses in Pacific Settlement. Science, 323:         
    479-483.