Re: [aroma.affymetrix] Changing the Aroma settings (RMA MedianPolish threshold ) in the latest version of the aroma.affymetrix
On Fri, Jun 5, 2015 at 8:09 AM, sanj 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 in
Re: [aroma.affymetrix] Changing the Aroma settings (RMA MedianPolish threshold ) in the latest version of the aroma.affymetrix
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.
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. 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/.
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.1R.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.0base64enc_0.1-3digest_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 >
> 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
> bu
Re: [aroma.affymetrix] Changing the Aroma settings (RMA MedianPolish threshold ) in the latest version of the aroma.affymetrix
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 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-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.
[aroma.affymetrix] Changing the Aroma settings (RMA MedianPolish threshold ) in the latest version of the aroma.affymetrix
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-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.
