Dear Gordon and List,
  
I have now found a way to show my plots regarding barley arrays, 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


       
---------------------------------

	[[alternative HTML version deleted]]

