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 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]
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 sood.s...@gmail.com javascript:
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,
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 sood.sanj...@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-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.
