Tips on Running CODEML on 17,000 genes?

[Iris Lawson-Ryan]

Hello, 

I am currently trying to run CODEML to find individual branch omega values and it seems to be taking a very long time to get through even one gene, let alone all 17,000 genes. 

When I start the analysis, it will get through to this section (see below, Species name changed for data embargo reasons) very quickly but then seems to stall there and not make any more progress. I have been running this gene for nearly 24 hours and it has not moved past this point. 

I will be calculating rtt omega after codeml is finished running so I need to run this as model = 1. 

Nei & Gojobori 1986. dN/dS (dN, dS) (Note: This matrix is not used in later ML. analysis. Use runmode = -2 for ML pairwise comparison.) species 1 species2 -1.0000 (0.0000 0.2854) species3 0.0164 (0.0037 0.2239) 0.0127 (0.0037 0.2896) species4 -1.0000 (0.0000 0.3099)-1.0000 (0.0000 0.3118) 0.0148 (0.0037 0.2474) species5 -1.0000 (0.0000 0.1998)-1.0000 (0.0000 0.2744) 0.0172 (0.0037 0.2139)-1.0000 (0.0000 0.2760) species6 -1.0000 (0.0000 0.2101)-1.0000 (0.0000 0.1912) 0.0211 (0.0037 0.1748)-1.0000 (0.0000 0.2022)-1.0000 (0.0000 0.2001) species7 0.0055 (0.0024 0.4472) 0.0050 (0.0024 0.4933) 0.0177 (0.0061 0.3478) 0.0048 (0.0024 0.5117) 0.0060 (0.0024 0.4068) 0.0074 (0.0024 0.3322) species8 -1.0000 (0.0000 0.2840)-1.0000 (0.0000 0.4340) 0.0106 (0.0037 0.3462)-1.0000 (0.0000 0.3708)-1.0000 (0.0000 0.3068)-1.0000 (0.0000 0.2959) 0.0053 (0.0024 0.4616) species9 -1.0000 (0.0000 0.2307)-1.0000 (0.0000 0.3089) 0.0123 (0.0037 0.3000)-1.0000 (0.0000 0.3225)-1.0000 (0.0000 0.2846)-1.0000 (0.0000 0.2005) 0.0069 (0.0024 0.3566)-1.0000 (0.0000 0.3308) species10 -1.0000 (0.0000 0.2189)-1.0000 (0.0000 0.2304) 0.0212 (0.0037 0.1737)-1.0000 (0.0000 0.2213)-1.0000 (0.0000 0.1891)-1.0000 (0.0000 0.1238) 0.0067 (0.0024 0.3664)-1.0000 (0.0000 0.2503)-1.0000 (0.0000 0.2194) species11 -1.0000 (0.0000 0.2504)-1.0000 (0.0000 0.2411) 0.0191 (0.0037 0.1933)-1.0000 (0.0000 0.2319)-1.0000 (0.0000 0.2092)-1.0000 (0.0000

Do you have any suggestions on how to speed this up? Since I have 17,000 genes to run through my pipeline I would like to be as efficient as possible!

Any thoughts would be helpful!

[Sandra AC]

Hi Iris,

Thanks for reaching out!

One possibility would be to run a phylogenetic software to infer the branch lengths of the tree topology you are fixing (e.g., IQ-TREE, RAxML-NG). In that way, CODEML will not spend time estimating these branch lengths and will just focus on estimating the other model parameters.

If you want to try this approach, you need to include "fix_blength = 2" in your control file. When setting variable "fix_blength" to "2", you are telling CODEML to use the branch lengths specified in the input tree file. Please note that the input tree file must be in PHYLIP format. E.g.: the first line must be a header where you specify the number of taxa and then the number of trees in the file; "703 1" would be a header for a 703-taxon tree ("1" refers to having only one tree in the file). Once you have inferred a tree topology with branch lengths with your preferred software, you then incorporate the Newick tree after the header (please remove other labels that may remain such as bootstrap support values; keep only the topology with the branch lengths). While fixing the branch lengths will reduce computational time, please note that the free-ratios model ("model = 1") is very parameter-rich, so it may still take a while to estimate the omega ratio for each branch. To avoid wasting computational resources, you may want to test this approach first with one gene. If this method works, you can then apply this to the rest of the genes you are analysing :)

Hope this helps!

All the best,
Sandy

[Sandra AC]

P.S.: Another possibility would be to run CODEML using model M0 ("NSsites = 0" + "model = 0") with "method = 1" to estimate the MLEs of the branch lengths (this run should be somewhat quick). Then, you can use the estimated branch lengths in your fixed tree topology as starting values ("fix_blength = 1") when running the branch model. You should keep the same settings for "CodonFreq" and "estFreq" in the control file in both runs. If you decided to try what I suggested in my previous message, you could check whether the branch lengths estimated by CODEML using this approach and those estimated by the phylogenetic software you chose are more or less similar.

All the best,
Sandy