Hi everyone, I am stuck with creating design matrix in LIMMA for my microarray experimental design. The experiment consists of 40 Agilent arrays. The subjects were treated with two stress conditions (S1 and S2). In each stress condition, the subjects were sampled in 5 timepoints (0, 6, 12, 24 and 48). Microarray experiment was performed in such a way that every fourth array of a particular timepoint, irrespective of stress treatment, is dye swapped. Basically, it is a time course experiment on the same subject with two different stresses. I want to to know which are the genes important/responsible for stress 1 or 2. The target file looks like below: 

The fact that two different pools were used means that you almost have two independent experiments. For a study like this, the standard limma method is to make two design matrices and combine them together:

design.S1 <- modelMatrix(targets[ 1:20,], ref="pool1")
design.S2 <- modelMatrix(targets[21:40,], ref="pool2")
design <- blockDiag(design.S1, design.S2)

I would also add an intercept term to allow for probe-specific dye effects, which tend to be important for Agilent arrays:

design <- cbind(Dye=1, design)

This is equivalent to Ryan's code, but a bit simpler.

Since your two stressors apparently have different reference samples (pool1 and pool2), I think the most straightforward way to implement this is by making a separate design matrix for each stressor and then combining them into a single block-diagonal matrix using Matrix::bdiag (correction: limma already provides a blockDiag function for exactly this purpose):

library(limma)
library(stringr)

# Indicate the reference sample for each row
targets$Ref <- str_replace(targets$Stress, "S", "pool")

# Split by reference and make the design matrix for each split
targets.split <- split(targets, targets$Ref)
designs <- lapply(targets.split, function(tgt) {
    modelMatrix(targets=droplevels(tgt), ref=tgt$Ref[1])
})

# Combine the separate design matrices into one
design <- as.matrix(do.call(blockDiag, designs))
## Ensure that the design has the rows in the same order as targets
design <- design[rownames(targets),]

## Define contrasts for an ANOVA test of all time points for each
## stressor
stressor1_contrasts <- sprintf("S1_Time%i - S1_Time0", c(6, 12, 24, 48))
stressor2_contrasts <- sprintf("S2_Time%i - S2_Time0", c(6, 12, 24, 48))
## Build the contrast matrix for each ANOVA test
stressor1_cmat <- makeContrasts(contrasts=stressor1_contrasts, levels=design)
stressor2_cmat <- makeContrasts(contrasts=stressor2_contrasts, levels=design)

From there, you can use the usual limma procedure of lmFit, contrasts.fit, eBayes, and topTable to get a table of genes for each ANOVA test.

If there are any additional variables to be taken into account, such as a subject effect, you'll need to reply with more information.

Woops -- provided a boneheaded answer in haste. Had a momentary lapse of reason regarding the nuts and bolts of two-color arrays.

The details of my original answer have been redacted to mitigate any damage this might cause a wary future traveler ...

Thanks to Ryan for the correction