[BioC] reading single channel Agilent data with limma [was arrayQualityMetrics doesn't work...]

[Wolfgang Huber]

Dear Gordon,

you are right, that use of the limma / read.maimages is rather odd. But 
from there Alex puts his data into an ExpressionSet 'esetPROC' and the 
arrayQualityMetrics error occurs with that.

Any analysis of these data (including arrayQualityMetrics) will only 
make sense after proper preprocessing, as you suggest, and this what 
Alex should do, and this is what the arrayQualityMetrics report should 
have told him to do.

Bottomline, my goal (to which we are close, but not there) is that 
arrayQualityMetrics should react gracefully even to the wildest 
instances of data, rather than stop with an error.

	Best wishes
	Wolfgang


Jun/10/12 5:12 AM, Gordon K Smyth scripsit::
> Hi Alex,
>
> I don't know arrayQualityMetrics, but you are using the limma package to
> read single-channel Agilent data in a way that I think might cause
> problems with down-stream analyses. Basically you're creating a
> two-color data object when your data is not actually of that type. This
> was a time when I suggested this sort of work-around as a stop-gap
> measure for some data problems, but hasn't been necessary for quite a
> few years.
>
> I'd also recommend that you do some background correction. If I
> understand your code correctly, I don't think it is currently making use
> of the background intensity column.
>
> There is a case study in the limma User's Guide that deals with single
> channel Agilent data. Could you please have a read of that for a cleaner
> way to read Agilent data?
>
> I don't know whether that will be enough to solve your
> arrayQualityMetrics problem, but perhaps it might.
>
> Best wishes
> Gordon
>
> ------------- original message -------------
> [BioC] arrayQualityMetrics() doesn't work for one-color non Affy arrays
> Alogmail2 at aol.com Alogmail2 at aol.com
> Fri Jun 8 09:39:21 CEST 2012
>
> Dear List,
>
> Could you share your experience with arrayQualityMetrics() for one-color
> non Affy arrays: it doesn't work for me (please see the code below).
>
> Thanks
>
> Alex Loguinov
>
> UC, Berkeley
>
>
>
>
>> options(error = recover, warn = 2)
>> options(bitmapType = "cairo")
>> .HaveDummy = !interactive()
>> if(.HaveDummy) pdf("dummy.pdf")
>
>> library("arrayQualityMetrics")
>
>> head(targets)
> FileName Treatment GErep Time Conc
> T0-Control-Cu_61_new_252961010035_2_4
> T0-Control-Cu_61_new_252961010035_2_4.txt C.t0.0 0 0 0
> T0-Control-Cu_62_new_252961010036_2_1
> T0-Control-Cu_62_new_252961010036_2_1.txt C.t0.0 0 0 0
> T0-Control-Cu_64_252961010031_2_2
> T0-Control-Cu_64_252961010031_2_2.txt C.t0.0 0 0 0
> T0-Control-Cu_65_new_252961010037_2_2
> T0-Control-Cu_65_new_252961010037_2_2.txt C.t0.0 0 0 0
> T04h-Contr_06_new_252961010037_2_4
> T04h-Contr_06_new_252961010037_2_4.txt C.t4.0 1 4 0
> T04h-Contr_10_new_252961010035_1_2
> T04h-Contr_10_new_252961010035_1_2.txt C.t4.0 1 4 0
>
>
>> ddaux = read.maimages(files = targets$FileName, source = "agilent",
> other.columns = list(IsFound = "gIsFound", IsWellAboveBG =
> "IsWellAboveBG",gIsPosAndSignif="gIsPosAndSignif",
> IsSaturated = "gIsSaturated", IsFeatNonUnifOF = "gIsFeatNonUnifOL",
> IsFeatPopnOL = "gIsFeatPopnOL", ChrCoord =
> "chr_coord",Row="Row",Column="Col"),
> columns = list(Rf = "gProcessedSignal", Gf = "gMeanSignal",
> Rb = "gBGMedianSignal", Gb = "gBGUsed"), verbose = T,
> sep = "\t", quote = "")
>
>
>> class(ddaux)
> [1] "RGList"
> attr(,"package")
> [1] "limma"
>> names(ddaux)
> [1] "R" "G" "Rb" "Gb" "targets" "genes" "source"
> "printer" "other"
>
>
> I could apply:
>>
>> class(ddaux$G)
> [1] "matrix"
>
>> all(rownames(targets)==colnames(ddaux$G))
> [1] TRUE
>
>> esetPROC = new("ExpressionSet", exprs = ddaux$G)
>
> But it results in errors:
>
>> arrayQualityMetrics(expressionset=esetPROC,outdir ="esetPROC",force =T)
>
> The directory 'esetPROC' has been created.
> Error: no function to return from, jumping to top level
>
> Enter a frame number, or 0 to exit
>
> 1: arrayQualityMetrics(expressionset = esetPROC, outdir = "esetPROC",
> force = T)
> 2: aqm.writereport(modules = m, arrayTable = x$pData, reporttitle =
> reporttitle, outdir = outdir)
> 3: reportModule(p = p, module = modules[[i]], currentIndex = currentIndex,
> arrayTable = arrayTableCompact, outdir = outdir)
> 4: makePlot(module)
> 5: print(_x at plot_ (mailto:x at plot) )
> 6: print.trellis(_x at plot_ (mailto:x at plot) )
> 7: printFunction(x, ...)
> 8: tryCatch(checkArgsAndCall(panel, pargs), error = function(e)
> panel.error(e))
> 9: tryCatchList(expr, classes, parentenv, handlers)
> 10: tryCatchOne(expr, names, parentenv, handlers[[1]])
> 11: doTryCatch(return(expr), name, parentenv, handler)
> 12: checkArgsAndCall(panel, pargs)
> 13: do.call(FUN, args)
> 14: function (x, y = NULL, subscripts, groups, panel.groups =
> "panel.xyplot", ..., col = "black", col.line = superpose.line$col,
> col.symbol =
> superpose.symb
> 15: .signalSimpleWarning("closing unused connection 5
> (Report_for_exampleSet/index.html)", quote(NULL))
> 16: withRestarts({
> 17: withOneRestart(expr, restarts[[1]])
> 18: doWithOneRestart(return(expr), restart)
>
> Selection: 0
>
>
> Error in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
> (converted from warning) Binning grid too coarse for current (small)
> bandwidth: consider increasing 'gridsize'
>
> Enter a frame number, or 0 to exit
>
> 1: arrayQualityMetrics(expressionset = esetPROC, outdir = "esetPROC",
> force = T)
> 2: aqm.writereport(modules = m, arrayTable = x$pData, reporttitle =
> reporttitle, outdir = outdir)
> 3: reportModule(p = p, module = modules[[i]], currentIndex = currentIndex,
> arrayTable = arrayTableCompact, outdir = outdir)
> 4: makePlot(module)
> 5: do.call(_x at plot_ (mailto:x at plot) , args = list())
> 6: function ()
> 7: meanSdPlot(x$M, cex.axis = 0.9, ylab = "Standard deviation of the
> intensities", xlab = "Rank(mean of intensities)")
> 8: meanSdPlot(x$M, cex.axis = 0.9, ylab = "Standard deviation of the
> intensities", xlab = "Rank(mean of intensities)")
> 9: smoothScatter(res$px, res$py, xlab = xlab, ylab = ylab, ...)
> 10: grDevices:::.smoothScatterCalcDensity(x, nbin, bandwidth)
> 11: KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, range.x
> = range.x)
> 12: warning("Binning grid too coarse for current (small) bandwidth:
> consider increasing 'gridsize'")
> 13: .signalSimpleWarning("Binning grid too coarse for current (small)
> bandwidth: consider increasing 'gridsize'", quote(KernSmooth::bkde2D(x,
> bandwidth = ba
> 14: withRestarts({
> 15: withOneRestart(expr, restarts[[1]])
> 16: doWithOneRestart(return(expr), restart)
>
> Selection: 0
>
>
>> sessionInfo()
> R version 2.14.2 (2012-02-29)
> Platform: i386-pc-mingw32/i386 (32-bit)
>
> locale:
> [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
> States.1252 LC_MONETARY=English_United States.1252
> [4] LC_NUMERIC=C LC_TIME=English_United
> States.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] CCl4_1.0.11 vsn_3.22.0
> arrayQualityMetrics_3.10.0 Agi4x44PreProcess_1.14.0 genefilter_1.36.0
> [6] annotate_1.32.3 AnnotationDbi_1.16.19 limma_3.10.3
> Biobase_2.14.0
>
> loaded via a namespace (and not attached):
> [1] affy_1.32.1 affyio_1.22.0 affyPLM_1.30.0
> beadarray_2.4.2 BiocInstaller_1.2.1 Biostrings_2.22.0
> [7] Cairo_1.5-1 cluster_1.14.2 colorspace_1.1-1
> DBI_0.2-5 grid_2.14.2 Hmisc_3.9-3
> [13] hwriter_1.3 IRanges_1.12.6 KernSmooth_2.23-7
> lattice_0.20-6 latticeExtra_0.6-19 plyr_1.7.1
> [19] preprocessCore_1.16.0 RColorBrewer_1.0-5 reshape2_1.2.1
> RSQLite_0.11.1 setRNG_2011.11-2 splines_2.14.2
> [25] stringr_0.6 survival_2.36-14 SVGAnnotation_0.9-0
> tools_2.14.2 XML_3.9-4.1 xtable_1.7-0
> [31] zlibbioc_1.0.1
>
> ______________________________________________________________________
> The information in this email is confidential and inte...{{dropped:19}}