Contents

1 Overview

In this document we illustrate some approaches to working with UK Biobank summary statistics. Be sure that that the python module ukbb_pan_ancestry has been installed where reticulate can find it. (We don’t use basilisk as of 12/24/2022 because of issues in the terra spark cluster.)

library(BiocHail)

If the above command indicates that BiocHail is not available, see the Installation section near the end of this document.

1.1 Initialization and description

1.1.1 Standalone

We have produced a representation of summary statistics for a sample of 9888 loci. This 5GB resource can be retrieved and cached with the following code:

hl = hail_init()
ss = get_ukbb_sumstat_10kloci_mt(hl) # can take about a minute to unzip 5GB
ss$count()   # but if a persistent MatrixTable is at the location given
## [[1]]
## [1] 9888
## 
## [[2]]
## [1] 7271
             # by env var HAIL_UKBB_SUMSTAT_10K_PATH it goes quickly

To get a description of available content, we need a python chunk:

ss$describe()

1.1.2 Terra

Here’s a basic description of the summary stats table, with code that works in terra.bio:

hl = bare_hail()
hl$init(idempotent=TRUE, spark_conf=list(
  'spark.hadoop.fs.gs.requester.pays.mode'= 'CUSTOM',
  'spark.hadoop.fs.gs.requester.pays.buckets'= 'ukb-diverse-pops-public',
  'spark.hadoop.fs.gs.requester.pays.project.id'= Sys.getenv("GOOGLE_PROJECT")))

We would need to use a python chunk to get the output, using gs:// storage references.

r.hl.read_matrix_table('gs://ukb-diverse-pops-public/sumstats_release/results_full.mt').describe()

1.2 Exploring the subset

1.2.1 Metadata on study phenotypes

As of Aug 22 the code will not compile in a vignette, but it all runs in the interpreter. The chunks are not evaluated but the code is in test_ukbb.R in tests/testthat.

We’ll collect the column metadata to start to understand details of annotation of phenotypes.

sscol = ss$cols()$collect() # OK because we are just working over column metadata
ss1 = sscol[[1]]
names(ss1)
ss1$get("phenocode")
ss1$get("description")

We’ve introduced a function that extracts selected fields for a given phenotype, that accommodates availability of results for specific populations.

multipop_df(ss1)

This can be iterated over all the elements of sscol to produce a comprehensive searchable table. Here’s a small illustration:

library(DT)
ndf = do.call(rbind, lapply(sscol[1:10], multipop_df))
datatable(ndf)

1.2.2 Metadata on loci

We’ll trim the material to be worked with by sampling both rows and columns. (2023.01.08: In future revisions we will be able to control the seed for random sampling.)

sss = ss$sample_rows(.01)$sample_cols(.01)
sss$count()

Row metadata are simple to collect:

rss = sss$rows()$collect()
rss1 = rss[[1]]
names(rss1)
names(rss1$locus)
rss1$locus$contig
sapply(c("contig", "position"), function(x) rss1$locus[[x]])

1.2.3 Summary statistics

The summary statistics themselves reside in entries of the MatrixTable. This can be expensive to collect and so filtering methods beyond random sampling must be mastered. But here is a basic view.

sse = sss$entries()$collect()
length(sse)
names(sse[[1]])
sse1 = sse[[1]]

The summary_stats component has the association p-values – log10 transformed?

length(sse1$summary_stats)
names(sse1$summary_stats[[1]])
sse1$summary_stats[[1]]$Pvalue

2 Exercises

2.1 Infrastructure

  • Define an interface to this subset of 10k loci that supports queries like
    • has disease x been studied in UK Biobank?
    • how many phenotypes have been studied in K populations, K=2, 3, …?
    • how consistent is the annotation – are numbers of controls and cases always recorded?
  • Do the UK Biobank portals produce information to resolve these questions?
    • if so, what are the API calls to obtain answers?
    • if not, what is missing from the portals to allow answers to be obtained?

2.2 Substantive

  • List the genes that are near the loci collected in the 10k random sample
  • (Hard) What is the most significant finding (position, phenotype, population) in the 10k subset?

3 Installing BiocHail

BiocHail should be installed as follows:

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("BiocHail")

SessionInfo

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.4.3    BiocHail_1.0.1   BiocStyle_2.28.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.7            utf8_1.2.3            generics_0.1.3       
##  [4] RSQLite_2.3.1         lattice_0.21-8        digest_0.6.33        
##  [7] magrittr_2.0.3        evaluate_0.21         grid_4.3.1           
## [10] bookdown_0.35         fastmap_1.1.1         blob_1.2.4           
## [13] jsonlite_1.8.7        Matrix_1.6-1          DBI_1.1.3            
## [16] BiocManager_1.30.22   httr_1.4.7            purrr_1.0.2          
## [19] fansi_1.0.4           scales_1.2.1          jquerylib_0.1.4      
## [22] cli_3.6.1             rlang_1.1.1           dbplyr_2.3.3         
## [25] basilisk.utils_1.12.1 munsell_0.5.0         bit64_4.0.5          
## [28] withr_2.5.0           cachem_1.0.8          yaml_2.3.7           
## [31] tools_4.3.1           dir.expiry_1.8.0      parallel_4.3.1       
## [34] memoise_2.0.1         dplyr_1.1.2           colorspace_2.1-0     
## [37] filelock_1.0.2        basilisk_1.12.1       BiocGenerics_0.46.0  
## [40] curl_5.0.2            reticulate_1.31       png_0.1-8            
## [43] vctrs_0.6.3           R6_2.5.1              magick_2.7.5         
## [46] BiocFileCache_2.8.0   lifecycle_1.0.3       bit_4.0.5            
## [49] pkgconfig_2.0.3       gtable_0.3.4          pillar_1.9.0         
## [52] bslib_0.5.1           Rcpp_1.0.11           glue_1.6.2           
## [55] highr_0.10            xfun_0.40             tibble_3.2.1         
## [58] tidyselect_1.2.0      knitr_1.43            htmltools_0.5.6      
## [61] rmarkdown_2.24        compiler_4.3.1