ML_genedata

Copyright Matthew Spencer, University of Liverpool, 2008.

PURPOSE

ML_genedata estimates mixture models for the gain and loss of gene families on a phylogenetic tree using maximum likelihood. Allows changes of evolutionary model within categories on specified edges. Calculates empirical Bayes posterior probabilities of each gene family being in each category. Also includes a simulator for the models that can be estimated.

INSTALLATION AND SYSTEM REQUIREMENTS

After unpacking, 

> make

should work on most linux/unix systems (uses gcc to compile). Tested on an Intel Xeon 3Ghz processor with 2 GB RAM running Debian Etch (kernel 2.6.18-6-686), with gcc 4:4.1.1-15. Requires the Gnu Scientific Library Version 1.8-2 or later (e.g. Debian packages libgsl0, libgsl0-dev).

USAGE

ML_genedata [-h] [-d] [-q] [-n] [-l] [-f] [-p pfile] [-c#] [-g#] [-s#] [-o:] [-m#] [infile] [outfile] [treefile]

-h: help

-d: print out genome data and unobservable patterns. Note that patterns that occur more than once will only be printed as a single column, but with a weight greater than 1. The row labelled 'Matches' reports the zero-indexed column of the printed data to which each gene family corresponds.

-q: quiet operation (default).

-n: noisy operation: settings, progress, etc.

-l: include gene family likelihoods and empirical Bayes posterior probabilities for each category and gene family in output file.

-f: fix edge lengths at input values.

-p pfile: read starting parameters (other than edge lengths) from file pfile with the following entries:

	Relative rates of gene gain and loss for each rate matrix (these will be rescaled to one expected event per unit time at equilibrium).
	Mixing probability for each category.
	Probability of gene family absence at root for each category.
	Shape parameter alpha for gamma distribution (has no effect unless -g used with more than one rate class).
	For example,
		1.2 0.1
		0.1 0.7
		0.25 0.75
		0.3 0.4
                0.5
	specifies relative gene gain and loss rates of 1.2 and 0.1 for the first rate matrix, relative gain and loss rates of 0.1 and 0.7 for the second rate matrix, mixing probabilities of 0.25 and 0.75 for two categories, and probabilities of gene absence at the root of 0.3 and 0.4 for these two categories, and shape parameter 0.5.

-c#: setting for conditioning on observable gene families. Options are 0 (all gene families are observable, default), 1 (condition on presence in at least one genome), 2 (condition on presence in at least two genomes), 3 (condition on presence in at least three genomes). Any gene family that appears in the data set but does not meet the observability criterion will have a weight of zero when calculating the likelihood and estimating parameters.

-g#: number of gamma rate categories (default 1).

-s#: constant by which to multiply all edge lengths when calculating likelihoods (default 1.0). Must be non-negative.

-m# : generate the specified number of simulated gene families, instead of using data in input file. The output from -m is a phylip-format file (readable as an input file by this program), with rows at the end containing the root state, the true category for each gene family, and the command used to generate the data. Can be used with -p, -c, -g and -s. Seeds the pseudorandom number generator with system time, so use a sleep command between reps if running multiple times from a script.

-o: setting for optimization method (parameters other than edge lengths). Options are S (simplex), Q (quasi-Newton, default) and N (no parameter optimization).

Input file format: a phylip format data matrix, arranged like an alignment (no input file with -m option). Taxon names must be no more than 10 characters and should not contain white space.  Number of taxa and gene families on the first row.  0 (absence) and 1 (presence) are the only symbols allowed.

The output is a text file containing:

-Version information

-The options used

-The tree topology with initial edge lengths, the rate matrix index for each category on each edge, and the initial parameter estimates.

-The log likelihood, gene family log likelihoods (if -l specified), and empirical Bayes posterior probabilities and category with maximum posterior probability (if -l specified and more than one category).  These are conditional on observability, and so are not meaningful for any gene families that appear in the data file but do not meet the specified observabilty criterion (-c). Such gene families are ignored when calculating the likelihood. The output includes a 0/1 vector where 0 indicates an unobservable and 1 an observable gene family.

-The estimated edge lengths and parameters (if these were not fixed). The parameters are given here as rate matrices, mixing probabilities, root probabilities of absence, and shape parameter.

-The gradient and Hessian for the parameters other than edge lengths. The parameters are given here as stationary probabilities of absence, mixing probabilities for all but the last category, root probabilities of absence, and shape parameter (the transformed values which are actually optimized are also given). The elements of the gradient for each parameter should be close to zero, and the Hessian should be negative definite, if an internal local optimum has been found. Standard errors for the parameters can be obtained from the square root of the diagonal elements of the inverse of the negative Hessian. For example, using Matlab:

Hess=[<copy and paste from output file>];
eig(Hess) %all eigenvalues should be negative
SE=sqrt(diag(inv(-Hess))) %standard errors

This analysis treats the edge lengths as known, which of course they aren't likely to be. However, if all the edge lengths were included, it's likely that the estimate of the Hessian would be unstable.

The gradient and Hessian are found for the transformed parameters, and then back-transformed using J'v for the gradient, and J'HJ for the Hessian, where J is the Jacobian (derivatives of transformed parameters with respect to untransformed parameters), v is the gradient with respect to transformed parameters, and H is the Hessian with respect to transformed parameters. The Jacobian hasn't been checked for cases with more than two major categories (where it won't be diagonal): use at your own risk, and report any problems to me.

