New upstream version 4.9h+dfsg (4a2cb6d7) · Commits · Debian Med / paml

doc/pamlHistory.txt

+41 −1

Original line number	Diff line number	Diff line
		@@ -8,6 +8,46 @@ discussion site
		https://groups.google.com/forum/#!forum/pamlsoftware.


		Version 4.9h, March 2018

		(*) mcmctree: gamma-Dirichlet versus conditional i.i.d. priors for
		rates for loci. Since 4.9d, the program and the documentation are
		inconsistent about the two priors, and which value (0 or 1) means
		which prior. I have now checked the program and the documentation to
		make sure that they are consistent:

		prior = 0: gamma-Dirichlet (dos Reis 2014). This is the default.
		prior = 1: conditional i.i.d. prior (Zhu et al. 2015).

		I believe these two are similar especially if the number of loci
		(partitions) is large, but no serious comparisons between the two
		priors have been published.

		Thanks to Adnan Moussalli for pointing out the errors.

		(*) codeml. It was discovered that the mechanistic amino acid
		substitution model implemented in Yang et al. (1998; see table 3),
		specified by seqtype = 2 model = 6, has been broken for a long time,
		since version 3.0 (2000) at least. Version 2.0 (1999) seems to be
		correct. This means that the model become broken soon since it was
		published. I have now fixed this.

		This model of amino acid substitution starts from a Markov chain for
		codons and then aggregate the states and merge the synonymous codons
		into one state (the coded amino acid). This is an approximate
		formulation since the process after state aggregation is not Markovian
		anymore.

		I have now added another codon-based amino acid substitution model
		that treats amino acids as ambiguities codons. The model is specified
		by seqtype = 2 model = 5. This is an exact formulation.


		(*) codeml. The number of categories in the BEB calculation under M2
		and M8 is unintentionally set to 4 rather than 10. I have changed
		this back to 10. The details of this calculation are in Yang et
		al. 2005 MBE.


		Version 4.9g, December 2017

		@@ -16,7 +56,7 @@ Version 4.9g, December 2017
		messages like "strange: f[ 5] = -0.0587063 very small." This bug was
		introduced in version 4.9b and affects versions 4.9b-f. A different
		bug was introduced in version 4.9f that causes the log likelihood
		function under the site model M8 (NSsites = 8) to calculated
		function under the site model M8 (NSsites = 8) to be calculated
		incorrectly. These are now fixed.

examples/mtCDNA/AAcodon.result.txt

0 → 100644

+112 −0

