Search
Question: Error from eBayes in limma package
0
gravatar for calen.p.ryan
3 months ago by
calen.p.ryan0 wrote:

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
In addition: Warning messages:
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   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[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:
[1] VennDiagram_1.6.17  futile.logger_1.4.3 gplots_3.0.1        readstata13_0.9.0  
[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  
ADD COMMENTlink modified 3 months ago by James W. MacDonald45k • written 3 months ago by calen.p.ryan0
1
gravatar for James W. MacDonald
3 months ago by
United States
James W. MacDonald45k wrote:

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.

ADD COMMENTlink written 3 months ago by James W. MacDonald45k

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

ADD REPLYlink modified 3 months ago • written 3 months ago by calen.p.ryan0

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.

ADD REPLYlink written 3 months ago by James W. MacDonald45k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 212 users visited in the last hour