RELAX does not reproduce its own estimates at default settings, and its output does not show it

Abstract

Selection-intensity estimates from RELAX are reported as a point value of K with a likelihood-ratio P. We report that, at default settings and on data of ordinary size, the program does not reproduce its own fits. Of 27 enzyme entries refitted under two optimiser configurations, none reproduced its log-likelihood to within 0.01 units; the median change was 103 units, the largest over 3,400, and four verdicts reversed. Eighty null orthologues reproduced none. A byte-identical command returned a distinct likelihood on every repetition, single-threaded, across three releases, and on alignments simulated under the fitted model, where 3.3 per cent of replicates reproduced. The documented random-number seed never reaches the generator when assigned on the command line, yet reads back as the value supplied. PAML localises the cause: its two-ratio model, without site classes, reproduced its log-likelihood for all 288 genes; its site-class models agreed for 27 to 67 per cent. The instability follows the mixture over sites, not the program. The output does not show it: 46 of 410 fits ended with a negative likelihood-ratio statistic, impossible under convergence, and 123 of 410 report a K re-estimated under a domain restriction rather than the unconstrained maximum. Of 234 published studies using RELAX, none reported a seed. Seeding while holding the thread count at one reproduced sixty of sixty runs on twenty genes under two releases; the seed alone reproduced none of five, and no documentation states the second condition. We recommend that fits be repeated and their dispersion published.