Hello,

I had a question about computationally singular matrices in DESeq and surrogate variables. After including all the SVs into my formula where I am analyzing paired samples (before/after treatment within a sample), I use the DESeq function and am returned with the error that my matrix is computationally singular. When I reduce the number of surrogate variables to 3, based on the n.sv() function using the "be" method, I no longer receive this error and am able to run the analysis. I want to know why the reduction of the number of SVs included in my design formula allows the dataset to run?

dds <- DESeqDataSetFromMatrix(countData = countdata, colData = phenotype, design = ~ study_id + SV1 + SV2 + SV3 + SV4 + SV5 + SV6 + SV7 + SV8 + SV9 + SV10 + SV11 + treatment)

I'm not getting a single output. it's a matrix of the pairwise study id's and their respective correlations. because i'm looking at before/after treatment within the same sample, these pairwise correlations are 1. Not all of them, just the ones within the same study id

I was expecting you would get a matrix of pairwise correlations between all the covariates in the design. There shouldn't be any covariates with correlation 1. Can you give an example?

(Intercept) 1 NA NA NA NA NA NA NA NA NA study

id121 NA 1.000000000 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.0303030303 studyid123 NA -0.030303030 1.00000000 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.0303030303 studyid124 NA -0.030303030 -0.03030303 1.00000000 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.0303030303 studyid125 NA -0.030303030 -0.03030303 -0.03030303 1.00000000 -0.03030303 -0.03030303 -0.03030303 -0.03030303 -0.0303030303 Warning message: In cor(model.matrix(design(dds), colData(dds))) : the standard deviation is zerothis is a portion of the matrix, and this is after the surrogate variables have been reduced to 3 (after using the n.sv(method = "be") in svaseq). i'm confused as to why this is happening. i suspect i may have set up my design formula incorrectly?

Oh, the diagonal is always 1, that's the correlation of a variable with itself. The concern would be off-diagonal high correlations.

So there is no concern with the output of this correlation? Were the 11 SV's just extremely correlated with the variation across the samples? Was it overcorrecting?

I don't know, but I would guess that some of those were too highly correlated with some other covariates.

thank you - it does seem that the first couple of SVs were highly correlated with cell type proportions.

Hi Michael,

I hope you don't mind if I ask you another question on this topic. I have a similar error where the vst function tells me my deseq object is computationally singular. The output from

`cor(model.matrix(design(dds), colData(dds)))`

shows that the correlation between some of my variables is 0.6 or -0.6. Therefore, I thought that that might be the problem. However, later I read this post: https://stats.stackexchange.com/questions/76488/error-system-is-computationally-singular-when-running-a-glm, where one of the answers suggests that there might be too big of a difference between the values of various variables in my design formula.The output of my metadata object looks like this:

As you can see, there is a pretty large difference between the libsize values and the sv1 values. The correlation between the libsize factors is 0.6, the correlation between the libsize factors and sv1 is also +0.6 or -0.6, which seems pretty high. I tried the suggestion from stackexchange myself and divided the values of the libsize variables by 10^6 so the libsize values and sv1 values would be closer to each other. Of course this does not influence the correlation, but it does solve my error! What's strange to me, is that using different factors to divide the values of libsize by (10^6 or 10^4), gives different outcomes in my heatmaps.

So my question is: how should I proceed? I can see that the libsize and sv1 values have a pretty large confounding effect on my dataset and therefore I want to add these variables to my design formula and then remove them with limma's removebatcheffect formula. However, I don't know how to transform the values of my libsize or sv1 variable correctly. Could you help me with this problem?

Center and scale numeric values, or at the least put them roughly on a scale between -10 and 10.

If you are just using

`vst`

, you can also use a design of ~1 which won't use the design at all.Note that current versions of DESeq2 also have a note about scaling numeric design variables, at object construction.

I centered and scaled the numeric values, but now I'm getting an error that the metadata is not full rank. I believe this is because some samples have approximately similar values for columns 2 and 3 but I don't think they are ever exactly the same. Is there a way to solve this error? I used the following code line for scaling and centering and a metadata example is below:

Check the correlation of the variables in your design.

E.g.:

In order to use

or

DDS has to be made:

However, creating DDS gives me the following error, so I cannot check the correlation between different variables.

I've tried to use the cor function on my

`centered_and_scaled_metadata`

but then I receive an error which says that my argument in cor must be numeric. To turn my metadata numeric, I must first delete my first column which contains characters as can be seen in the previous post. Now I can finally produce a numeric metadata with`cor(centered_and_scaled_metadata)`

which gives me the following:So there is no collinearity between variables and I do not understand why I'm receiving the full rank error. Could you help?

I'd recommend fixing the class of these columns ahead of time. Set the columns that are supposed to be numeric to numeric, likewise for the factor columns.

E.g., if you have a data.frame

`x`

, at the top of the script, set the proper classes:Do this before you start running DESeq2 or other diagnostics.

I'd recommend fixing the class of these columns ahead of time. Set the columns that are supposed to be numeric to numeric, likewise for the factor columns.

E.g., if you have a data.frame

`x`

, at the top of the script, set the proper classes:Do this before you start running DESeq2 or other diagnostics.

Hi Michael I'd love to hear your input on this issue

I would love some feedback on the error if you have the time