[BioC] paired design, LIMMA
Gordon Smyth
smyth at wehi.EDU.AU
Mon Apr 23 09:37:10 CEST 2007
>[BioC] paired design, LIMMA
>Lev Soinov lev_embl1 at yahoo.co.uk
>Fri Apr 20 12:37:40 CEST 2007
>
>Dear List,
>
> I am learning about a simple paired design in LIMMA and am
> playing with a small dataset of 6 Affymetrix barley chips, 3
> treated and 3 untreated. I have some problems with interpreting the
> results and would be grateful for any comments/suggestions.
> Experiment: sample pairs (treated & untreated) were prepared in
> three biological replicates, using the same protocol (same
> treatments, etc.) but separately from each other (in different times).
> For all genes with negative fold changes, adj. p values for
> moderated t statistics appear to be higher than 0.1 (the smallest
> adj. p value among down-regulated genes is 0.139). Besides, only
> two down-regulated probes have abs(logFC)>log2(1.5).
>
> Questions:
> 1. From your experience, is the fact that among significantly
> regulated genes I only get up-regulated ones an indication of a
> problem with the data (log_intensity plots and boxplots did not
> flag up any significant problems)?
Well, if this is biologically infeasible, then it would seem to
indicate a problem.
> 2. With moderated t-statistics I am getting no significantly
> down-regulated genes but with ordinary t-statistics I get more than
> 4000 down-regulated probes with adj. p <0.05. Is this common?
Ordinary t-tests typically throw up a lot of spuriously DE genes,
which have very small standard deviations, low fold changes and low
expression levels.
The difference here between moderated and ordinary t-test suggests to
me that all the apparently down-regulated probes are in the lower
expression range. This does suggest to me a problem with the data. A
fitted model MA-plot might throw some light on the situation.
Best wishes
Gordon
> I also wonder why the difference between adj. p values for
> moderated and ordinary t statistics is so huge, i.e. moderated adj.
> p values for down-regulated genes are all higher than 0.1, while
> ~4000 down-regulated probes have ordinary adj.p<0.05.
>
> With kind regards,
> Lev.
> Bioinformatician, UK.
>
> I am using the following code (as described in the LIMMA user
> guide, p.40, 8.3 Paired Samples):
> memory.limit(size = 2048)
> data<-ReadAffy(widget=TRUE)
> sampleNames(data)
> temp<-rma(data)
> targets <- readTargets("Targets.txt")
> Pair <- factor(targets$Pair)
> Treat <- factor(targets$Treatment, levels=c("A","B"))
> design <- model.matrix(~Pair+Treat)
> fit_pair <- lmFit(temp, design)
> fit_pair <- eBayes(fit_pair)
> topTable(fit_pair, coef="TreatB")
>
> My targets file is:
> FileName Pair Treatment
> Bar18 1 A
> Bar19 1 B
> Bar20 2 A
> Bar21 2 B
> Bar22 3 A
> Bar23 3 B
>
>
>Ordinary t statistics for paired test were calculated using:
>
> >tstat.ord <- fit_pair$coef / fit_pair$stdev.unscaled / fit_pair$sigma
>
> >p.value.ord <- 2 * pt( abs(tstat.ord), df=fit_pair$df.residual,
> lower.tail=FALSE)
>
> >pvalue.ord.adj <- p.adjust(p.value.ord, method="fdr")
>
> My data and session info:
> >data
> AffyBatch object
> size of arrays=712x712 features (23770 kb)
> cdf=Barley1 (22840 affyids)
> number of samples=6
> number of genes=22840
> annotation=barley1
> > sessionInfo()
> R version 2.4.1 (2006-12-18)
> i386-pc-mingw32
> attached base packages:
> [1] "tcltk" "tools" "stats" "graphics" "grDevices"
> "utils" "datasets" "methods" "base"
> other attached packages:
> barley1cdf tkWidgets DynDoc
> widgetTools limma affy affyio Biobase
> "1.14.0" "1.12.1" "1.12.0" "1.10.0" "2.9.8"
> "1.12.2" "1.2.0" "1.12.2"
>
> There are NO significantly down-regulated genes, see the TreatB
> column below:
> > results<-decideTests(fit_pair)
> > summary(results)
> (Intercept) Pair2 Pair3 TreatB
> -1 0 407 161 0
> 0 4 21977 22453 22819
> 1 22836 456 226 21
>
> topTable (probe IDs are removed for brevity):
> ID logFC AveExpr t P.Value
> adj.P.Val B
> 1 1.4 5.9 15.2 0.0000002 0.003 5.3
> 2 6.1 6.7 14.6 0.0000003 0.003 5.1
> 3 3.1 6.4 14.3 0.0000003 0.003 5.0
> 4 1.2 7.6 12.6 0.0000010 0.005 4.6
> 5 3.0 7.7 12.4 0.0000011 0.005 4.5
> 6 1.5 4.6 11.7 0.0000017 0.007 4.3
> 7 1.9 6.7 11.1 0.0000026 0.007 4.0
> 8 1.0 4.6 11.1 0.0000026 0.007 4.0
> 9 3.3 5.2 11.0 0.0000029 0.007 4.0
> 10 1.8 7.5 10.1 0.0000055 0.013 3.6
> 11 1.5 6.4 9.5 0.0000089 0.018 3.3
> 12 1.1 8.2 9.4 0.0000097 0.018 3.2
> 13 3.5 6.4 8.4 0.0000223 0.039 2.7
> 14 1.0 4.4 8.3 0.0000256 0.040 2.6
> 15 1.1 6.4 8.3 0.0000260 0.040 2.6
> 16 1.8 4.6 8.1 0.0000290 0.041 2.5
> 17 0.9 9.4 8.0 0.0000345 0.044 2.4
> 18 0.9 6.9 7.9 0.0000349 0.044 2.3
> 19 4.5 7.2 7.9 0.0000372 0.045 2.3
> 20 1.0 7.4 7.8 0.0000397 0.045 2.3
> 21 0.8 7.0 7.7 0.0000437 0.048 2.2
> 22 2.1 3.4 7.6 0.0000484 0.050 2.1
> 23 0.9 6.5 7.2 0.0000701 0.070 1.8
> 24 2.6 6.0 7.2 0.0000734 0.070 1.8
> 25 0.6 9.0 7.1 0.0000781 0.071 1.7
> 26 1.2 5.1 7.1 0.0000808 0.071 1.7
> 27 1.3 6.5 7.0 0.0000867 0.072 1.7
> 28 1.5 4.8 7.0 0.0000891 0.072 1.6
> 29 3.8 5.2 6.9 0.0000945 0.072 1.6
> 30 2.3 6.0 6.9 0.0000959 0.072 1.6
More information about the Bioconductor
mailing list