Dear Prof. Smyth and limma experts,

I would like to clarify if my design and logic in applying limma is correct. I have a timecourse where the blood was taken from participants before vaccination (Day0), after first vaccination (Day1) and after second vaccination eight weeks later (Day57). These participants were subsequently challenged with the infectious pathogen and monitored.

Some individuals were protected while others developed the disease (and treated). Our aim is to look for genes correlated with the protection status.

If we subset the dataset for each timepoint and analyzed for protection status, only a couple of genes are significant. However, when we analyze all time points together with adjustment for timepoint, we obtain many more significant genes. The phenotype data looks something like this:

     SubjectID   Timepoint   Status              Age  Sex  Batch

     1001           Day0          NotProtected   20     M     B1

     1001           Day1          NotProtected   20     M     B2

     1001           Day57        NotProtected   20     M     B3

     1015           Day0          Protected         35     F     B2

     1015           Day1          Protected         35     F     B3

     1015           Day57        Protected         35     F     B1

      ...               ...                ...                     ...      ...    ....

And here is the model matrix and limma call for analyzing all time points jointly.

   mm <- model.matrix( ~ -1 + Age + Sex + Batch + Status + Timepoint, data=pheno )

   dupCor <- duplicateCorrelation( EXPRS, design=mm, block=pheno$SubjectID )$consensus 

   fit <- lmFit( EXPRS, design=mm, block=pheno$Subject, correlation=dupCor )

I am slightly concerned that we are violating some basic assumption here. Or is the blocking for subject and duplicateCorrelation sufficient? Any advice is greatly appreciated. Many thanks.

Regards, Adai

<font color="#0782c1">--</font>

Adaikalavan Ramasamy

Senior Leadership Fellow in Bioinformatics

Head of the Transcriptomics Core Facility

The Jenner Institute, University of Oxford

Roosevelt Drive, Oxford OX3 7DQ

Email: adaikalavan.ramasamy@ndm.ox.ac.uk

Office: 01865 617 100

Mob: 07906 308 465

http://www.jenner.ac.uk/transcriptomics-facility

The Age factor in your design seems a bit strange to me. If this is a numeric vector, the corresponding coefficient would be modelling genes where the (log-)expression is linearly proportional to the age of the subject. This seems unrealistic, e.g., would you really expect to see genes that double (or halve) expression in 50 year olds compared to 25 year olds?

In any case, the age effect is implicitly accounted for when you block on SubjectID. This means that there's no need to specify Age in the design if you're not explicitly interested in it. The same could be argued for Sex, though the specification of this factor in design is more sensible and could be retained (at the cost of residual degrees of freedom for variance estimation).

To address your main concern, blocking on SubjectID with duplicateCorrelation has some (not unreasonable) limitations. For example, it requires homogeneity in the subject effects, i.e., you shouldn't have systematic differences between subjects or an outlier subject. duplicateCorrelation also estimates a common correlation across all genes and uses that consensus value in generalized least squares. This is necessary to improve precision, but it also introduces some (acceptable) bias when the correlations are different between genes. The ideal strategy would consider subject-specific effects on a gene-by-gene basis. This could be done by putting SubjectID in the design, but the (interesting) coefficient for Status would then become unestimable. Thus, blocking with duplicateCorrelation is probably the best approach in this situation.

As an aside, one should avoid subsetting the data prior to a limma analysis. This will reduce the amount of information available to estimate the dispersion, resulting in greater uncertainty and a decrease in the number of DE calls.

I am also using this approach for each individual having more than one sample, each sample with a different condition.  One problem I have run into is that duplicateCorrelation takes a very long time to run, several days even.  Is there anyway to expedite the process?  It seems like it should be an "embarrassingly parallel" process and could easily threaded or run on a cluster.