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]]