Original line number	Diff line number	Diff line
		codeml results for fitting codon-based amino acid models.
		the data are 12 mt proteins from 7 apes, the small dataset used in yang et al. 1998, table 2 & 3.


		* seqfile = mtCDNApri.nuc * sequence data filename
		seqfile = mtCDNApri.aa * sequence data filename
		treefile = mtCDNApri.trees * tree file name

		(((human, (chimpanzee, bonobo)), gorilla), (orangutan, Sumatran), gibbon);

		********************
		CODONML
		Fequal (seqtype = 1 model = 0 CodonFreq = 0)
		TREE # 1: (((1, (2, 3)), 4), (5, 6), 7); MP score: 4299
		lnL(ntime: 11 np: 13): -31744.953377 +0.000000
		8..9 9..10 10..1 10..11 11..2 11..3 9..4 8..12 12..5 12..6 8..7
		0.120080 0.058139 0.175107 0.098210 0.072451 0.060123 0.206933 0.224103 0.115862 0.114939 0.408898 7.417478 0.112776
		SEs for parameters:
		0.009345 0.006777 0.009287 0.007425 0.005568 0.005205 0.010416 0.011610 0.008006 0.007893 0.015467 0.310582 0.003409

		Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).

		tree length = 1.654845

		(((1: 0.175107, (2: 0.072451, 3: 0.060123): 0.098210): 0.058139, 4: 0.206933): 0.120080, (5: 0.115862, 6: 0.114939): 0.224103, 7: 0.408898);

		(((human: 0.175107, (chimpanzee: 0.072451, bonobo: 0.060123): 0.098210): 0.058139, gorilla: 0.206933): 0.120080, (orangutan: 0.115862, Sumatran: 0.114939): 0.224103, gibbon: 0.408898);

		kappa (ts/tv) = 7.41748
		omega (dN/dS) = 0.11278


		********************
		CODONML
		F3x4 (seqtype = 1 model = 0 CodonFreq = 2)
		TREE # 1: (((1, (2, 3)), 4), (5, 6), 7); MP score: 4299
		check convergence..
		lnL(ntime: 11 np: 13): -29967.856102 +0.000000
		8..9 9..10 10..1 10..11 11..2 11..3 9..4 8..12 12..5 12..6 8..7
		0.261697 0.100524 0.242894 0.120349 0.078171 0.066639 0.286667 0.506947 0.158670 0.142655 0.829442 14.249448 0.041084
		SEs for parameters:
		0.023672 0.013801 0.014988 0.011874 0.006748 0.006435 0.018655 0.032072 0.013309 0.012734 0.043077 0.638686 0.001380

		Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).

		tree length = 2.794654

		(((1: 0.242894, (2: 0.078171, 3: 0.066639): 0.120349): 0.100524, 4: 0.286667): 0.261697, (5: 0.158670, 6: 0.142655): 0.506947, 7: 0.829442);

		(((human: 0.242894, (chimpanzee: 0.078171, bonobo: 0.066639): 0.120349): 0.100524, gorilla: 0.286667): 0.261697, (orangutan: 0.158670, Sumatran: 0.142655): 0.506947, gibbon: 0.829442);

		kappa (ts/tv) = 14.24945
		omega (dN/dS) = 0.04108


		******************************************************
		AAML JTT
		(seqtype = 2 model = 3)
		TREE # 1: (((1, (2, 3)), 4), (5, 6), 7); MP score: 895
		lnL(ntime: 11 np: 11): -14717.981418 +0.000000
		8..9 9..10 10..1 10..11 11..2 11..3 9..4 8..12 12..5 12..6 8..7
		0.025579 0.009654 0.022079 0.012328 0.011335 0.010281 0.026980 0.052187 0.026525 0.019651 0.062350
		SEs for parameters:
		0.003197 0.002018 0.002777 0.002097 0.001932 0.001865 0.003193 0.004412 0.003071 0.002679 0.004757

		tree length = 0.278951

		(((1: 0.022079, (2: 0.011335, 3: 0.010281): 0.012328): 0.009654, 4: 0.026980): 0.025579, (5: 0.026525, 6: 0.019651): 0.052187, 7: 0.062350);

		(((human: 0.022079, (chimpanzee: 0.011335, bonobo: 0.010281): 0.012328): 0.009654, gorilla: 0.026980): 0.025579, (orangutan: 0.026525, Sumatran: 0.019651): 0.052187, gibbon: 0.062350);



		*********************
		AA-CODONML, Old model FromCodon
		(seqtype = 2 model = 6)
		TREE # 1: (((1, (2, 3)), 4), (5, 6), 7); MP score: 895
		lnL(ntime: 11 np: 12): -14718.224885 +0.000000
		8..9 9..10 10..1 10..11 11..2 11..3 9..4 8..12 12..5 12..6 8..7
		0.027189 0.010346 0.023115 0.012595 0.011879 0.010512 0.028689 0.056919 0.026419 0.021490 0.070764 9.156815
		SEs for parameters:
		0.003186 0.002070 0.002831 0.002126 0.001994 0.001832 0.003110 0.004354 0.003100 0.002754 0.004703 0.710971

		tree length = 0.299918

		(((1: 0.023115, (2: 0.011879, 3: 0.010512): 0.012595): 0.010346, 4: 0.028689): 0.027189, (5: 0.026419, 6: 0.021490): 0.056919, 7: 0.070764);

		(((human: 0.023115, (chimpanzee: 0.011879, bonobo: 0.010512): 0.012595): 0.010346, gorilla: 0.028689): 0.027189, (orangutan: 0.026419, Sumatran: 0.021490): 0.056919, gibbon: 0.070764);

		kappa (ts/tv) = 9.15682


		*********************
		AA-CODONML, new model FromCodon0
		(seqtype = 2 model = 5)
		TREE # 1: (((1, (2, 3)), 4), (5, 6), 7); MP score: -1
		lnL(ntime: 11 np: 13): -14707.663779 +0.000000
		8..9 9..10 10..1 10..11 11..2 11..3 9..4 8..12 12..5 12..6 8..7
		0.636398 0.245832 0.527920 0.290247 0.270174 0.240427 0.658706 1.320827 0.602244 0.493482 1.621211 9.246897 0.031208
		SEs for parameters:
		0.321763 0.132591 0.271862 0.153691 0.142901 0.127706 0.332545 0.664754 0.316022 0.256229 0.821216 0.716647 0.016303

		Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).

		tree length = 6.907467

		(((1: 0.527920, (2: 0.270174, 3: 0.240427): 0.290247): 0.245832, 4: 0.658706): 0.636398, (5: 0.602244, 6: 0.493482): 1.320827, 7: 1.621211);

		(((human: 0.527920, (chimpanzee: 0.270174, bonobo: 0.240427): 0.290247): 0.245832, gorilla: 0.658706): 0.636398, (orangutan: 0.602244, Sumatran: 0.493482): 1.320827, gibbon: 1.621211);

		kappa (ts/tv) = 9.24690
		omega (dN/dS) = 0.03121

