coseq plotting error: no method for coercing this S4 class to a vector
1
0
Entering edit mode
minyaaa9058 ▴ 20
@minyaaa9058-22760
Last seen 2.2 years ago

Hi!

I'm using coseq to cluster my RNAseq data processed with Tximport and then DESeq. I have followed the post from https://support.bioconductor.org/p/118586/#118671 to modify my code so everything was running smoothly. However, when I tried to plot my coseq result, I got an error message and I couldn't figure out how to fix it (sorry I'm still a beginner in R). My code was:

# my input data dds.raw was generated using Tximport: 
# dds.raw <- DESeqDataSetFromTximport(txi.kallisto, colData=sampleTable, design=~condition)
>dds <- DESeq(dds.raw, test = "LRT", reduced = ~1)
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 132 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

> res <- results(dds, independentFiltering=FALSE)

> run <- coseq(counts(dds, norm=TRUE)[which(res$padj < 0.05),], 
+              normFactors = "DESeq", transformation = "logclr", 
+              model = "kmeans", parallel = TRUE, meanFilterCutoff = 5, K=2:30)
****************************************
coseq analysis: kmeans approach & logclr transformation
K = 2 to 30 
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running K = 2 ...
Running K = 3 ...
Running K = 4 ...
Running K = 5 ...
...
Running K = 28 ...
Running K = 29 ...
Running K = 30 ...

> run
An object of class coseqResults
 4382 features by 32 samples. 
 Models fit: K = 2 ... 30
 Chosen clustering model: K = 16

> plot(run, graphs="boxplots")
Error in as.vector(x, mode = "numeric") : 
  no method for coercing this S4 class to a vector

Thank you very much!

coseq RNAseq clustering • 395 views
ADD COMMENT
0
Entering edit mode

Hi!

Thanks for your question. I haven't been able to reproduce your error on my end, so I could use a little more detail from you to figure this out.

Specifically, here's a fully reproducible example that runs with no error for me using some of the tximportData kallisto format example data (I'm guessing that's the pipeline you used given the commented-out first line of your code above?):

library(tximportData)
library(tximport)
library(readr)
library(DESeq2)
library(coseq)

## Read in some default tximport data, run DESeq2
dir <- system.file("extdata", package = "tximportData")
samples <- read.table(file.path(dir, "samples.txt"), header = TRUE)
files <- file.path(dir, "kallisto_boot", samples$run, "abundance.h5")
names(files) <- paste0("sample", 1:6)
txi.kallisto <- tximport(files, type = "kallisto", txOut = TRUE)
sampleTable <- data.frame(condition = factor(rep(c("A", "B"), each = 3)))
rownames(sampleTable) <- colnames(txi.kallisto$counts)
dds.raw <- DESeqDataSetFromTximport(txi.kallisto, 
                                    colData=sampleTable, 
                                    design=~condition)
dds <- DESeq(dds.raw, test = "LRT", reduced = ~1)
res <- results(dds, independentFiltering=FALSE)

run <- coseq(counts(dds, norm=TRUE),  
             normFactors = "DESeq", transformation = "logclr", 
             model = "kmeans", parallel = FALSE, 
             meanFilterCutoff = 5, K=2:15)

plot(run, graphs="boxplots")

Otherwise, the only thing I've changed with respect to your code is that I removed the filter on padj < 5% in coseq given that this is an artificial example and there are not very many DE genes. Could you run that and let me know whether it works for you? Also, the full output of sessionInfo() would be useful to know what versions of packages you're working with.

If you do still see an error, it would also be immensely helpful to me if you could send a fully reproducible example for me to test out on my end.

Thanks! Andrea (coseq developer)

ADD REPLY
0
Entering edit mode

Thank you for your quick response!

I was able to make the plot following your code, and even with adding [which(res$padj < 0.05),]

The only difference between today’s run and yesterday’s (which produced error) is that yesterday I was directly loading a .RData, which is saved from dds.raw <- DESeqDataSetFromTximport(txi.kallisto, colData=sampleTable, design=~condition) But today I went through the whole process without loading the data, and I think it’s very likely that problem. I’m sorry I should have tried that yesterday before I asked about the error.

ADD REPLY
0
Entering edit mode

Ok, thanks for letting me know! Don't hesitate to submit another question if you run into other issues.

ADD REPLY
0
Entering edit mode
andrea.rau ▴ 80
@andrearau-7032
Last seen 5 weeks ago
INRAE / Jouy en Josas, France

It seems the question was linked to an issue that arose when loading an existing DESeqDataSet object from an .RData, but the problem disappeared when rerunning code from start to finish.

ADD COMMENT

Login before adding your answer.

Traffic: 237 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6