Question: lesser fold change value using RMA background correction
2
gravatar for Paul
12 months ago by
Paul80
India
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)

 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?

ADD COMMENTlink modified 12 months ago • written 12 months ago by Paul80
1
d1<-ReadAffy()

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

ADD REPLYlink written 12 months ago by Kevin Blighe48k

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

ADD REPLYlink written 12 months ago by Paul80
2

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 12 months ago • written 12 months ago by Kevin Blighe48k

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 12 months ago • written 12 months ago by Paul80
1

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 12 months ago • written 12 months ago by Kevin Blighe48k

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 12 months ago • written 12 months 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")
design.mat

correct this if I am wrong

ADD REPLYlink written 12 months ago by Paul80
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: 887 users visited in the last hour