examples/mtCDNA/codeml.AAcodon.ctl

0 → 100644

+69 −0

Original line number	Diff line number	Diff line
		* seqfile = mtCDNApri.nuc * sequence data filename
		seqfile = mtCDNApri.aa * sequence data filename
		treefile = mtCDNApri.trees * tree file name

		outfile = mlc * main result file name

		noisy = 3 * 0,1,2,3,9: how much rubbish on the screen
		verbose = 1 * 0: concise; 1: detailed, 2: too much
		runmode = 0 * 0: user tree; 1: semi-automatic; 2: automatic
		* 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise

		seqtype = 2 * 1:codons; 2:AAs; 3:codons-->AAs
		CodonFreq = 0 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
		* 4:F1x4MG, 5:F3x4MG, 6:FMutSel0, 7:FMutSel

		* estFreq = 0
		* hkyREV = 0 * 0: HKY-like; 1: GTR(REV)-like

		* ndata = 1 * number of data sets or loci
		* bootstrap = 0 * generate bootstrap data sets

		model = 5
		* models for codons:
		* 0:one, 1:b, 2:2 or more dN/dS ratios for branches, 6:FromCodon
		* models for AAs or codon-translated AAs:
		* 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
		* 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)

		NSsites = 0 * 23 24 25 26 * 23 24 25 26 * 0:one w; 1:NearlyNeutral; 2:PositiveSelection; 3:discrete;
		* 4:freqs; 5:gamma; 6:2gamma; 7:beta; 8:beta&w+; 9:betaγ
		* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
		* 13:3normal>0;
		* 22:M2a_Old(M2a_rel);
		* 23:Tgamma; 24:Tinvgamma; 25:Tgamma+1; 26:Tinvgamma+1.

		clock = 0 * 0:no clock, 1:global clock; 2:local clock
		aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
		aaRatefile = jones.dat * for aa seqs under model = 3 (empirical+F)
		* dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own

		icode = 1 * 0:universal code; 1:mammalian mt; 2-10:see below
		Mgene = 0
		* codon: 0:rates, 1:separate; 2:diff pi, 3:diff kappa, 4:all diff
		* AA: 0:rates, 1:separate

		* NShmm = 0 * 1: hidden Markov model

		fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
		kappa = 3 * initial or fixed kappa
		fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
		omega = 1.5 * initial or fIf yoixed omega, for codons or codon-based AAs

		fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
		alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
		Malpha = 0 * different alphas for genes
		ncatG = 10 * # of categories in dG of NSsites models

		getSE = 1 * 0: don't want them, 1: want S.E.s of estimates
		RateAncestor = 1 * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)
		Small_Diff = 1e-8
		* cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)?
		* fix_blength = 1 * 0: ignore, -1: random, 1: initial, 2: fixed
		* method = 0 * Optimization method 0: simultaneous; 1: one branch a time

		* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
		* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt.,
		* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt.,
		* 10: blepharisma nu., 11: Yang's regularized code
		* These codes correspond to transl_table 1 to 11 of GenBank.

examples/mtCDNA/mtCDNApri.trees

+3 −0

Original line number	Diff line number	Diff line
		7 1
		(((human, (chimpanzee, bonobo)), gorilla), (orangutan, Sumatran), gibbon);

		// end of file

		(((human: 0.254974, (chimpanzee: 0.078779, bonobo: 0.067071): 0.120425): 0.104179, gorilla: 0.302845): 0.289306, (orangutan: 0.168115, Sumatran: 0.141455): 0.573030, gibbon: 0.930806); [unrooted tree, M0 F3x4 branch lengths]

src/baseml.c

+3 −3

Original line number	Diff line number	Diff line
		@@ -83,7 +83,7 @@ struct TREEN {
		} nodes, gnodes, nodes_t[2 NS - 1];


		/* for sptree.nodes[].fossil: lower, upper, bounds, gamma, inverse-gamma */
		/* for stree.nodes[].fossil: lower, upper, bounds, gamma, inverse-gamma */
		enum { LOWER_F = 1, UPPER_F, BOUND_F } FOSSIL_FLAGS;
		char *fossils[] = { " ", "L", "U", "B" };

		@@ -95,7 +95,7 @@ struct SPECIESTREE {
		double age, pfossil[7]; /* lower and upper bounds or alpha & beta */
		double lnrates; / log rates for loci */
		} nodes[2 * NS - 1];
		} sptree;
		} stree;
		/* all trees are binary & rooted, with ancestors unknown. */

		struct DATA { /* locus-specific data and tree information */
		@@ -203,7 +203,7 @@ int main(int argc, char *argv[])
		error2("oom blengths0");
		}

		SetMapAmbiguity();
		SetMapAmbiguity(com.seqtype, 0);

		/* AllPatterns(fout); */