%calculateGWASpowerR2 Computes power and PGS R-squared in a meta-GWAS. % % Author: R. de Vlaming % Date: May 2, 2016. % % [dP,dR2] = calculateGWASpowerR2(iC,vN,vH2,mCGR,iM,iS) returns meta-GWAS % power per trait-associated SNP and expected PGS R-squared in a hold-out % sample with SNP-weightes based on the meta-analysis results. % % Meta-analysis of GWAS results from iC studies with sample size per % study specified in vector vN (iC-dimensional) and SNP-heritability per % study (including hold-out sample) in vector vH2 ((iC+1)-dimensional). % Cross-study genetic correlations are specified in (iC+1)-by-(iC+1) % matrix mCGR. The effective number of causal SNPs is given by iM and the % total effective number of SNPs by iS. % % Example: Compute (i) power of a meta-analysis of two studies with 40k % and 60k individuals and SNP-heritability of 60% and 40% respectively, % for a trait with 1k independent causal SNPs, with a genetic correlation % of 0.9 between the meta-analysis studies, assuming 100k independent % SNPs in total in the genome, and (ii) expected PGS R-square in the % hold-out sample, with heritability 50% and a genetic correlation of 0.8 % with study 1 and 0.7 with study 2, using SNP-weights from the % meta-analysis to construct the PGS. % % mCGR = [1.0 0.9 0.8 ;... % 0.9 1.0 0.7 ;... % 0.8 0.7 1.0]; % [dP,dR2] = calculateGWASpowerR2(2,[40000 60000],[0.6 0.4 0.5],... % mCGR,1000,100000); % function [dP,dR2] = calculateGWASpowerR2(iC,vN,vH2,mCGR,iM,iS) % determine threshold for genome-wide significant Z-score dPval = 5*(10^(-8)); dZhit = norminv(dPval/2,0,1); % determine sample-size of the meta-analysis iN = sum(vN(1:iC)); % compute lowest eigenvalue of the CGR matrix dMinEigCGR = min(eig(mCGR)); if dMinEigCGR > 0 % if CGR matrix is positive definite: continue % retrieve Cholesky decomposition of CGR matrix mGamma = chol(mCGR)'; % set Var(Z-stat) of meta and numerator corr(PGS,Y) in hold-out dVarZmeta = 0; dNumeratorR = 0; % for each genetic factor in the meta-analysis for iI = 1:iC % set the standard-error contributed to the Z-statistic dThisSE = 0; % for each study affected by this genetic factor for iJ = iI:iC % update the standard-error contributed to the Z-statistic dThisSE = dThisSE + ... (mGamma(iJ,iI)*vN(iJ)*sqrt(vH2(iJ)/(iM-vH2(iJ)))); % update numerator of expected R-squared according to relation % hold-out study and current meta study with this factor dNumeratorR = dNumeratorR + ... (mGamma(iC+1,iI)*mGamma(iJ,iI)*... vN(iJ)*sqrt(vH2(iJ)/(iM-vH2(iJ)))); end % add contribution to the Var(Z-stat) of this factor dThisVar = dThisSE^2; dVarZmeta = dVarZmeta + dThisVar; end % comute final R-squared numerator and denominator, and Var(Z-stat) dNumeratorR2 = dNumeratorR^2; dDenominatorR2 = dVarZmeta + ((iN*iS)/iM); dVarZmeta = 1 + (dVarZmeta/iN); % compute power resulting from Var(Z-stat) dP = 1 - erf(abs(dZhit)/sqrt(2*(dVarZmeta))); % compute PGS R-squared for the hold-out sample dR2 = vH2(iC+1)*dNumeratorR2/dDenominatorR2; else % if CGR matrix is not positive definite: terminate disp(['ERROR: CGR matrix is not positive definite.']); % return NaNs dP = NaN; dR2 = NaN; end end