[BioC] Please advice

[Naomi Altman]

Dear Mali,
I have a bit more time now.  I think I answered your first question 
in my last email.  Limma does not directly do a "main effects" and 
interactions type of F-test, although you can can force the issue by 
using the F.p.value and several different sets of contrasts.

The second approach is not valid.  The t-test denominator will be 
incorrect because you did not account for batch and treatment.

Best wishes,
Naomi


Regarding

At 07:05 AM 3/7/2011, you wrote:
>Dear Naomi
>First I would like to apologize for writing you directly.
>My name is Mali Salmon, I am a researcher in Tel Hashomer hospital in Israel.
>I am currently analyzing microarray data using limma package, and I 
>would like to have some advice from a statistician who has an 
>experience with microarray analysis using Limma and bioconductor.
>I was searching the mailing list, and saw your name as one who 
>assist many people with their analysis.
>I have a simple factorial design which I build the design and 
>contrast matrix for, and I just want to be sure that what I did is 
>correct, I would be very grateful is you could give me your opinion.
>
>I have two factors, one is strain with three levels (C,D,S), and a 
>treatment factor with two levels: treated (DSR) and untreated (NO), 
>I also have batch effect which I included in the matrix
>
> > targets
>
>    FileName Treatment Strain batch
>1     NO2_C        NO      C     1
>2    DSR2_C       DSR      C     1
>3     NO2_D        NO      D     1
>4    DSR2_D       DSR      D     1
>5     NO2_S        NO      S     1
>6    DSR2_S       DSR      S     1
>7     NO3_C        NO      C     2
>8    DSR3_C       DSR      C     2
>9     NO3_D        NO      D     2
>10   DSR3_D       DSR      D     2
>11    NO3_S        NO      S     2
>12   DSR3_S       DSR      S     2
>
>I am trying to follow example 8.7 in the user guide (factorial design)
>In order to build the design matrix this is what I did:
>
>TS <- paste(targets$Treatment, targets$Strain, sep=".")
>TS <- factor(TS, unique(TS))
>design<-model.matrix(~0+TS)
>colnames(design) <- levels(TS)
>batch<-c(0,0,0,0,0,0,1,1,1,1,1,1)
>design<-cbind(design,batch)
>
> > design
>    NO.C DSR.C NO.D DSR.D NO.S DSR.S batch
>1     1     0    0     0    0     0     0
>2     0     1    0     0    0     0     0
>3     0     0    1     0    0     0     0
>4     0     0    0     1    0     0     0
>5     0     0    0     0    1     0     0
>6     0     0    0     0    0     1     0
>7     1     0    0     0    0     0     1
>8     0     1    0     0    0     0     1
>9     0     0    1     0    0     0     1
>10    0     0    0     1    0     0     1
>11    0     0    0     0    1     0     1
>12    0     0    0     0    0     1     1
>
>
>Now the questions I would like to answer are:
>1. what is the main effect of the strain (here I want to ignore the 
>treatment, just to find the difference between different strains)
>2. what is the main effect of the treatment (genes that change their 
>expression because of treatment, again I want to ignore the strain here).
>3. genes that are respond to treatment in each strain + interaction
>
>Below is how I created the contrast.matrix, the first three 
>contrasts are to answer the first question, contrast number 4 is to 
>answer question 2, and the last 6 contrasts to answer question 3
>
>cont.matrix <- makeContrasts(
>      DvsC=(DSR.D+NO.D)-(DSR.C+NO.C),
>
>      SvsC=(DSR.S+NO.S)-(DSR.C+NO.C),
>      DvsS=(DSR.D+NO.D)-(DSR.S+NO.S),
>      DSRvNO=(DSR.C+DSR.D+DSR.S)-(NO.C+NO.D+NO.S),
>      C.DSRvsNO=DSR.C-NO.C,
>      D.DSRvsNO=DSR.D-NO.D,
>      S.DSRvsNO=DSR.S-NO.S,
>      DiffDvsC=(DSR.D-NO.D)-(DSR.C-NO.C),
>      DiffSvsC=(DSR.S-NO.S)-(DSR.C-NO.C),
>      DiffSvsD=(DSR.S-NO.S)-(DSR.D-NO.D),
>      levels=design)
>fit2 <- contrasts.fit(fit, cont.matrix)
>fit2 <- eBayes(fit2)
>
>Is this matrix correct?
>The batch is included in the design matrix, is that mean that the 
>batch will be removed?
>
>One of the comparisons I am interested with is the different between 
>strains (first question). Would you suggest to do a separate 
>"several group" analysis (example 8.6 in limma user guide) instead 
>of the one indicated above?
>The analysis for answer the first question would be:
>
>design.Strain <- model.matrix(~ 0+factor(c(1,1,2,2,3,3,1,1,2,2,3,3)))
>colnames(design.Strain) <- c("C", "D", "S")
>#correct for batch effect
>batch<-c(0,0,0,0,0,0,1,1,1,1,1,1)
>design.Strain<-cbind(design.Strain,batch)
>fit.Strain <- lmFit(selDataMatrix, design.Strain)
>contrast.matrix.Strain <- makeContrasts(D-C, S-D, S-C, levels=design.Strain)
>fit2.Strain <- contrasts.fit(fit.Strain, contrast.matrix.Strain)
>fit2.Strain <- eBayes(fit2.Strain)
>
>I have tried the two approaches, and I got different list of 
>differential expressed genes, so which one is preferred?
>
>Looking forward to your reply, and apologising for bothering you
>Thank you
>Mali
>
>