Question: Why my limma DEG result in top Table all without significant meaning,showing logFC <1.5 or logFC>-1.5 ?
0
gravatar for BioLite
7 months ago by
BioLite20
BioLite20 wrote:

Hi, all I'm a newbie, and I'm dealing with Agilent single channel microarray when I filtering DEGs using limma, I find all my results showing logFC>-1.5 or logFC<1.5, without significant expression. I have been searching reason for a long time, but I still don't know what my faults are, I really need some help, thanks! My code as following:

treatment<-targets$Treatment
Level<-c("neurologically healthy control","Parkinson disease")
Group<-factor(treatment,levels = Level)
Design<-model.matrix(~0+Group)
colnames(Design) <- c("Control","PD")
fita<-lmFit(exprs(eSet),Design)
cont.matrix <- makeContrasts(Diff=PD-Control, levels=Design)
fitb<- contrasts.fit(fita, cont.matrix)
fitc<- eBayes(fitb)

##result
Result<-topTable(fitc,adjust="BH",number = Inf,p.value = 0.05)
dt_result<-decideTests(fitc,lfc =1.5,p.value=0.05,method="separate")
summary(dt_result)range

(Result$logFC)
[1] -1.471461  1.192352
toptable deg limma ;logfc • 493 views
ADD COMMENTlink modified 7 months ago • written 7 months ago by BioLite20
1

One issue is here:

dt_result<-decideTests(fitc,lfc = log2(1.5),p.value=0.05,method="separate")

The parameter, lfc, expects a value that is already on the log (base 2) scale. So, the value that you are passing as the cut-off is actually:

log2(1.5) = 0.5849625

You need to use:

dt_result <- decideTests(fitc, lfc = 1.5, p.value=0.05, method="separate")

Kevin

ADD REPLYlink modified 7 months ago • written 7 months ago by Kevin Blighe41k

Yes Kevin, thank for your answer. But my key question still alive. Once I set parameter lfc=1.5, summary(dt_result) showing Diff Down 0 NotSig 38822 Up 0 That's because of my topTable results with no changes, with nothing significant expression. I really don't know how to get rid of this situation. Could you help me, Kevin? Very thanks!

ADD REPLYlink written 7 months ago by BioLite20

You are in good hands, now that Devon has added a comment. Please follow up with both him and Fabio (Fabio is another great guy)

ADD REPLYlink written 7 months ago by Kevin Blighe41k

Kevin, thank you very much! Devon and Fabio both experienced than me, I will learn modestly from them. Thanks again, wish you have a good day!

ADD REPLYlink modified 7 months ago • written 7 months ago by BioLite20

Here is my complete code, really hope to get some help.

rm(list = ls())
setwd("E:/Practice/Aglient/PD")
library(Agi4x44PreProcess)
library(hgug4112a.db)
library(limma)

##loading data
targets<-read.targets("targets.txt")
raw_data<-read.AgilentFE(targets = targets,makePLOT = FALSE)

##Background Correction and Normalization between arrays
data_NORM<-BGandNorm(raw_data,
                     BGmethod="half",
                     NORMmethod="quantile",
                     foreground="MeanSignal",
                     background="BGMedianSignal",
                     offset=50,makePLOTpre=FALSE,
                     makePLOTpost=FALSE)

##filtering probe
ddFILT<-filter.probes(data_NORM,
                      control=TRUE,
                      wellaboveBG=TRUE,
                      isfound=TRUE,
                      wellaboveNEG=TRUE,
                      sat=TRUE,
                      PopnOL=TRUE,
                      NonUnifOL=TRUE,
                      nas=TRUE,
                      limWellAbove=75,
                      limISF=75,
                      limNEG=75,
                      limSAT=75,
                      limPopnOL=75,
                      limNonUnifOL=75,
                      limNAS=100,
                      makePLOT=TRUE,
                      annotation.package="hgug4112a.db",
                      flag.counts=TRUE,targets)

##summarizing
data_filter<-summarize.probe(ddFILT, makePLOT=FALSE, targets)
dim(data_filter)

