Merge "examples-peter" with "master"

.gitignore

@@ -16,7 +16,6 @@ src/lightspeed*
 src/lapack*
 src/progress.m
 src/compute_pve2.m
-src/rss_bvsr_mdp.m
 src/rss_bvsr_cap.m
 src/propose_gamma_cap.m
 src/rss_bslmm_mwt.m

examples/example1.m

@@ -1,10 +1,15 @@
 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
-example_data = matfile('example1.mat');
+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 
@@ -16,31 +21,39 @@
 
 fprintf('Data set is loaded ... \n');
 
-% check model assumption
+% Check model assumption.
 chat = sqrt((betahat(:).^2) ./ (Nsnp(:).*(se(:).^2) + betahat(:).^2));
-fprintf('Look at the five-number summary of log10 sample phenotype-genotype correlations: \n')
-disp(prctile(log10(chat), 0:25:100));
+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
+% 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
+
+% 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);
+  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('example1_rssbvsr.mat', 'betasam', 'gammasam', 'hsam', 'logpisam', 'pvesam', 'Naccept', 'runtime', '-v7.3');
+
+% 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');
 
@@ -50,20 +63,26 @@
 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);
+[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);
+  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('example1_rssbslmm.mat', 'bsam', 'zsam', 'lpsam', 'hsam', 'rsam', 'pvesam', 'Naccept', 'runtime', '-v7.3');
+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');
 
@@ -74,20 +93,26 @@
 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);
+[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);
+  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('example1_rssash.mat', 'bsam', 'zsam', 'wsam', 'lsam', 'pvesam', 'Naccept', 'runtime', '-v7.3');
+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');
 

examples/example2.m

@@ -1,10 +1,15 @@
 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
-example_data = matfile('example2.mat');
+fprintf('Loading data.\n');
+example_data = matfile(cat(2,working_dir,'/','example2.mat'));
 
 betahat = example_data.betahat; % single-SNP effect size estimates
 se 	= example_data.se; 	% standard error of betahat
@@ -16,32 +21,41 @@
 Nthin = 9e1;
 
 % load the three types of ld matrix
-genotype_data = matfile('genotype2.mat');
+genotype_data = matfile(cat(2,working_dir,'/','genotype2.mat'));
 
 % fit rss-bvsr model when R is the sample correlation in the panel
 cohort_R = genotype_data.cohort_R;
 tic;
-[betasam, gammasam, hsam, logpisam, Naccept] = rss_bvsr(betahat, se, cohort_R, Nsnp, Ndraw, Nburn, Nthin);
+[betasam, gammasam, hsam, logpisam, Naccept] = ...
+  rss_bvsr(betahat, se, cohort_R, Nsnp, Ndraw, Nburn, Nthin);
 runtime = toc;
-save('example2_rss_c.mat', 'betasam', 'gammasam', 'hsam', 'logpisam', 'Naccept', 'runtime', '-v7.3');
+save(cat(2,working_dir,'/','example2_rss_c.mat'),...
+     'betasam', 'gammasam', 'hsam', 'logpisam', 'Naccept',...
+     'runtime', '-v7.3');
 clearvars betasam gammasam hsam logpisam Naccept runtime;
 fprintf('RSS-C is done ... \n');
 
 % fit rss-bvsr model when R is the sample correlation in the panel
 shrink_R = genotype_data.shrink_R;
 tic;
-[betasam, gammasam, hsam, logpisam, Naccept] = rss_bvsr(betahat, se, shrink_R, Nsnp, Ndraw, Nburn, Nthin);
+[betasam, gammasam, hsam, logpisam, Naccept] = ...
+  rss_bvsr(betahat, se, shrink_R, Nsnp, Ndraw, Nburn, Nthin);
 runtime = toc;
-save('example2_rss_s.mat', 'betasam', 'gammasam', 'hsam', 'logpisam', 'Naccept', 'runtime', '-v7.3');
+save(cat(2,working_dir,'/','example2_rss_s.mat'),...
+     'betasam', 'gammasam', 'hsam', 'logpisam', 'Naccept',...
+     'runtime', '-v7.3');
 clearvars betasam gammasam hsam logpisam Naccept runtime;
 fprintf('RSS is done ... \n');
 
 % fit rss-bvsr model when R is the sample correlation in the panel
 panel_R = genotype_data.panel_R;
 tic;
-[betasam, gammasam, hsam, logpisam, Naccept] = rss_bvsr(betahat, se, panel_R, Nsnp, Ndraw, Nburn, Nthin);
+[betasam, gammasam, hsam, logpisam, Naccept] = ...
+  rss_bvsr(betahat, se, panel_R, Nsnp, Ndraw, Nburn, Nthin);
 runtime = toc;
-save('example2_rss_p.mat', 'betasam', 'gammasam', 'hsam', 'logpisam', 'Naccept', 'runtime', '-v7.3');
+save(cat(2,working_dir,'/','example2_rss_p.mat'),...
+     'betasam', 'gammasam', 'hsam', 'logpisam', ...
+     'Naccept', 'runtime', '-v7.3');
 clearvars betasam gammasam hsam logpisam Naccept runtime;
 fprintf('RSS-P is done ... \n');
 

