Question: lesser fold change value using RMA background correction
gravatar for Paul
2.1 years ago by
Paul80 wrote:

I was trying to do background correction using RMA for my GEO dataset, since the dataset was produced using Affymetrix array (GPL570). And preprocessing using limma package in R.

I got this tutorial for normalization From Data import to Normalization in Microarray Analysis using in R (Part I)

But unable to know, how to process this normalized output and give it as an input to identify the differentially expressed genes using limma. I have written the following code, but it returns 0 deg (topTable)



    data.rma <- expresso(d1,bgcorrect.method="rma",normalize.method="quantiles",pmcorrect.method="pmonly",summary.method="medianpolish")
    eset <- exprs(data.rma)

    sample <- factor(rep(c("Case","Cont"), each = 10))
    design.mat <- model.matrix(~0+sample)
    colnames(design.mat) <- levels(sample)

    fit <- lmFit(eset,design.mat)
    fit3 <- eBayes(fit)
    deg <- topTable(fit3, coef = 2, p.value = 0.05, adjust.method = 'BH', number = nrow(eset), lfc > 2)

When I remove lfc criteria, I get genes with very less fold change value, none are above 2. On the contrary, when I perform mas5 background correction I get high fold change values.. How do I get the actual fold change values?

ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by Paul80

You are not reading in any data (?) Which GEO dataset are you using?

ADD REPLYlink written 2.1 years ago by Kevin Blighe65k

I have set the path before d1<-ReadAffy(). d1 is not empty. I am using GSE4757 dataset

ADD REPLYlink written 2.1 years ago by Paul80

Okay. What is the lowest P value obtained when you just run:

topTable(fit3, coef = 'Diff', p.value = 1, adjust.method = 'fdr', lfc = 0, number = nrow(eset.rma))

Also, the argument that you pass to lfc is expected to already be on the log (base 2) scale, i.e., you don't have to log the value yourself with lfc = log2(1.5) (it would just be lfc = 1.5)

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by Kevin Blighe65k

I have modified the code a bit, I am not getting a proper fold change values.. using RMA my fold change is going down and using Mas5 my fold change values are going high. How do I correct the fold change values?

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by Paul80

Can you actually show the fold-changes obtained by the different normalisation strategies?, like, plot the sorted / ordered fold changes with plot(). Can you also show the boxplot of your normalised data?

Go here to learn how to share images: How to add images to a Biostars post

Also, you should consider using oligo, and not affy. For some Affymetrix arrays, oligo is actually the better option.

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by Kevin Blighe65k

I couldn't plot the whole logFC, but following is the head(genes) in descending order of logFC

                logFC   AveExpr             t     P.Value       adj.P.Val   B
223500_at   1.338822221 9.704255821 4.121405486 0.000474727 0.107099415 -0.000301553
231783_at   1.26555711  7.575587468 2.809189993 0.010424161 0.194320113 -2.682451156
205924_at   1.229768314 6.736765908 3.507544701 0.002062693 0.130621425 -1.280355215
204540_at   1.203772435 9.295652818 3.406487382 0.002619724 0.135774909 -1.488344843
231391_at   1.198685321 5.519248145 1.592223094 0.126076553 0.383783955 -4.736587099
AFFX-BioB-3_at  1.182337133 9.389058914 1.621105002 0.119714959 0.376757921 -4.696859358
218952_at   1.173692029 8.701987814 2.904299409 0.008404878 0.183376674 -2.497512205
237053_at   1.159737875 7.195353692 4.556380728 0.00016674  0.091335433 0.909233566
234946_at   1.153359901 6.23010351  3.489931506 0.002150626 0.131163846 -1.316691485
225969_at   1.148201276 7.386474415 4.94689063  6.55E-05    0.091335433 1.717141363
234024_at   1.136922184 6.326235237 5.112017647 4.42E-05    0.091335433 2.054745916
220551_at   1.130914731 6.877391666 2.52829159  0.01942516  0.227296977 -3.213009976
ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by Paul80

I think my design matrix is not proper, it should be like

design.mat <- cbind(c(1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1))
colnames(design.mat) <- c("Case","Cont")

correct this if I am wrong

ADD REPLYlink written 2.1 years ago by Paul80
Please log in to add an answer.


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