Could you help me with interpretation of codeml results?

[František Zedek]

Hi PAML fans,

I have duplicated gene that has two functional domains, A and B. I used site models to test for positive selection. When I analyzed full-length gene, LRT was not significant either for M1 vs M2 or M7 vs M8 (F3x4, F61 was the same), but it was not zero and BEB identified site(s) in "domain B" that I had suspected to evolve adaptively. When I analyzed the domain B separately, I got significant LRT for both M1/M2 and M7/M8 with F3x4 option, but only M7/M8 was significant with F61 option. BEB identified suspected sites with high posterior probability. Sum of dS over all branches was always around 0.7. 

I also ran branch model for "domain B" only. I added an outgroup. Since I had no apriori hypothesis about which of the paralogs should evolve under positive selection, I ran free-ratio model. It identified several branches under positive selection (under both F3x4 and F61) and LRT remained significant even after Bonferroni correction.

Taken the results of site and branch models into account, can I conclude that my gene of interest (or the domain B at least) has been evolving under positive selection and that the selection pressure was rather weak?

Thank you for any any advice!

Franta

[Ziheng]

there is evidence for positive selection, but the evidence is not extremely strong and as a result you get significant results only in some analyses but not in all of them. Nite that in all cases there is evidence, but the evidence reached statistical significance only in some analyses.
if you have other sources of information, for example, structural or experimental results, you can try to integrate them.
Or if you can add more sequences, perhaps the power of the test may improve.
The test is known to lack power, so that is one interpretation.
Note that this is about the strength of the evidence for positive selection, not the strength of positive selection. The latter is reflected in the estimates of p2 and w2, the proportion of sites under positive selection and also the dn/ds ratio.
Ziheng

[František Zedek]

Dear prof Yang,

many thanks for your reply. Unfortunately, so far I have no more sequences to add to improve the power of the test. However, the sites that have been identified by BEB lie in the loop region that directly interacts with DNA (the gene is centromeric histone H3), which is consistent with previous experiments and structural analyes.

I am not sure how can I decide whether the selection is strong or weak. I mean, is there any justifiable treshold? In my analyses full-length protein, the model M8 showed, p1=0.02 and w=6.39 under F3x4 and p1=0.16 and w=2.58 under F61. When just the histone fold domain is analyzed p1=0.04 and w=5.56 under F3x4 and p1=0.08 and w=2.98 under F61. Very similar results were obtained with p2 and w2 in M2a model. It seems to me as weak positive selection, but I have no justification for it except a feeling that the numbers are low. Would you be so kind and show me a direction or perhaps just recommend some papers to read on the subject of selection strength?

Best wishes,

Franta

[Ziheng]

the tests are designed to test whether there is positive selection, not to measure the strength of positive selection. i think typically the first question is more important, put in another way, it is hard just to decide whether there is positive selection, not to mention deciding on how strong the positive selection is.
one tentative measure we used was p_s * w_s, the proportion of sites under positive selection multiplied by the w ratio for the site class. Those two parameters are expected to be strongly negatively correlated, and that is certainly the case in your data as you sometimes get a large p and small w and sometimes a smaller p and larger w.
ziheng