Question: frma normalization function not working on ExpressionFeatureSet object
0
gravatar for asalimih
6 weeks ago by
asalimih20
asalimih20 wrote:

Hello,
I have an ExpressionFeatureSet object which is created by ae2bioc function from ArrayExpress package. Now i want to perform a Frozen Robust Multiarray Analysis using frma function from frma package but i get this error:

object must be of class AffyBatch, ExonFeatureSet, or GeneFeatureSet.

as i wanted to change some of srdf files of ArrayExpress manually I downloaded all raw and annotation data and used getAE function following by ae2bioc function (instead of main ArrayExpress function). and as a result it provided me an ExpressionFeatureSet object instead of AffyBatch. So my question is:

  • can i convert ExpressionFeatureSet to any of those classes or is there any other way to solve my problem?

any help would be appreciated.

Here is my Code:

AEset = getAE("E-MEXP-950", path = "Data/E-MEXP-950" ,type='full',local = T,sourcedir = "Data/E-MEXP-950")
rawset= ae2bioc(mageFiles = AEset)
A = frma(rawset$`A-AFFY-33`)
B = frma(rawset$`A-AFFY-34`)
ADD COMMENTlink modified 6 weeks ago by Kevin Blighe50k • written 6 weeks ago by asalimih20

Please show some code that allows reproduction. [ Please read before posting a question ] -- How To Ask A Good Question

ADD REPLYlink written 6 weeks ago by ATpoint24k
1

I edited my post and added my code.

ADD REPLYlink written 6 weeks ago by asalimih20
2
gravatar for Kevin Blighe
6 weeks ago by
Kevin Blighe50k
Kevin Blighe50k wrote:

I can reproduce the error by running the following code:

require(ArrayExpress)
AEset <- getAE("E-MEXP-950", type='full')
rawset <- ae2bioc(mageFiles = AEset)

require(frma)
A <- rawset$`A-AFFY-33`
frma(A)

Error in frma(A) : 
  object must be of class AffyBatch, ExonFeatureSet, or GeneFeatureSet.

There is a way around this issue, though.

You should have the CEL files on your local drive. Even with my command, above, the raw files are downloaded automatically and their paths are stored here:

as.character(A@protocolData@data@.Data[[1]])
 [1] "C:/Users/kevin/Documents/PN5538_U133A_1_1_1.CEL"        "C:/Users/kevin/Documents/PN72_U133A_1_1_1.CEL"         
 [3] "C:/Users/kevin/Documents/kiel1_U133A_1_1_1.CEL"         "C:/Users/kevin/Documents/PDN15Normal_U133A_1_1_1.CEL"  
 [5] "C:/Users/kevin/Documents/PDN16Normal_U133A_1_1_1.CEL"   "C:/Users/kevin/Documents/PDN17Normal_U133A_1_1_1.CEL"  
 [7] "C:/Users/kevin/Documents/PDN8_3_U133A_1_1_1.CEL"        "C:/Users/kevin/Documents/PKN13Normal_U133A_1_1_1.CEL"  
 [9] "C:/Users/kevin/Documents/PKN15Normal_U133A_1_1_1.CEL"   "C:/Users/kevin/Documents/PKN16Normal_U133A_1_1_1.CEL"  
[11] "C:/Users/kevin/Documents/PKN9Normal_U133A_1_1_1.CEL"    "C:/Users/kevin/Documents/pt10_U113A_1_1_1.CEL"         
[13] "C:/Users/kevin/Documents/33B_U133A_1_1_1.CEL"           "C:/Users/kevin/Documents/35B_U133A_1_1_1.CEL"          
[15] "C:/Users/kevin/Documents/39B_U133A_1_1_1.CEL"           "C:/Users/kevin/Documents/PT55_U133A_1_1_1.CEL"         
[17] "C:/Users/kevin/Documents/56A_U133A_1_1_1.CEL"           "C:/Users/kevin/Documents/PDT121-2Tumor_U133A_1_1_1.CEL"
[19] "C:/Users/kevin/Documents/PDT14Tumor_U133A_1_1_1.CEL"    "C:/Users/kevin/Documents/PD51_3_U133A_1_1_1.CEL"       
[21] "C:/Users/kevin/Documents/PDT81-2_U133A_1_1_1.CEL"       "C:/Users/kevin/Documents/PDT91-2_U133A_1_1_1.CEL"      
[23] "C:/Users/kevin/Documents/PKT4_1Tumor_U133A_1_1_1.CEL"   "C:/Users/kevin/Documents/PKT5Tumor_U133A_1_1_1.CEL"    
[25] "C:/Users/kevin/Documents/PKT9Tumor_U133A_1_1_1.CEL"

