Dear List,
   
  Sorry to bother you all with this.
  I am continuing my thread on Rat arrays (previous messages are below).
  I created MAplots (for each comparison in fit2_q, using plotMA(fit2_q) ) and a QC report for my arrays http://photo.all-faqs.info/index2.html
   
  I know that among you there are many Affy QC gurus! Would you say that arrays from the 4th treatment (arrays 11-15 in the QC report) should be discarded?
   
  Thank you for your help,
  Lev. 


Lev Soinov <lev_embl1@yahoo.co.uk> wrote:    Dear List,
   
  In addition to my previous e-mail (below):
  I am also getting very strange results such as: for coef=6 the control probe AFFX_Rat_GAPDH_5_at (number 31093 on the Rat230_2 array) has M=1.25, adj.P=3.02E-07 and B=10.44
  There are some other control probes showing the same  behaviour.
  I am lost in what might happen here.
   
  Looking forward to hearing your comments,
  Lev.
  P.S.
  My data:
  > data
AffyBatch object
size of arrays=834x834 features (157597 kb)
cdf=Rat230_2 (31099 affyids)
number of samples=29
number of genes=31099
annotation=rat2302
  

Lev Soinov <lev_embl1@yahoo.co.uk> wrote:
  Dear List,

I am analysing 25 rat affy arrays: 5 treatments, 5 rats each.
I am using the following code:

memory.limit(size = 2048)
data<-ReadAffy()

temp_q = rma(data)
design <- model.matrix(~0+factor(c(1,1,1,1,1,2,2,2,2,2,4,4,4,4,4,5,5,5,5,5,6,6,6,6,6)))
colnames(design) <- c("group1", "group2","group4","group5","group6")
contrast.matrix <- makeContrasts(group2-group1, group6-group5, group6-group4, group6-group2, group4-group2, group5-group1, levels=design)
fit_q <- lmFit(temp_q, design)
fit2_q <- contrasts.fit(fit_q, contrast.matrix)
fit2_q <- eBayes(fit2_q)

I am getting the following results:
> results<-decideTests(fit2_q, method="separate", lfc=0, p=0.05)
> summary(results)
group2 - group1 group6 - group5 group6 - group4 group6 - group2 group4 - group2 group5 - group1
-1 2366 102 10677 3836 7596 3728
0 26550 30848 11118 22437 15761 22113
1 2183 149 9304 4826 7742 5258

The number of differentially expressed genes is enormous and is significantly different than if I run lmFit separately for each comparison, for example:

design_small <- model.matrix(~0+factor(c(1,1,1,1,1,2,2,2,2,2)))
colnames(design_small) <- c("group1", "group2")
contrast.matrix_small <- makeContrasts(group2-group1, levels=design_small)
fit <- lmFit(exprs(temp_q)[,1:10], design_small)
fit2 <- contrasts.fit(fit, contrast.matrix_small)
fit2 <- eBayes(fit2)

> results<-decideTests(fit2)
> summary(results)
group2 - group1
-1 3518
0 24100
1 3481

topTable gives different rankings of genes for these two cases:
> topTable(fit2_q,coef=1, adjust="BH")
ID logFC AveExpr t P.Value adj.P.Val B
3743 1371194_at 3.951198 7.965214 20.72545 1.678042e-17 5.218543e-13 28.96959
1852 1369303_at 3.945072 8.300152 19.03215 1.304008e-16 1.478600e-12 27.20154
908 1368359_a_at 1.696257 11.202488 18.96096 1.426348e-16 1.478600e-12 27.12316
13780 1381233_at 1.501010 8.518978 17.75213 6.849205e-16 5.325085e-12 25.73860
5278 1372729_at 3.521604 6.794043 17.55966 8.867461e-16 5.515383e-12 25.50839
16984 1384448_at 1.824844 8.203814 16.57946 3.433390e-15 1.636902e-11 24.29159
1324 1368775_at 2.949045 6.712478 16.52965 3.684464e-15 1.636902e-11 24.22771
19739 1387671_at 1.685727 5.202906 16.00705 7.810413e-15 3.036201e-11 23.54500
728 1368179_at 0.950436 8.586395 15.68615 1.251490e-14 4.324455e-11 23.11423
21641 1389573_at 1.239271 8.752593 15.47619 1.711045e-14 5.321179e-11 22.82747

> topTable(fit2, coef=1, adjust="BH")
ID logFC t P.Value adj.P.Val B
908 1368359_a_at 1.696257 38.43022 2.425595e-14 7.543358e-10 20.90599
1324 1368775_at 2.949045 35.15163 7.313681e-14 1.137241e-09 20.25415
1852 1369303_at 3.945072 29.91043 5.369030e-13 5.565715e-09 18.94744
16984 1384448_at 1.824844 27.83420 1.301827e-12 1.012138e-08 18.31582
7124 1374575_at 1.519371 25.14516 4.532823e-12 2.819325e-08 17.37644
6614 1374065_at 1.666809 24.35419 6.705340e-12 3.475490e-08 17.07025
5278 1372729_at 3.521604 22.84314 1.467501e-11 6.317391e-08 16.44241
6430 1373881_at 1.568560 22.65308 1.625105e-11 6.317391e-08 16.35918
21641 1389573_at 1.239271 22.18231 2.099743e-11 7.255545e-08 16.14868
3743 1371194_at 3.951198 21.69581 2.751449e-11 8.166672e-08 15.92444

QCs show that the data is OK. Am I making some silly mistake?
Could you please help?
Lev.
P.S.
> sessionInfo()
R version 2.4.1 (2006-12-18) 
i386-pc-mingw32 
locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom.1252;LC_MONETARY=English_United Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
attached base packages:
[1] "tools" "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base" 
other attached packages:
rat2302cdf limma affy affyio Biobase 
"1.14.0" "2.9.8" "1.12.2" "1.2.0" "1.12.2" 


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

[[alternative HTML version deleted]]

_______________________________________________
Bioconductor mailing list
Bioconductor@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor

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


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

	[[alternative HTML version deleted]]

