Hello, I am having a lot of difficulty setting up my design matrix. I have been searching extensively online for an answer but I think I lack some of the background to understand how to do this since my experiment is fairly complicated. I think I have a basic grasp of setting up model matrices in simple experiments but contrast coding is still a complete mystery to me.

This is a RIP-Seq experiment where I pulled down a protein a of interest with an antibody and then extracted the attached RNA which was then sequenced. This was performed in either wild-type (WT) tissue or protein knock-out (KO) where the KO is assumed to represent some kind of noise or background for the experiment. WT or KO tissue is from animals injected with various drugs.

My basic design is as follows:

I have tried setting up my model matrix a few different ways using model.matrix() but I don't always see all the coefficients I am interested in. The kind of interaction I am interested in would be either Control_WT:Control_KO:Drug2_WT:Drug2_KO or Drug1_WT:Drug1_KO:Drug2_WT:Drug2_KO. In both of these cases, what I am really interested in is the WT conditions fold-change (i.e. Control over Drug2) but I want the KO noise to be modeled in as well. This 4-way interaction seems completely wrong and way too complicated. How is such modeling normally done? Limma refuses to use such columns if I try ("coefficients not estimable").

This brings me to a second question. Say I just had one coefficient ("Drug1_WT" for example). What exactly is limma giving back to me after lmfit() > eBates() > topTable()$AveExpr? This condition isn't being compared to anything is it? Or would it be compared to any other coefficients that come before it if I used an additive model to set up the regression equation? The limma reference seems to suggest that the "AveExpr" column is the fold change compared to all other columns in the regression equation. If that is the case, then [limma documentation does not say this, see Gordon Smyth's answer below] should my coefficients really be Control_WT:Control_KO and then Drug2_WT:Drug2_KO? Then when I perform the fit and use topTable() to find the Control_WT:Control_KO coefficient, the AveExpr column would reflect the fold change relative to Drug2_WT:Drug2_KO?

Sorry if this is really confusing. Any guidance on how to best set up my model matrix (and/or contrast matrix if it's even needed) would be greatly appreciated.

The best way to proceed in a situation like this is to just compute the mean expression for all groups that you care about and then do the contrasts by hand. This is actually pretty simple - it's just algebra - so don't over think things.

You care about groups like Drug2_KO, which is the set of KO animals that got drug 2, right? So you can set up your model matrix like

grps <- paste(df$Condition, df$Genotype, sep = "_")

## where df is the name of your data.frame containing the group information

design <- model.matrix(~ 0 + grps)

## this will now compute the mean of each group

colnames(design) <- gsub("grps", "", colnames(design))

## get rid of extraneous cruft in the column names

contrast <- makeContrast(interaction1 = (Control_WT - Control_KO) - (Drug2_WT - Drug2_KO), levels = design)

This contrast will now calculate (literally!) (Control_WT - Control_KO) - (Drug2_WT - Drug2_KO), which is the interaction between KO and WT according to drug treatment. Or conversely, the interaction between drug and control according to genotype. Or in English, it finds genes that react to the drug2 treatment differently, depending on if the animals are WT or KO.

You can add any number of different interactions, to find genes that react to the different drugs depending on the genotype status by emulating what I have already done.

Does that help?

Your experimental design is a pretty standard one and is dealt with in Section 9.5 of the limma User's Guide. James has recommended that you follow the approach described in Section 9.5.2.

I am puzzled by your remarks about AveExpr. The limma User's Guide (page 60) says that

"The AveExpr column gives the average log2-expression level for that gene across all the arrays and channels in the experiment."

help("topTable") says the same thing:

"AveExpr: average log2-expression for the probe over all arrays and channels"

Is that not clear? I can't see how either of these statements could possibly suggest a fold change compared to anything, let alone to "all the other columns" as you've interpreted it as.