Hi Michael and other DESeq2 users,

My dilemma is setting the best analysis for this design. 

I have been trying the following code:

dds <-DESeqDataSetFromMatrix(countData=counttable, colData=colData, design= ~group)

dds <- DESeq (dds)

This works but I want to take account the patient so making a paired design. 

So I did this:

design(dds) <- ~ group + patient

then rerun DESeq again,

I got the model matrix not full rank error. I would like to do a paired analysis where I have to compare CondA_7 of a specific subject to CondA_0 of the same patient, CondA_63 of the same patient vs CondA_0 of the same patient then compare them as a group afterwards versus other groups, taking account the time and the condition for each.So I tried this:

mm1 <- model.matrix(~ group + group:patient, colData)

mm1
idx <- which(colSums(mm1 == 0) == nrow(mm1))
mm1 <- mm1[,-idx]
mm1
mm0 <- model.matrix(~ group, colData)
colnames(mm1)

dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds, modelMatrix=mm1)

error: inv(): matrix seems singular
Error: inv(): matrix seems singular

dds <- estimateDispersionsFit(dds)

found already estimated fitted dispersions, removing these

dds <- estimateDispersionsMAP(dds, modelMatrix=mm1)

found already estimated dispersions, removing these
using supplied model matrix
dds <- nbinomLRT(dds, full=mm1, reduced=mm0)

using supplied model matrix
found results columns, replacing these
3 rows did not converge in beta, labelled in mcols(object)$fullBetaConv. Use larger maxit argument with nbinomLRT

Did I run it correctly? What am I missing? I have search on most paired analysis and tried most of them but I think I am missing a great point to make this work. Im currently using DESeq2 v. 1.12.3.

Here is the colData:

colData

Thank you!

I'd first recommend that you have three separate variables for use in the design. You already have 'time' and 'patient' and likewise create 'condition' with levels A,B,C,D.

Secondly, you need to make sure all three of these are factor. Time should be a factor with three levels: 0, 7, 63.

Then you can follow the advice in the vignette in Section 3.12 “Model matrix not full rank"

It begins with: 

"Consider an experiment with grouped individuals, where we seek to test the group-specific effect of a treatment, while controlling for individual effects."

The vignette describes group, individual and condition. For your experiment, these are A,B,C,D (what you have called condition), patient, and time. You will have to create a nested variable for the patient variable as in the vignette example. Do you have equal number of patients in each condition?

Hi Mike,

I checked all and they are factors. Yes they also do have equal number of patients in each condition which is 5 patients, 3 time points.

This is what I did:

condition=factor(colData$condition)

time=factor(colData$time)
patient=factor(colData$patient)

> class (time)
[1] "factor"
> class(condition)
[1] "factor"
> class(patient)
[1] "factor"

Then I ran them all again for DESeqDataSetFromMatrix with this design.

dds<-DESeqDataSetFromMatrix(countData=counttable,
                           colData=colData,
                           design=~time+condition+patient)

Still didn't work,

Then I tried similar to the vignette,

colData$ind.n<-factor(rep(rep(1:2,each=3),10))
model.matrix(~group+group:ind.n+group:condition, colData)

I think I figured out, I will just make a nested design at the end. Thanks.

Another question:

Am I correct that log2Foldchange >1 means that the output given by DESeq2 is a fold change of 2. 

Is there also a way of getting a ratio of fold change in DESeq wherein after doing the contrast of a group you like, I can also compare the ratio of two fold changes among each other and get DEgenes?

CondADay7vsCondADay0<-results(dds, contrast=c('group','CondADay7', 'CondADay0'))

CondDDay7vsCondDDay0<-results(dds, contrast=c('group','CondDDay7', 'CondDDay0'))

then

CondADay7 vs CondDay0 / CondDDay7 vs CondDDay0

A colleague of mine suggested to add another column called nested.patient then the design will be nested.patient+group, is this correct?Can you explain why not design= ~group+nested.patient?