Last updated: 2022-05-16

Checks: 7 0

Knit directory: logistic-susie-gsea/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). 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(20220105) 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 814cab6. 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:    .DS_Store
    Ignored:    .RData
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/figure/
    Ignored:    cache/
    Ignored:    library/
    Ignored:    renv/library/
    Ignored:    renv/staging/
    Ignored:    staging/

Untracked files:
    Untracked:  .ipynb_checkpoints/
    Untracked:  _targets.R
    Untracked:  _targets.html
    Untracked:  _targets.md
    Untracked:  _targets/
    Untracked:  additive.l5.gonr.aggregate.scores
    Untracked:  analysis/alpha_ash_v_point_normal.Rmd
    Untracked:  analysis/de_droplet_noshrink.Rmd
    Untracked:  analysis/de_droplet_noshrink_logistic_susie.Rmd
    Untracked:  analysis/fetal_reference_cellid_gsea.Rmd
    Untracked:  analysis/fixed_intercept.Rmd
    Untracked:  analysis/gsea_made_simple.Rmd
    Untracked:  analysis/iDEA_examples.Rmd
    Untracked:  analysis/latent_gene_list.Rmd
    Untracked:  analysis/linear_method_failure_modes.Rmd
    Untracked:  analysis/linear_regression_failure_regime.Rmd
    Untracked:  analysis/linear_v_logistic_pbmc.Rmd
    Untracked:  analysis/logistic_susie_veb_boost_vs_vb.Rmd
    Untracked:  analysis/logistic_susie_vis.Rmd
    Untracked:  analysis/logistic_variational_bound.Rmd
    Untracked:  analysis/logsitic_susie_template.Rmd
    Untracked:  analysis/pcb_scratch.Rmd
    Untracked:  analysis/redundancy_is_a_problem.Rmd
    Untracked:  analysis/references.bib
    Untracked:  analysis/roadmap.Rmd
    Untracked:  analysis/simulations.Rmd
    Untracked:  analysis/simulations_l1.Rmd
    Untracked:  analysis/tccm_vs_logistic_susie.Rmd
    Untracked:  analysis/template.Rmd
    Untracked:  analysis/test.Rmd
    Untracked:  build_site.R
    Untracked:  code/html_tables.R
    Untracked:  code/latent_logistic_susie.R
    Untracked:  code/logistic_susie_data_driver.R
    Untracked:  code/marginal_sumstat_gsea_collapsed.R
    Untracked:  code/point_normal.R
    Untracked:  code/sumstat_gsea.py
    Untracked:  code/susie_gsea_queries.R
    Untracked:  data/adipose_2yr_topsnp.txt
    Untracked:  data/de-droplet/
    Untracked:  data/deng/
    Untracked:  data/fetal_reference_cellid_gene_sets.RData
    Untracked:  data/human_chimp_eb/
    Untracked:  data/pbmc-purified/
    Untracked:  data/wenhe_baboon_diet/
    Untracked:  docs.zip
    Untracked:  export/
    Untracked:  index.md
    Untracked:  l1.sim.aggregate.scores
    Untracked:  linear_v_logistic.png
    Untracked:  presentations/
    Untracked:  references.bib
    Untracked:  simulation_targets/

