Question: PAML codeML 4.8 wont give maximum likelihood nor omega estimates
0
gravatar for mtd66
12 months ago by
mtd660
mtd660 wrote:

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γ
               * 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 • 406 views
ADD COMMENTlink written 12 months ago by mtd660

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.

ADD REPLYlink written 12 months ago by Brice Sarver3.5k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1563 users visited in the last hour