function [pID,pN] = FDR(p,q)
% FORMAT pt = FDR(p,q)
% 
% p   - vector of p-values
% q   - False Discovery Rate level
%
% pID - p-value threshold based on independence or positive dependence
% pN  - Nonparametric p-value threshold
%______________________________________________________________________________
% @(#)FDR.m	1.3 Tom Nichols 02/01/18
% Uppdated to be used in DeltaProt 2.0 Steinar Thorvaldsen Oct. 7, 2009.

p = sort(abs(p(:))); %in DeltaProt the direction of change is often stored in the sign of P-value
%p = p(~isnan(p)); %remove NaNs from a vector p.
V = length(p);
I = (1:V)';

cVID = 1;
cVN = sum(1./(1:V));

pID = p(max(find(p<=I/V*q/cVID)));
pN = p(max(find(p<=I/V*q/cVN)));