[BioC] paired design, LIMMA

[Gordon Smyth]

Dear Lev,

Just one final comment from me. You original question related to the 
difference between ordinary and moderated t-tests. The MA-plot that 
you give here would suggest that the the ordinary t-statistic result 
of 4000 down-regulated genes is completely spurious. The lack of any 
sizeable down fold changes would seem to show that the moderated 
t-test is correct for your data. Do you not agree?

Best wishes
Gordon

At 01:11 AM 4/05/2007, Lev Soinov wrote:
>Dear Gordon and List,
>
>I have now found a way to show my plots regarding barley arrays, 
><http://photo.all-faqs.info/>http://photo.all-faqs.info/.
>They are:
>- plotMA created with plotMA(fit_pair, array=4)
>- boxplots and a density plot for RMA normalised and raw data
>- MVA plots for RMA normalised data, created with 
>mva.pairs(exprs(temp), log.it=FALSE)
>
>Could you, please have a look at them and see whether you find 
>anything strange?
>
>My scripts for the analysis of the data and previous e-mails are below.
>
>With kind regards,
>Lev.
>
>
>
> > Thank you,
> > Lev.
> >
> > Gordon Smyth wrote:
> >
>[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