Mixed-model in Limma

[KamilMaliszArdigen]

Description of feature

We would like to extend this pipeline logic to allow processing of data and model described bellow. Can you suggest how we should approach this topic?

Consider the following experiment frame: Patient	Condition	Tissue
1	Diseased	A
1	Diseased	B
2	Diseased	A
2	Diseased	B
3	Diseased	A
3	Diseased	B
4	Normal	A
4	Normal	B
5	Normal	A
5	Normal	B
6	Normal	A
6	Normal	B

[pinin4fjords]

Just to be super clear, can you illustrate the differences in how you would construct the models for this please?

[mzenczak]

Just to be super clear, can you illustrate the differences in how you would construct the models for this please?

@pinin4fjords if I understood you correctly here is the answer to your question:

Following the steps for analyzing mixed models, we first combine the experimental factors into a single factor, as shown in the table below (serving as the sample.sheet). This newly created factor becomes our contrast variable, which we then use in the model construction. See the code snippet below for details:

Patient	Condition	Tissue	Condition.Tissue
1	Diseased	A	Diseased.A
1	Diseased	B	Diseased.B
2	Diseased	A	Diseased.A
2	Diseased	B	Diseased.B
3	Diseased	A	Diseased.A
3	Diseased	B	Diseased.B
4	Normal	A	Normal.A
4	Normal	B	Normal.B
5	Normal	A	Normal.A
5	Normal	B	Normal.B
6	Normal	A	Normal.A
6	Normal	B	Normal.B

DGE <- DGEList(count.table)
DGE <- calcNormFactors(DGE)
contrast_variable <- "Condition.Tissue"
model <- paste("~ 0 +", contrast_variable)
design <- model.matrix(as.formula(model), data = sample.sheet)
voom <- voom(counts = DGE, design = design)
corfit <- duplicateCorrelation(object = voom, design = design, block = Patient)
fit <- lmFit(object = voom, design = design, block = Patient, correlation = corfit$consensus.correlation)

For example, to define a contrast that tests for deferentially expressed genes between diseased tissue A and normal tissue A, you could construct a contrast like this:

contrast <- "Condition.Tissue_Diseased.A_vs_Normal.A = Diseased.A - Normal.A"
contrast_matrix <- makeContrasts(contrasts = contrast, levels = design)
fit <- contrasts.fit(fit, contrast_matrix)
fit <- eBayes(fit)
results <- topTable(fit)

Hope this answers your question. :)

[pinin4fjords]

IIUC the contrast can be dealt with by supplying the workflow with a sample sheet containing the composite variable, and building contrasts on that variable in the normal way.

After that the only missing thing is the duplicateCorrelation treatment, which could be PR'd to the limma module as a new feature.

[grst]

There's also the DREAM method from the bioconductor package variancePartition which is an extension of limma for linear mixed effects models:

https://gabrielhoffman.github.io/variancePartition/articles/dream.html

It allows to specify random effects using an extension of Wilkinson formulas that is also adopted by other R packages and might be more convenient than what's described in the limma vignette:

form <- ~ Disease + (1 | Individual)