Stochmap version 1.0: stochastic mapping on phylogenetic trees
Matthew Spencer,
School of Environmental Sciences,
University of Liverpool,
Liverpool,
L69 3GP,
UK.

m.spencer@liverpool.ac.uk

LICENSE INFORMATION

 Copyright 2009 Matthew Spencer, Simon Whelan

 This file is part of stochmap.

    stochmap is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    stochmap is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with stochmap.  If not, see <http://www.gnu.org/licenses/>.

 **********************************

 This program uses substantial amounts of code from MrBayes 3.1.2:
 mbmath_sub.c (just added a few things from bayes.c, command.c, model.c, mb.c to mbmath.c, to avoid having to compile the whole of MrBayes)
 mb.h
 mbmath_sub.h
 globals.h

 **********************************
 This program uses my_getopt, for which the license is given below:

my_getopt - a command-line argument parser
Copyright 1997-2001, Benjamin Sittler

Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without
restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.

*******************************************/
ABOUT STOCHMAP

Stochmap is designed to produce estimates of the ways in which a data set departs from what would be expected under a given evolutionary model. When used with software for visualizing trees and protein structures, it can help to assess the appropriateness of a model, and to identify unusual evolutionary events. It is described in Spencer et al (unpublished).


USING STOCHMAP

To compile under Linux:

cd source
make

To run:

./stochmap [infile] [outfile] [lfile]

infile is an input file containing information on transition probability matrices and partial likelihoods for each site on each edge, generated by Leaphy (Whelan, 2007). An example is given below.

outfile is an output file containing information on conditional and prior expectations of numbers of events, calculated using the stochastic mapping methods described in Minin and Suchard (2008). 

The first four lines contain information on the command used to generate the file. 

The next block of data contains the following information, in tab-separated columns:

The branch (numbered as in infile)
The process (should always be zero in the current version)
The site (numbered as in infile)
The conditional expected number of labelled events for the given site on the given branch
The prior expected number of labelled events for the given site on the given branch
Prior variance of expected number of events: unlikely to be useful because it is correct only for a tree with one branch, and increasingly overestimates variability as branches get further from the leaves.
A z-score for the given site on the given branch (conditional-prior)/sqrt(variance). This is unlikely to be useful. Instead, the parametric bootstrap approach described in Spencer et al. (unpublished) should be used to assess the magnitude of departures from prior expectation for branches and sites.

The next block contains the following tab-separated columns:

Branch (numbered as in infle)
The conditional expected number of labelled events on the branch, summed over sites
The prior expected number of labelled events on the branch, summed over sites
A branch zscore (should be ignored: unlikely to be useful)
The branch length

The next block gives site information as tab-separated columns:

The site number
The conditional expected number of labelled events at the site, summed over branches
The prior expected number of labelled events at the site, summed over branches
A site z-score (should be ignored, unlikely to be useful)

Finally, an overall z-score (should be ignored, unlikely to be useful).

lfile is a text file indicating the kinds of events that should be counted. It should consist of a tab-separated matrix of 1s and 0s, the same size as the instantaneous rate matrix for the evolutionary model. Events that should be counted are coded 1, and other events as 0. The diagonal elements should be 0. Examples are given below for counting changes of rate class and substitutions in a covarion model. The matrix should be symmetrical.

EXAMPLES

To analyze the EF-1alpha alignment described in Spencer et al. (unpublished):

The following will give quick visualizations of differences between conditional and prior numbers of events. Matlab (The Mathworks, Inc. Natick, MA) and PyMol (Schrodinger, 2008) are needed. In the paper, we also used parametric bootstraps to get p-values for these differences, but this requires the Leaphy software (Whelan, 2007) and takes much longer.

cd source

./stochmap ../examples/EF1A.partL ../examples/EF1A.shift.out ../examples/L100shift.txt

./stochmap ../examples/EF1A.partL ../examples/EF1A.subs.out ../examples/L100subs.txt

EF1A.partL was generated by Leaphy (Whelan, 2007) from the micro* alignment from Inagaki et al. (2004: alignment supplied by Andrew Roger), under the covarion model described in Spencer et al. (unpublished), using the M+F tree topology (EF1alpha_Feb11_fullnames.tree) from Inagaki et al. (2004).

The first command gets the differences between conditional and prior expectations of changes of rate class. The second command gets the differences between conditional and prior expectations of substitutions. The lfiles contain 100 x 100 matrices, because we are fitting an amino acid covarion model with 5 rate classes.

To generate data files for colouring trees and structures by the difference between conditional and prior expected numbers of changes of rate class: 

cd matlab
matlab&

In Matlab (The Mathworks, Inc., Natick, MA, tested with R2009b):

EF1alpha_example_script


This generates the following files in the examples directory:

EF1alpha_shiftcolour.dendro can be opened in Dendroscope (Huson et al. 2007), which will display the tree topology with edges coloured by the difference between conditional and prior expected numbers of changes of rate class. The colour map is similar to the 'cool' colour map defined in Matlab, with pink for edges having higher conditional than expected events, and cyan for lower. The topology of the tree is read from EF1alpha_Feb11_fullnames.dendro in the examples directory. The matlab function colortreeedges2.m, called by EF1alph_example_script.m, uses the edge weights in the .dendro file to match up data from the .out file with edges on the tree. Thus the required input .dendro files should be generated from Dendroscope by opening the tree file (e.g. EF1alpha_Feb11_fullnames.tree), selecting View|Show edge weights, and then File|Save as (save as a Dendroscope file).

siteshiftminus.txt contains values of conditional-prior expected changes of rate class where those values are negative, and 99 where they are non-negative. siteshiftplus.txt is similar, but for positive values. To visualize these on the structure (Andersen et al 2001, pdb file 1IJF.pdb):

pymol&

Then on the pymol command linem, type

@EF1alpha.pml

to run a Python script that will colour sites with positive values pink, sites with negative values cyan, and the EF-1beta molecule yellow.

REFERENCES

Andersen, G. R., Valente, L., Pedersen, L., Kinzy, T. G., and Nyborg, J. (2001). Crystal structures of nucleotide exchange intermediates in the eEF1A-eEF1B complex. Nature Structural Biology, 8(6):531-534.

Huson, D.H., Richter, D.C., Rausch, C., Dezulian, T., Franz, M., and Rupp, R. (2007). Dendroscope: an interactive viewer for large phylogenetic trees. BMC Bioinformatics, 8:460.

Inagaki, Y., Susko, E., Fast, N.M. and Roger, A.J. (2004). Covarion shifts cause a long-branch attraction artifact that unites microsporidia and archaebacteria in EF-1alpha phlogenies. Molecular Biology and Evolution 21(7):1340-1349.

Minin, V. N. and Suchard, M. A. (2008). Fast, accurate and simulation-free stochastic mapping. Philosophical Transactions of the Royal Society of London Series B, 364:3985-3995.

Schrodinger, LLC (2008). The PyMOL molecular graphics system, version 1.1r2pre.

Spencer, M., Blackburne, B., and Whelan, S. (unpublished). Identifying changes in evolutionary processes using stochastic mapping and covarion models.

Whelan, S. (2007). New approaches to phylogenetic tree search and their application to large numbers of protein alignments. Systematic Biology 56(5):727-740.
