Hi Everybody,

I have a question regarding limma design matrix. I searched for "limma and design" but could not find an answer for my question. I will be happy if you could send me the link if a similar question has been answered previously!

I am re-analyzing a data set from GEO (GSE74615) to compare with some published microarray data and an RNA-Seq that we have.

I generated the target file as bellow:

FileName    Cy3    Cy5
US12302316_252665511011_S01_GE2-v5_95_Feb07_1_1.txt    microglia_WT    microglia_AD
US12302316_252665511011_S01_GE2-v5_95_Feb07_1_2.txt    microglia_AD    microglia_WT
US12302316_252665511011_S01_GE2-v5_95_Feb07_1_3.txt    microglia_WT    microglia_AD
US12302316_252665511011_S01_GE2-v5_95_Feb07_1_4.txt    microglia_AD    microglia_WT
US12302316_252665511015_S01_GE2-v5_95_Feb07_1_1.txt    microglia_WT    microglia_AD
US12302316_252665511015_S01_GE2-v5_95_Feb07_1_2.txt    microglia_AD    microglia_WT
US12302316_252665511015_S01_GE2-v5_95_Feb07_1_3.txt    microglia_AD    microglia_WT
US12302316_252665511016_S02_GE2-v5_95_Feb07_1_1.txt    astrocytes_AD    astrocytes_WT
US12302316_252665511016_S02_GE2-v5_95_Feb07_1_2.txt    astrocytes_WT    astrocytes_AD
US12302316_252665511016_S02_GE2-v5_95_Feb07_1_3.txt    astrocytes_WT    astrocytes_AD
US12302316_252665511016_S02_GE2-v5_95_Feb07_1_4.txt    astrocytes_AD    astrocytes_WT

What I want to do for  DE gene analysis is to compare the followings:

1) microglia_AD vs microglia_WT

2) astrocytes_AD vs astrocytes_WT

3) astrocytes_AD vs microglia_AD

4) astrocytes_WT vs microglia_WT

I did the first 2 comparisons with the script below

library("limma")
wtFun <- function(x) as.numeric(x$ControlType == 0)
targets <- readTargets()
microglia <- targets[!(targets$Cy3 %in% "astrocytes_AD") & !(targets$Cy5 %in% "astrocytes_AD"),]
astrocytes <- targets[!(targets$Cy3 %in% "microglia_AD") & !(targets$Cy5 %in% "microglia_AD"),]
#1) microglia_AD vs microglia_WT
RG <- read.maimages(microglia, source = "agilent", annotation = "ProbeName", wt.fun = wtFun)
RG <- backgroundCorrect(RG, method = "normexp", offset = 50)
MA <- normalizeBetweenArrays(RG, method = "quantile")
MA <- avereps(MA, ID = MA$genes$ProbeName)
design <- modelMatrix(microglia, ref = "microglia_WT")
fit <- lmFit(MA, design)
fit <- eBayes(fit, trend = TRUE)
allResults <- topTable(fit, n = Inf)

#2) astrocytes_AD vs astrocytes_WT
RG <- read.maimages(astrocytes, source = "agilent", annotation = "ProbeName", wt.fun = wtFun)
RG <- backgroundCorrect(RG, method = "normexp", offset = 50)
MA <- normalizeBetweenArrays(RG, method = "quantile")
MA <- avereps(MA, ID = MA$genes$ProbeName)
design <- modelMatrix(astrocytes, ref = "astrocytes_WT")
fit <- lmFit(MA, design)
fit <- eBayes(fit, trend = TRUE)
allResults <- topTable(fit, n = Inf)

is there an easier way to design a design matrix that can do all the comparisons?

thank you very much in advance,

ilyas.

The experimental design does not include any direct comparisons between astrocytes and microglia, so you cannot make comparisons 3 and 4 in a conventional two-color microarray analysis. To make the first two comparisons, I would do it like this:​

library("limma")
targets <- readTargets("targets.txt")
RG <- read.maimages(targets, source = "agilent")
keep <- RG$genes$ControlType==0
RG2 <- RG[keep,]
RG2 <- backgroundCorrect(RG2, method = "normexp", offset = 16)
MA <- normalizeWithinArrays(RG2, method = "loess")
dm <- modelMatrix(targets[1:7,], ref="microglia_WT")
da <- modelMatrix(targets[8:11,], ref="astrocytes_WT")
design <- cbind(Dye=1, blockDiag(dm,da))
keep <- rowMeans(MA$A) > 6
fit <- lmFit(MA[keep,], design)
fit <- eBayes(fit, trend=TRUE, robust=TRUE)
topTable(fit, coef="microglia_AD")
topTable(fit, coef="astrocytes_AD")

Note: I made a correction to the design matrix in the above code 14 hours after my original post.

To make all four comparisons that you want to do, it is necessary to undertake a separate channel analysis instead:

targetsC <- targetsA2C(targets)
group <- factor(targetsC$Target)
group <- relevel(group,ref="microglia_WT")
designsc <- model.matrix(~0+group)
colnames(designsc) <- levels(group)
MA <- normalizeWithinArrays(RG2, method = "loess")
MA <- normalizeBetweenArrays(MA, method="Aquantile")
keep <- rowMeans(MA$A) > 6
isc <- intraspotCorrelation(MA[keep,], designsc)
fitsc <- lmscFit(MA[keep,], designsc, isc$consensus)
cont.matrix <- makeContrasts(
   mAD.mWT = microglia_AD - microglia_WT,
   aAD.aWT = astrocytes_AD - astrocytes_WT,
   aAD.mAD = astrocytes_AD - microglia_AD,
   aWT.mWT = astrocytes_WT - microglia_WT,
   levels = designsc
)
fit2 <- contrasts.fit(fitsc, cont.matrix)
fit2 <- eBayes(fit2, trend=TRUE, robust=TRUE)
topTable(fit2, coef="mAD.mWT")
topTable(fit2, coef="aAD.aWT")
topTable(fit2, coef="aAD.mAD")
topTable(fit2, coef="aWT.mWT")

Note that the separate channel analysis is more powerful and will generally give smaller p-values than the conventional analysis for the same comparisons.