ukbrapR (phonetically: 'U-K-B-wrapper') is an R package for working in the UK Biobank Research Analysis Platform (RAP). The aim is to make it quicker, easier, and more reproducible.
Since
v0.2.0ukbrapR works best on a "normal" cluster using RStudio and raw data from the table-exporter. Old Spark functions are still available but are not updated.
Wrapped server icon by DALL-E
In the DNAnexus Tools menu launch Posit Workbench and start an RStudio environment.
# install current version
remotes::install_github("lcpilling/ukbrapR")There are three main groups of functions:
- 🧬 Genetics: extract genotypes from Bulk data, create polygenic score
- 📋 Diagnoses: ascertain from health records and self-reported illness data, determine date first diagnosed
- 🛠️ Utilities: label UK Biobank data fields, upload/download files from RAP, and pull phenotypes from Spark
Bulk imputed genotypes and variant calls from Whole Genome Sequencing are available and can be easily accessed in an RStudio instance.
extract_variants() by default uses bgenix and plink to subset the imputed BGEN files and read it quickly and easily into R.
The only required input is a data frame (or path to file) containing "rsid" and "chr" variables. See function documentation for further details/options (including the available make_imputed_bed() and load_bed() internal functions).
varlist <- data.frame(rsid=c("rs1800562","rs429358"), chr=c(6,19))
imputed_genotypes <- extract_variants(varlist)
#> ~10 seconds
dim(imputed_genotypes)
#> [1] 487409 3By setting option source="dragen" the function will instead use tabix and plink to subset the DRAGEN WGS pVCF files. This requires "pos" in the input data frame (build 38). It is much slower than the imputed version, so unless you actually need the WGS calls it is not necessary.
varlist_b38 <- data.frame(rsid=c("rs1800562","rs429358"), chr=c(6,19), pos=c(26092913,44908684))
dragen_genotypes <- extract_variants(varlist_b38, source="dragen")
#> ~3 minutes (takes ~90 seconds per pVCF file)The highlight of developing this feature was naming the internal function make_dragen_bed() 🐉 🛏️
create_pgs() takes a data frame containing a list of variant associations with a trait and creates a weighted allele score using plink. By default it uses the imputed genotypes.
The only required input is a data frame (or path to file) containing rsid, chr, pos, effect_allele, other_allele, beta. For DRAGEN pos is build 38.
# weights from GWAS of liver cirrhosis (Innes 2020 Gastroenterology doi:10.1053/j.gastro.2020.06.014)
varlist_pgs <- readr::read_tsv(system.file("files", "pgs_liver_cirrhosis.txt", package="ukbrapR"))
head(varlist_pgs)
#> rsID CHR POS effect_allele other_allele effect_weight locus_name
#> <chr> <dbl> <dbl> <chr> <chr> <dbl> <chr>
#> 1 rs2642438 1 220796686 A G -0.177 MARC1
#> 2 rs11925835 3 56831417 T C -0.235 ARHGEF3
#> 3 rs72613567 4 87310241 TA T -0.166 HSD17B13
#> 4 rs2954038 8 125495147 C A 0.16 TRIB1
#> 5 rs11065384 12 120985482 T C 0.275 HNF1A
#> 6 rs28929474 14 94378610 T C 0.561 SERPINA1
liver_pgs <- create_pgs(
in_file=varlist_pgs, # can be a data frame or file path
out_file="liver_cirrhosis.imputed.pgs", # {optional} prefix for created .bed and .tsv files
pgs_name="liver_cirrhosis_pgs") # {optional} variable name
#> → Extracting 9 variants from 8 imputed files
#> ✔ PGS created! See file liver_cirrhosis.imputed.pgs.tsv
#> ~1 minute
summary(liver_pgs$liver_cirrhosis_pgs)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.00000 0.06006 0.08200 0.08589 0.10722 0.26639Diagnosis of conditions in UK Biobank participants come from multiple data sources. {ukbrapR} makes it fast and easy to ascertain diagnoses from multiple UK Biobank data sources in the DNAnexus Research Analysis Platform (RAP). Follow the below steps. See the website article for more details.
This only needs to happen once per project. Run export_tables() to submit the table-exporter jobs to save the required files to the RAP persistent storage. ~10Gb of text files are created, costing ~£0.15 per month to store.
For a given set of diagnostic codes get the participant Electronic Medical Records (EMR) and self-reported illess data. Returns a list containing up to 6 data frames: the subset of the clinical files with matched codes.
Codes need to be provided as a data frame with two fields: vocab_id and code. Valid code vocabularies are:
ICD10(for searching HES diagnoses, cause of death, and cancer registry)ICD9(for searching older HES diagnosis data)Read2andCTV3(for GP clinical events)OPCS3andOPCS4(for HES operations)ukb_cancerandukb_noncancer(for self-reported illness at UK Biobank assessments - all instances will be searched)
# example diagnostic codes for CKD
codes_df_ckd <- ukbrapR:::codes_df_ckd
head(codes_df_ckd)
#> condition vocab_id code
#> 1 ckd ICD10 N18.3
#> 2 ckd ICD10 N18.4
#> 3 ckd ICD10 N18.5
#> ...
# get diagnosis data - returns list of data frames (one per source)
diagnosis_list <- get_diagnoses(codes_df_ckd)
#> 7 ICD10 codes, 40 Read2 codes, 37 CTV3 codes
#> ~2 minutes
# N records for each source
nrow(diagnosis_list$gp_clinical) # 29,083
nrow(diagnosis_list$hesin_diag) # 206,390
nrow(diagnosis_list$death_cause) # 1,962Identify the date first diagnosed for each participant from any of datasets searched with get_diagnoses() (cause of death, HES diagnoses, GP clinical, cancer registry, HES operations, and self-reported illness fields).
Also included are:
- a
srcfield indicating the source of the date of first diagnosis. - a
binfield indicating the cases [1] and controls [0]. This relies on a small number of baseline fields also exported. Thedffield for the controls is the date of censoring (currently 30 October 2022). - a
bin_prevfield indicating whether the case was before the UK Biobank baseline assessment
# for each participant, get Date First diagnosed with the condition
# {optional} add a prefix to the variable names with "prefix"
diagnosis_df <- get_df(diagnosis_list, prefix="ckd")
#> ~2 seconds
# how many cases ascertained?
table(diagnosis_df$ckd_bin)
#> 0 1
#> 470334 31935
# source of earliest diagnosis date
table(diagnosis_df$ckd_src)
#> death gp hes selfrep_i0 selfrep_i1 selfrep_i2 selfrep_i3
#> 224 12394 19310 85 16 63 3
# date of diagnosis for prevalent cases (i.e., before UK Biobank baseline assessment)
summary(diagnosis_df$ckd_df[ diagnosis_df$ckd_bin_prev == 1 ])
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> "1958-01-01" "2006-06-21" "2007-01-12" "2006-06-24" "2007-11-19" "2010-06-16" The default get_df() behaviour is to use all available codes. However the most time-efficient way to get multiple conditions is to run get_diagnoses() once for all codes for the conditions you wish to ascertain, then get the "date first diagnosed" for each condition separately. In the codes data frame you just need a field indicating the condition name, that will become the variable prefixes.
# combine haemochromatosis and CKD codes together
# each contain there columns: condition, vocab_id, and code
# where `condition` is either "hh" or "ckd" and will become the variable prefix
codes_df_combined = rbind(ukbrapR:::codes_df_hh, ukbrapR:::codes_df_ckd)
# get diagnosis data - returns list of data frames (one per source)
diagnosis_list <- get_diagnoses(codes_df_combined)
# for each participant, get Date First diagnosed with the condition
diagnosis_df = get_df(diagnosis_list, group_by="condition")
# each condition has full set of output
table(diagnosis_df$hh_bin)
#> 0 1
#> 500254 2015
table(diagnosis_df$ckd_bin)
#> 0 1
#> 470334 31935 In the above example we also included a UK Biobank self-reported illness code for haemochromatosis, that was also ascertained (the Date First is run on each condition separately, they do not all need to have the same data sources).
- Label UK Biobank data fields with
label_ukb_fields() - Upload/download files between worker and RAP with
upload_to_rap()anddownload_from_rap() - Pull phenotypes from Spark instance with
get_rap_phenos()
Please report any bugs or issues, and feel free to suggest changes as pull requests. Alternatively, feel free to contact me via e-mail [email protected]