Hi everyone,

I have a bit of a problem with DESeq2 which I have already seen described and answered before, but I think I’m missing some key information to correctly tackle it. To summarise the experiment, I have a set of samples from patients that either respond well to a certain treatment (R) or don’t respond well (NR), and samples were obtained from each patient at up to three different (t0, t1 and t2). Each sample obtained was subdivided into cells positive for a certain cell surface marker (pos) and negative for that marker (neg). Right now, I am only interested in comparing gene expression between two discrete groups: between responders and non-responders, only at t0 and of positive cells, at t0 of negative cells, at t1 of positive cells, and at t1 of negative cells. For this, I created a column which nests all these features into a single definition, with these groups being R_pos_t0, NR_pos_t0, R_neg_t0, NR_neg_t0, R_pos_t1 and NR_pos_t1, R_neg_t1 and NR_neg_t1 (as well as R_pos_t2 and R_pos_t2 which won’t be compared to non-responders since there are no t2 non-responder samples) so it should be quite straight forward. However, I now have quite a few more variables I wish to add in order to control for these variables that might affect response. These variables are: percentage positivity for another marker (pdl1), sex, age, smoking habits (smoke), ecog, stage of disease, histology of disease, previous treatment (chemo), recist, and type of current treatment (treatment).

I understand that the correct model design would be the following:

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition + condition:pdl1 + condition:sex + condition:age + condition:smoke + condition:ecog + condition:stage + condition:histology + condition:chemo + condition:recist + condition:treatment)

However, I get to a problem which I have already seen described, which is that the model cannot be used because the resulting model is not full rank since it produces columns with zero values due to all the interactions. I have read from responses here that it is possible to create a custom model matrix and feed it to DESeq2 after eliminating zero columns, which I assume would be done like this:

mm = model.matrix(~ condition + condition:pdl1 + condition:sex + condition:age + condition:smoke + condition:ecog + condition:stage + condition:histology + condition:chemo + condition:recist + condition:treatment, colData(dds))

mm <- mm[, colSums(mm == 0) != nrow(mm)] 

dds = DESeq(dds, full=mm, betaPrior=FALSE)

This would allow to first create the custom matrix with all interactions, then eliminate all zero columns from said matrix, and then feed the matrix to DESeq2. However, what I am missing is how do I get to this point. How can I create a custom model matrix that requires colData information from a dds object that in turn hasn’t been generated due to the non zero values in the columns? This is the step that I am missing, creating the model matrix if DESeqDaraSetFromMatrix has not been able to generate a dds object. It’s probably something very silly that I am missing but it is something that I cannot find or that I’m not looking for properly.

Thank you very much in advance for any responses!

Kind regards Jesús

How can I create a custom model matrix that requires colData information from a dds object that in turn hasn’t been generated due to the non zero values in the columns?

Two options:

1) you can use a data.frame instead of colData(se). E.g.: model.matrix(~var1 + var2, data=df)

2) you can specify ~1 to create a dds then access colData(dds), then re-assign the design with design(dds) <- ...