I have a problem in model formula for design matrix. The data looks like:

I would like to have a design matrix that allows me to do contrast in the paired/blocking samples like:

Disease1_P0vsDisease1_T0

Disease1_P0vs_Disease2_P0

How should I set the model formula? Thanks for your help!

For your first contrast, I would do something like:

Group <- paste0(Disease, Treatment, Time)
design <- model.matrix(~ 0 + factor(ID) + Group)
design <- design[,!grepl("P0", colnames(design))] # get to full rank.

The first three coefficients represent the patient effects and are not interesting. The next three coefficients represent the log-fold changes of P2, T0 and T2 over P0 for disease 1, "averaged" over patients 1 and 3. The last coefficient represents the log-fold change of T0 over P0 for disease 2.

So, the first contrast is fairly straightforward; just drop coefficient 5 (i.e., Groupdisease1T0).

Your second contrast is much harder as a comparison between disease states involves a comparison between patients. This means we can't block on ID in the design matrix, as your patients are nested within diseases; thus any ID blocking factor will absorb the disease effect. The only solution is to take only one sample from each patient (i.e., the P0 sample) and fit a separate model to those samples. Otherwise, the correlations between samples from the same patient will compromise error control in edgeR.

The alternative is to switch to voom-limma and use duplicateCorrelation() instead.

Hi Aaron, thanks for your reply.

ID is sample ID (patient identifier). Do you need additional info?

Thanks