xwas demo
1. Setup
We’ll assume you have compiled the xwas library from latest Github source. If you are using our provided Docker container then you can assume this is true. The two commands below will install in your R environment if necessary.
> library(devtools)
> install_github('chiragjp/xwas')If you performed a git checkout of the code itself and ran your own build then the default build puts the library into ~/.R_libs/xwas/ and assumes your R library is set accordingly. The command below outputs the search path for R packages.
.libPaths()## [1] "/Users/npho/.R_libs"
## [2] "/Library/Frameworks/R.framework/Versions/3.2/Resources/library"Now you can load the xwas package and are good to go.
library(xwas)## Loading required package: parallel## Loading required package: survey## Loading required package: grid##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchart## Loading required package: survival## Warning: package 'survival' was built under R version 3.2.5##
## Attaching package: 'xwas'## The following object is masked from 'package:stats':
##
## qqplotFurther validate it was successfully loaded by checking the search() output.
search()## [1] ".GlobalEnv" "package:xwas" "package:survival"
## [4] "package:survey" "package:grid" "package:parallel"
## [7] "package:stats" "package:graphics" "package:grDevices"
## [10] "package:utils" "package:datasets" "package:methods"
## [13] "Autoloads" "package:base"2. Package Survey
The xwas library comes with some data, load included NHANES datasets for downstream analyses.
data(nhanes)
Figure 1: the two data.frames when you load the nhanes data set are MainTable and VarDescription.
The data set is divided among two data.frames as diagrammed in Figure 1. The MainTable has all the individuals represented by a row and columns with the co-variates. The VarDescription table has the description of what the co-variates are from the columns in MainTable.
3. Perform Regression Analysis
Let’s look at telomere length and see what environental exposures could influence it.
First, what is the outcome variable codified as? Let’s refer to the VarDescription table.
VarDescription[grep("telomere", VarDescription$var_desc, ignore.case=T), ]## tab_name tab_desc series module
## 31731 telo_a Telomere Mean and Standard Deviation 1999-2000 laboratory
## 31741 telo_b Telomere Mean and Standard Deviation 2001-2002 laboratory
## var var_desc analyzable binary_ref_group
## 31731 TELOMEAN Mean Telomere Length 1 NA
## 31741 TELOMEAN Mean Telomere Length 1 NA
## comment_var version_date is_comment is_weight is_questionnaire
## 31731 <NA> 2014-11-19 0 0 0
## 31741 <NA> 2014-11-19 0 0 0
## is_ecological is_binary is_ordinal categorical_ref_group
## 31731 0 0 0 NA
## 31741 0 0 0 NA
## categorical_levels N age_range_low age_range_high age_median
## 31731 <NA> 7827 19 85 47
## 31741 <NA> 7827 19 85 47
## female_N female_pct white_N white_pct mexican_N mexican_pct black_N
## 31731 4056 0.5182062 3965 0.5065798 1876 0.2396831 1333
## 31741 4056 0.5182062 3965 0.5065798 1876 0.2396831 1333
## black_pct other_hispanic_N other_hispanic_pct other_N other_pct
## 31731 0.1703079 417 0.05327712 236 0.03015204
## 31741 0.1703079 417 0.05327712 236 0.03015204
## series_num category sub_category
## 31731 1 aging <NA>
## 31741 2 aging <NA>Mean telomere length is codified as TELOMEAN in the MainTable by looking at the VarDescription$var column. From the VarDescription$series column we can also see that the 1999-2000 and 2001-2002 cohorts are associated with this measurement.
Now lets conduct an EWAS on TELOMEAN, focusing only on variables present in these two cohorts where the mean telomere length (TELOMEAN) was measured.
VarDescription.telo <- subset(VarDescription, series == '1999-2000' | series == '2001-2002')Quick check on what this looks like by tabulating the data.frame.
sort(table(VarDescription.telo$category))##
## aging cotinine cognitive functioning
## 2 2 4
## sexual behavior immunization social support
## 4 6 6
## alcohol use blood pressure street drug
## 8 8 8
## polyflourochemicals phytoestrogens dioxins
## 11 12 13
## diakyl hormone smoking family
## 14 14 16
## housing furans phthalates
## 18 19 19
## nutrients viral infection physical fitness
## 26 26 30
## hydrocarbons body measures blood
## 37 38 40
## heavy metals occupation volatile compounds
## 42 42 51
## bacterial infection smoking behavior pcbs
## 53 55 59
## pesticides disease biochemistry
## 64 80 93
## supplement use pharmaceutical food component recall
## 128 252 256pie(sort(table(VarDescription.telo$category)))
The top two categories are “food component recall” and “pharmaceutical” with 256 and 252 entries, respectively. Let us focus our EWAS on the environmental exposures.
useTheseCategories <- c('dioxins', 'furans', 'diakyl', 'heavy metals', 'hormone', 'hydrocarbons', 'pcbs', 'pesticides', 'physical fitness', 'phytoestrogens', 'polyflourochemicals', 'smoking behavior', 'smoking family', 'volatile compounds')
useTheseCategories <- sort(useTheseCategories) # store pre-sorted for easier processing
VarDescription.telo <- VarDescription.telo[VarDescription.telo$category %in% useTheseCategories, ]Another inspection will confirm that we’re looking at much more tractable set of data.
pie(sort(table(VarDescription.telo$category)))
sort(table(VarDescription.telo$category))##
## polyflourochemicals phytoestrogens dioxins
## 11 12 13
## diakyl hormone smoking family
## 14 14 16
## furans physical fitness hydrocarbons
## 19 30 37
## heavy metals volatile compounds smoking behavior
## 42 51 55
## pcbs pesticides
## 59 64The top two categories are now “pesticides” and “pcbs” with 64 and 59 entries, respectively. We’re close to performing the actual EWAS but will now define specific co-variates to adjust for in each regression. The vector is defined below as adjustfor and then the combined with unique co-variates in our cohort as well as the dependent variable of mean telomere length we are interested in studying (TELOMEAN).
We are also going to derive a new variable RIDAGEYR2 which is the age-squared to test for non-linear effects associated with age.
MainTable$RIDAGEYR2 <- MainTable$RIDAGEYR^2The newData variable will contain all the pertinent records for analysis. However, the data is right-skewed so we both scale and log transform the data into a new object.
adjustfor <- c('RIDAGEYR', 'RIDAGEYR2', 'female', 'mexican', 'black', 'other_eth', 'other_hispanic', 'INDFMPIR')
variablesToKeep <- c(unique(VarDescription.telo$var), 'TELOMEAN', adjustfor, 'SDMVPSU', 'SDMVSTRA', 'WTMEC2YR', 'WTMEC4YR')
newData <- subset(MainTable, SDDSRVYR <= 2)[, variablesToKeep] ## only keep those that are
newLogData <- cbind(newData[, is.categorical(newData, lower=0)], apply(newData[, !is.categorical(newData, lower=0)], 2, function(x) { scale(log(x+1e-10)) })) # log transform the dataThe xwas library also has integrated support for the survey package to perform survey-weighted analysis. You will have to define your own survey::svydesign objects but it’s taken as an argument to most xwas functions. We define our study in the dsn variable then perform an EWAS both with and with-out the survey design for comparison.
dsn <- svydesign(ids=~SDMVPSU, strata=~SDMVSTRA, weight=~WTMEC4YR, nest=T, data=subset(newLogData, WTMEC4YR > 0)) # create the design object
telo.norm <- xwas(data=subset(newLogData, WTMEC4YR > 0), depvar="TELOMEAN", adjvars=adjustfor, verbose=F) # ewas WITHOUT survey design## <simpleError in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...): 0 (non-NA) cases>
## <simpleError in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...): 0 (non-NA) cases>telo.surv <- xwas(data=subset(newLogData, WTMEC4YR > 0), depvar="TELOMEAN", adjvars=adjustfor, design=dsn, verbose=F) # ewas WITH survey design## Warning in glm.fit(x = structure(numeric(0), .Dim = c(0L, 10L), .Dimnames =
## list(: no observations informative at iteration 1## Warning: glm.fit: algorithm did not converge## <simpleError in glm.fit(x = structure(numeric(0), .Dim = c(0L, 10L), .Dimnames = list( NULL, c("(Intercept)", "I(scale(SMD100MN))", "RIDAGEYR", "RIDAGEYR2", "female", "mexican", "black", "other_eth", "other_hispanic", "INDFMPIR")), assign = 0:9), y = structure(numeric(0), class = "AsIs", "`scaled:center`" = 0.00244357177919902, "`scaled:scale`" = 0.999429581414316), weights = numeric(0), start = NULL, etastart = NULL, mustart = NULL, offset = NULL, family = structure(list(family = "gaussian", link = "identity", linkfun = function (mu) mu, linkinv = function (eta) eta, variance = function (mu) rep.int(1, length(mu)), dev.resids = function (y, mu, wt) wt * ((y - mu)^2), aic = function (y, n, mu, wt, dev) { nobs <- length(y) nobs * (log(dev/nobs * 2 * pi) + 1) + 2 - sum(log(wt)) }, mu.eta = function (eta) rep.int(1, length(eta)), initialize = expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$link == "inverse" && any(y == 0)) || (family$link == "log" && any(y <= 0)))) stop("cannot find valid starting values: please specify some") mustart <- y }), validmu = function (mu) TRUE, valideta = function (eta) TRUE), .Names = c("family", "link", "linkfun", "linkinv", "variance", "dev.resids", "aic", "mu.eta", "initialize", "validmu", "valideta"), class = "family"), control = structure(list( epsilon = 1e-08, maxit = 25, trace = FALSE), .Names = c("epsilon", "maxit", "trace")), intercept = TRUE): object 'fit' not found>## Warning in glm.fit(x = structure(numeric(0), .Dim = c(0L, 10L), .Dimnames =
## list(: no observations informative at iteration 1
## Warning in glm.fit(x = structure(numeric(0), .Dim = c(0L, 10L), .Dimnames =
## list(: glm.fit: algorithm did not converge## <simpleError in glm.fit(x = structure(numeric(0), .Dim = c(0L, 10L), .Dimnames = list( NULL, c("(Intercept)", "I(scale(SMD100FL))", "RIDAGEYR", "RIDAGEYR2", "female", "mexican", "black", "other_eth", "other_hispanic", "INDFMPIR")), assign = 0:9), y = structure(numeric(0), class = "AsIs", "`scaled:center`" = 0.00244357177919902, "`scaled:scale`" = 0.999429581414316), weights = numeric(0), start = NULL, etastart = NULL, mustart = NULL, offset = NULL, family = structure(list(family = "gaussian", link = "identity", linkfun = function (mu) mu, linkinv = function (eta) eta, variance = function (mu) rep.int(1, length(mu)), dev.resids = function (y, mu, wt) wt * ((y - mu)^2), aic = function (y, n, mu, wt, dev) { nobs <- length(y) nobs * (log(dev/nobs * 2 * pi) + 1) + 2 - sum(log(wt)) }, mu.eta = function (eta) rep.int(1, length(eta)), initialize = expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$link == "inverse" && any(y == 0)) || (family$link == "log" && any(y <= 0)))) stop("cannot find valid starting values: please specify some") mustart <- y }), validmu = function (mu) TRUE, valideta = function (eta) TRUE), .Names = c("family", "link", "linkfun", "linkinv", "variance", "dev.resids", "aic", "mu.eta", "initialize", "validmu", "valideta"), class = "family"), control = structure(list( epsilon = 1e-08, maxit = 25, trace = FALSE), .Names = c("epsilon", "maxit", "trace")), intercept = TRUE): object 'fit' not found>As you can see the EWAS with survey-weighted analysis takes longer to perform but it is technically more accurate as you can see within the visualizations below.
4. Visualization
Now that we have performed the heavy lifting of EWAS analysis we can use the packaged visualization functions, specifically manhattan, volcano, and qqplot to gather insights into the relationship of environmental exposures to mean telomere length (TELOMEAN).
par(mfrow=c(1,1), xpd=TRUE) # setup to plot outside of area if necessaryFirst look at a Manhattan plot of the telomere length without survey-weighted analysis.
manhattan(telo.norm)
Since we looked at the co-variates by sub-categories of environmental exposure it also makes sense to show this visually. If we provided a mapping of variables to categories the Manhattan plot can show this.
groups <- subset(VarDescription, category %in% useTheseCategories)[, c("category", "var")]
groups <- unique(groups)Now by passing in this additional information we can show a more visually clear relationship. Let’s look at the telomere study without survey-weighted analysis.
manhattan(telo.norm, group.by=groups, legend=T, xlab="", ylim=c(0,14), nlabels=5)
Now we can compare this to the study done with survey-weighted analysis. You can see that there is less noise in the results and the significant co-variates are clearer.
manhattan(telo.surv, group.by=groups, legend=T, xlab="", ylim=c(0,14), nlabels=5)
Given the results, let’s convert to a more readable table and codify all the variables to their descriptions. We’ll just print out the top 10 results by significance.
top <- xwas.to.df(telo.surv) # convert xwas results to a data.frame
# create a table with the labels and descriptions for each variable used
labels <- subset(VarDescription, category %in% useTheseCategories)[, c("var", "var_desc", "category")]
labels <- labels[!duplicated(labels$var), ]
rownames(labels) <- labels$var
# match variables from the xwas output to the labels object
top <- merge(top, labels, all.x=TRUE, by="row.names")
rownames(top) <- top$Row.names
top <- top[, -which(colnames(top) %in% c("Row.names", "var"))]
top <- top[order(top$q), ]
# show the top 10 results with p-values and descriptions
head(top, n=10)## p x q
## LBXBCD 0.0001500260 -0.06771447 0.03675636
## LBX138 0.0001689830 0.14720785 0.04123185
## LBXHXC 0.0003406962 0.14354116 0.08278917
## LBX180 0.0004093170 0.16330759 0.09905470
## WTMEC4YR 0.0004678815 0.12921712 0.11275945
## CVDVOMAX 0.0005148082 0.09742866 0.12355396
## LBX099 0.0006243116 0.11235434 0.14921046
## LBX153 0.0008250694 0.15831092 0.19636652
## URXACE 0.0009448780 0.07508905 0.22393607
## LBX170 0.0010012468 0.11929804 0.23629425
## var_desc category
## LBXBCD Cadmium (ug/L) heavy metals
## LBX138 PCB138 & 158 (ng/g) pcbs
## LBXHXC 3,3,4,4,5,5-hxcb (fg/g) pcbs
## LBX180 PCB180 (ng/g) pcbs
## WTMEC4YR <NA> <NA>
## CVDVOMAX Predicted VO2max (ml/kg/min) physical fitness
## LBX099 PCB99 (ng/g) pcbs
## LBX153 PCB153 (ng/g) pcbs
## URXACE Acetochlor mercapturate (ug/L) result pesticides
## LBX170 PCB170 (ng/g) pcbsAnother xwas visualization is the volcano plot that allows you to quickly discern both significance and effect size.
volcano(telo.norm, xlim=c(-.2,.2), xlab="Association Size")
Notice how the visualizations change with the survey-weighted analysis results.
volcano(telo.surv, xlim=c(-.2,.2), xlab="Association Size")
We can also compare our output above to Figure 2 from the paper we’re replicating.
Figure 2: this is grabbed from Figure 2 from the paper we’re replicating for comparison, link above.
An xwas visualization to determine significance is the qqplot. Any deviation from the line y=x indicates true signal in the associations.
qqplot(telo.norm, xlim=c(0,14), ylim=c(0,14))
Now with the survey-weighted data you can see that the signal:noise ratio is much more distinct.
qqplot(telo.surv, xlim=c(0,14), ylim=c(0,14))
Additional visualizations and other features coming soon in future releases.