[BioC] Model Design

Sander [guest] guest at bioconductor.org
Thu Aug 21 14:57:22 CEST 2014


Dear edgeR maintainers,

I have some troubles setting up the correct design for my DE experiment. I now how my experiment design should look like and that is like this: 
> design
      subject1 subject2 subject3 subject4
S7309        1        0        0        0
S7310        1        0        0        0
S7311        0        1        0        0
S7312        0        1        0        0
S7313        0        0        1        0
S7314        0        0        1        0
S7315        0        0        0        1
S7316        0        0        0        1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$subject
[1] "contr.treatment"

After that i make my contrasts:
> makeContrasts(subject1,subject2,subject3,subject4, levels = design)
          Contrasts
Levels     subject1 subject2 subject3 subject4
  subject1        1        0        0        0
  subject2        0        1        0        0
  subject3        0        0        1        0
  subject4        0        0        0        1

This all looks right since i need to compare the subjects with each other. But when i run my analysis function which looks like this:

library( edgeR );
library( ggplot2 );
library( reshape );
library( FactoMineR );

analyse <- function( counts, design, contrast, name, style ) {
  counts <- counts[ rowSums( counts, na.rm = TRUE ) > 0, ];
  y <- DGEList( counts = counts, genes = rownames( counts ) );
  y <- calcNormFactors( y );
  y <- estimateGLMCommonDisp( y, design );
  y <- estimateGLMTrendedDisp( y, design, df = 5 );
  y <- estimateGLMTagwiseDisp( y, design );
  
  fit <- glmFit( y, design );
  lrt <- glmLRT( fit, contrast = contrast );
  de <- decideTestsDGE( lrt, p = 0.05, adjust = "BH" );
  cpmY <- cpm( y );
  
  daf <- designAsFactor( design );
  orderedDesign <- design[ order( daf, names( daf ) ), ];
  tab <- data.frame(
    row.names = rownames( cpmY ),
    genes = rownames( cpmY ),
    de = de,
    cpmY[ ,order( daf, names( daf ) ) ]
  );
  
  aRepTab <- topTags( lrt, n = nrow( counts ) )$table;
  aRepTab$rank <- 1:nrow( counts );
#  repTab <- tab[ match( aRepTab$genes, rownames( tab ) ), ];
  
  repTab <- merge( aRepTab, tab, by = "genes", sort = FALSE );
  repTab <- repTab[ order( repTab$rank ), ];
#  data.frame(
#    row.names = rownames( aRepTab ),
#    aRepTab,
#    tab[ match( aRepTab$genes, tab$genes ), ]
#  )
  
  list(
    name = name,
    y = y,
    fit = fit,
    lrt = lrt,
    de = de,
    tab = tab,
    style = style,
    repTab = repTab,
    orderedDesign = orderedDesign
  );
}

I got the following error: 

Error in mglmLevenberg(y, design = design, dispersion = dispersion, offset = offset,  : 
  BLAS/LAPACK routine 'DGEMM ' gave error code -13 
5 mglmLevenberg(y, design = design, dispersion = dispersion, offset = offset, 
    weights = weights, coef.start = start, maxit = 250) 
4 glmFit.default(glmfit$counts, design = design0, offset = glmfit$offset, 
    weights = glmfit$weights, dispersion = glmfit$dispersion, 
    prior.count = 0) 
3 glmFit(glmfit$counts, design = design0, offset = glmfit$offset, 
    weights = glmfit$weights, dispersion = glmfit$dispersion, 
    prior.count = 0) 
2 glmLRT(fit, contrast = contrast) at diffexpr.R#15
1 analyse(counts, design, contrast, countsId, style) 

I tried to use different models but i cannot succeed to avoid the error for this comparison. (other comparisons do succeed) Any hints will be very appreciated.

Thanks in advance!

Best regards,

Sander

 -- output of sessionInfo(): 

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] FactoMineR_1.26 reshape_0.8.5   ggplot2_1.0.0   reshape2_1.4    edgeR_3.6.7     limma_3.20.8   

loaded via a namespace (and not attached):
 [1] car_2.0-20           cluster_1.15.2       colorspace_1.2-4     digest_0.6.4         grid_3.1.0           gtable_0.1.2        
 [7] htmltools_0.2.4      lattice_0.20-29      leaps_2.9            MASS_7.3-33          munsell_0.4.2        nnet_7.3-8          
[13] plyr_1.8.1           proto_0.3-10         Rcpp_0.11.2          rmarkdown_0.2.49     scales_0.2.4         scatterplot3d_0.3-35
[19] stringr_0.6.2        tools_3.1.0          yaml_2.1.13 

--
Sent via the guest posting facility at bioconductor.org.



More information about the Bioconductor mailing list