Hi Rosemarie,

I can offer only a few quick suggestions from my own experiences with PAML.

- 4 species is on the smaller scale. Can you mine other plant databases like phytozome or GenBank to build a larger tree?

- Your alignment looks wonky; lots of indels, and indels in every species. The latter should be removed, and collapsed into a single concatenated sequence. See here from PAML's FAQ: "Alignment gaps are more difficult. PAML does not have any methods of dealing with them

properly. The two options are (1) to remove them, which you could do by manually removing

them or by choosing cleandata = 1; and (2) to treat them as ambiguity (undetermined)

nucleotides or amino acids. The latter is the default behavior if gaps are in the data. Neither is

ideal. One obvious effect is that both strategies under-estimate sequence divergences, while the

effects on other analyses might not be so clear. Personally I think sites at which most sequences

have data except for one or two sequences should perhaps be kept while sites at which all

sequences except one or two have alignment gaps had better be removed."

- On the topic of sequence alignment, are you sure you are comparing the same sequence across your taxa? Could there be a paralogous sequence in there? This could cause elevated substitution rates.

- Have you tried aligning with the protein sequence, and then back translating with the mRNA sequence to fill in the aa?

I generally think the data is fine. I was able to translate your codon alignments and based on blastp, they are all homologous. Of course, deeper divergences will generally yield more alignment uncertainty, but that's just the way it is. If they are paralogs, well you might want to investigate this. The synonymous evolutionary distance of Anoph is quite large compared to the rest of the tree, you might want to check this sequence specifically.

Adam raises some very good points concerning more sequence data. More data will give you more power, especially for testing individual sites. If you do not have >10 sequences, you just don't have much power for detecting specific sites under positive selection (Murrell et al, 2010 - considering site classes are considered fixed effects). You might want to consider testing if a proportion of sites are under selection on specific branches, this is the Nielsen and Yang Branch-Site test and the paml manual gives instructions on how to configure the ctl file for this test.

I've attached a zip file of 5 different branch tests to show how to constrain dN/dS on different branches to generate nested models allowing application of the likelihood ratio test (LRT). You should be able to see that allowing Albo to have a separate dN/dS estimate (relaxed in this case) from all other branches is the preferred model in this set of tests. The LRT is calculated by comparing this model (H2) to the null model. For a nice step by step tutorial on these perform these branch tests and others, see Joe Bielawski's PAML tutorial at: http://myweb.dal.ca/js551958/PAML_lab/lab.html

PAML does use ML to estimate dN/dS. It can do everything that is appropriate for the scope of this data. HYPHY can implement some interesting tests and has some computational benefits, but PAML can be more user friendly for implementing these standard tests.

I don't know what the tests are in the other papers you mentioned, but based on the distribution of sites, these are the site specific tests for positive selection using the models Yang named M0, M1a, M2a, M3, M7, and M8). This is covered in Bielawski's tutorial, but remember, more data means more power for these tests. I would not be surprised if you did not detect any sites under positive selection with just 4 taxa. You'll probably want to try and get more sequences from genbank. Otherwise, do it anyway and see what happens!

I hope this helps!

HI Adam and George in response to your questions:

- 4 species is on the smaller scale. If you do not have >10 sequences, you just don't have much power for detecting specific sites under positive selection (Murrell et al, 2010 - considering site classes are considered fixed effects). 

I've mined Genbank, and the genus (despite having a lot of species in it) only has these two as having any of the immunity genes I'm interested in.

- You might want to consider testing if a proportion of sites are under selection on specific branches, this is the Nielsen and Yang Branch-Site test and the paml manual gives instructions on how to configure the ctl file for this test.

I'll have a look through the paml manual and the nice documentation that is around. I'm very good at following instructions.

- Your alignment looks wonky; lots of indels, and indels in every species. The latter should be removed, and collapsed into a single concatenated sequence.

I've actually deleted a lot of sequence data already, which was even worse in terms of how many indels there were. I guess I left these in as past experience gave me more accurate tree distances with PhyML.

When I'm only working with 4 sequences, it's hard to just keep the sites that have 3+ sequences on offer.

- On the topic of sequence alignment, are you sure you are comparing the same sequence across your taxa? Could there be a paralogous sequence in there? This could cause elevated substitution rates.

- Have you tried aligning with the protein sequence, and then back translating with the mRNA sequence to fill in the aa?

I translated the original mRNA sequences, aligned at the aa level using Se-Al, and then went back to nt. I'm certain that I am comparing the same sequence across the taxa, but it is possible they are too different (they are all mosquito sequences, but they differ in how well the genus has been probed into).

I've attached a zip file of 5 different branch tests to show how to constrain dN/dS on different branches to generate nested models allowing application of the likelihood ratio test (LRT).

