This is a test version of Biostars. For the public version, visit https://www.biostars.org.
lesser fold change value using RMA background correction

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)

 library("limma")
    library("affy")
    library("gcrma")

    setwd('E://GSE...7/')
    d1<-ReadAffy()

    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)
    design.mat

    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?

microarray deg limma rma normalization
• 3,108 views • link
written by Paul Paul 
d1<-ReadAffy()

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

• 1 views • link
written by Kevin Blighe Kevin Blighe

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

• 1 views • link
written by Paul Paul

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)

• 1 views • link
written by Kevin Blighe Kevin Blighe

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?

• 1 views • link
written by Paul Paul

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.

• 1 views • link
written by Kevin Blighe Kevin Blighe

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
• 1 views • link
written by Paul Paul

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")
design.mat

correct this if I am wrong

• 1 views • link
written by Paul Paul

0 answers

No answers yet.

Log in to answer this question.

More posts like this