Entering edit mode
Hi I have done a DESeq2 DGE analysis and it showed an error when I was using ashr shrinkage
using 'ashr' for LFC shrinkage. If used in published research, please cite:
Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
https://doi.org/10.1093/biostatistics/kxw041
warning: solve(): system is singular (rcond: 4.52254e-18); attempting approx solution
warning: solve(): system is singular (rcond: 7.84057e-18); attempting approx solution
warning: solve(): system is singular (rcond: 3.78635e-19); attempting approx solution
warning: solve(): system is singular (rcond: 6.26207e-19); attempting approx solution
warning: solve(): system is singular (rcond: 7.78489e-19); attempting approx solution
warning: solve(): system is singular (rcond: 1.14601e-18); attempting approx solution
I assume that stems from co-linearity of some of the factors? Is there any fix? May I ask if the result is trustworthy? And if I am going to use it is there anything I should be careful of?
The code I used:
st$sv <- svobj$sv #st is the sample table. svobj$sv is the sva variables estimated for batch effect correction
ddsTxi.1 <- DESeqDataSetFromTximport(txi, #txi is the tximport from kalisto of a cancer patient cohort
colData = st,
design = ~ 0 + Final.Label + sv) #Final label is the subtype of the cancer of each patient
# Pre-filtering. Filter out samples with TPM below 2 in less than 4% of the sample (there are some rare subtypes)
library(genefilter)
abundance = txi$abundance
abs_crit = pOverA(0.04, 2)
abs_filter <- filterfun(abs_crit)
keep <- genefilter(abundance, abs_filter)
dds.1 <- ddsTxi.1[keep,]
dds.1 <- DESeq(dds.1)
res.1 <- lfcShrink(dds.1, contrast = c(-1/12,-1/12,-1/12,-1/12,1,-1/12,-1/12,-1/12,-1/12,-1/12,-1/12,-1/12,-1/12,0), type = "ashr")
Edit: Thank you Michael for pointing towards the solution here: Warnings when using "ashr" in lfcshrinkage for DESeq2
Hi Michael. Thanks for the heads-up. I actually read some of the posts here and on Biostars regarding to the same error, of which most of them you have replied. It seems like most people tried to fix it by filtering out more genes with low counts but the problem persists if I interpreted what they replied correctly. Do you recall what kind of fix it is? Maybe it is just a few google search away but I have the wrong keywords.
I searched ashr and did a time sort:
https://support.bioconductor.org/post/search/?query=Ashr&order=date
It’s the most recent after this one.