Hi Sergei,

I am running AnalyzeCodonData.bf using the MG94W9 codon model with local parameters and then running post_sns.bf on the results.  I believe that there is a bug in the dN and dS branch length calculations.  If an Esyn branch value is in units of synonymous substitutions per nucleotide site (as reported in the output), then we should be able to convert this to a dS branch length by multiplying Esyn by the number of nucleotide sites and then dividing by the number of synonymous sites.  Right?  And similarly for dN values.

Esyn [syn subs/nt site] x num nt sites / num syn sites  = dS [syn subs/syn site]

Ensyn [nsyn subs/nt site] x num nt sites / num nsyn sites = dN [nsyn subs/nsyn site]


However, the dS and dN values are 3 times smaller than expected based on these formulas (sample output posted below).  I originally thought that the Esyn and Ensyn values were actually in units of subs per codon rather than subs per nt site as reported in the post_sns.bf output.  But comparison with PAML dN and dS results suggests that HYPHY's dN and dS values are being divided by three incorrectly.  Can you shed some insight?  


Thanks!
Jeff



Here's the sample output:
-------------------------------
Total nucleotide sites :2286
Synonymous  sites      :527.448
Nonsynonymous  sites   :1664.93

E[Syn subs/nucleotide site] tree:
     (spec1:0.391003,spec2:0.584268,spec3:0.355384)

E[Non-syn subs/nucleotide site] tree:
     (spec1:0.132946,spec2:0.0815103,spec3:0.0653267)

dS tree:
     (spec1:0.564879,spec2:0.844087,spec3:0.513421)

dN tree:
     (spec1:0.0608464,spec2:0.0373055,spec3:0.0298986)

Dear Jeff,

You are quite right; there was an incorrect factor of 1/3 incorporated into the output! Thanks so much for pointing this out. I'll put the fix in the next build.
In the meantime, change lines 87 and 88 in post_sns.bf to read:

[code]
sSites  = 3*filteredData.sites/sSites;
nsSites = 3*filteredData.sites/nsSites;
[/code]

Cheers,
Sergei

Great, thanks for confirming, and glad I could help!

Jeff