I recently performed REL on an alignment and found that the non-positive discrete model had lower AIC than the unconstrained model. I ran both constrained and unconstrained versions of four-category "nonsynonymous only" and "dual variable" models. The nonsynonymous, constrained model performed best.

I then ran the alignment through FEL and found two sites under significant (p < 0.05) positive selection. These two sites (and one other) were also found positively selected (Bayes Factor > 50) in the marginals of the non-positive discrete model.

1. How could positive selection be inferred for any site in a model in which the highest dN is constrained to be 1?

2. I originally chose REL because of its ability to detect weak positive selection. It appears it can miss selection when it operates on very few sites, in that it "glosses over" the rate categories assigned to those sites. Is this a fair interpretation? My goal is to assess the strength of selection in different genes, and I'm unsure how to handle a REL that returns the non-positive discrete model as best fitting and FEL (and especially the same REL) returning positively selected sites.

Background: Alignments have 70-150 taxa, >200 sites.

Thanks!

Sarah

Dear Sarah,

1. This could be a bug in how the posterior probabilities are tabulated.
2. Your interpretation is correct: rate 'smoothing' is a big issue with REL models.

PARRIS (implemented on Datamonkey, and also in HyPhy, under Selection/Recombination) is more robust than REL with non-positive parameterization, because it is more numerically stable and better parameterized.

For comparing distribution of dN/dS between genes, have you tried dNdSdistributioncomparison.bf under Codon Selection Analyses?

Cheers,
Sergei

Hi, Sergei,

Thanks.

dNdSdistributioncomparison.bf is what I'm using for the second round of tests (i.e., to compare strengths of selection after I figure out what's positively selected using REL).

I'm looking at PARRIS. My genes don't recombine--would its advantage be that it fits discrete classes of omega instead of dN (it's not obvious how that would make the inference stronger), or is it something more subtle?

Thanks again.

Sarah

Hi Sarah,

PARRIS is a more natural parameterization for enforcing dN<=dS, because instead of parameterizing dN and dS separately (like REL does), it directly controls dN/dS ratios.

In terms of comparing selection between genes, dNdSdistributioncomparison.bf makes more sense than REL, because it specifically tackles the statistical comparison (e.g. is the proportion of sites with dN>dS is equal between two genes), whereas in REL it is secondary and ad hoc (e.g. gene 1 has more selected sites than gene 2).

HTH,
Sergei

Hi, Sergei,

Thanks again. I am running dNdSdistributioncomparison on the alignments that show positive selection to see whether they show interesting differences in the strength of positive selection.

I plan to do PARRIS on all the other alignments to confirm that they really don't show evidence of positive selection.

I'm thinking there might be a bug in the NPD nonsynonymous-variation-only model, and/or I don't understand the model's assumptions or output. The dual model, as expected, returns rate distributions with max dN/dS = 1.00. In the four alignments I've tested so far, the nonsynonymous NPD models return rate distributions containing at least one rate >>1. I though these models constrained the max dN (beta) to one, and that dS (alpha) is also one.

Background: I'm using a bash script that executes NPD and IDD models with (syn=4,nonsyn=4) categories.

I really appreciate your help. (I'm writing a dissertation chapter and trying to sort through results.)

Sarah