##create an expressionSet
eSet<-build.eset(data_filter, targets, makePLOT=FALSE,annotation.package="hgug4112a.db")

##DEG filter
treatment<-targets$Treatment
Level<-c("neurologically healthy control","Parkinson disease")
Group<-factor(treatment,levels = Level)
Design<-model.matrix(~0+Group)
colnames(Design) <- c("Control","PD")
fita<-lmFit(exprs(eSet),Design)
cont.matrix <- makeContrasts(Diff=PD-Control, levels=Design)
fitb<- contrasts.fit(fita, cont.matrix)
fitc<- eBayes(fitb)

##DEG result
Result<-topTable(fitc,adjust="BH",number = Inf,p.value = 0.05,sort.by = "p")
dt_result<-decideTests(fitc,lfc = 1.5,p.value=0.05,method="separate")
summary(dt_result)
ADD REPLYlink modified 7 months ago by genomax65k • written 7 months ago by BioLite20

Does Result indicate you have any DE genes?

ADD REPLYlink written 7 months ago by Devon Ryan89k

Nothing. I just checked my code again, finding some wrongs. So sad, I always make mistakes. Let me correct it and then discuss with you.

Devon, thanks for your caring! Have a good day!

BTW, my data from this article,https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002794. I am trying to replicate the Agilent part results (the first part) of this article. As a newbie, It's apparently I failed completely.

Best wish for you again.

ADD REPLYlink modified 7 months ago • written 7 months ago by BioLite20

If there really are DE genes in that dataset (from the paper I guess one would expect so) then it's likely that there's either an outlier sample that needs excluding or some of the sample labels are swapped. Make a PCA plot and do some clustering and see if that indicates either of those.

ADD REPLYlink written 7 months ago by Devon Ryan89k

So sad, I always make mistakes.

Do not worry - we all make mistakes. At least you are honest about it, whereas others would prefer to lie about their mistakes. You are a great person.

ADD REPLYlink written 7 months ago by Kevin Blighe41k

How many replicates do you have per condition?

Do you see very high variation between replicates?

Did you try doing the analysis with a toy data set (probably provided by the package developer) to see if using that you got significant results?

Are there additional software packages you may try to see if the result is due to some problem in the analysis with this specific package?

ADD REPLYlink modified 7 months ago • written 7 months ago by Fabio Marroni2.1k

Yes, Fabio. Q1, 27 PD, and 26 control samples;

Q2, my analysis data from this published article,https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002794. The bad samples have been removed by the authors, I get the cleaned data from ArrayExpress,https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-812/

Q3, not yet.

Q4, I made this analysis using limma at first and got the same results.

I reviewed my practice article again, thinking there may be some faults in my limma lmFit object and model design. I will correct them firstly.

Thanks for your precious time. Fabio, wish you have a great day!

ADD REPLYlink written 7 months ago by BioLite20

I added code markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:

101010 Button

ADD REPLYlink written 7 months ago by WouterDeCoster38k

In addition, cross-posted: https://support.bioconductor.org/p/112875/

ADD REPLYlink written 7 months ago by Kevin Blighe41k

Yes, thanks for your sweet tip, I will notice my code layout next time.

ADD REPLYlink written 7 months ago by BioLite20
1
gravatar for Gordon Smyth
7 months ago by
Gordon Smyth750
Australia
Gordon Smyth750 wrote:

The Agi4x44PreProcess package was removed from Bioconductor many years ago. You must be using very old versions of R and Bioconductor. I would certainly not advise you to preprocess the Agilent microarray data in the complicated way that you have. Why not use more up-to-date software and a pipeline as recommended in the limma User's Guide?

I analysed this data myself and had no difficulty finding DE genes. See my answer on the Bioconductor Support site https://support.bioconductor.org/p/112875/

ADD COMMENTlink modified 7 months ago • written 7 months ago by Gordon Smyth750

Dear Gordon,

Actually, I worked with limma at first. However, I couldn't understand the probe filter part in the user guide.

Following your tips, I will analysis with limma again.

Thanks!

ADD REPLYlink modified 7 months ago • written 7 months ago by BioLite20
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: 1337 users visited in the last hour