So, you just need to read those guys / gals into an AffyBatch object via the affy package:

files <- as.character(A@protocolData@data@.Data[[1]])
myAffyObject <- affy::ReadAffy(filenames = files)

Then, you can use frma():

myAffyObjectNormalised <- frma::frma(myAffyObject)
myAffyObjectNormalised
ExpressionSet (storageMode: lockedEnvironment)
assayData: 22283 features, 25 samples 
  element names: exprs, se.exprs 
protocolData: none
phenoData
  sampleNames: PN5538_U133A_1_1_1.CEL PN72_U133A_1_1_1.CEL ... PKT9Tumor_U133A_1_1_1.CEL (25 total)
  varLabels: sample
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu133a

It should automatically detect the annotation that is required. When you run frma(), you may have to install a new package, hgu133afrmavecs, and re-run the command (I did).

Kevin

ADD COMMENTlink modified 6 weeks ago • written 6 weeks ago by Kevin Blighe50k

Thank you very much. yes it works but there is another problem. if you try your code on A-AFFY-34 you'll face an error saying hgu133bfrmavecs package is required. and actually we have no package named hgu133bfrmavecs .

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by asalimih20

That is unfortunate. I may take a look later, however, you could also try to email the package maintainer: https://bioconductor.org/packages/release/bioc/html/frma.html

ADD REPLYlink written 6 weeks ago by Kevin Blighe50k

I took another look but I am not sure how to circumvent that second issue. If you look at the hgu133afrmavecs package data, though, it seems to contain experimentally-derived values (I had planned to create a custom one for HGU133b, but now not sure):

require(hgu133afrmavecs)
data(hgu133afrmavecs)
str(hgu133afrmavecs)
List of 6
 $ normVec        : Named num [1:247965] 15.3 15.5 15.6 15.7 15.7 ...
  ..- attr(*, "names")= chr [1:247965] "129340" "213420" "396671" "82246" ...
 $ probeVec       : Named num [1:247965] -0.239 0.736 1.332 2.768 3.005 ...
  ..- attr(*, "names")= chr [1:247965] "129340" "213420" "396671" "82246" ...
 $ probeVarWithin : Named num [1:247965] 0.161 0.124 0.148 0.127 0.128 ...
  ..- attr(*, "names")= chr [1:247965] "129340" "213420" "396671" "82246" ...
 $ probeVarBetween: Named num [1:247965] 1.199 0.439 0.754 0.78 0.948 ...
  ..- attr(*, "names")= chr [1:247965] "129340" "213420" "396671" "82246" ...
 $ probesetSD     : Named num [1:22283] 1.008 0.619 0.568 0.455 0.334 ...
  ..- attr(*, "names")= chr [1:22283] "1007_s_at" "1053_at" "117_at" "121_at" ...
 $ medianSE       : Named num [1:22283] 0.196 0.171 0.184 0.235 0.219 ...
  ..- attr(*, "names")= chr [1:22283] "1007_s_at" "1053_at" "117_at" "121_at" ...

Unless you really need to use frma(), just use affy::rma(), which is AOK for both of these microarrays.

ADD REPLYlink written 6 weeks ago by Kevin Blighe50k

Thank you for your time and help. Here is the package maintainer's reply:

Unfortunately, I never created the vectors for the hgu133b platform. You could create your own using the frmaTools package, but that's fairly involved.

for now I'll just use rma. but I'll try frmaTools in future too.

ADD REPLYlink written 6 weeks ago by asalimih20

Great - thanks for the feedback. I took a quick look and it seems possible to create your own vectors. Take a look at this reproducible example:

require(ArrayExpress)
AEset <- getAE("E-MEXP-950", type='full')
rawset <- ae2bioc(mageFiles = AEset)

B <- rawset$`A-AFFY-34`
files <- as.character(B@protocolData@data@.Data[[1]])
myAffyObject <- affy::ReadAffy(filenames = files)

require(frmaTools)
hgu133bfrmavecs <- makeVectorsAffyBatch(
  files,
  batch.id = rep('batch1', length(files)))

require(frma)
myAffyObjectNormalised <- frma(
  myAffyObject,
  input.vecs = hgu133bfrmavecs)

The only problem is that I run into an error eventually when running makeVectorsAffyBatch(), relating to the medianSE vector. Let me know if you also run into an error.

ADD REPLYlink written 6 weeks ago by Kevin Blighe50k

yes i could reproduce the error :

Beginning Median SE Calculation ... 
Error in if (any(w < 0)) { : missing value where TRUE/FALSE needed
ADD REPLYlink written 6 weeks ago by asalimih20
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2047 users visited in the last hour