On Fri, Jun 5, 2015 at 8:09 AM, sanj <sood.sanj...@gmail.com> wrote: > Hi Henrik, > > Thanks for replying. > > We are not completely sure if its the issue with the version of > aroma.affymetrix or not. > > We were trying to use the function mufColumn from FIRMAGene library on HTA > platform with Brain array gene level cdf. However it didn't work and we got > the following error: > > Error in mufC(x) : NA/NaN/Inf in foreign function call (arg 1) > > On exploring further we realised we had NaN in our residual object due to > which mufColumn was throwing an error. It seems like aroma.affymetrix > settings are creating issues in our analysis.
This looks very much related to thread 'mufColumns on genes with more than 500 probes across the gene' (https://groups.google.com/forum/#!topic/aroma-affymetrix/JaEz5h71yqg). Are you related? > > Since it is HTA platform so we have greater number of probes for some of the > units(genes with more then 500 probes/cells). It seems like for few of the > genes which have number of probes greater then 500 it does not estimate any > probe level data due to which residual object has NaNs in it (it is always > the last probe in these genes that has NaN residue). What we are not clear > is why does it work for some genes with more then 500 probes and not for > others? Do you have any idea about it? > > Also when we changed the aroma settings as follows: > > setOption(aromaSettings, "models/RmaPlm/medianPolishThreshold", c(3500,6)) > > It works fine and we are able to compute mufScores. Ok, so setting options works; it's all about the effect of setting this particular option. > Though we think we have > solved the error we are still not clear about what is exactly the reason for > this dependency on settings. Also we are not very clear about the > explanation given about aroma settings on the this link > http://www.aroma-project.org/settings/. So, you increasing 'models/RmaPlm/medianPolishThreshold' to c(3500,6) you're asking to only use median polish to step in for probesets where there are more than 3,500 probes. The default is that this happens for probesets with more than 500 probes. It looks like by increasing the option to a large enough threshold, you are simply avoiding using median polish at all. Internally, stats::medpolish(y, trace.iter=FALSE) and more importantly, it is using the default na.rm=FALSE. Now, if you pass NA to medpolish(..., na.rm=FALSE) it should throw error: Error in if (converged) break : missing value where TRUE/FALSE needed So, I'd expect that you'd see that also when calling: > fit(plmTr, verbose=verbose) I'm also confused why you get NA probe signals in the first place. I recommend that you try to trace where they come from, e.g. do you have them in your original CEL files, are introduced somehow in one of the preprocessing steps? Are the NAs across all arrays, or do they originate from a single array? The extractMatrix() method can be useful for this. /Henrik > > Following is the code: > > library(affxparser) > library(aroma.affymetrix) > library(FIRMAGene) > library(fdrtool) > library("limma") > library("statmod") > > verbose <- Arguments$getVerbose(-30); timestampOn(verbose) > > > setOption(aromaSettings, "models/RmaPlm/medianPolishThreshold", c(3500,6)) > setOption(aromaSettings, "memory/ram", 4.0) > > getOption(aromaSettings) > > #Provide the name of folder in rawData with CEL files > name <- 'TWIN' > > chipType <- 'hta20_Hs_ENSG_19' > > cdf <- AffymetrixCdfFile$byChipType(chipType, tags="binary") > print(cdf) > > cs <- AffymetrixCelSet$byName("TWIN", cdf=cdf,verbose=verbose) > setCdf(cs,cdf) > > # Background correction > bc <- RmaBackgroundCorrection(cs) > print(bc) > csBC <- process(bc,verbose=verbose) > > # Quantile normalization > qn <- QuantileNormalization(csBC, typesToUpdate="pm") > csN <- process(qn, verbose=verbose) > > # convert to unique cell > cdfu <- getUniqueCdf(cdf,verbose=-80) > csNU <- convertToUnique(csN,verbose=-20) > > # fit probe level model > plmTr <- ExonRmaPlm(csN, mergeGroups=TRUE) > print(plmTr) > fit(plmTr, verbose=verbose) > > # get Gene-level expression intensities > cesTr <- getChipEffectSet(plmTr) > trFit <- extractDataFrame(cesTr, addNames=TRUE) > > trFit <- data.frame(trFit$unitName, trFit[,6:ncol(trFit)]) > colnames(trFit)[colnames(trFit) == 'trFit.unitName'] <- 'gene' > > ###Calculate Residuals > rs <- calculateResidualSet(plmTr, verbose=verbose) > rsu <- readUnits(rs, verbose = verbose) > > nm <- getNames(cs) > cls <- 1:length(nm) > > cat("Extracting standardized residuals.\n") > > rsu1 <- lapply(rsu,FUN=function(u,cls) { > r <- log2(u[[1]]$eps) > m <- mad(r) > calcMeans(r/m,cls) > },cls=cls) > nProbe <- sapply(rsu,FUN=function(u) nrow(u[[1]]$eps)) > > rm(rsu) > gco <- gc() > > > upLimit <- 3485 # the length of unit with maximum number of probes > w <- (nProbe >= minProbes) & (nProbe <= upLimit) > nProbe <- nProbe[w] > > mufScores <- t(sapply( rsu1[w], FUN=function(u) mufColumns(u))) > > The session info is: > > R version 3.1.2 (2014-10-31) > Platform: x86_64-apple-darwin10.8.0 (64-bit) > > locale: > [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] stringr_1.0.0 statmod_1.4.21 limma_3.22.7 > fdrtool_1.2.14 FIRMAGene_0.9.8 aroma.light_2.2.1 > aroma.affymetrix_2.13.1 aroma.core_2.13.0 R.devices_2.13.0 > R.filesets_2.7.1 R.utils_2.0.2 > [12] R.oo_1.19.0 R.methodsS3_1.7.0 affxparser_1.38.0 > > loaded via a namespace (and not attached): > [1] aroma.apd_0.6.0 base64enc_0.1-3 digest_0.6.8 DNAcopy_1.40.0 > magrittr_1.5 matrixStats_0.14.0 PSCBS_0.44.0 R.cache_0.11.0 > R.huge_0.9.0 R.rsp_0.20.0 stringi_0.4-1 tools_3.1.2 > > Regards, > S > > > On Thursday, June 4, 2015 at 5:31:46 PM UTC+1, Henrik Bengtsson wrote: >> >> Hi, >> >> could you please share the exact calls you are using and elaborate a >> bit more what makes you think it's not working? That'll help me help >> you. >> >> Henrik >> >> On Thu, Jun 4, 2015 at 6:40 AM, sanjana <sood.s...@gmail.com> wrote: >> > Hi, >> > >> > Our group is working with HTA platform and using aroma.affymetrix for >> > the >> > analysis. We intend to change the RMA median polish threshold in aroma >> > settings. We are able to achieve this in version 2.13.1 but somehow not >> > in >> > version 2.13.2? >> > >> > Has anyone else experienced the same thing? Is it possible there is a >> > bug >> > in the 2.13.2 version? >> > >> > Thanks, >> > >> > -- >> > -- >> > 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 with website http://www.aroma-project.org/. >> > To post to this group, send email to aroma-af...@googlegroups.com >> > To unsubscribe and other options, go to >> > http://www.aroma-project.org/forum/ >> > >> > --- >> > You received this message because you are subscribed to the Google >> > Groups >> > "aroma.affymetrix" group. >> > To unsubscribe from this group and stop receiving emails from it, send >> > an >> > email to aroma-affymetr...@googlegroups.com. >> > For more options, visit https://groups.google.com/d/optout. > > -- > -- > 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 with website http://www.aroma-project.org/. > To post to this group, send email to aroma-affymetrix@googlegroups.com > To unsubscribe and other options, go to http://www.aroma-project.org/forum/ > > --- > You received this message because you are subscribed to the Google Groups > "aroma.affymetrix" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to aroma-affymetrix+unsubscr...@googlegroups.com. > For more options, visit https://groups.google.com/d/optout. -- -- 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 with website http://www.aroma-project.org/. To post to this group, send email to aroma-affymetrix@googlegroups.com To unsubscribe and other options, go to http://www.aroma-project.org/forum/ --- You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group. To unsubscribe from this group and stop receiving emails from it, send an email to aroma-affymetrix+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
