Hi there,

I'd like some guidance on how best to carry out differential expression analysis using the limma pipeline with my data set. Below is representative of my targets data.frame:

Sample   Patient   Tissue   Time   Treatment   Class
  1        1      Diseased   0      Topical      A
  2        1      Diseased   2      Topical      A
  3        1      Diseased   4      Topical      A
  4        1      Healthy    0      Topical      A
  5        1      Healthy    2      Topical      A
  6        1      Healthy    4      Topical      A
  7        2      Diseased   0      Systemic     B
  8        2      Diseased   2      Systemic     B
  9        2      Diseased   4      Systemic     B
  10       2      Healthy    0      Systemic     B
  11       2      Healthy    2      Systemic     B
  12       2      Healthy    4      Systemic     B

I have data from patients undergoing either topical or systemic treatment, with samples taken at three time points: 0 (baseline) = before starting treatment 2 = two weeks after starting treatment 4 = four weeks after starting treatment

Samples were also taken from two tissue types in each patient = healthy and diseased.

Each patient also has either class A disease or class B disease.

There are a number of questions that I'd like to answer here. I'll state my intended approach and then hopefully I can get some pointers if I have the wrong idea.

(1) What is the difference between class A disease and class B disease at week 0?

Here, I'm guessing a simple two group comparison would be best. I.e. subset my data by grabbing just the week 0 samples and then do a two group comparison.

(2) What is the effect of each treatment type in each tissue type over time?

Here, section 9.7 on multi-level experiments in the limma manual seems most relevant. I'm thinking I should analyse each treatment type separately by subsetting my data, and then proceed like so:

> Treat <- factor(paste(targets$Tissue,targets$Time,sep="."))
> design <- model.matrix(~0+Treat)
> colnames(design) <- levels(Treat)

> corfit <- duplicateCorrelation(eset,design,block=targets$Patient)
> corfit$consensus

> fit <- lmFit(eset,design,block=targets$Patient,correlation=corfit$consensus)

> cm <- makeContrasts(
+ Diseased_twoweekchange = Diseased.2-Diseased.0,
+ Diseased_fourweekchange = Diseased.4-Diseased.0,
+ Healthy_twoweekchange = Healthy.2-Healthy.0,
+ Healthy_fourweekchange = Healthy.4-Healthy.0,
+ levels=design)

(3) Is the effect of each treatment type in each tissue type over time different in Class A compared to Class B?

This is the question I'm most unsure about. Should I proceed like I did in question 2 but use a three way combined factor (Tissue.Time.Class) and then compare the DEG lists after?

Thank you.

You can subset the data to answer different questions, but you will miss out on a lot of the power of limma by doing that. It is usually better to analyse all the samples at once:

• create a single Treat factor by combining Tissue, Time, Treatment and Class
• estimate the intra-patient correlation using duplicateCorrelation
• use makeContrasts to make any comparison(s) you need.

This gives an integrated analysis that allows you to answer all your questions in a unified way.