Hi folks,

We are using pseudobulk DESeq2 (DESeq2_1.26.0) approach to perform DGE testing on single cell datasets.

We are interested in assessing condition induced differences, while controlling for other variables such as RIN, PMI, Age and tissue (shown in metadata).

   Age RIN PMI  Tissue  Condition
1   45 9.0  16      A   Disease
2   71 8.8  35      A   Control
3   47 9.3  25      A   Disease
4   38 7.3  14      A   Disease
5   58 7.4  21      A   Disease
6   85 8.8  26      A   Control
7   56 7.3  18      A   Disease
8   53 9.2  27      A   Disease
9   64 7.6  20      A   Control
10  51 9.4  26      A   Disease
11  54 9.3  15      A   Disease
12  67 7.8  23      A   Control
13  37 9.2  17      A   Disease
14  43 7.9  32      A   Disease
15  38 7.1  27      A   Control
16  43 8.2  32      A   Disease
17  37 8.6  32      A   Control
18  57 8.8  11      A   Control
19  42 8.2  21      A   Control
20  47 8.1  24      A   Control
21  53 9.2  27      A   Disease
22  67 7.2  23      B   Control
23  58 8.2  21      B   Disease
24  38 7.5  27      B   Control
25  67 7.2  12      B   Disease
26  64 7.3  20      B   Control
27  58 8.2  29      B   Disease
28  91 7.0  25      C   Control
29  89 7.4  34      C   Control
30  64 7.2  20      C   Control
31  81 6.9  10      C   Disease
32  74 6.8  29      C   Disease
33  84 7.6  28      C   Disease

However, when using:

dds <- DESeqDataSetFromMatrix(pseudobulk, sample_table, ~Age + RIN +PMI+tissue+disease)

following warning is shown:

  converting counts to integer mode
     the design formula contains one or more numeric variables with integer values,
     specifying a model with increasing fold change for higher values.
     did you mean for this to be a factor? if so, first convert
     this variable to a factor using the factor() function
     the design formula contains one or more numeric variables that have mean or
     standard deviation larger than 5 (an arbitrary threshold to trigger this message).
     it is generally a good idea to center and scale numeric variables in the design
     to improve GLM convergence.
    Warning message:
    In DESeqDataSet(se, design = design, ignoreRank) :
     some variables in design formula are characters, converting to factors

Going further with this approach, we were able to obtain some results

dds  <- DESeq(dds, test = "LRT", reduced = ~Age+ RIN +PMI+tissue, sfType = "poscounts", useT = T, minmu = 1e-6, minReplicatesForReplace=Inf)

As a second approach, in order to avoid above warning, we converted numeric variables to factors using factor() and call the same function

 dds <- DESeqDataSetFromMatrix(pseudobulk, sample_table, ~Age + RIN +PMI+tissue+disease)

following error is thrown:

 Error in checkFullRank(modelMatrix) : 
     the model matrix is not full rank, so the model cannot be fit as specified.
     One or more variables or interaction terms in the design formula are linear
     combinations of the others and must be removed.

In the first approach where some variables were numeric, the warning was given but it worked and we obtained some results.

In the second approach it didn’t work at all. Now we are not quite sure, why there is no linear combination when dealing with numerics, but it appears when converting them to factors.

So the question is, is it right, to convert the numeric values here as factors and if so, how we can check if there is really a linear combination? As we are dealing with more then 30 replicates, it is not obvious to check, if there are some, as most of the numeric variables are highly diverse between the replicates.

The second question is, is it ok to regress so many different variables to test only condition induced differences?

Thanks a lot for your help!

When you convert the numeric data to factors, there will definitely be an issue, as reported by DESeq2. The problem is that, once converted to factors, each individual number is treated as its own category / level, and the model becomes impossible to fit; hence, the checkFullRank errror.

As per the initial warning message, you could try to scale these variables via the scale() function, which will, per variable, center the variable (subtract the mean), and then divide by the standard deviation - this is essentially a Z transform, and it makes the distribution of each variable more amenable to a regression.

A question back to you, though: do you need these extra covariates in the design formula? Does there exist preliminary or previously published evidence that they are likely confounders of gene expression in your disease / condition under study? I ask because I frequently see people filling up their design formula with many parameters (and I used to do it all of the time, too), when the best solution is to keep the model relatively simple and mitigate biasing factors via a solid study design.

Kevin