Dear Michael and DESeq2 users, I started using DESeq2 three weeks ago. Even I was having some normal problems for a newby person my experience was really positive and I was and I am detecting more DEG than with DESeq(1). However suddenly I am encountering a new problem that wasn't coming up previously: > resSig <- res[ res$padj < .1, ] Error in normalizeSingleBracketSubscript(i, x, byrow = TRUE, exact = FALSE) : subscript contains NAs For me this is pretty weird because I've been using the exactly same RScript and datasets for the last 3 weeks and I wasn't having this error before. My RScript is: #Count matrix input Cele_MRvsER_old = read.csv (file.choose(), header=TRUE, row.names=1) CeleDesign <- data.frame( row.names = colnames(Cele_MRvsER_old), condition = factor(c("MR", "MR", "ER", "ER", "ER", "ER"))) dds <- DESeqDataSetFromMatrix(countData = Cele_MRvsER_old, colData = CeleDesign, design = ~ condition) dds #Est size factor = normalize for library size dds<- estimateSizeFactors(dds) ddsLocal <- estimateDispersions(dds, fitType="local") ddsLocal <- nbinomWaldTest(ddsLocal) #plot dispersion plotDispEsts(ddsLocal) #Differential expression analysis resultsNames(ddsLocal) res <- results(ddsLocal, name= "condition_MR_vs_ER") res <- res[order(res$padj),] head(res) plotMA(ddsLocal) sum(res$padj < .1, na.rm=TRUE) #filter for upregulated and downregulated genes resSig <- res[ res$padj < .1, ] head(resSig) write.table(as.data.frame(resSig [ order(resSig$padj), ]), file="BRvsMR_old_MostSigGenes.txt") write.csv(as.data.frame(resSig[ order( resSig$log2FoldChange, -resSig$baseMean ), ]), file="BRvsMR_old_DownRegulated.csv") write.csv(as.data.frame(resSig[ order( -resSig$log2FoldChange, -resSig$baseMean ),]), file="BRvsMR_old_UpRegulated.csv") #up-regulated head(resSig[order(resSig$log2FoldChange,-resSig$baseMean),]) #down-regulated head(resSig[order(-resSig$log2FoldChange,-resSig$baseMean),]) #list most sig genes head( resSig[ order(resSig$padj), ] ) And my sessionInfo(): > sessionInfo() R version 3.0.2 (2013-09-25) Platform: x86_64-apple-darwin10.8.0 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] DESeq2_1.2.0 RcppArmadillo_0.3.920.1 Rcpp_0.10.5 GenomicRanges_1.14.1 [5] XVector_0.2.0 IRanges_1.20.0 BiocGenerics_0.8.0 loaded via a namespace (and not attached): [1] annotate_1.40.0 AnnotationDbi_1.24.0 Biobase_2.22.0 DBI_0.2-7 [5] genefilter_1.44.0 grid_3.0.2 lattice_0.20-24 locfit_1.5-9.1 [9] RColorBrewer_1.0-5 RSQLite_0.11.4 splines_3.0.2 stats4_3.0.2 [13] survival_2.37-4 tools_3.0.2 XML_3.95-0.2 xtable_1.7-1* * Any idea about why am I having this error now but not in the past and about how to fix it? Thank you!!! [[alternative HTML version deleted]]

Hi, On Fri, Oct 25, 2013 at 10:08 AM, Alicia R. P?rez-Porro <alicia.r.perezporro at="" gmail.com=""> wrote: > Dear Michael and DESeq2 users, > > I started using DESeq2 three weeks ago. Even I was having some normal > problems for a newby person my experience was really positive and I was and > I am detecting more DEG than with DESeq(1). > > However suddenly I am encountering a new problem that wasn't coming up > previously: > >> resSig <- res[ res$padj < .1, ] > Error in normalizeSingleBracketSubscript(i, x, byrow = TRUE, exact = FALSE) > : > subscript contains NAs resSig is a DataFrame, right? Does changing your code to this do the trick: R> resSig <- subset(res, padj < 0.1) Using `subset` instead of `[ ... ]` handles the case when there are NA's in the column you are filtering against. Otherwise you could do: R> resSig <- res[!is.na(res$padj) & res$padj < 0.1, ] > For me this is pretty weird because I've been using the exactly same > RScript and datasets for the last 3 weeks and I wasn't having this error > before. Wierd, indeed ... perhaps you recently upgraded to the new bioc2.13 universe? -steve -- Steve Lianoglou Computational Biologist Bioinformatics and Computational Biology Genentech

