PAML codeML 4.8 wont give maximum likelihood nor omega estimates
0
0
Entering edit mode
21 months ago
mtd66 • 0

I am using codeML 4.8 from PAML and while I can get it to run without giving me error messages it only gives me a matrix with what I assume contains dn/ds ratios between individuals.

Here is my code.ctl file:

  seqfile = temp2  * sequence data filename
treefile = testtree      * tree structure file name
outfile = mlc           * main result file name

noisy = 0  * 0,1,2,3,9: how much rubbish on the screen
verbose = 2  * 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 = 1  * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2  * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table

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

model = 2
* models for codons:
* 0:one, 1:b, 2:2 or more dN/dS ratios for branches
* 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  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
* 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
* 13:3normal>0

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

fix_kappa = 0  * 1: kappa fixed, 0: kappa to be estimated
kappa = 2  * initial or fixed kappa
fix_omega = 0  * 1: omega or omega_1 fixed, 0: estimate
omega = 0.5 * initial or fixed 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 = 8  * # of categories in dG of NSsites models

getSE = 1  * 0: don't want them, 1: want S.E.s of estimates
RateAncestor = 0  * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)

Small_Diff = .45e-6
cleandata = 0  * remove sites with ambiguity data (1:yes, 0:no)?
*  fix_blength = 0  * 0: ignore, -1: random, 1: initial, 2: fixed
method = 0  * Optimization method 0: simultaneous; 1: one branch a time


my "testtree" file looks like this:

((zzzzz1mAF),(zzzzz3mAF))


and my sequence looks like:

>zzzzz1mAF
ATGCATTCCCGTGAAAGTTCAGGGGGAGGCTGGGAAGGGGGCCCTCGCAAGTACGCCTAC
AATCGCGAGGGAGCCCGGGAACGAGAGATGCACCTCTATCTTCCATACTGGGGGTTTTCG
GATTCCGAGAAACATGCCCAGCAATCGATTGGAACCCATCCGAATGTGTCACGCGACAAT
GTCTTGACAGTCTTTGCGCAACTTGCTGCACTACGAATGGGCGCTCAAAGAGCCCTGATA
TCGCTTTTTGATAAGACCATGCAACATGTCGTAGCCGAATCAACCCCGGGACTAAGCTTG
>zzzzz3mAF
ATGCATTCCCGTGAAAGTTCAGGGGGAGGCTGGGAAGGGGGCCCTCGCAAGTACGCCTAC
AATCCCGAGGGAGCCCGGGAACGAGAGATGCACCTCTATCTTCCATACTGGGGGTTTTCG
GATTCCGAGAAACATGCCCAGCAATCGATTGGAACCCATCCGAATGTGTCACGCGACAAT
GTCTTGACAGTCTTTGCGCAACTTGCTGCACTACGAATGGGTGCTCAAAGAGCCCTGATA
TCACTTTTTGATAAGACCATGCAACATGTCGTTGCCGAAGCAACCCCGGGACTAAGCTTG


For the sake of simplicity I have narrowed this down to two sequences in two groups. In reality I have three groupings with a 10-15 individuals in each. No matter what I try I cant get it to give me more information than described above.

paml alignment dn/ds • 668 views
0
Entering edit mode

1) Do you have enough variability in your sequences? 2) Are you providing enough sequences at an appropriate tree height to expect to recover anything? 3) You're allowing more than one ratio among branches (model = 2) yet assuming a single rate among sites (NSsites = 0). Is this what you're trying to do? If you're looking for omega estimates for sites under positive selection, do you want to vary your site models and select among them (e.g., using a likelihood ratio test for nested models)? This is outlined in the manual.