[BioC] lmfit zero genes fit model

[Urska Cvek]

Hello BioC,

I am using the limma package and the lmfit function, as described in
section 8.7 here:
http://bioconductor.org/packages/1.9/bioc/vignettes/limma/inst/doc/usersguide.pdf.
I have used this same package before, on a different data set, and it
worked nice, so I am thinking there's something up with my data that I
don't understand completely.

I have 8 Affy chips, for a factorial design with Deletion and Strain
as the two factors:
Filename    Deletion    Strain
Ht-2002.cel    h            t2
Ht-2003.cel    h            t2
Hn2-2002.cel  h           n2
Hn2-2003.cel  h           n2
Lt-2002.cel     l            t2
Lt-2003.cel     l            t2
Ln2-2002.cel   l           n2
Ln2-2003.cel   l           n2

readAffy function is used to read in the AffyBatch object "a" using
the above PhenoData.
Since these experiments were done as replicate pairs, I first
normalize the data by pairs.

normalize(a[,1:2]), and repeate for the other three

normalized.a = merge(tmp1,tmp2)
normalized.a = merge(normalized.a, tmp3)
normalized.a = merge(normalized.a, tmp4)

x=rma(normalized.a, normalize=FALSE)

I look at the boxplots for the raw intensities, and expression set "x"
intensities, and the second set looks much better.

TS <- paste(pd$Deletion, pd$Strain, sep=".")    # extract all
combinations in one vector
TS

# fit a model with a coefficient for each of the four factor combinations
# and then extract the comparisons of interest as contrasts
TS <- factor(TS, levels=c("h.t2", "h.n2", "l.t2", "l.n2"))
design <- model.matrix(~0+TS)

> design
  h.t2 h.n2 low.t2 low.n2
1         1       0        0      0
2         1       0        0      0
3         0       1        0      0
4         0       1        0      0
5         0       0        1      0
6         0       0        1      0
7         0       0        0      1
8         0       0        0      1

colnames(design) <- levels(TS)
fit <- lmFit(x, design)

cont.matrix <- makeContrasts(
  H.n2vst2=h.t2-h.n2, L.n2vst2=l.t2-l.n2,
  Diff=(h.n2-h.t2)-(l.n2-l.t2),
  levels=design)

> cont.matrix
            Contrasts
Levels      H.n2vsi2 L.n2vsi2 Diff
  h.t2              1        0       -1
  h.n2             -1       0        1
  low.t2           0        1        1
  low.n2          0       -1       -1

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

# display the counts
results <- decideTests(fit2)
vennDiagram(results)

This vennDiagram is empty (has zeros for counts). There are no genes
that would be differentially expressed in other of the two linear
models, which I find very surprising. I know that it's possible that
not many genes would be differentially expressed, but I don't
understand why there are no genes that would fit the model at all.
I have tried just using the rma on the whole set of 8 chips as well,
without the separate normalization step, but this does not yield much
luck either. I also added the fourth replicate to the experiment,
making it a total of 12 chips, but that chip had such different
distribution that I removed it from the analysis.

What should I do to be able to answer the above questions? Is the
caused by my data and I cannot find such genes?

Thank you in advance, your help is appreciated.

U.C.