Yes! This what was I had done before. Phew. I did look at that tutorial a while back in my searching (I distinctly remember reading about 'more beer time' and despairing) - I just didn't know that the results I was getting were reasonable, or even answering one of the questions I was asking.

I also found that Albo did have a different omega to the rest, but I wasn't sure if I had asked enough questions of the surrounding branches (where the branch in question has a different dn/ds to all of the other branches, including the one it is most closely related to).

In one of the other cases, the df ended up being 0, which I assumed couldn't have happened. On other documentation though, the df is said to always be 1. I thought it would depend on how many cases you were comparing.

HYPHY can implement some interesting tests and has some computational benefits, but PAML can be more user friendly for implementing these standard tests.

I might leave HYPHY well enough alone then. Arguing with just one program is enough for me.

I don't know what the tests are in the other papers you mentioned, but based on the distribution of sites, these are the site specific tests for positive selection using the models Yang named M0, M1a, M2a, M3, M7, and M8).

Yes, these were the site specific tests suggested. I didn't realise however that different tree files were needed in each case (duh, of course they have to be different MLs)). 

Thanks so much for your help. I was so relieved to see that someone had answered my question. I have a meeting with my supervisor this afternoon - perhaps she'll have a good idea for fishing more sequence data from other contacts in the field.

Hmm immunity genes. These seem like they would be fast evolving; almost like an arms race between disease and ways to combat the disease.

How are you certain that they are orthologous? Your particular gene wouldn't be part of a larger gene family? You could include other homologous/paralogous members of this gene family in your tree, and other gene members in other closely related genera. I might be wrong, but there are a lot of mosquito genomes sequenced and annotated? You could run OrthoMCL to find putative families that house your gene(s) of interest, and paralogs. You might have something similar to what is observed in primates and ordorant receptors. See: http://www.pnas.org/content/106/50/21247.full.pdf+html.

With regards to df=0, some models are NOT comparable or nested. Specifically, the LRT is performed between (H1 and null), (H2 and null), (H3 and null), (H4 and H2), (H5 and H4). This keeps the df =1 for each test. We only need to compare H4 to H2, since H2 is our preferred 2 free-parameter model and the best interpretation. For example, comparing H4 to H1 would be misleading since it would give a significant p-value by LRT, not because the 3 free parameters are preferred, but only because Anoph is allowed to vary in this case. Comparing H4 to H2 only makes it clear that giving Anoph its own omega is what improves the likelihood. 

You can set up tests where df > 1, but since it sounds like you will be looking at many genes, try to keep your tests with df =1 so they are all approximately distributed as Chi2(1) and apply some multiple testing correction to the distribution of tests (Generally correcting for the FDR (Benjamini-Hochberg) or positive FDR (Storres)).

With regards to the site-specific models (M0 - M8), please refer to Joe's tutorial. He explains things far better than I could hope to and I suspect better than most.

If your supervisor can help you get more data, that is the best scenario. Otherwise, it sounds like your data is in order and you have an interesting question working with immunity genes; just set up the tests that are biologically interesting, follow the stats, and see if your story emerges.

Using the tool to plot transitions and transversions vs F84 distances in each pairwise comparison, it is clear that your data is far beyond saturation of silent sites with mutations.  The cactus gene/protein you have aligned here is also found in other insects and in vertebrates too.  I made an alignment of the catus protein (it is not named cactus in the vertebrates) and built a quick neighbor-joining tree from the protein sequences, just to put the mosquito data in perspective.  Mosquitoes shared a common ancestor with Drosophila longer ago than mammals shared a common ancestor with birds, according to this tree (or the insect gene/protein evolved at a faster rate, but branch lengths being nearly equal for all lineages argues against that).

My sequences have the GenBank accession numbers in the labels, in case you want to go get the DNA sequences instead of the protein sequences.

I am confused by you referring to cactus as an immunity gene/protein because it seems to be a dorsal/ventral development regulator in most of the literature about it.  There does only seem to be one single copy of this gene in each organism's genome, my BLAST searching did not for example hit several related genes/proteins in any one insect or vertebrate genome.  

The random effects likelihood (REL) test in HYPHY/datamonkey.org is equivalent to the Nielson and Yang model implemented in PAML.  I have used both and find that HYPHY is a lot more user friendly (e.g., no need to setup Ctl files for null and alternative hypothesis and then compare them using likelihood ratio tests).  I suppose it also depends on whether you are using the GUI or command line version of PAML.

Thanks everyone! My supervisor ended up changing the project, and so I am working on Drosophila genes now. There is such a lovely lot of them to choose from, it's making working with Paml a lot easier. It couldn't have been done without your help.