[BioC] Problems in limma package

[Gordon Smyth]

At 10:19 PM 30/03/2004, Binita Dutta wrote:
>Dear All,
>
>I have done cDNA microarray (dye swap) experiments  (two repeats)  with 
>samples from wild type and mutant  mice(experiment is very similar to 
>Swirl experiment given in tutorials). I tried normalising data with 
>"limma" package with following commands:
>library(limma)
>Warning message:
>package limma was built under R version 1.8.1
>RG<-read.maimages(files=c("binita1.txt","binita2.txt","binita3.txt","binita4.txt"), 
>columns=list(Rf = "CH1_NBC_INT", Rb = "CH1_SPOT_BKGD", Gf = "CH2_NBC_INT", 
>Gb = "CH2_SPOT_BKGD", Name="GENE_DESCRIPTION",verbose=TRUE,sep=" 
>\t",  quote = "\"'", dec = "."))
>  RG$genes<-read.table("binitaSample2.txt", sep="\t",header=TRUE,  quote = 
> "\"'",fill=TRUE)
>layout=list(ngrid.r=2,ngrid.c=12,nspot.r=45,nspot.c=21)
>  RG$printer<-layout
>  RG<-backgroundCorrect(RG, method="minimum")
>  MA<-normalizeWithinArrays(RG,layout=RG$printer)
>plotMA(MA)
>As expected, I get MA plot
>
>However, when i try normalising  the data between the arrays with 
>following commands:
>MA<-normalizeBetweenArrays(MA)
>  plotMA(MA)
>The graph flips  (if you want i can send graphs)
>
>i.e X axis is reversed

This doesn't make any sense to me, I don't think there is any way that the 
X axis could reverse. You would have to provide some evidence that 
something is wrong.

>design<-c(-1,1,-1,1)
>
>fit<-lmFit(MA,design)
>fit<-eBayes(fit)
>  qqt(fit$t,df=fit$df.prior+fit$df.residual,pch=16,cex=0.1)
>top<-topTable(fit,number=22680,adjust="fdr")
>ord<-order(fit$lods,decreasing=TRUE)
>top30<-ord[1:30]
>plot(A,M,pch=16,cex=0.1)
>text(A[top30],M[top30],labels=MA$genes[top30,"SPOT.LABEL"],cex=0.8,col="blue")
>I get following errors:
>
>1) P.Values which i obtain is 1 or above 1, i tried adjusting P.Value with 
>"Holms" etc but get the same result. However, the same experiment when i 
>analyse thorugh other progams, the P.Values are very less (less than 
>0.001) for differentially expressed genes.

The obvious explanation is that you have assigned treatments incorrectly, 
e.g., your design matrix is wrong. Are you saying that *all* of your 
p.values are equal to 1?

Are you claiming that you have p.values above 1? As far as I know, that 
cannot occur.

>3) i have problems in subseting topTables also.
>subset<-subset(topTable,P.Value<0.01,select=MA$genes)
>Error in subset.default(topTable, P.Value < 0.01, select = MA$genes) :
>         Object "P.Value" not found

You are trying to subset a function rather than a data.frame! The *output* 
from topTable() would be a data.frame. In any case, have you considered 
using topTable with a smaller 'number' to do what you want?

>2)on MA plot i want to highlight the top 30 genens with "SPOT LABEL" which 
>is there on MA$genes, but the program picks up  some number which does not 
>corresonds to the SPOT.LABEL.

This is most likely because your SPOT.LABEL is a factor rather than a 
character vector. Try setting

   MA$genes$SPOT.LABEL <- as.character(MA$genes$SPOT.LABEL)

Gordon

>I have shown here, top 10 genes
>SPOT.LABEL CIU CLONE.ID     M         A     t P.Value     B
>4459        4459 MM6    25179  1.91  1.627766 11.18       1 -2.82
>14466      14466 MM6    14674  1.62  0.695737 10.35       1 -2.84
>17413      17413 MM6    17720  2.48 -1.036538 10.08       1 -2.85
>21140      21140 MM6    16316  1.57  0.000431  9.80       1 -2.85
>18070      18070 MM6     9920  1.61  1.888210  9.61       1 -2.86
>776          776 MM6    27323  1.72  1.873411  9.45       1 -2.86
>12014      12014 MM6    25084  1.74  1.226841  9.43       1 -2.87
>19423      19423 MM6    20962 -1.44  4.009106 -9.22       1 -2.87
>17258      17258 MM6    13529  2.31  0.289187  9.06       1 -2.88
>21565      21565 MM6    27496  1.60  1.794250  8.78       1 -2.89
>
>Help in this regards will be highly appreciated.
>Sincerely yours,
>Binita
>
>==============================
>
>Binita Dutta, PhD
>MicroArray Facility(MAF)
>UZ Gasthuisberg
>Onderwijs en Navorsing
>Herestraat 49
>3000 Leuven
>Belgium
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at stat.math.ethz.ch
>https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor