[BioC] Very low P-values in limma

[Mayer, Claus-Dieter]

Dear Paul!

What is the corfit$consensus value? If it is close to 1, using limma with duplicateCorrelation or averaging (or indeed just using any of the duplicates) should give you very similar results.
If however corfit$consensus is close to 0, limma would indeed treat the duplicates similar to independent replicates. Obviously a corfit$consensus value close to 0 would be quite worrying as that would indicate a very high technical variability (or some subtle mistake in your code). Anyway, checking that value might help to solve the puzzle.

Claus

> -----Original Message-----
> From: bioconductor-bounces at stat.math.ethz.ch [mailto:bioconductor-
> bounces at stat.math.ethz.ch] On Behalf Of Paul Geeleher
> Sent: 23 October 2009 15:13
> To: Gordon K Smyth
> Cc: Bioconductor mailing list
> Subject: Re: [BioC] Very low P-values in limma
>
> Dear list,
>
> I have a query, I've rerun my analysis by averaging out the within
> array duplicates using the avedups() function instead of
> duplicateCorrelation().
>
> aveMat <- avedups(mat at hx, ndups=4, spacing=1, weights=NULL)
> fit <- lmFit(aveMat, design)
> ebayes <- eBayes(fit)
> topTable(ebayes, coef = 2, adjust = "BH", n = 10)
>
> When I do this my most significant p-values drop from 10^-7 to .06!
>
> This seems dubious? It seems like to get values as low as 10^-7 the
> duplicateCorrelation() function must be treating the within array
> replicates in a similar way it would separate samples, this seems
> suspicious to me...
>
> -Paul.
>
>
>
> On Fri, Oct 23, 2009 at 3:45 AM, Gordon K Smyth <smyth at wehi.edu.au> wrote:
> > Dear Paul,
> >
> > I'm not quite sure why you're so shocked about getting significant
> results.
> >  Usually people complain to me when they don't get significance :).
> >
> > limma is able to obtain much higher significance levels than a t-test
> > because it (i) leverages information from the within-array replicates
> and
> > (ii) borrows information across probes.  These two approaches in
> combination
> > increase the effective degrees of freedom dramatically.
> >
> > We would hardly go to so much time and trouble to develop new
> statistical
> > methods if they didn't improve on ordinary t-tests.
> >
> > If I were you, I'd be looking at the sizes of the fold changes, and
> whether
> > the results seem biologically sensible and agree with known information,
> and
> > so on.
> >
> > Best wishes
> > Gordon
> >
> >
> > ---------- original message -----------
> > [BioC] Very low P-values in limma
> > Paul Geeleher paulgeeleher at gmail.com
> > Thu Oct 22 11:34:01 CEST 2009
> >
> > The value of df.prior is 3.208318. Hopefully somebody can help me here
> > because I'm still really at a loss on why these values are so low.
> >
> > Here's the script. I'm fairly sure it conforms to the documentation:
> >
> > library(affy)
> >
> > setwd('/home/paul/PhD_Stuff/miRNA_guilleme_project/HERArrays/')
> >
> > # create a color pallete for boxplots
> > cols <- brewer.pal(12, "Set3")
> >
> > # This defines the column name of the mean Cy3 foreground intensites
> > #Cy3 <- "F532 Median"
> > # background corrected value
> > Cy3 <- "F532 Median - B532"
> >
> > # This defines the column name of the mean Cy3 background intensites
> > Cy3b <- "B532 Mean"
> >
> > # Read the targets file (see limmaUsersGuide for info on how to create
> this)
> > targets <- readTargets("targets.csv")
> >
> > # read values from gpr file
> > # because there is no red channel, read Rb & R to be the same values as
> G
> > and Gb
> > # this is a type of hack that allows the function to work
> > RG <- read.maimages( targets$FileName,
> > source="genepix",columns=list(R=Cy3,G=Cy3, Rb=Cy3b, Gb=Cy3b))
> >
> > # remove the extraneous values red channel values
> > RG$R <- NULL
> > RG$Rb <- NULL
> >
> > # this line of code will remove any negative values which cause errors
> > further on
> > RG$G[RG$G<0] <- 0
> >
> > # create a pData for the data just read
> > # this indicates which population each array belongs to
> > # a are "her-" and b are "her+" (almost certain)
> > pData <- data.frame(population = c('a', 'a', 'a', 'a', 'b', 'b', 'b'))
> > rownames(pData) <-  RG$targets$FileName
> >
> > # create design matrix
> > design <- model.matrix(~factor(pData$population))
> >
> > # In my .gpr files all miRNAs contain the string "-miR-" or "-let-" in
> their
> > name
> > # so the grep function can be used to extract all of these, removing
> > # all control signals and printing buffers etc.
> > # You need to check your .gpr files to find which signals you should
> > extract.
> > miRs <- c(grep("-miR-", RG$genes$Name), grep("-let-", RG$genes$Name))
> > RG.final <- RG[miRs, ]
> >
> > # load vsn library, this contains the normalization functions
> > library('vsn')
> >
> > # Do VSN normalization and output as vns object
> > mat <- vsnMatrix(RG.final$G)
> >
> > # insert rownames into matrix with normalized data
> > # this will mean that the gene names will appear in our final output
> > rownames(mat at hx) <- RG.final$genes$Name
> >
> > # my .gpr files contain 4 "within-array replicates" of each miRNA.
> > # We need to make full use of this information by calculating the
> duplicate
> > correlation
> > # in order to use the duplicateCorrelation() function below,
> > # we need to make sure that the replicates of each gene appear in
> sequence
> > in this matrix,
> > # so we sort the normalized values so the replicate groups of 4 miRs
> appear
> > in sequence
> > mat at hx <- mat at hx[order(rownames(mat at hx)), ]
> >
> > # calculate duplicate correlation between the 4 replicates on each array
> > corfit <- duplicateCorrelation(mat, design, ndups=4)
> >
> > # fit the linear model, including info on duplicates and correlation
> > fit <- lmFit(mat, design, ndups=4, correlation=corfit$consensus)
> >
> > # calculate values using ebayes
> > ebayes <- eBayes(fit)
> >
> > # output a list of top differnetially expressed genes...
> > topTable(ebayes, coef = 2, adjust = "BH", n = 10)
> >
> >
> > -Paul.
> >
> >
> > On Thu, Oct 22, 2009 at 12:09 AM, Wei Shi <shi at wehi.edu.au> wrote:
> >
> >> Dear Paul:
> >>
> >>  The low p-values you have got are not surprising to me. I have got
> >
> > even
> >>
> >> lower FDR adjusted p-values than yours. This just means you got
> >> significantly differentially expressed genes. But on the other hand, I
> >
> > did
> >>
> >> see much higher adjusted p-values in some of my microarray analyses. In
> >
> > that
> >>
> >> case, I will explore the data in more depth such as looking at batch
> >
> > effect
> >>
> >> etc.
> >>
> >> Cheers,
> >> Wei
> >>
> >>
> >> Paul Geeleher wrote:
> >>
> >>> Hi folks, I'm analyzing microRNA data using limma and I'm wondering
> >
> > about
> >>>
> >>> the validity of the p-values I'm getting out. Its a simple 'Group A Vs
> >>> Group
> >>> B' experimental design. 4 arrays in one group, 3 in the other and 4
> >>> duplicate spots for each miRNA on each array.
> >>>
> >>> The lowest adjusted p-values in the differential expression analysis
> >
> > are
> >>>
> >>> in
> >>> the region of 10^-7.
> >>>
> >>> Its been pointed out to me that plugging the point values from each
> >
> > sample
> >>>
> >>> into a regular t-test you get p=0.008, which then also needs to take
> >
> > the
> >>>
> >>> multiple test hit. Can anybody explain why limma is giving me such
> >
> > lower
> >>>
> >>> values and if they are valid?
> >>>
> >>> I can provide more information if required.
> >>>
> >>> Thanks,
> >>>
> >>> Paul.
> >>>
> >>>
> >>>
> >>
> >
> >
> > --
> > Paul Geeleher
> > School of Mathematics, Statistics and Applied Mathematics
> > National University of Ireland
> > Galway
> > Ireland
> >
> >
>
>
>
> --
> Paul Geeleher
> School of Mathematics, Statistics and Applied Mathematics
> National University of Ireland
> Galway
> Ireland
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor


The University of Aberdeen is a charity registered in Scotland, No SC013683.