[BioC] Model Design

[James W. MacDonald]

Hi Sander,

On Aug 21, 2014 8:59 AM, "Sander [guest]" <guest at bioconductor.org> wrote:
>
> 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)

This is where you have made a mistake. Toy should be subtracting one
subject from another, depending on the contrasts you care about.

The resulting contrasts matrix should have both positive and negative
values, not just zeros and ones.

See the edgeR users guide or limma users guide for examples.

Best,

Jim

>           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.
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor

	[[alternative HTML version deleted]]