Hi Biostars

I'm doing positive selection analysis on a relatively large set of orthologous genes using PAML's implementation of the branch-site model. I've got the results of my likelihood ratio test (LRT) for each gene, and according to this (1) and other papers, the branch-site model follows a 1:1 mixture of chi^2 and point mass 0 distribution, so the critical values for 5% and 1% are 2.71 and 5.41 respectively.

My problem is that I need to perform a Benjamini-Hochberg correction, and for that I need the specific p-value for each gene, not just to decide whether it passes either critical value.

How would I go about calculating the actual p-value from each LRT result and the number of degrees of freedom? I can write a script in python/perl (maybe R if necessary) if I know how to actually calculate the p-values. Any help would be greatly appreciated!

(1) Statistical properties of the branch-site test of positive selection http://www.ncbi.nlm.nih.gov/pubmed/21087944