Tree file format: a Newick tree with initial edge lengths. Subsequent lines (one per category of gene families) must contain the same tree topology, with the indices of the rate matrix used on that edge in that category replacing the edge lengths. The indices should be non-negative integers, and you should use all the integers between zero and the largest index. For example, use indices {0,1,2} rather than {1,2,3} or {0,1,3}. There is no ordering associated with the indices: rate matrix 1 does not imply higher rates than rate matrix 0. All rate matrices are scaled so that there is one expected event per unit time at equilibrium, but they will differ in the relative rates of gene gains and losses. The number of rate matrices does not have to be the same as the number of categories, and a rate matrix may appear in multiple categories.


EXAMPLES

1. Estimate log likelihood for a model with two categories of gene families, one with the same rates everywhere, and the other with different rates in edges leading only to parasites. Assume edge lengths are fixed (-f) at those in the input tree file (examples/Fedges.tree), and other parameters are fixed (-oN) at the values given in Fparms.txt (-p examples/Fparms.txt). Assume gene families are observable if they appear in at least three genomes (-c3). Use gamma rate variation with four rate classes (-g4). Noisy operation (-n). Input data are examples/COG_bacteria_input.phy. Output to F.out:

./ML_genedata -c3 -n -g4 -oN -f -p examples/Fparms.txt examples/COG_bacteria_input.phy F.out examples/Fedges.tree

2. Generate a simulated data set using parameters in Fparms.txt (-p examples/Fparms.txt) and edge lengths in examples/Fedges.tree. Two categories of gene families, one with the same rates everywhere, the other with different rates on edges leading only to parasites. Gene families are observable if they appear in at least three genomes (-c3). Four gamma rate classes (-g4). 3944 observable gene families (-m3944). Output to simF.txt: 

./ML_genedata -c3 -g4 -m3944 -p examples/Fparms.txt simF.txt examples/Fedges.tree

3. Estimate log likelihood for the same model as example 1, but optimize both edge lengths (initial values in examples/start.tree) and other parameters (initial values in examples/startparms.txt), and write information on gene family log likelihoods and posterior probabilities (-l). Output to Ffit.out:

./ML_genedata -c3 -l -n -g4 -p examples/startparms.txt examples/COG_bacteria_input.phy Ffit.out examples/start.tree

CONTACT

Matthew Spencer, School of Biological Sciences, University of Liverpool, Crown Street, Liverpool, L69 7ZB, UK (m.spencer@liverpool.ac.uk)

BUGS/TO DO/NOTES

11/7/08 Processor time wraps around and will give the wrong answer in long runs (wall clock time is OK).

11/7/08 -oQ (BFGS quasi-Newton optimizer) is much slower than -oS (simplex optimizer). This may be because the BFGS in GSL is inefficient. Newer GSL versions have bfgs2, which is supposed to be a much better implementation (but this hasn't made it into Debian stable yet).

11/7/08 Possible speedups include using Brent's method instead of golden section to optimize edges, and using analytical derivatives where these are available.

7/8/08 Edges of exactly zero can cause problems: a pair of edges of length zero with different character states will give zero likelihood. Set the minimum edge length to 1e-6 to avoid these.

8/8/08 If tree file and genome data file have different numbers of genomes, it will seg fault.

5/9/08 Optimizer sometimes takes bad steps that give NaN log likelihood. This seems to happen only with very large numbers of gene families, but detailed cause not yet investigated.

9/9/08 Compiling and linking in Cygwin with GSL 1.11-1 or Kubuntu Gutsy Gibbon with GSL 1.9-3 gives slightly different numerical results to Debian Etch with GSL 1.8-2. With all parameters fixed, the log likelihoods are identical (GSL used oly in the gamma distribution function). With no gamma rates and only edge length optimization, the results are identical (GSL not used at all). With parameter optimization, there are small differences (GSL optimizer used repeatedly). This strongly suggests dependence on the GSL version or BLAS library and accumulating differences in optimization. The differences are probably not large enough to matter (typically first decimal place of log likelihood, second decimal place of parameter estimates).

ACKNOWLEDGMENTS

Parts of the code are derived from fastdnaml by Gary J. Olsen, and seq-gen by Andrew Rambaut and Nick Grassly. See the start of ML_genedata.c for details.

This work was funded by BBSRC grant BB/E019005/1 to MS. Part of this work was carried out during the phylogenetics programme at the Isaac Newton Institute for Mathematical Sciences, Cambridge. We are grateful to David Bryant, Miklos Csuros, Brian Golding, Weilong Hao, Tal Pupko, Ed Susko and especially Ofir Cohen for ideas, test data sets and help with code development.