Unstaged changes:
    Modified:   .gitignore
    Deleted:    _simulation_targets.R
    Modified:   _targets.Rmd
    Modified:   _targets.yaml
    Modified:   analysis/alpha_for_single_cell.Rmd
    Modified:   analysis/baboon_diet.Rmd
    Modified:   analysis/example_pbmc.Rmd
    Modified:   analysis/gseabenchmark_tcga.Rmd
    Modified:   analysis/human_chimp_eb_de_example.Rmd
    Modified:   analysis/linear_v_logistic_start.Rmd
    Modified:   analysis/single_cell_pbmc.Rmd
    Modified:   analysis/single_cell_pbmc_l1.Rmd
    Deleted:    analysis/summary_stat_gsea_univariate_simulations.Rmd
    Modified:   analysis/the_big_geneset.Rmd
    Modified:   code/fit_baselines.R
    Modified:   code/fit_logistic_susie.R
    Modified:   code/fit_mr_ash.R
    Modified:   code/fit_susie.R
    Deleted:    code/logistic_susie_vb.R
    Deleted:    code/logistic_susie_veb_boost.R
    Modified:   code/marginal_sumstat_gsea.R
    Deleted:    code/simulate_gene_lists.R
    Modified:   code/tccm_ebnm.R
    Modified:   code/utils.R
    Deleted:    target_components/factories.R
    Deleted:    target_components/methods.R

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 (analysis/getting_started.Rmd) and HTML (docs/getting_started.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 814cab6 karltayeb 2022-05-16 wflow_publish(“analysis/getting_started.Rmd”)
html 03c133f karltayeb 2022-05-16 Build site.
Rmd a8038f9 karltayeb 2022-05-16 wflow_publish(“analysis/getting_started.Rmd”)

Introduction

Here we demonstrate the basics of how to use the gseasusie package.

devtools::install_github('karltayeb/gseasusie')
library(gseasusie)
library(tidyverse)

Loading gene sets

gseasusie ships with lots of gene set databases! There are other packages that curate geneset databases WebGestaltR, msigdbr, and stephenslab/pathways All I did here was write a few functions to load them into a standard format.

Each gene set list contains three elements: * X: a gene x gene set indicator matrx * geneSet: a two column data.frame (geneSet and gene) mapping genes to gene sets * geneSetDes: a written description of what the gene set does

gseasusie::load_gene_sets lets you load multiple gene sets. Note: the first time you try to load gene sets it might take a while because we need to build the indicator matrices (and the implimentation is far from optimized). After your fit it onces it will cache to a folder ./cache/resources/ relative to your working directory/project directory

eventually, the package will ship with some data

genesets <- gseasusie::load_gene_sets(c('c2', 'all_msigbd'))
loading gene set from msigdbr: C2
genesets$c2$geneSet %>% head()
# A tibble: 6 × 2
  geneSet gene  
  <chr>   <chr> 
1 M1423   79026 
2 M1423   214   
3 M1423   91369 
4 M1423   8289  
5 M1423   594   
6 M1423   146556
genesets$c2$geneSetDes %>% head()
# A tibble: 6 × 4
  geneSet gs_cat gs_subcat description                                          
  <chr>   <chr>  <chr>     <chr>                                                
1 M1423   C2     CGP       Genes down-regulated in AtT20 cells (pituitary cance…
2 M1458   C2     CGP       Genes up-regulated in AtT20 cells (pituitary cancer)…
3 M1481   C2     CGP       Genes down-regulated in GH3 cells (pituitary cancer)…
4 M1439   C2     CGP       Genes up-regulated in GH3 cells (pituitary cancer) a…
5 M2509   C2     CGP       Major ELAVL4 [GeneID=1996] associated mRNAs encoding…
6 M2096   C2     CGP       Genes down-regulated in the RCC4 cells (renal cell c…
dim(genesets$c2$X)
Loading required package: Matrix

Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':

    expand, pack, unpack
NULL

How to format data

Managing multiple experiment

If you have multiple experiments that you want to run enrichment for, I’ve found it useful to format as follows. data is a list of data frames, one for each experiment. The dataframes can have whatever information you’d like, but they must have the following columns: * ENTREZID: the gene sets above use ENTREZID, so map your gene names to this! * threshold.on: this might be a p-value, adjusted pvalue, effect size, etc * and beta:: right now this is only used for sign information, so if you don’t care specifically about up or down regulated genes just make a dummy column with 1s.

Eventually, beta and se columns will be important for the summary stat enrichment analaysis.

source('code/load_data.R')
data <- load_sc_pbmc_deseq2()
Loading required package: DESeq2
Loading required package: S4Vectors
Warning: package 'S4Vectors' was built under R version 4.1.3
Loading required package: methods
Loading required package: stats

Attaching package: 'stats'
The following objects are masked from 'package:Matrix':

    cov2cor, toeplitz, update
The following objects are masked from 'package:dplyr':

    filter, lag
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: graphics

Attaching package: 'graphics'
The following object is masked from 'package:stats4':

    plot
The following object is masked from 'package:Matrix':

    image

Attaching package: 'BiocGenerics'
The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs
The following objects are masked from 'package:dplyr':

    combine, intersect, setdiff, union
The following objects are masked from 'package:base':

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which.max, which.min

Attaching package: 'S4Vectors'
The following objects are masked from 'package:Matrix':

    expand, unname
The following objects are masked from 'package:dplyr':

    first, rename
The following object is masked from 'package:tidyr':

    expand
The following objects are masked from 'package:base':

    expand.grid, I, unname
Loading required package: IRanges

Attaching package: 'IRanges'
The following objects are masked from 'package:dplyr':

    collapse, desc, slice
The following object is masked from 'package:purrr':

    reduce
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: 'matrixStats'
The following object is masked from 'package:dplyr':

    count

Attaching package: 'MatrixGenerics'
The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Attaching package: 'Biobase'
The following object is masked from 'package:MatrixGenerics':

    rowMedians
The following objects are masked from 'package:matrixStats':

    anyMissing, rowMedians
'select()' returned 1:many mapping between keys and columns
data$`CD19+ B` %>% head()
          ENSEMBL  ENTREZID     baseMean log2FoldChange     lfcSE       pvalue
1 ENSG00000237683      <NA> 2.327879e-03     0.29159244 0.3257749 8.206061e-06
2 ENSG00000228463    728481 9.995371e-05     1.05065835 1.0336767 4.722355e-02
3 ENSG00000228327      <NA> 2.304559e-03    -0.47955339 0.4029420 5.056276e-02
4 ENSG00000237491 105378580 9.536208e-04     0.03026135 0.3774938 8.526449e-01
5 ENSG00000225880     79854 9.567947e-03     0.68820333 0.1604846 5.072878e-11
6 ENSG00000230368    284593 5.068608e-04     1.78293901 0.5604303 2.269509e-04
          padj        beta        se threshold.on
1 1.781359e-05  0.29159244 0.3257749 8.206061e-06
2 6.266937e-02  1.05065835 1.0336767 4.722355e-02
3 6.659266e-02 -0.47955339 0.4029420 5.056276e-02
4 8.561610e-01  0.03026135 0.3774938 8.526449e-01
5 1.568993e-10  0.68820333 0.1604846 5.072878e-11
6 4.258453e-04  1.78293901 0.5604303 2.269509e-04

Binarizing the data

At the end of the day, you’re going to need a binary gene x gene set matrix X and a binary gene list y. However you want to get those is fine.

One benefit of organizing your data as in the previous section, gseasusie ships with some helpful functions to format/prepare the data. gseasusie::prep_binary_data will binarize the data by thresholding on threshold.on.

db <- 'c2'  # name of gene set database to use in `genesets`
experiment = 'CD19+ B'  # name of experiment to use in `data`
thresh = 1e-4  # threshold for binarizing the data
bin.data <- gseasusie::prep_binary_data(genesets[[db]], data[[experiment]], thresh)

X <- bin.data$X
y <- bin.data$y

Fitting different enrichment models/methods

Now we can fit our enrichment models. NOTE: The marginal regressions are implemented in python (basilisk will spin up a conda environment with necessary dependencies, which will take some time the first time you run it, but will run quickly after!)

# fit logistic susie
logistic.fit <- gseasusie::fit_logistic_susie_veb_boost(X, y, L=20)
fitting logistic susie (via VEB.Boost
ELBO: -9062.304
9.019 sec elapsed
# fit linear susie
# (all of the functions that work with logistic susie fits should work with regular susie fits)
linear.fit <- susieR::susie(X, y)

# compute odds ratios, and pvalues under hypergeometric (one-sided) and fishers exact (two-sided) tests
ora <- gseasusie::fit_ora(X, y)
computin ORA statistics...
4.092 sec elapsed
# like ORA but performed implimented as a univariate logistic regression
marginal_regression <- gseasusie::fit_marginal_regression_jax(X, y)
fitting marginal logistic regressions...
34.585 sec elapsed
# redo the logistic regression ORA, conditional on enrichments found by logistic susie
# in practice, we introduce logistic susie predictions as an offset in the new model
residual_regression <- gseasusie::fit_residual_regression_jax(X, y, logistic.fit)
fitting residual logistic regressions...
29.936 sec elapsed
res = list(
  fit=logistic.fit,
  ora=ora,
  marginal_reg = marginal_regression,
  residual_reg = residual_regression
)

Visualizations

We can visualize our results with a volcano plot (more plots coming) The large circles highlight gene sets in a 95% credible set. The color indicates which SuSiE component)

gseasusie::enrichment_volcano(logistic.fit, ora)
Joining, by = "geneSet"
Joining, by = c("geneSet", "component")
Joining, by = c("geneSet", "component")
Joining, by = "geneSet"

If we account for the predictions made by (logitic) SuSiE, we can see that a lot (but not all of) the enrichment signal has been accounted for.

gseasusie::residual_enrichment_histogram(marginal_regression, residual_regression)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Interactive tables

We can produce an interactive table to explore the enrichment results:

gseasusie::interactive_table(logistic.fit, ora)
Joining, by = "geneSet"
Joining, by = c("geneSet", "component")
Joining, by = c("geneSet", "component")
Joining, by = "geneSet"
Joining, by = "geneSet"

sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] graphics stats4   stats    methods  utils    base    

other attached packages:
 [1] DESeq2_1.34.0               SummarizedExperiment_1.24.0
 [3] Biobase_2.54.0              MatrixGenerics_1.6.0       
 [5] matrixStats_0.62.0          GenomicRanges_1.46.1       
 [7] GenomeInfoDb_1.30.1         IRanges_2.28.0             
 [9] S4Vectors_0.32.4            BiocGenerics_0.40.0        
[11] Matrix_1.4-1                forcats_0.5.1              
[13] stringr_1.4.0               dplyr_1.0.9                
[15] purrr_0.3.4                 readr_2.1.2                
[17] tidyr_1.2.0                 tibble_3.1.7               
[19] ggplot2_3.3.6               tidyverse_1.3.1            
[21] gseasusie_0.0.0.9000       

loaded via a namespace (and not attached):
  [1] readxl_1.4.0           backports_1.4.1        workflowr_1.7.0       
  [4] plyr_1.8.7             splines_4.1.2          crosstalk_1.2.0       
  [7] BiocParallel_1.28.3    digest_0.6.29          htmltools_0.5.2       
 [10] fansi_1.0.3            magrittr_2.0.3         memoise_2.0.1         
 [13] tzdb_0.3.0             Biostrings_2.62.0      annotate_1.72.0       
 [16] modelr_0.1.8           colorspace_2.0-3       blob_1.2.3            
 [19] rvest_1.0.2            haven_2.5.0            xfun_0.30             
 [22] crayon_1.5.1           RCurl_1.98-1.6         jsonlite_1.8.0        
 [25] genefilter_1.76.0      survival_3.3-1         glue_1.6.2            
 [28] gtable_0.3.0           zlibbioc_1.40.0        XVector_0.34.0        
 [31] DelayedArray_0.20.0    RcppZiggurat_0.1.6     scales_1.2.0          
 [34] mvtnorm_1.1-3          DBI_1.1.2              data.tree_1.0.0       
 [37] Rcpp_1.0.8.3           xtable_1.8-4           reticulate_1.25       
 [40] bit_4.0.4              htmlwidgets_1.5.4      httr_1.4.3            
 [43] dir.expiry_1.2.0       RColorBrewer_1.1-3     ellipsis_0.3.2        
 [46] pkgconfig_2.0.3        reshape_0.8.9          XML_3.99-0.9          
 [49] farver_2.1.0           sass_0.4.1             dbplyr_2.1.1          
 [52] locfit_1.5-9.5         utf8_1.2.2             here_1.0.1            
 [55] tidyselect_1.1.2       labeling_0.4.2         rlang_1.0.2           
 [58] later_1.3.0            AnnotationDbi_1.56.2   reactR_0.4.4          
 [61] munsell_0.5.0          cellranger_1.1.0       tools_4.1.2           
 [64] cachem_1.0.6           cli_3.3.0              generics_0.1.2        
 [67] RSQLite_2.2.13         broom_0.8.0            evaluate_0.15         
 [70] fastmap_1.1.0          yaml_2.3.5             grDevices_4.1.2       
 [73] org.Hs.eg.db_3.14.0    knitr_1.39             bit64_4.0.5           
 [76] fs_1.5.2               KEGGREST_1.34.0        reactable_0.2.3       
 [79] whisker_0.4            tictoc_1.0.1           xml2_1.3.3            
 [82] compiler_4.1.2         rstudioapi_0.13        filelock_1.0.2        
 [85] susieR_0.11.92         png_0.1-7              VEB.Boost_0.0.0.9037  
 [88] spatstat.utils_2.3-0   reprex_2.0.1           geneplotter_1.72.0    
 [91] bslib_0.3.1            stringi_1.7.6          highr_0.9             
 [94] basilisk.utils_1.6.0   lattice_0.20-45        vctrs_0.4.1           
 [97] pillar_1.7.0           lifecycle_1.0.1        BiocManager_1.30.17   
[100] jquerylib_0.1.4        bitops_1.0-7           irlba_2.3.5           
[103] httpuv_1.6.5           R6_2.5.1               promises_1.2.0.1      
[106] renv_0.15.4            mr.ash.alpha_0.1-42    assertthat_0.2.1      
[109] rprojroot_2.0.3        withr_2.5.0            datasets_4.1.2        
[112] GenomeInfoDbData_1.2.7 parallel_4.1.2         hms_1.1.1             
[115] grid_4.1.2             basilisk_1.6.0         rmarkdown_2.14        
[118] Rfast_2.0.6            emulator_1.2-21        git2r_0.30.1          
[121] mixsqp_0.3-43          lubridate_1.8.0