nested samples in design DEseq2

Question

I'm running a DEG analysis with Deseq2 by specifing the following design that includes 12 samples from 6 pts in 4 different conditions:

coldata <- data.frame("sample"=c(
  "T58", "T86", "T87",
  "D58", "D86", "D87",
  "C62", "C80", "C84",
  "C62", "C80", "C84"
))
coldata$Pts <- c(
  "m58", "m86", "m87",
  "m58", "m86", "m87",
  "c62", "c80", "c84",
  "c62", "c80", "c84"
)
coldata$class <- c(
  "Treated", "Treated", "Treated",
  "Distal", "Distal", "Distal",
  "CLEFT", "CLEFT", "CLEFT",
  "CRIGHT", "CRIGHT", "CRIGHT"
)

I would model the differences among the conditions by considering that the same individuals are nested in different conditions

I tried the following code

ddsMat <- DESeqDataSetFromMatrix(
  countData = as.matrix(matrix),
  colData = coldata,
  design = ~ class+Pts
)

but this produces the famous error

"the model matrix is not full rank..."

If I reduce the previous dataset to

coldata2 <- as.data.frame(c(
  "T58", "T86", "T87",
  "D58", "D86", "D87"
))
coldata2$Pts <- c(
  "m58", "m86", "m87",
  "m58", "m86", "m87"
)
coldata2$class <- c(
  "Treated", "Treated", "Treated",
  "Distal", "Distal", "Distal"
)

then everything goes well and I'm able to create a DESeq object with the following design:

ddsMat <- DESeqDataSetFromMatrix(
  countData = as.matrix(matrix[,1:6]),
  colData = coldata2,
  design = ~ class+Pts
)

Please can you provide to me some suggestion on how to solve the problem?

Thank you

Answer 1

It is actually simple. The first two conditions are nested with the Pt's starting with m and the other two conditions are nested with those PT's starting with c. That means you cannot run this design. The only solution is to subset the data into two chunks and analyze separately.

Answer 2

The DESeq2 documentation explains this error as follows:

While most experimental designs run easily using design formula, some design formulas can cause problems and result in the DESeq function returning an error with the text: “the model matrix is not full rank, so the model cannot be fit as specified.” There are two main reasons for this problem: either one or more columns in the model matrix are linear combinations of other columns, or there are levels of factors or combinations of levels of multiple factors which are missing samples.

In your case, I think the error is popping up because there are too many defined values in the "Pts" variable. DESeq2 doesn't require you to separately define replicates as a variable; it assumes that different data columns are taken from different biological replicates.

Getting the models right is tricky, and needs a deeper discussion about what the goal of analysis should be. My guess is that the following should work, but only if you're not interested in comparing expression results from different patients:

coldata$Pts <- c(
  "m", "m", "m",
  "m", "m", "m",
  "c", "c", "c",
  "c", "c", "c"
)
coldata$class <- c(
  "Treated", "Treated", "Treated",
  "Distal",  "Distal",  "Distal",
  "CLEFT",   "CLEFT",   "CLEFT",
  "CRIGHT",  "CRIGHT",  "CRIGHT"
)

This may not be what you want to do; other ways to deal with this are explained in the DESeq2 documentation on the Bioconductor website:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#model-matrix-not-full-rank

Answer 3

check if Pts and class are of factor class or not by using is.factor(coldata$Pts) and is.factor(coldata$class)

if it is not then convert it into factor by

coldata$Pts = factor(coldata$Pts)

coldata$class = factor(coldata$class)

Then try to design your matrix (A design matrix is a matrix of values of the grouping variable). I think this would work as it works with

design <- model.matrix(~0+coldata$Pts+coldata$class) in limma package.

You can also check in this RNA-seq analysis workflow of DESeq2 where they also used nested samples for their design matrix.

Answer 4

Your original was,

ddsMat <- DESeqDataSetFromMatrix(countData = as.matrix(matrix), # colData = coldata, design = ~ class+Pts)

What happens when you use

ddsMat <- DESeqDataSetFromMatrix(countData = as.matrix(matrix), colData = coldata, design = ~ class+Pts)

Thus # is removed. The suspicion is # is causing the code to fall over. # is used for comments in R so the presence of this within a line of code will cause the following code colData = coldata, design = ~ class+Pts to not be interpreted because it can't be seen.