-
Notifications
You must be signed in to change notification settings - Fork 7
/
GCcorrect_logpvec.m
39 lines (29 loc) · 1.27 KB
/
GCcorrect_logpvec.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
function [logpvec_gc, sig0] = GCcorrect_logpvec(logpvec,ivec0,use_standard_gc,pruneidx)
% GCCORRECT_LOGPVEC Perform a correction for genomic control on -log(p)
%
% Reference: Devlin, B., and Kathryn Roeder. "Genomic control for
% association studies." Biometrics 55.4 (1999): 997-1004.
if ~islogical(ivec0), error('ivec0 must be a logical vector of control elements'); end
if ~exist('use_standard_gc', 'var'), use_standard_gc = false; end
if ~exist('pruneidx', 'var'), pruneidx = []; end
if isempty(pruneidx), pruneidx = true(size(logpvec)); end
randprune_n = size(pruneidx, 2);
sig0 = nan(1, randprune_n);
for randprune_i = 1:randprune_n
% Get median
logpvec0 = logpvec(ivec0 & pruneidx(:, randprune_i));
zvec0 = abs(norminv(1/2*10 .^ -logpvec0));
if use_standard_gc
sig0(randprune_i) = sqrt(median(zvec0(~isnan(zvec0)).^2) ./ chi2inv(0.5,1));
else
fracvec = 1 - logspace(log10(0.5), -2, 50);
prc = prctile(zvec0.^2, 100*fracvec);
prc = reshape(prc, 1, length(prc)); % octave-fix: convert prc nx1 to 1xn
sig0vec = sqrt(prc ./ chi2inv(fracvec, 1));
sig0(randprune_i) = median(sig0vec);
end
end
% Normalize by median
zvec = abs(norminv(1/2*10 .^ -logpvec));
logpvec_gc = -log10(2*normcdf(-abs(zvec/median(sig0))));
end