A detailed look at some of the SuSiE fits¶
In comparison with CAVIAR follow-ups; workflow implemented in this notebook.
In [1]:
%revisions -s -n 10
Conclusions¶
- For these identified disagreement between CAVIAR in Example 1 and 2, if we set
estimate_prior_variance=TRUEinsusieR::susiewe remove the differences. - CAVIAR and DAP are in better agreement using default parameters.
Example 1¶
In [2]:
%preview /home/gaow/Documents/GIT/LargeFiles/JointLCL/AS_output/fastqtl_qqnorm_ASintron_RNAseqGeuvadis_YangVCF_100Kb/PVE_Prior/CAVIAR_follow_up/chr11_63987798_63988488_clu_2010.CAVIAR.png
We trace back to its input data-set and see what is going on.
In [3]:
ls /home/gaow/Documents/GIT/LargeFiles/JointLCL/AS_output/fastqtl_qqnorm_ASintron_RNAseqGeuvadis_YangVCF_100Kb/PVE_Prior/SuSiE_CS_*/chr11_63987798_63988488_clu_2010*
In [4]:
fn = "/home/gaow/Documents/GIT/LargeFiles/JointLCL/AS_output/fastqtl_qqnorm_ASintron_RNAseqGeuvadis_YangVCF_100Kb/PVE_Prior/SuSiE_CS_2/chr11_63987798_63988488_clu_2010.rds"
susie = readRDS(fn)
dat = readRDS(susie$input)[[susie$idx]]
saveRDS(dat, "~/chr11_63987798_63988488_clu_2010_data.rds")
In [5]:
top_idx = which.max(abs(dat$z_score))
b = rep(0, length(dat$z_score))
b[top_idx] = 1
In [6]:
dat$pve
Fit SuSiE with prior variance set to PVE estimate¶
In [7]:
library(susieR)
set.seed(1)
fitted = susie(dat$X, dat$y, L=5,
estimate_residual_variance=TRUE,
prior_variance=dat$pve/(1-dat$pve),
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
In [8]:
susie_pplot(pip, res = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
To plot the fitting process:
In [9]:
susie_iterplot(fitted, 3, 'example_1')
In [10]:
%preview example_1.gif
In [11]:
fitted$niter
In [12]:
fitted$elbo
However, if I set L=2 as I did for CAVIAR (-c 2), the result is more consistent with CAVIAR.
In [13]:
set.seed(1)
fitted = susie(dat$X, dat$y, L=2,
estimate_residual_variance=TRUE,
prior_variance=dat$pve/(1-dat$pve),
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
susie_pplot(pip, res = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
In [14]:
fitted$niter
In [15]:
fitted$elbo
Takes longer to converge but the ELBO is bigger for L=2.
In [16]:
susie_iterplot(fitted, 3, 'example_1_l2')
In [17]:
%preview example_1_l2.gif
Analysis with CAVIAR -c 3¶
In [18]:
%cd ~/tmp/14-Jul-2018
In [28]:
library(abind)
mm_regression = function(X, Y, Z=NULL) {
if (!is.null(Z)) {
Z = as.matrix(Z)
}
reg = lapply(seq_len(ncol(Y)), function (i) simplify2array(susieR:::univariate_regression(X, Y[,i], Z)))
reg = do.call(abind, c(reg, list(along=0)))
# return array: out[1,,] is betahat, out[2,,] is shat
return(aperm(reg, c(3,2,1)))
}
sumstats = mm_regression(as.matrix(dat$X), as.matrix(dat$y))
dat$Y = matrix(dat$y, length(dat$y), 1)
saveRDS(list(data=dat, sumstats=sumstats), 'N2.with_sumstats.rds')
In [7]:
export PATH=~/GIT/github/mvarbvs/dsc/modules/linux:$PATH
Rscript ~/GIT/software/susieR/inst/code/caviar.R input=\"N2.with_sumstats.rds\" output=\"N2finemapping.CAVIAR\" args=\"-c\ 3\"
In [11]:
caviar = readRDS("N2finemapping.CAVIAR.rds")
snp = caviar[[1]]$snp
pip = snp[order(as.numeric(snp$snp)),]$snp_prob
susie_pplot(pip, dtype='PIP', b=b, main = 'CAVIAR')
The result is less "clean", but still CAVIAR is able to capture the top signal.
Analysis with DAP-g¶
In [22]:
export PATH=~/GIT/github/mvarbvs/dsc/modules/linux:$PATH
python ~/GIT/software/susieR/inst/code/dap-g.py N2.with_sumstats.rds N2finemapping.DAP -ld_control 0.20 --all
In [23]:
dap = readRDS("N2finemapping.DAP.rds")
snp = dap[[1]]$snp
pip = snp[order(as.numeric(snp$snp)),]$snp_prob
susie_pplot(pip, dtype='PIP', b=b, main = 'DAP')
DAP seems to work fine in this case.
SuSiE with L 10¶
In [25]:
set.seed(1)
fitted = susie(dat$X, dat$y, L=10,
estimate_residual_variance=TRUE,
prior_variance=dat$pve/(1-dat$pve),
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
susie_pplot(pip, res = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
Interestingly, 1 CS is picked up, that does not include the signal.
SuSiE with L = 5 but estimated prior variance¶
In [19]:
set.seed(1)
fitted = susie(dat$X, dat$y, L=5,
estimate_residual_variance=TRUE,
estimate_prior_variance=TRUE,
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
susie_pplot(pip, res = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
We now get the same output as CAVIAR.
Example 2¶
In [20]:
%preview /home/gaow/Documents/GIT/LargeFiles/JointLCL/AS_output/fastqtl_qqnorm_ASintron_RNAseqGeuvadis_YangVCF_100Kb/PVE_Prior/CAVIAR_follow_up/chr10_114186117_114186976_clu_3277.CAVIAR.png
In [21]:
ls /home/gaow/Documents/GIT/LargeFiles/JointLCL/AS_output/fastqtl_qqnorm_ASintron_RNAseqGeuvadis_YangVCF_100Kb/PVE_Prior/SuSiE_CS_*/chr10_114186117_114186976_clu_3277*
In [22]:
fn = "/home/gaow/Documents/GIT/LargeFiles/JointLCL/AS_output/fastqtl_qqnorm_ASintron_RNAseqGeuvadis_YangVCF_100Kb/PVE_Prior/SuSiE_CS_2/chr10_114186117_114186976_clu_3277.rds"
susie = readRDS(fn)
dat = readRDS(susie$input)[[susie$idx]]
saveRDS(dat, "~/chr10_114186117_114186976_clu_3277_data.rds")
In [23]:
top_idx = which.max(abs(dat$z_score))
b = rep(0, length(dat$z_score))
b[top_idx] = 1
In [24]:
dat$pve
In [26]:
library(susieR)
set.seed(1)
fitted = susie(dat$X, dat$y, L=5,
estimate_residual_variance=TRUE,
prior_variance=dat$pve/(1-dat$pve),
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
susie_pplot(pip, res = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
To track iterations:
In [22]:
susie_iterplot(fitted, 3, 'example_2')
In [23]:
%preview example_2.gif
In [20]:
fitted$niter
In [21]:
fitted$elbo
If I set L=1:
In [30]:
library(susieR)
set.seed(1)
fitted = susie(dat$X, dat$y, L=1,
estimate_residual_variance=TRUE,
prior_variance=dat$pve/(1-dat$pve),
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
susie_pplot(pip, fitted = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
In [31]:
fitted$elbo
The ELBO here is in fact better.
In [33]:
susie_iterplot(fitted, 3, 'example_2_l1')
In [34]:
%preview example_2_l1.gif
CAVIAR and DAP other parameters¶
CAVIAR is more robust to parameter choices here.
In [29]:
sumstats = mm_regression(as.matrix(dat$X), as.matrix(dat$y))
dat$Y = matrix(dat$y, length(dat$y), 1)
saveRDS(list(data=dat, sumstats=sumstats), 'N2.with_sumstats.rds')
In [14]:
export PATH=~/GIT/github/mvarbvs/dsc/modules/linux:$PATH
Rscript ~/GIT/software/susieR/inst/code/caviar.R input=\"N2.with_sumstats.rds\" output=\"N2finemapping.CAVIAR\" args=\"-c\ 3\"
In [16]:
caviar = readRDS("N2finemapping.CAVIAR.rds")
snp = caviar[[1]]$snp
pip = snp[order(as.numeric(snp$snp)),]$snp_prob
susie_pplot(pip, dtype='PIP', b=b, main = 'CAVIAR')
In [30]:
export PATH=~/GIT/github/mvarbvs/dsc/modules/linux:$PATH
python ~/GIT/software/susieR/inst/code/dap-g.py N2.with_sumstats.rds N2finemapping.DAP -ld_control 0.20 --all
In [32]:
dap = readRDS("N2finemapping.DAP.rds")
snp = dap[[1]]$snp
pip = snp[order(as.numeric(snp$snp)),]$snp_prob
susie_pplot(pip, dtype='PIP', b=b, main = 'DAP')
Here DAP also does not think the top z-score is strong enough to be boosted.
SuSiE with prior estimated¶
In [33]:
set.seed(1)
fitted = susie(dat$X, dat$y, L=5,
estimate_residual_variance=TRUE,
estimate_prior_variance=TRUE,
tol=1e-3, track_fit=TRUE)
sets = susie_get_CS(fitted,
coverage = 0.95,
X = dat$X,
min_abs_corr = 0.4)
pip = susie_get_PIP(fitted, sets$cs_index)
susie_pplot(pip, res = fitted, CS = sets, max_cs = 0.4, dtype='PIP',b=b)
Now SuSiE gives exactly same result as CAVIAR.
In [36]:
fitted$sa2
In [ ]:
Copyright © 2016-2020 Gao Wang et al at Stephens Lab, University of Chicago