Hi,
to outrule that you get outdated cached results from runs with an old
version of aroma.affymetrix, you could always clear your
~/R.cache/aroma.affymetrix/<chip type>/ directory where <chip type>
specifies the particular chip type you are working with. Even better
is to delete all of ~/R.cache/aroma.affymetrix/. FYI, I try my best
to protect against this when code is updated, but bugs always slip
through. When you delete this directory, several things will have to
be reprocessed so it might take some time.
Cheers
/Henrik
On Fri, Sep 5, 2008 at 11:08 AM, Mark Robinson <[EMAIL PROTECTED]> wrote:
>
>
> Hi again Dick.
>
> I wasn't successful in reproducing the problem you have. I downloaded
> some MoGene data from Affy and processed it with:
>
> .... your code to get 'idx.MoGene' ...
>
> d<-doEverything("affytissue","MoGene-1_0-st-v1",doExon=F,getExpression =
> TRUE, units = idx.MoGene,verbose=-20)
>
> .. and it all ran no problem ...
>
>> d$d[1:5,1:2]
> MouseTP_Heart_01_mGENE MouseTP_Heart_02_mGENE
> 10338001 -Inf -Inf
> 10338003 -Inf -Inf
> 10338004 9.667838 9.427953
> 10338017 13.039836 12.586355
> 10338025 9.113231 8.620904
>
> ... so the units you excluded are represented by -Inf in the data table.
>
> For what its worth, below is my 'sessionInfo()', but there are only minor
> package differences to yours.
>
> Next step would be to run your code with a fresh R session (maybe set
> verbose=-20 to give more output) and give the result of a 'traceback()'.
>
> Cheers,
> Mark
>
>
>
>
> R version 2.7.1 (2008-06-23)
> i386-apple-darwin8.10.1
>
> locale:
> en_AU.UTF-8/en_AU.UTF-8/C/C/en_AU.UTF-8/en_AU.UTF-8
>
> attached base packages:
> [1] tools stats graphics grDevices utils datasets methods
> [8] base
>
> other attached packages:
> [1] preprocessCore_1.2.0 affyio_1.8.0 Biobase_2.1.7
> [4] aroma.affymetrix_0.9.4 aroma.apd_0.1.3 R.huge_0.1.6
> [7] affy_1.18.2 affxparser_1.12.2 aroma.core_0.9.4
> [10] sfit_0.1.5 aroma.light_1.8.1 digest_0.3.1
> [13] matrixStats_0.1.3 R.rsp_0.3.4 R.cache_0.1.7
> [16] R.utils_1.0.4 R.oo_1.4.5 R.methodsS3_1.0.3
>>
>
>
>
>
>>
>> Hi Mark,
>>
>> I did add those two arguments so as to have a little more flexibility.
>> Maybe my passing units to fit() is no longer correct. This code did
>> work ok several months ago :-)
>>
>> Maybe fit() modifies the plm object incorrectly now when I pass units
>> to it? Or, maybe the getChipEffectSet(plm) call does something
>> strange now? I added the units argument to have the fit() step not
>> take several hours (remove the probes with so many replicates, I
>> think). This was something I did after posting the long compute times
>> issue to the aroma.affymetrix group and someone responded with the
>> explanation about units and such.
>>
>>
>> doEverything<-function(name,chip,checkChipType=FALSE,
>> doExon=FALSE,doCore=FALSE,force=FALSE,doMerge=FALSE,doResiduals=FALSE,
>>
>> doWeights=FALSE,getExpression=FALSE,verbose=Arguments$getVerbose(-8),units=NULL,
>> doNorm=TRUE) {
>>
>> cs <- AffymetrixCelSet$fromName(name, chipType=chip,
>> checkChipType=checkChipType)
>> if (doCore)
>> chip<-paste(chip,",core",sep="")
>> cdf <- AffymetrixCdfFile$fromChipType(chip)
>> setCdf(cs, cdf)
>> csN <- cs
>> if(doNorm){
>> bc <- RmaBackgroundCorrection(cs)
>> csBC <- process(bc,verbose=verbose,force=force)
>> setCdf(csBC, cdf)
>> qn <- QuantileNormalization(csBC, typesToUpdate="pm")
>> csN <- process(qn, verbose=verbose,force=force) #time required
>> setCdf(csN, cdf)
>> }
>> if (doExon)
>> plm <- ExonRmaPlm(csN,mergeGroups=doMerge)
>> else
>> plm <- RmaPlm(csN)
>> fit(plm, verbose=verbose,units=units) #time required
>> qam <- QualityAssessmentModel(plm)
>> chp<-getChipEffectSet(plm)
>> d<-NULL; ws<-NULL; rs<-NULL
>> if (getExpression) {
>> if (!doExon | doMerge)
>> d<-getExprs(chp,doExon=F)
>> else
>> d<-getExprs(chp,doExon=T)
>> }
>> if (doResiduals)
>> rs<-calculateResiduals(plm,verbose=verbose,force=force)
>> if (doWeights)
>> ws<-calculateWeights(plm,verbose=verbose,force=force)
>>
>> list(name=name,chip=chip,d=d,cs=csN,cdf=cdf,plm=plm,qam=qam,chp=chp,res=rs,wts=ws)
>>
>> }
>>
>> ces <- chp
>> getExprs <- function(ces,doExon=FALSE, ...) {
>> cdf <- getCdf(ces);
>> theta <- extractMatrix(ces, ..., returnUgcMap=TRUE);
>> ugcMap <- attr(theta, "unitGroupCellMap");
>> if (doExon) {
>> x<-getFile(ces,1)
>> y<-readUnits(x)
>> cn<-names(unlist(y,recursive=F))
>> rownames(theta) <-cn
>> } else {
>> un<-getUnitNames(cdf, ugcMap[,"unit"]);
>> rownames(theta) <- un
>> }
>> log2(theta)
>>
>> }
>>
>> Again, I am very grateful for your help.
>>
>> Cheers,
>> Dick
>>
>>
>> On Thu, Sep 4, 2008 at 12:53 PM, Mark Robinson <[EMAIL PROTECTED]>
>> wrote:
>>>
>>>
>>> Hi Dick.
>>>
>>> Can you send the 'doEverything' function you are using. The one I have
>>> (the one posted a few months ago) doesn't even have arguments 'doNorm'
>>> and
>>> 'units' as you do. Is this something you've added ?
>>>
>>> Mark
>>>
>>>
>>>>
>>>> Hi Mark,
>>>>
>>>> Thanks for your help on this. I tried the cdf file you are using (I
>>>> was already using this version, but, just in case I was doing
>>>> something wrong I downloaded the one you suggested).
>>>>
>>>> I still get the same error as before:
>>>>
>>>> Error in list("doEverything("ChinMouseST_08.07.31",
>>>> "MoGene-1_0-st-v1", getExpression = TRUE, units = idx.MoGene, doNorm =
>>>> TRUE)" = <environment>, :
>>>>
>>>> [2008-09-04 10:11:27] Exception: Argument 'units' is of unknown
>>>> structure.
>>>> at throw(Exception(...))
>>>> at throw.default("Argument 'units' is of unknown structure.")
>>>> at throw("Argument 'units' is of unknown structure.")
>>>> at getCellMap.ChipEffectFile(this, units = units, ..., verbose =
>>>> less(verbose))
>>>> at getCellMap(this, units = units, ..., verbose = less(verbose))
>>>> at withCallingHandlers(expr, warning = function(w)
>>>> invokeRestart("muffleWarning"))
>>>> at suppressWarnings({
>>>> at getDataFlat.ChipEffectFile(cf, units = ugcMap, fields = field,
>>>> verbose = less(verbose))
>>>> at getDataFlat(cf, units = ugcMap, fields = field, verbose =
>>>> less(verbose))
>>>> at extractMatrix.ParameterCelSet(ces, ..., returnUgcMap = TRUE,
>>>> field = "theta")
>>>> at NextMethod("extractMatrix", this, ..., field = field)
>>>> at extractMatrix.ChipEffectSet(ces, ..., returnUgcMap = TRUE)
>>>> at extractMatrix(ces, ..., returnUgcMap = TRUE)
>>>> at getExprs(chp, doExon = F)
>>>> at doEverything("ChinMouseST_08.07.31", "MoGene-1
>>>>
>>>> and my sessionInfo() is:
>>>> R version 2.7.1 (2008-06-23)
>>>> x86_64-redhat-linux-gnu
>>>>
>>>> locale:
>>>> LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
>>>>
>>>> attached base packages:
>>>> [1] tools stats graphics grDevices datasets utils methods
>>>> base
>>>>
>>>> other attached packages:
>>>> [1] preprocessCore_1.2.1 affyio_1.8.1 Biobase_2.0.1
>>>> aroma.affymetrix_0.9.4 aroma.apd_0.1.3 R.huge_0.1.6
>>>> [7] affy_1.18.2 affxparser_1.12.2 aroma.core_0.9.4
>>>> sfit_0.1.5 aroma.light_1.8.1 digest_0.3.1
>>>> [13] matrixStats_0.1.3 R.rsp_0.3.4 R.cache_0.1.7
>>>> R.utils_1.0.4 R.oo_1.4.5 R.methodsS3_1.0.3
>>>>
>>>>
>>>> Is it possible the way I'm defining units is now wrong? Here is my
>>>> code for that:
>>>>
>>>> cdf <- AffymetrixCdfFile$fromChipType("MoGene-1_0-st-v1",tags="r3")
>>>> pathname <- getPathname(cdf);
>>>> unitGroupSizes <- readCdfNbrOfCellsPerUnitGroup(pathname);
>>>> names(unitGroupSizes) <- NULL;
>>>> unitGroupSizes <- lapply(unitGroupSizes, unname);
>>>> nbrOfCellsPerUnit <- sapply(unitGroupSizes, sum);
>>>> sizeDistr <- table(nbrOfCellsPerUnit);
>>>> print(sizeDistr);
>>>> #cat(sprintf("%6s: %7d", names(sizeDistr), sizeDistr), sep="\n");
>>>> idx <- which(unitGroupSizes>100)
>>>> names(idx) <- NULL
>>>> idx.MoGene <- setdiff(1:35512,idx)
>>>> length(idx.MoGene)
>>>>
>>>> I must of gotten this bit of code from someone else, as I'm sure I
>>>> wouldn't have thought this up on my own.
>>>>
>>>> I'd be very grateful for any help or pointers you might have.
>>>>
>>>> Thanks much,
>>>> Dick
>>>>
>>>> On Tue, Sep 2, 2008 at 8:33 AM, Mark Robinson <[EMAIL PROTECTED]>
>>>> wrote:
>>>>>
>>>>> Hi Dick.
>>>>>
>>>>> Sorry for the slow reply on this, I've been traveling for the past few
>>>>> days ...
>>>>>
>>>>> The code you ran works ok for me (in terms of getting plotRle/Nuse
>>>>> output) with data from the HuGene chip that I have. What CDF are you
>>>>> using for MoGene?
>>>>>
>>>>> Are you using the one from:
>>>>> http://groups.google.com/group/aroma-affymetrix/web/mogene-1-0-st-v1 ?
>>>>>
>>>>> If not, give that one a try. If so, I'll have to dig a little
>>>>> further.
>>>>>
>>>>> Let me know.
>>>>>
>>>>> Cheers,
>>>>> Mark
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On 30/08/2008, at 3:59 AM, [EMAIL PROTECTED] wrote:
>>>>>
>>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I am getting an error from code that used to work and was wondering
>>>>>> if
>>>>>> anyone might have a suggestion.
>>>>>>
>>>>>> I am using the doEverything function (code from Mark Robinson,
>>>>>> several
>>>>>> months ago), which works great with the following arguments:
>>>>>>
>>>>>> e.chinMus2 <- doEverything("ChinMouseST_08.07.31", "MoGene-1_0-
>>>>>> st-
>>>>>> v1", getExpression=FALSE, units=idx.MoGene, doNorm=TRUE)
>>>>>>
>>>>>> If I then try plotNuse:
>>>>>>
>>>>>> plotNuse(e.chinMus2$qam, ylim=c(0.9,1.2), pch=".", main="NUSE
>>>>>> ChinMouseST_08.07.31",xaxt="n",xlab=NULL)
>>>>>>
>>>>>> Error in list("plotNuse(e.chinMus2$qam, ylim = c(0.9, 1.2), pch =
>>>>>> ".",
>>>>>> main = "NUSE ChinMouseST_08.07.31", xaxt = "n", xlab = NULL)" =
>>>>>> <environment>, :
>>>>>>
>>>>>> [2008-08-29 10:36:00] Exception: Argument 'units' is of unknown
>>>>>> structure.
>>>>>> at throw(Exception(...))
>>>>>> at throw.default("Argument 'units' is of unknown structure.")
>>>>>> at throw("Argument 'units' is of unknown structure.")
>>>>>> at getCellMap.ChipEffectFile(this, units = units, ..., verbose =
>>>>>> less(verbose))
>>>>>> at getCellMap(this, units = units, ..., verbose = less(verbose))
>>>>>> at withCallingHandlers(expr, warning = function(w)
>>>>>> invokeRestart("muffleWarning"))
>>>>>> at suppressWarnings({
>>>>>> at getDataFlat.ChipEffectFile(this, units = ugcMap, fields = field,
>>>>>> verbose = less(verbose))
>>>>>> at getDataFlat(this, units = ugcMap, fields = field, verbose =
>>>>>> less(verbose))
>>>>>> at extractMatrix.ParameterCelFile(avg, field = "theta", units =
>>>>>> ugcMap, verbose = less(verbose, 5))
>>>>>> at NextMethod("extractMatrix", this, ..., field = field)
>>>>>> at extractMatrix.ChipEffectFile(avg, field = "theta", units =
>>>>>> ugcMap, verbose = less(verbose, 5))
>>>>>> at extractMatrix(avg, field = "theta", units = ugcMap, verbose =
>>>>>> less(v
>>>>>>
>>>>>> This is the same error I get when I try doEverything with
>>>>>> getExpression=TRUE, which calls from within doEverything:
>>>>>>
>>>>>> d<-getExprs(chp,doExon=F)
>>>>>>
>>>>>> and getExprs calls extractMatrix which is where the error gets
>>>>>> started
>>>>>> (I think):
>>>>>>
>>>>>> getExprs <- function(ces,doExon=FALSE, ...) {
>>>>>> cdf <- getCdf(ces);
>>>>>> theta <- extractMatrix(ces, ..., returnUgcMap=TRUE);
>>>>>> ugcMap <- attr(theta, "unitGroupCellMap");
>>>>>> if (doExon) {
>>>>>> x<-getFile(ces,1)
>>>>>> y<-readUnits(x)
>>>>>> cn<-names(unlist(y,recursive=F))
>>>>>> rownames(theta) <-cn
>>>>>> } else {
>>>>>> un<-getUnitNames(cdf, ugcMap[,"unit"]);
>>>>>> rownames(theta) <- un
>>>>>> }
>>>>>> log2(theta)
>>>>>>
>>>>>> }
>>>>>>
>>>>>> The units argument I pass to doEverything is what I do to get around
>>>>>> the long compute times (a problem I posted last Dec 13 and got lots
>>>>>> of
>>>>>> great help):
>>>>>>
>>>>>> cdf <- AffymetrixCdfFile$fromChipType(chip="MoGene-1_0-st-v1")
>>>>>> pathname <- getPathname(cdf);
>>>>>> unitGroupSizes <- readCdfNbrOfCellsPerUnitGroup(pathname);
>>>>>> names(unitGroupSizes) <- NULL;
>>>>>> unitGroupSizes <- lapply(unitGroupSizes, unname);
>>>>>> nbrOfCellsPerUnit <- sapply(unitGroupSizes, sum);
>>>>>> sizeDistr <- table(nbrOfCellsPerUnit);
>>>>>> print(sizeDistr);
>>>>>> idx <- which(unitGroupSizes>100)
>>>>>> names(idx) <- NULL
>>>>>> idx.MoGene <- setdiff(1:35512,idx)
>>>>>> length(idx.MoGene) # 35456
>>>>>>
>>>>>> I think all my packages are up to date:
>>>>>> sessionInfo()
>>>>>> R version 2.7.1 (2008-06-23)
>>>>>> x86_64-redhat-linux-gnu
>>>>>>
>>>>>> locale:
>>>>>> LC_CTYPE
>>>>>> =
>>>>>> en_US
>>>>>> .UTF
>>>>>> -8
>>>>>> ;LC_NUMERIC
>>>>>> =
>>>>>> C
>>>>>> ;LC_TIME
>>>>>> =
>>>>>> en_US
>>>>>> .UTF
>>>>>> -8
>>>>>> ;LC_COLLATE
>>>>>> =
>>>>>> en_US
>>>>>> .UTF
>>>>>> -8
>>>>>> ;LC_MONETARY
>>>>>> =
>>>>>> C
>>>>>> ;LC_MESSAGES
>>>>>> =
>>>>>> en_US
>>>>>> .UTF
>>>>>> -8
>>>>>> ;LC_PAPER
>>>>>> =
>>>>>> en_US
>>>>>> .UTF
>>>>>> -8
>>>>>> ;LC_NAME
>>>>>> =
>>>>>> C
>>>>>> ;LC_ADDRESS
>>>>>> =C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
>>>>>>
>>>>>> attached base packages:
>>>>>> [1] tools stats graphics grDevices datasets utils
>>>>>> methods base
>>>>>>
>>>>>> other attached packages:
>>>>>> [1] preprocessCore_1.2.1 affyio_1.8.1
>>>>>> Biobase_2.0.1 aroma.affymetrix_0.9.4 aroma.apd_0.1.3
>>>>>> R.huge_0.1.6
>>>>>> [7] affy_1.18.2 affxparser_1.12.2
>>>>>> aroma.core_0.9.4 sfit_0.1.5 aroma.light_1.8.1
>>>>>> digest_0.3.1
>>>>>> [13] matrixStats_0.1.3 R.rsp_0.3.4
>>>>>> R.cache_0.1.7 R.utils_1.0.4 R.oo_1.4.5
>>>>>> R.methodsS3_1.0.3
>>>>>>
>>>>>> I also tried passing the units argument to getExprs()
>>>>>>
>>>>>> d<-getExprs(chp,doExon=F, units=idx.MoGene)
>>>>>>
>>>>>> but that didn't help, same error (Exception: Argument 'units' is of
>>>>>> unknown structure).
>>>>>>
>>>>>> I assume I don't have the right version of some package or the code
>>>>>> that used to work now doesn't because I've done something wrong.
>>>>>> Maybe there is a more recent version of doEverything() and getExprs()
>>>>>> that doesn't have this problem?
>>>>>>
>>>>>> Ultimately I would like to use the QC calls such as plotNuse() and
>>>>>> plotRle() that I've found very useful in preprocessing affy 3' chips.
>>>>>>
>>>>>> I'd very much appreciate any help, pointers, suggestions anyone might
>>>>>> have.
>>>>>>
>>>>>> Thanks very much,
>>>>>> Dick
>>>>>>
>>>>>>
>>>>>> >
>>>>>
>>>>> ------------------------------
>>>>> Mark Robinson
>>>>> Epigenetics Laboratory, Garvan
>>>>> Bioinformatics Division, WEHI
>>>>> e: [EMAIL PROTECTED]
>>>>> e: [EMAIL PROTECTED]
>>>>> p: +61 (0)3 9345 2628
>>>>> f: +61 (0)3 9347 0852
>>>>> ------------------------------
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> >
>>>>>
>>>>
>>>> >
>>>>
>>>
>>>
>>>
>>> >
>>>
>>
>> >
>>
>
>
>
> >
>
--~--~---------~--~----~------------~-------~--~----~
When reporting problems on aroma.affymetrix, make sure 1) to run the latest
version of the package, 2) to report the output of sessionInfo() and
traceback(), and 3) to post a complete code example.
You received this message because you are subscribed to the Google Groups
"aroma.affymetrix" group.
To post to this group, send email to aroma-affymetrix@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/aroma-affymetrix?hl=en
-~----------~----~----~----~------~----~------~--~---
