Replacing aspects of the DESeqDataSet

[Rivalyn]

Upon analysing my RIPseq data with DESeq2, I noticed I had a large number of genes for which the fold change was read as "NA". This is because many of the genes have no reads in them. I wanted to change all the reads in my DESeqDataSet with a value of "0" to a value of "1" in order to avoid divide-by-zero errors resulting in "NA" fold changes. I thought it would work something like this:

Trimmed_BIP_dds[Trimmed_BIP_dds$assay == 0] <- 1

But that gives me this error, probably because DESeqDataSets can't be indexed like normal data frames:

Error in Trimmed_BIP_dds[Trimmed_BIP_dds$assay == 0] <- 1 : 
  object of type 'S4' is not subsettable

The help file on SummarizedExperiment instances (of which DESeqDataSet is a subtype) says this:

x[i,j], x[i,j] <- value:
Create or replace a subset of x. i, j can be numeric, logical, character, or missing. value must be a SummarizedExperiment instance with dimensions, dimension names, and assay elements consistent with the subset x[i,j] being replaced.

My problem is that I don't really understand what that means. I thought I should do something like this:

Trimmed_BIP_dds$assay[==0] <-1

But that doesn't work either. I also tried extracting the counts with the assay function and then changing them inside the matrix:

Trimmed_BIP_matrix <- assay(Trimmed_BIP_dds)
Trimmed_BIP_matrix[Trimmed_BIP_matrix == 0] <- 1

But then I couldn't figure out how to get the new matrix back into the DESeqDataSet object.

Can anyone help me with this?

[dpryan]

Given a DeseqDataSet named dds:

IDX <- which(assay(dds) == 0)
assay(dds)[IDX] <- 1

Having said that, if the genes don't have any reads for any of the samples then a fold-change of NA would seem completely appropriate. Further, I really wouldn't recommend mucking with raw counts like this, you'll inevitably change a few of the results (so don't blame me if you shoot yourself in the foot).

[Rivalyn]

A lot of our read counts look like this:

Control = 0
Treated = tons

In this case it seems to me that we're discarding perfectly useful data due to a simple math error. I see your point, though.

Thanks for the help, I'll give it a whirl once I get back to work tomorrow.