Hi,
I am having problems fitting a model (paired design with multiple factors) using limma, I have tried different ways of specifying the design and contrast matrix, but I stuck in the middle somewhere, so I post them here, maybe someone have already worked it out and know how to properly fit this:

I have data from 10 animals, sampled at two different time points (D0 and D2), animals have different phenotypes (LS and PS).
I want to do paired tests for

Differential expression between D0-D2 in LS
Differential expression between D0-D2 in PS
Differential expression between LS and PS

See detailed information below, I will appreciate if someone can give some suggestion to what I have did.

kind regards,
Tinghua


Here is my targets frame:

AffyCell   Pig           Type        Day
C01          P01          LS            D0
C02          P01          LS            D2
C03          P02          LS            D0
C04          P02          LS            D2
C05          P03          LS            D0
C06          P03          LS            D2
C07          P04          LS            D0
C08          P04          LS            D2
C09          P05          PS            D0
C10          P05          PS            D2
C11          P06          PS            D0
C12          P06          PS            D2
C13          P07          PS            D0
C14          P07          PS            D2
C15          P08          PS            D0
C16          P08          PS            D2
C17          P09          PS            D0
C18          P09          PS            D2
C19          P010        PS            D0
C20          P010        PS            D2

I want to know:

1. Which genes respond to stimulation in LS,
2. Which genes respond to stimulation in PS, and
3. Which genes respond differently in LS compared to PS.

I have tested three kinds of design: paired, factorial, factorial and paired.
The "factorial" works fine, but this is not what I want. My experiment is paired experiment, I mean the Affychips, for example C01 and C02 are from the sample animal of different time point.

Basically I stuck here when I trying to fit the paired model:

design <- model.matrix(~pig+type)
# colnames(design) <- c(levels(pig), levels(type.day))
fit <- lmFit(eset, design)

This line "fit <- lmFit(eset, design)" throw out an error like this:

> fit <- lmFit(eset, design)
Coefficients not estimable: typePS
Warning message:
Partial NA coefficients for 24123 probe(s)

> fit1 <- eBayes(fit)
Warning message:
In ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim) :
  Estimation of var.prior failed - set to default value

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

The following is the coding

library(affy)
library("arrayQualityMetrics")
library(limma)

# read affy data

cel.dir = "Cells"
exp.dir = "exp.expression"
qc.dir = "Qc"

eset <- justRMA(celfile.path=cel.dir)
exp <- exprs(eset)
write.table(exp, file = exp.dir, sep = "\t")

arrayQualityMetrics(expressionset=eset, outdir=qc.dir, force = TRUE)

# limma

targets <- read.csv("targets.csv")
targets

pig <- factor(targets$Pig)
pig

type <- factor(targets$Type)
type

day <- factor(targets$Day)
day

type.day <- factor(paste(targets$Type, targets$Day, sep="."))
type.day

# paired

design <- model.matrix(~pig+type)
# colnames(design) <- c(levels(pig), levels(type.day))
fit <- lmFit(eset, design)

fit1 <- eBayes(fit)
topTable(fit, coef="typePS")

results <- decideTests(fit1)
vennDiagram(results)

# factorial

design <- model.matrix(~0+type.day)
# colnames(design) <- c(levels(type.day))
fit <- lmFit(eset, design)


cont.matrix <- makeContrasts(
                LS.D2vD0=LS.D2-LS.D0,
                PS.D2vD0=PS.D2-PS.D0,
                Diff.PSvLS=(PS.D2-PS.D0)-(LS.D2-LS.D0),
                levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

results <- decideTests(fit2)
vennDiagram(results)


# factorial and paired

design <- model.matrix(~pig+type.day)
# colnames(design) <- c(levels(pig), levels(type.day))
fit <- lmFit(eset, design)


cont.matrix <- makeContrasts(
                LS.D2vD0=LS.D2-LS.D0,
                PS.D2vD0=PS.D2-PS.D0,
                Diff.PSvLS=(PS.D2-PS.D0)-(LS.D2-LS.D0),
                levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

results <- decideTests(fit2)
vennDiagram(results)


	[[alternative HTML version deleted]]

