Analysts may wish to further adjust data for biological covariates.
We provide an additional function,
recompute_derived_biomarkers() to recompute all composite
biomarkers and ratios from 107 non-derived biomarkers, which is useful
for ensuring data consistency when adjusting for unwanted biological
variation. A companion function,
recompute_derived_biomarker_qc_flags() will aggregate the
QC flags for the biomarkers underlying each composite biomarker and
ratio.
Note these functions assume the data has been returned to absolute
units after adjusting for covariates. For example the ratio of two
biomarkers A and B is computed as A/B, which may not be true if the two
biomarkers are on different scales (e.g. regression residuals) after
adjustment.
A worked example of recomputing derived biomarkers after adjusting
for age, sex, and BMI follows:
library(ukbnmr)
library(data.table) # for fast reading and writing of csv files using fread() and fwrite()
library(MASS) # Robust linear regression
# Load exported field data saved by the Table Exporter tool on the RAP
exported <- fread("path/to/exported_ukbiobank_phenotype_data.csv")
# First we need to remove the effects of technical variation before removing the
# biological covariates
processed <- remove_technical_variation(exported)
tech_qc <- processed$biomarkers
# Write out all the component results as above
fwrite(tech_qc, file="path/to/nmr_biomarker_data.csv")
fwrite(processed$biomarker_qc_flags, file="path/to/nmr_biomarker_qc_flags.csv")
fwrite(processed$sample_processing, file="path/to/nmr_sample_qc_flags.csv")
fwrite(processed$log_offset, file="path/to/nmr_biomarker_log_offset.csv")
fwrite(processed$outlier_plate_detection, file="path/to/outlier_plate_info.csv")
# Second, we can now create an additional dataset that has been adjusted for
# age, sex, and BMI.
# First, we need to extract the relevant fields from UK biobank and format as
# above, assuming the relevant fields for age (field #21003), sex (field #31),
# and BMI (field #21001) have also been saved by the Table Exporter tool in
# path/to/exported_ukbiobank_phenotype_data.csv
baseline_covar <- exported[, .(eid, age=p21003_i0, sex=p31_i0, bmi=p21001_i0)]
repeat_covar <- exported[, .(eid, age=p21003_i1, sex=p31_i1, bmi=p21001_i1)]
covar <- rbind(idcol="visit_index", "0"=baseline_covar, "1"=repeat_covar)
covar[, visit_index := as.integer(visit_index)]
# Combine covariates with processed NMR biomarker data. Filter the processed
# biomarkers only to the non-derived set, as we will plan to rederive the sums
# and ratios after adjusting for age sex and bmi
non_derived <- ukbnmr::nmr_info[Type == "Non-derived", Biomarker]
non_derived <- intersect(non_derived, names(nmr)) # In case some biomarkers are not in the data
bio_qc <- tech_qc[, .SD, .SDcols=c("eid", "visit_index", non_derived)]
bio_qc <- covar[bio_qc, on = .(eid, visit_index), nomatch=0]
# Transform to long format so we can operate on all biomarkers in bulk
bio_qc <- melt(bio_qc, id.vars=c("eid", "visit_index", "age", "sex", "bmi"),
variable.name="biomarker")
# Log transform biomarkers prior to adjustment. To do so, we need to add a small
# offset to any measurements equal to 0. This code is copied directly from the
# internals of the remove_technical_variation() function
log_offset <- bioqc[!is.na(value),
.(Minimum=min(value), Minimum.Non.Zero=min(value[value != 0])),
by=biomarker]
log_offset[, Log.Offset := ifelse(Minimum == 0, Minimum.Non.Zero / 2, 0)]
bio_qc[log_offset, on = .(biomarker), log_value := log(value + Log.Offset)]
# Now adjust for age, sex, and (log) BMI using robust linear regression from the MASS package
bio_qc[, adjusted := rlm(biomarker ~ age + factor(sex) + log(bmi))$residuals, by=biomarker]
# The adjusted residuals are normally distributed with mean 0. We can rescale to
# absolute units by subtracting the mean of the original log concentrations, then
# undoing the log transformation. This code is copied directly from the internals
# of the remove_technical_variation() function.
bio_qc[, adjusted := adjusted + as.vector(coef(rlm(value ~ 1)))[1], by=biomarker]
bio_qc[, adjusted := exp(adjusted)]
bio_qc[log_offset, on = .(biomarker), adjusted := adjusted - Log.Offset] # remove log offset
# Some values that were 0 are now < 0, apply small right shift for these biomarkers
# (the shift should be very small, essentially numeric error. This is worth double
# checking!). Again, this code is copied directly from the internals of the
# remove_technical_variation() function.
shift <- bio_qc[, .(Right.Shift=-pmin(0, min(adjusted))), by=biomarker]
log_offset <- log_offset[shift, on = .(biomarker)]
bio_qc[log_offset, on = .(biomarker), adjusted := adjusted + Right.Shift]
# Cast back to wide format so that each biomarker has its own column
bioqc <- dcast(bio_qc, eid + visit_index ~ biomarker, value.var="adjusted")
# Now we recompute the composite biomarkers and derived ratios and save the result.
bio_qc <- recompute_derived_biomarkers(bio_qc)
fwrite(bio_qc, file="path/to/nmr_biomarkers_adjusted_for_covariates.csv")
# You may also want to aggregate and save the quality control flags for each sample from
# the biomarkers underlying each derived biomarker or ratio, adding them as additional
# columns to the input data (see help("recompute_derived_biomarker_qc_flags")).
biomarker_qc_flags <- recompute_derived_biomarker_qc_flags(nmr)
fwrite(biomarker_qc_flags, file="path/to/biomarker_qc_flags.csv")
