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

[Gordon K Smyth]

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 intend...{{dropped:4}}