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 097406c. 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/

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
Rmd 097406c Xiang Zhu 2020-06-24 wflow_publish(“rmd/example_5b.Rmd”)

Overview

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.

Step-by-step illustration

Fitting the baseline model with example5_null.m

Since 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.

sbatch example5_null.sbatch

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 midway427

Note 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.mat

Fitting the enrichment model with example5_gsea.m

Unlike 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:

>> run example5_gsea.m

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  double

There 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.

Prioritizing pathway genes with example5_gene.m

Now 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     1

As 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:

>> run example5_gene.m

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  double

Similar 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] 36

We 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

More examples

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/.


devtools::session_info()
─ 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