rss/examples/example1.m at b97507335a7af6916ab7ea61917883b1bca7d8e9 · stephenslab/rss · GitHub

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
clear;

% Set this to the directory containing example1.mat and where the output
% files will be stored.
working_dir = '.';

% add search paths
addpath(genpath('../src'));

% load summary-level data
fprintf('Loading data.\n');
example_data = matfile(cat(2,working_dir,'/','example1.mat'));

R   = example_data.R; 	% population LD matrix
bwd = example_data.bwd; % bandwidth of R
BR  = example_data.BR; 	% banded storage of R

betahat = example_data.betahat; % single-SNP effect size estimates
se 	= example_data.se; 	% standard error of betahat
Nsnp 	= example_data.Nsnp; 	% sample size of each SNP

fprintf('Data set is loaded ... \n');

% Check model assumption.
chat = sqrt((betahat(:).^2) ./ (Nsnp(:).*(se(:).^2) + betahat(:).^2));
fprintf('Look at the five-number summary of log10 sample\n');
fprintf('phenotype-genotype correlations: \n')
disp(percentile(log10(chat),0:0.25:1));

% Fit rss-bvsr model.
fprintf('Start RSS-BVSR analysis ... \n');
Ndraw = 2e6;
Nburn = 2e5;
Nthin = 9e1;
tic;

% 1. Simulate posterior samples via mcmc.
[betasam, gammasam, hsam, logpisam, Naccept] = ...
    rss_bvsr(betahat, se, R, Nsnp, Ndraw, Nburn, Nthin);

% 2. Compute the posterior samples of PVE.
matrix_type 	= 1;
M 		= length(hsam);
pvesam 		= zeros(M,1);
progress_bar 	= progress('init','start PVE calculation');
for i = 1:M
  pvesam(i) = compute_pve(betasam(i,:), betahat, se, Nsnp, bwd, BR,...
                          matrix_type);
  progress_bar	= progress(progress_bar, i/M);
end
runtime = toc;

% 3. Save output.
save(cat(2,working_dir,'/','example1_rssbvsr.mat'),...
     'betasam', 'gammasam', 'hsam', 'logpisam', 'pvesam', ...
     'Naccept', 'runtime', '-v7.3');
clearvars betasam gammasam hsam logpisam pvesam Naccept runtime;
fprintf('RSS-BVSR analysis is done ... \n');

% fit rss-bslmm model
fprintf('Start RSS-BSLMM analysis ... \n');
Ndraw = 2e6;
Nburn = 2e5;
Nthin = 9e1;
tic;

% 1. simulate posterior samples via mcmc
[bsam, zsam, lpsam, hsam, rsam, Naccept] = ...
    rss_bslmm(betahat, se, R, Nsnp, Ndraw, Nburn, Nthin);

% 2. compute the posterior samples of pve
matrix_type 	= 1;
M 		= length(hsam);
pvesam 		= zeros(M,1);
progress_bar 	= progress('init','start PVE calculation');
for i = 1:M
  pvesam(i) = compute_pve(bsam(i,:), betahat, se, Nsnp, bwd, BR, matrix_type);
  progress_bar = progress(progress_bar, i/M);
end
runtime = toc;

% 3. save output
save(cat(2,working_dir,'/','example1_rssbslmm.mat'),...
     'bsam', 'zsam', 'lpsam', 'hsam', 'rsam', 'pvesam',...
     'Naccept', 'runtime', '-v7.3');
clearvars bsam zsam lpsam hsam rsam pvesam Naccept runtime;
fprintf('RSS-BSLMM analysis is done ... \n');

% fit rss-ash model
fprintf('Start RSS-ASH analysis ... \n');
Ndraw = 5e7;
Nburn = 1e7;
Nthin = 1e3;
sigma_beta = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3];
tic;

% 1. simulate posterior samples via mcmc
[bsam, zsam, wsam, lsam, Naccept] =...
    rss_ash(betahat, se, R, Nsnp, sigma_beta, Ndraw, Nburn, Nthin);

% 2. compute the posterior samples of pve
matrix_type 	= 1;
M 		= length(lsam);
pvesam 		= zeros(M,1);
progress_bar 	= progress('init','start PVE calculation');
for i = 1:M
  pvesam(i) = compute_pve(bsam(i,:), betahat, se, Nsnp, bwd, BR, matrix_type);
  progress_bar	= progress(progress_bar, i/M);
end
runtime = toc;

% 3. save output
save(cat(2,working_dir,'/','example1_rssash.mat'),...
     'bsam', 'zsam', 'wsam', 'lsam', 'pvesam', 'Naccept', ...
     'runtime', '-v7.3');
clearvars bsam zsam wsam lsam pvesam Naccept runtime;
fprintf('RSS-ASH analysis is done ... \n');

exit;