Last updated: 2020-06-24
Checks: 7 0
Knit directory: rss/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200623) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 1e806af. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rproj.user/
Ignored: .spelling
Ignored: examples/example5/.Rhistory
Ignored: examples/example5/Aseg_chr16.mat
Ignored: examples/example5/example5_simulated_data.mat
Ignored: examples/example5/example5_simulated_results.mat
Ignored: examples/example5/ibd2015_path2641_genes_results.mat
Untracked files:
Untracked: docs_old/
Unstaged changes:
Modified: rmd/_site.yml
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (rmd/example_5b.Rmd) and HTML (docs/example_5b.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 5db4b3e | Xiang Zhu | 2020-06-24 | Build site. |
| Rmd | 097406c | Xiang Zhu | 2020-06-24 | wflow_publish(“rmd/example_5b.Rmd”) |
This is Part B of Example 5, which illustrates how to perform enrichment analysis of GWAS summary statistics based on variational Bayes (VB) inference of RSS-BVSR model. This part describes an end-to-end enrichment analysis of inflammatory bowel disease (IBD) GWAS summary statistics (Liu et al, 2015) and a gene set named “IL23-mediated signaling events” (Pathway Commons 2, PID, 37 genes) using RSS. This part illustrates the actual data analyses performed in Zhu and Stephens (2018).
To reproduce results of Example 5 Part B, please use scripts in the directory example5, and follow the step-by-step guide below. Before running any script in example5, please install the VB subroutines of RSS. Please find installation instructions here.
Since a genome-wide enrichment analysis is conducted here, this part is more complicated than Example 5 Part A, which is based on a relatively small simulated dataset. It is advisable to go through Example 5 Part A before diving into this real data example.
Note that the working directory here is assumed to be rss/examples/example5. Please modify certain scripts accordingly if a different directory is used.
example5_null.mSince we have to perform Bayesian multiple regression analyses of 1.1 million common SNPs multiple times when fitting the baseline model, we need to use the memory-efficient and parallel implementation of RSS rss_varbvsr_bigmem_squarem.m, which requires MATLAB Parallel Computing Toolbox. If you do not have this toolbox available, please skip this section, and use our result files for the following enrichment analysis instead (see Step 4).
Step 1. Download the input data file ibd2015_sumstat.mat, which contains the GWAS summary statistics and LD matrix estimates. This file is large (17G) because it has a LD matrix of 1.1 million common SNPs. Please contact me (xiangzhu[at]uchicago.edu) if you have trouble accessing this file.
Before proceeding to the next step, let’s look at the contents of ibd2015_sumstat.mat.
>> sumstat = matfile('ibd2015_sumstat.mat');
>> sumstat
sumstat =
matlab.io.MatFile
Properties:
Properties.Source: 'ibd2015_sumstat.mat'
Properties.Writable: false
SiRiS: [22x1 cell]
betahat: [22x1 cell]
chr: [22x1 cell]
pos: [22x1 cell]
se: [22x1 cell]Here GWAS summary statistics and LD estimates are stored as cell arrays. For each Chromosome j,
betahat{j,1} stores single-SNP effect size estimates of all SNPs on Chromosome j;se{j,1} stores standard errors of betahat{j, 1};chr{j,1} and pos{j, 1} store physical positions of these SNPs (GRCh37 build);SiRiS{j,1} stores a sparse matrix, defined as repmat((1./se),1,p) .* R .* repmat((1./se)',p,1), where R is the estimated LD matrix of these p SNPs.Step 2. Specify several dataset-specific variables that are required by null_template.m, a template script that fits a baseline model.
% specify trait-specific information
trait_name = 'ibd2015';
sample_size = (12882+21770); % cases: 12,882; controls: 21,770
% specify grid of hyper-parameters
h_rv = 0.3;
theta0_rv = (-2.9:0.025:-2.85);
% specify stage of analysis
stage = 'step1';
% specify random start and algorithm
myseed = 459;
method = 'squarem';Note that here the grid of hyper-parameters is set minimum for illustration purpose. In practice we are using much larger grid of hyper-parameters (because we seldomly have a sensible guess for these parameters and we should learn them from data); please see Supplementary Tables 6-7 of Zhu and Stephens (2018) for more details.
Step 3. Fit baseline models across a grid of hyper-parameters in parallel. We write a simple sbatch script, example5_null.sbatch, and submit it to a cluster where Slurm has been installed.
The script example5_null.sbatch makes three copies of example5_null.m (because there are three sets of (h, theta0) in the grid specified in Step 2), and then assigns one copy to each (h, theta0) setting.
After the sbatch submission, these three jobs run in three different nodes simultaneously.
xiangzhu@midway-login1: squeue -u xiangzhu
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
29451148_3 sandyb ibd2015 xiangzhu R 1:29 1 midway249
29451148_1 sandyb ibd2015 xiangzhu R 2:42 1 midway426
29451148_2 sandyb ibd2015 xiangzhu R 2:42 1 midway427Note that example5_null.sbatch is designed for the computing clusters at University of Chicago Research Computing Center. You may need to modify this script if you plan to run it on a different environment.
Step 4. Aggregate the baseline results for enrichment analyses. As soon as the three jobs in Step 3 are completed, the baseline results are saved as mat files (version 7.3).
In case you cannot fit the baseline models in your own computing environment, you can download the following baseline result files from https://projects.rcc.uchicago.edu/mstephens/rss_wiki/example5/.
xiangzhu@midway-login1: ls ibd2015_null_*.mat
ibd2015_null_h_30_theta0_285_seed_459_squarem_step1.mat
ibd2015_null_h_30_theta0_288_seed_459_squarem_step1.mat
ibd2015_null_h_30_theta0_290_seed_459_squarem_step1.matexample5_gsea.mUnlike fitting the baseline model, fitting the enrichment model does not require any specialized toolbox or any high-performance computing cluster. With the baseline result files listed above in place (i.e. downloaded and saved in the working directory rss/examples/example5), everyone should be able to run this part in a modern desktop by simply typing the following line in a MATLAB console:
Since we use additional pathway information when fitting the enrichment model, here we need to download and add a new data file ibd2015_sumstat_path2641.mat to the working directory rss/examples/example5. Let’s look at the content of this file first.
>> sumstat = matfile('ibd2015_sumstat_path2641.mat');
>> sumstat
sumstat =
matlab.io.MatFile
Properties:
Properties.Source: 'ibd2015_sumstat_path2641.mat'
Properties.Writable: false
SiRiS: [2990x2990 double]
betahat: [2990x1 double]
chr: [2990x1 double]
pos: [2990x1 double]
se: [2990x1 double]
snps: [2990x1 double]This new file contains GWAS summary statistics and LD estimates of SNPs that are “inside” a biological pathway named “IL23-mediated signaling events” (Pathway Commons 2, PID, 37 genes). As in Zhu and Stephens (2018), here we annotate a SNP as “inside a pathway” if this SNP is within 100 kb of transcribed region of any member gene in this pathway.
The results of fitting the enrichment model are saved as ibd2015_gsea_seed_459_path2641_squarem.mat.
>> load ibd2015_gsea_seed_459_path2641_squarem.mat
>> whos
Name Size Bytes Class Attributes
alpha 2990x3x101 7247760 double
h 1x1 8 double
log10bf 1x1 8 double
logw0 3x1 24 double
logw1 3x101 2424 double
method 1x7 14 char
mu 2990x3x101 7247760 double
myseed 1x1 8 double
runtime 1x1 8 double
s 2990x3x101 7247760 double
snps 2990x1 23920 double
theta 101x1 808 double
theta0 3x1 24 doubleThere are three groups of output variables in the result file.
The first group consists of user-specified quantities in example5_gsea.m: method (model fitting algorithm), myseed (seed of random number generator), h & theta0 & theta (grid of hyper-parameters), and snps (indices of “inside-pathway” SNPs in the genome).
The second group consists of estimated variational parameters alpha & mu & s under the enrichment model, where alpha(:,i,j) & mu(:,i,j) & s(:,i,j) correspond to estimation under the hyper-parameter setting [h, theta0(i), theta(j)]. Note that these variational parameter estimates can be further used to prioritize genes within an enriched pathway; please see the next section.
The third group consists of estimated variational lower bounds under the baseline (logw0) and enrichment (logw1) model, and the log 10 enrichment Bayes factor (log10bf).
The log 10 enrichment Bayes factor is 22.1355, suggesting strong enrichment of genetic associations within “IL23-mediated signaling events” for IBD. This seems consistent with the important role of IL23 in immune response.
example5_gene.mNow we use the baseline and enrichment model fitting results above to prioritize 37 member genes within the target pathway. Similar to Example 5 Part A, we need to download and save a gene information file ibd2015_path2641_genes.mat in the working directory rss/examples/example5 first.
>> load ibd2015_path2641_genes.mat
>> whos
Name Size Bytes Class Attributes
Aseg 1081481x37 160000304 double sparse
chr 1081481x1 4325924 int32
gene_chr 37x1 296 double
gene_id 37x1 4494 cell
gene_start 37x1 296 double
gene_stop 37x1 296 double
pos 1081481x1 4325924 int32
>> whos Aseg
Name Size Bytes Class Attributes
Aseg 1081481x37 160000304 double sparse
>> full(unique(Aseg(:)))'
ans =
0 1As shown above, there are 1,081,481 genome-wide SNPs and 37 pathway member genes (both based on GRCh Build 37) in this simulated dataset. Here we assign a SNP j to a gene g (i.e. Aseg(j,g)==1) if and only if SNP j is within 100 kb of transcribed region of gene g (i.e. [gene_start-100e3 gene_stop+100e3]).
With ibd2015_path2641_genes.mat in place, we can then prioritize pathway genes by simply typing the following line in a MATLAB console:
The results of prioritizing pathway genes are saved as ibd2015_path2641_genes_results.mat.
>> load ibd2015_path2641_genes_results.mat
>> whos
Name Size Bytes Class Attributes
b_p1 37x1 296 double
b_p2 37x1 296 double
e_p1 37x1 296 double
e_p2 37x1 296 double
gene_chr 37x1 296 double
gene_id 37x1 4494 cell
gene_start 37x1 296 double
gene_stop 37x1 296 doubleSimilar to Example 5 Part A, [b/e]_p_[1/2] denotes the posterior probability of each locus (gene with 100 kb window) containing at lease 1/2 trait-associated SNPs under the baseline/enrichment model.
We can easily load the gene-level results in R as follows:
> suppressPackageStartupMessages(library(R.matlab))
> suppressPackageStartupMessages(library(dplyr))
>
> res <- R.matlab::readMat("ibd2015_path2641_genes_results.mat")
>
> df <- data.frame(id=unlist(res$gene.id), chr=res$gene.chr, start=res$gene.start, stop=res$gene.stop, b_p1=res$b.p1, e_p1=res$e.p1, b_p2=res$b.p2, e_p2=res$e.p2)
> names(df) <- c("gene","chr","start","stop","b_p1","e_p1","b_p2","e_p2")
> df <- dplyr::arrange(df, -e_p1, -b_p1)
>
> sum(df$e_p1 >= df$b_p1)
[1] 35
> sum(df$e_p2 >= df$b_p2)
[1] 36We can see that the enrichment model produces stronger gene-level signals for almost all 37 pathway genes than the baseline model:
| gene | chr | start | stop | b_p1 | e_p1 | b_p2 | e_p2 |
|---|---|---|---|---|---|---|---|
| IL23R | 1 | 67632169 | 67725662 | 1.0000 | 1.0000 | 1.0000 | 1.0000 |
| IL12B | 5 | 158741791 | 158757481 | 1.0000 | 1.0000 | 0.9994 | 1.0000 |
| IL19 | 1 | 206972215 | 207016326 | 1.0000 | 1.0000 | 0.4685 | 0.9933 |
| JAK2 | 9 | 4985245 | 5128183 | 1.0000 | 1.0000 | 0.7649 | 0.9995 |
| IFNG | 12 | 68548550 | 68553521 | 1.0000 | 1.0000 | 0.2420 | 0.9234 |
| STAT3 | 17 | 40465343 | 40540513 | 0.9999 | 1.0000 | 0.6887 | 0.9944 |
| CCL2 | 17 | 32582296 | 32584222 | 0.9994 | 1.0000 | 0.0941 | 0.8375 |
| IL18RAP | 2 | 103035254 | 103069025 | 0.9994 | 1.0000 | 0.3131 | 0.9374 |
| IL18R1 | 2 | 102979097 | 103015218 | 0.9995 | 1.0000 | 0.3393 | 0.9346 |
| TYK2 | 19 | 10461204 | 10491248 | 0.9629 | 0.9999 | 0.1699 | 0.8993 |
| STAT5A | 17 | 40439565 | 40463961 | 0.9998 | 0.9997 | 0.0618 | 0.6792 |
| IL2 | 4 | 123372625 | 123377650 | 0.2659 | 0.9989 | 0.0370 | 0.9074 |
| STAT4 | 2 | 191894306 | 192015925 | 0.7813 | 0.9976 | 0.1240 | 0.8868 |
| STAT1 | 2 | 191833762 | 191878976 | 0.7642 | 0.9953 | 0.0735 | 0.7793 |
| IL6 | 7 | 22766766 | 22771621 | 0.2801 | 0.9776 | 0.0364 | 0.7868 |
| NFKBIA | 14 | 35870716 | 35873960 | 0.2497 | 0.9298 | 0.0318 | 0.7264 |
| CXCL9 | 4 | 76922623 | 76928641 | 0.4166 | 0.9189 | 0.0991 | 0.7151 |
| PIK3R1 | 5 | 67511584 | 67597649 | 0.1645 | 0.9125 | 0.0136 | 0.6976 |
| PIK3CA | 3 | 178866311 | 178952500 | 0.1889 | 0.9094 | 0.0185 | 0.6923 |
| IL18 | 11 | 112013976 | 112034840 | 0.1543 | 0.8683 | 0.0105 | 0.5506 |
| IL24 | 1 | 207070788 | 207077484 | 0.1595 | 0.8586 | 0.0121 | 0.5677 |
| IL17F | 6 | 52101484 | 52109298 | 0.1160 | 0.8515 | 0.0069 | 0.5720 |
| IL17A | 6 | 52051185 | 52055436 | 0.1155 | 0.8504 | 0.0068 | 0.5697 |
| MPO | 17 | 56347217 | 56358296 | 0.1085 | 0.8483 | 0.0059 | 0.5614 |
| NFKB1 | 4 | 103422486 | 103538459 | 0.1073 | 0.8290 | 0.0059 | 0.5313 |
| RELA | 11 | 65421067 | 65430443 | 0.1103 | 0.8263 | 0.0061 | 0.5195 |
| IL1B | 2 | 113587337 | 113594356 | 0.1193 | 0.8263 | 0.0071 | 0.5224 |
| IL12RB1 | 19 | 18170371 | 18197697 | 0.1286 | 0.8109 | 0.0078 | 0.4811 |
| NOS2 | 17 | 26083792 | 26127555 | 0.1178 | 0.7861 | 0.0069 | 0.4544 |
| ALOX12B | 17 | 7975954 | 7991021 | 0.0773 | 0.7716 | 0.0030 | 0.4345 |
| CD3E | 11 | 118175295 | 118186890 | 0.0946 | 0.7692 | 0.0044 | 0.4299 |
| SOCS3 | 17 | 76352859 | 76356158 | 0.1541 | 0.7264 | 0.0115 | 0.3666 |
| CXCL1 | 4 | 74735109 | 74736959 | 0.0813 | 0.7103 | 0.0033 | 0.3517 |
| ITGA3 | 17 | 48133340 | 48167849 | 0.0577 | 0.6682 | 0.0017 | 0.3035 |
| IL23A | 12 | 56732663 | 56734194 | 0.0370 | 0.5185 | 0.0007 | 0.1672 |
| CD4 | 12 | 6898638 | 6929976 | 0.0337 | 0.4747 | 0.0006 | 0.1366 |
| TNF | 6 | 31543350 | 31546112 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
The enrichment analyses of 31 complex traits and 4,026 gene sets reported in Zhu and Stephens (2018) are essentially 124,806 “replications” of the example above (with larger grids on hyper-parameters). Our full analysis results are publicly available at https://xiangzhu.github.io/rss-gsea/.
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.0.1 (2020-06-06)
os macOS Catalina 10.15.5
system x86_64, darwin17.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/Los_Angeles
date 2020-06-24
─ Packages ───────────────────────────────────────────────────────────────────
package * version date lib source
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.0)
backports 1.1.8 2020-06-17 [1] CRAN (R 4.0.0)
callr 3.4.3 2020-03-28 [1] CRAN (R 4.0.0)
cli 2.0.2 2020-02-28 [1] CRAN (R 4.0.0)
crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.0)
desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.0)
devtools 2.3.0 2020-04-10 [1] CRAN (R 4.0.0)
digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.0)
ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.0)
evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.0)
fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.0)
fs 1.4.1 2020-04-04 [1] CRAN (R 4.0.0)
git2r 0.27.1 2020-05-03 [1] CRAN (R 4.0.0)
glue 1.4.1 2020-05-13 [1] CRAN (R 4.0.0)
htmltools 0.5.0 2020-06-16 [1] CRAN (R 4.0.0)
httpuv 1.5.4 2020-06-06 [1] CRAN (R 4.0.0)
knitr 1.29 2020-06-23 [1] CRAN (R 4.0.0)
later 1.1.0.1 2020-06-05 [1] CRAN (R 4.0.0)
magrittr 1.5 2014-11-22 [1] CRAN (R 4.0.0)
memoise 1.1.0 2017-04-21 [1] CRAN (R 4.0.0)
pkgbuild 1.0.8 2020-05-07 [1] CRAN (R 4.0.0)
pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.0)
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.0)
processx 3.4.2 2020-02-09 [1] CRAN (R 4.0.0)
promises 1.1.1 2020-06-09 [1] CRAN (R 4.0.0)
ps 1.3.3 2020-05-08 [1] CRAN (R 4.0.0)
R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.0)
Rcpp 1.0.4.6 2020-04-09 [1] CRAN (R 4.0.0)
remotes 2.1.1 2020-02-15 [1] CRAN (R 4.0.0)
rlang 0.4.6 2020-05-02 [1] CRAN (R 4.0.0)
rmarkdown 2.3 2020-06-18 [1] CRAN (R 4.0.0)
rprojroot 1.3-2 2018-01-03 [1] CRAN (R 4.0.0)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.0)
stringi 1.4.6 2020-02-17 [1] CRAN (R 4.0.0)
stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.0)
testthat 2.3.2 2020-03-02 [1] CRAN (R 4.0.0)
usethis 1.6.1 2020-04-29 [1] CRAN (R 4.0.0)
whisker 0.4 2019-08-28 [1] CRAN (R 4.0.0)
withr 2.2.0 2020-04-20 [1] CRAN (R 4.0.0)
workflowr * 1.6.2 2020-04-30 [1] CRAN (R 4.0.0)
xfun 0.15 2020-06-21 [1] CRAN (R 4.0.0)
yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.0)
[1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library