examples/example3.m

@@ -1,10 +1,15 @@
 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
-example_data = matfile('example1.mat');
+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 
@@ -22,40 +27,49 @@
 
 abs_diff = abs(se_1 - se_2);
 fprintf('Look at the five-number summary of the absolute difference between two definitions of SE: \n')
-disp(prctile(log10(abs_diff), 0:25:100)); % require stat toolbox
+disp(percentile(log10(abs_diff),0:0.25:1));
+
+% fit rss-bvsr model with the first definition of se
+fprintf('Start RSS-BVSR analysis on the first definition of SE... \n');
 
 % mcmc info
 Ndraw = 2e6;
 Nburn = 2e5;
 Nthin = 9e1;
 
-% fit rss-bvsr model with the first definition of se
-fprintf('Start RSS-BVSR analysis on the first definition of SE... \n');
-
 tic;
-% 1. simulate posterior samples via mcmc
-[betasam, gammasam, hsam, logpisam, Naccept] = rss_bvsr(betahat, se_1, R, Nsnp, Ndraw, Nburn, Nthin);
-% 2. compute the posterior samples of pve
+
+% 1. Simulate posterior samples via mcmc.
+[betasam, gammasam, hsam, logpisam, Naccept] = ...
+    rss_bvsr(betahat, se_1, 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_1, Nsnp, bwd, BR, matrix_type);
-	progress_bar 	= progress(progress_bar, i/M);
+  pvesam(i) 	= compute_pve(betasam(i,:), betahat, se_1, Nsnp, bwd, BR, matrix_type);
+  progress_bar 	= progress(progress_bar, i/M);
 end
 runtime = toc;
-% 3. save output
-save('example3_se1.mat', 'betasam', 'gammasam', 'hsam', 'logpisam', 'pvesam', 'Naccept', 'runtime', '-v7.3');
+
+% 3. Save output.
+save(cat(2,working_dir,'/','example3_se1.mat'),...
+     'betasam', 'gammasam', 'hsam', 'logpisam', 'pvesam',...
+     'Naccept', 'runtime', '-v7.3');
 clearvars betasam gammasam hsam logpisam pvesam Naccept runtime;
 
 % fit rss-bvsr model with the second definition of se
 fprintf('Start RSS-BVSR analysis on the second definition of SE... \n');
 
 tic;
-% 1. simulate posterior samples via mcmc
-[betasam, gammasam, hsam, logpisam, Naccept] = rss_bvsr(betahat, se_2, R, Nsnp, Ndraw, Nburn, Nthin);
-% 2. compute the posterior samples of pve
+
+% 1. Simulate posterior samples via mcmc.
+[betasam, gammasam, hsam, logpisam, Naccept] = ...
+    rss_bvsr(betahat, se_2, R, Nsnp, Ndraw, Nburn, Nthin);
+
+% 2. Compute the posterior samples of PVE.
 matrix_type     = 1;
 M               = length(hsam);
 pvesam          = zeros(M,1);
@@ -65,7 +79,9 @@
         progress_bar    = progress(progress_bar, i/M);
 end
 runtime = toc;
-% 3. save output
-save('example3_se2.mat', 'betasam', 'gammasam', 'hsam', 'logpisam', 'pvesam', 'Naccept', 'runtime', '-v7.3');
-clearvars betasam gammasam hsam logpisam pvesam Naccept runtime;
 
+% 3. Save output.
+save(cat(2,working_dir,'/','example3_se2.mat'),...
+     'betasam', 'gammasam', 'hsam', 'logpisam', 'pvesam',...
+     'Naccept', 'runtime', '-v7.3');
+clearvars betasam gammasam hsam logpisam pvesam Naccept runtime;

examples/example3_plots.m

@@ -0,0 +1,36 @@
+% Run this script after example3.m. Note: This script only works in MATLAB
+% 2014a or later since "histogram" function was introduced in that version
+% of MATLAB.
+se1 = matfile('example3_se1.mat');
+se2 = matfile('example3_se2.mat');
+
+figure(1);
+subplot(2,2,1);
+scatter(mean(se1.betasam)', mean(se2.betasam)','b','*');
+% Avoid using Statistics Toolbox:
+hold on
+plot([-1 2],[-1 2],'-','Color','red');
+hold off
+% ozline = refline([1 0]);
+% ozline.Color = 'r';
+xlabel('using se_1','Interpreter','none');
+ylabel('using se_2','Interpreter','none');
+title('Posterior mean of \beta');
+
+subplot(2,2,2);
+h1=histogram(se1.hsam,'Normalization','pdf'); 
+hold on; 
+h2= histogram(se2.hsam,'Normalization','pdf');
+title('Posterior density of h');
+
+subplot(2,2,3);
+h1=histogram(se1.logpisam,'Normalization','pdf'); 
+hold on; 
+h2= histogram(se2.logpisam,'Normalization','pdf');
+title('Posterior density of log(\pi)');
+
+subplot(2,2,4);
+h1=histogram(se1.pvesam,'Normalization','pdf'); 
+hold on; 
+h2=histogram(se2.pvesam,'Normalization','pdf');
+title('Posterior density of PVE');