hi Alicia, The error is new because you are now using the v1.2 of DESeq2 which came with the Bioconductor 2.13 release this month. Maybe the system you are working on automatically upgraded? The error you see is saying that you cannot subset the results DataFrame with a vector that contains NAs. In the newest release of DESeq2 we have included automatic independent filtering to optimize the number of genes which will have FDR less than a given threshold. The genes which are filtered out get an adjusted p-value of NA (this behavior can be turned off by setting independentFiltering=FALSE to the results() function). But your script should work if you instead use the line: resSig <- res[ which(res$padj < .1), ] Mike On Fri, Oct 25, 2013 at 1:08 PM, Alicia R. Pérez-Porro < alicia.r.perezporro@gmail.com> wrote: > Dear Michael and DESeq2 users, > > I started using DESeq2 three weeks ago. Even I was having some normal > problems for a newby person my experience was really positive and I was and > I am detecting more DEG than with DESeq(1). > > However suddenly I am encountering a new problem that wasn't coming up > previously: > > > resSig <- res[ res$padj < .1, ] > Error in normalizeSingleBracketSubscript(i, x, byrow = TRUE, exact = > FALSE) : > subscript contains NAs > > For me this is pretty weird because I've been using the exactly same > RScript and datasets for the last 3 weeks and I wasn't having this error > before. > > My RScript is: > > #Count matrix input > Cele_MRvsER_old = read.csv (file.choose(), header=TRUE, row.names=1) > CeleDesign <- data.frame( > row.names = colnames(Cele_MRvsER_old), > condition = factor(c("MR", "MR", "ER", "ER", "ER", "ER"))) > dds <- DESeqDataSetFromMatrix(countData = Cele_MRvsER_old, > colData = CeleDesign, > design = ~ condition) > dds > > #Est size factor = normalize for library size > dds<- estimateSizeFactors(dds) > ddsLocal <- estimateDispersions(dds, fitType="local") > ddsLocal <- nbinomWaldTest(ddsLocal) > #plot dispersion > plotDispEsts(ddsLocal) > > > #Differential expression analysis > resultsNames(ddsLocal) > res <- results(ddsLocal, name= "condition_MR_vs_ER") > res <- res[order(res$padj),] > head(res) > plotMA(ddsLocal) > sum(res$padj < .1, na.rm=TRUE) > > > #filter for upregulated and downregulated genes > resSig <- res[ res$padj < .1, ] > head(resSig) > write.table(as.data.frame(resSig [ order(resSig$padj), ]), > file="BRvsMR_old_MostSigGenes.txt") > write.csv(as.data.frame(resSig[ order( resSig$log2FoldChange, > -resSig$baseMean ), ]), > file="BRvsMR_old_DownRegulated.csv") > write.csv(as.data.frame(resSig[ order( -resSig$log2FoldChange, > -resSig$baseMean ),]), > file="BRvsMR_old_UpRegulated.csv") > > #up-regulated > head(resSig[order(resSig$log2FoldChange,-resSig$baseMean),]) > #down-regulated > head(resSig[order(-resSig$log2FoldChange,-resSig$baseMean),]) > > #list most sig genes > head( resSig[ order(resSig$padj), ] ) > > And my sessionInfo(): > > > sessionInfo() > R version 3.0.2 (2013-09-25) > Platform: x86_64-apple-darwin10.8.0 (64-bit) > > locale: > [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] parallel stats graphics grDevices utils datasets methods > base > > other attached packages: > [1] DESeq2_1.2.0 RcppArmadillo_0.3.920.1 Rcpp_0.10.5 > GenomicRanges_1.14.1 > [5] XVector_0.2.0 IRanges_1.20.0 BiocGenerics_0.8.0 > > loaded via a namespace (and not attached): > [1] annotate_1.40.0 AnnotationDbi_1.24.0 Biobase_2.22.0 > DBI_0.2-7 > [5] genefilter_1.44.0 grid_3.0.2 lattice_0.20-24 > locfit_1.5-9.1 > [9] RColorBrewer_1.0-5 RSQLite_0.11.4 splines_3.0.2 > stats4_3.0.2 > [13] survival_2.37-4 tools_3.0.2 XML_3.95-0.2 > xtable_1.7-1* * > > > Any idea about why am I having this error now but not in the past and > about how to fix it? > > Thank you!!! > [[alternative HTML version deleted]]