Error from eBayes in limma package
1
0
Entering edit mode
@calenpryan-11107
Last seen 4.1 years ago

Hi experts,

I am trying to compare methylation for people with four combinations of childhood/adulthood socioeconomic status: 160 = LOW/LOW, 90 = LOW/HIGH, 90 = HIGH/LOW, 160 = HIGH/LOW.

I want to look for differences in DNAm between these groups.

My design matrix appears as follows:

      (Intercept) SEA_grphilo SEA_grplohi SEA_grplolo sex1
20867           1           0           0           0    1
20373           1           0           0           1    1
21987           1           0           0           0    1
20685           1           0           0           0    1
20827           1           0           0           0    1
22421           1           0           0           1    1

I run lmFit() no problem:

> fit.SEA_grp_chg <-lmFit(matrix, design, na.action = na.omit)

I run eBayes to calculate the variance:

> fit.SEA_grp_chg <-eBayes(fit.SEA_grp_chg, design)

But I get the following error:

Error in rep.int(0, ntarget) : invalid 'times' value
1: In if (ntarget < 1) return(NA) :
the condition has length > 1 and only the first element will be used
2: In 1:ntarget :
numerical expression has 2465 elements: only the first used
3: In 1:ntarget :
numerical expression has 2465 elements: only the first used

I don't know if this is something to be concerned about, and if/how it might be affecting my outcome, and can't seem to find anything about it.

---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago)

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] grid      parallel  stats     graphics  grDevices utils     datasets  methods
[9] base

other attached packages:
[5] Biobase_2.34.0      BiocGenerics_0.20.0 limma_3.30.13

loaded via a namespace (and not attached):
[1] Rcpp_0.12.11         knitr_1.15.1         magrittr_1.5         munsell_0.4.3
[5] colorspace_1.3-2     R6_2.2.0             plyr_1.8.4           dplyr_0.5.0
[9] caTools_1.17.1       tools_3.3.2          gtable_0.2.0         KernSmooth_2.23-15
[13] DBI_0.6              lambda.r_1.1.9       gtools_3.5.0         yaml_2.1.14
[17] lazyeval_0.2.0       assertthat_0.1       tibble_1.2           ggplot2_2.2.1
[21] futile.options_1.0.0 bitops_1.0-6         gdata_2.18.0         scales_0.4.1  
eBayes limma error • 981 views
1
Entering edit mode
@james-w-macdonald-5106
Last seen 59 minutes ago
United States

Why are you including your design matrix? The help page for eBayes says:

Usage:

ebayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4),
trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
eBayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4),
trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
treat(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))

Arguments:

fit: an 'MArrayLM' fitted model object produced by 'lmFit' or
'contrasts.fit'.  For 'ebayes' only, 'fit' can alternatively
be an unclassed list produced by 'lm.series', 'gls.series' or
'mrlm' containing components 'coefficients',
'stdev.unscaled', 'sigma' and 'df.residual'.

proportion: numeric value between 0 and 1, assumed proportion of genes
which are differentially expressed

stdev.coef.lim: numeric vector of length 2, assumed lower and upper
limits for the standard deviation of log2-fold-changes for
differentially expressed genes

trend: logical, should an intensity-trend be allowed for the prior
variance? Default is that the prior variance is constant.

robust: logical, should the estimation of 'df.prior' and 'var.prior'
be robustified against outlier sample variances?

winsor.tail.p: numeric vector of length 1 or 2, giving left and right
tail proportions of 'x' to Winsorize. Used only when
'robust=TRUE'.

lfc: the minimum log2-fold-change that is considered
scientifically meaningful

There is no mention of a design matrix there, and since you are using positional arguments, eBayes thinks your design matrix is the 'proportion' argument.

0
Entering edit mode

Thanks James, you just pointed out an obvious error in my code I was totally staring at yet still missing - now I feel silly!

0
Entering edit mode

By blocking on sex, you are accounting for female-specific differences. So one way you could interpret that contrast is that you are comparing the hilo and lohi samples after adjusting for sex. But this implies that any change in expression due to sex is simply a consistent shift in the expression level.

In other words, the 'Female' coefficient in this model is computing the mean difference between Female and Male subjects, and is using that to adjust the female observations down (or up) to the same level as the male subjects.

This assumes that for all the genes the only difference between males and females is that the gene expression for one sex is consistently shifted by a semi-constant amount, regardless of whatever hihi, hilo, etc are measuring.

It is always possible that a given gene reacts differently to changes in SEA_grph in females than in males, in which case you need to fit an interaction term between sex and whatever the SEA_grph thing is measuring, in order to capture any changes in gene expression between the SEA_grph levels that are sex-dependent.