Bioinformatics Stack Exchange is a question and answer site for researchers, developers, students, teachers, and end users interested in bioinformatics. It only takes a minute to sign up.
Sign up to join this community
I obtained a matrix of RNA-seq count data that has been normalized by DESeq2's median of ratio method.
I know that DESeq2 wants to take in un-normalized counts, but I do not have access to those data.
How do I best proceed here if I want to perform DEG analysis using DESeq2?
I know I can always start from .fastq files, but that would be so much extra work.
I don't think I can un-normalize the data because I don't know the scaling factors.
Any suggestions?
You can read in the normalized count table and don't normalize the data, but my advice here is not to do that. The rounding of the normalized matrix introduces some noise, but I think the larger issue is how sure are you that the table you are working with, is exactly a count table of normalized counts from DESeq2 ?
It really doesn't take much time to align and get count table again, rather than going back and forth wondering whether you did the normalization correctly.
I can show below with an example how it will look like, first we make an example dataset and obtain the normalized counts:
library(DESeq2)
set.seed(111)
sz = runif(6,min=0.5,max=1.5)
x = makeExampleDESeqDataSet(sizeFactors=sz,m=6)
x = estimateSizeFactors(x)
ncounts = counts(x,normalize=TRUE)
x = estimateDispersions(x)
Now suppose we feed this into a new DESeq2 object, but before that you need to do is check most of the size factors are close to 1:
estimateSizeFactorsForMatrix(round(ncounts))
sample1 sample2 sample3 sample4 sample5 sample6
1.039179 1.045991 1.048276 1.045872 1.032420 1.052006
If you see something above that suggest a deviation from 1, something is wrong from the normalization and do not even think about using it.
Below you can use the normalized counts and you can see the difference is not huge:
dds = DESeqDataSetFromMatrix(round(ncounts),colData(x),~condition)
sizeFactors(dds) = 1
dds = estimateDispersions(dds)
Compare the new count table:
head(counts(dds,normalize=TRUE))
sample1 sample2 sample3 sample4 sample5 sample6
gene1 0 2 0 1 0 0
gene2 9 8 16 6 1 2
gene3 23 40 8 38 11 41
gene4 0 0 0 6 8 0
gene5 3 2 0 5 1 2
gene6 0 1 3 19 4 9
And your input:
head(ncounts)
sample1 sample2 sample3 sample4 sample5 sample6
gene1 0.000000 2.2658464 0.000000 0.9310803 0.000000 0.000000
gene2 8.693480 8.3081033 16.251453 5.5864817 1.121781 2.086619
gene3 23.472395 40.0299524 7.584012 38.1742919 11.217808 40.689076
gene4 0.000000 0.0000000 0.000000 5.5864817 7.852465 0.000000
gene5 3.477392 1.5105642 0.000000 4.6554014 1.121781 2.086619
gene6 0.000000 0.7552821 3.250291 18.6216058 4.487123 9.389787
We can check the dispersion, for lowly expressed genes you see some difference, but on the whole difference
with(elementMetadata(x),plot(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#FFC93C80"))
with(elementMetadata(dds),points(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#DBF6E980"))
Again, I would say I agree with @ATpoint comments that if you can, start from the fastq again.. don't work with data you have doubts about.
You might be able to cheat by giving DESeq a custom size factor matrix where all the values are 1. Then the software will feel like it's done the normalization step, but the 'normalized' values will be what you want.
