This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Error in `[.data.frame`(res.df, , cols) : undefined columns selected when exporting glimma

Hi,

I am trying to export my RNAseq results with glimma and I am following this workshop (https://bioconductor.github.io/BiocWorkshops/rna-seq-data-analysis-with-deseq2.html#building-reports) unfortuneatly I am geting an error as per the title.

Please would someone be able to educate me as to why this is occuring? My code below:

#Glimma
library(org.Hs.eg.db)
res.s$symbol <- mapIds(org.Hs.eg.db,keys = rownames(res.s),column = "SYMBOL",keytype="ENSEMBL")
table(rownames(res.s) %in% keys(org.Hs.eg.db, "ENSEMBL"))

library("Glimma")
tmp <- "Z:/R.DESeq2"
status <- as.numeric(res.s$padj < .1)
anno <- data.frame(GeneID=rownames(res.s), symbol=res.s$symbol)
glMDPlot(res.sLFC, status=status, counts=counts(dds.s,normalized=TRUE),
         groups=dds.s$comparison, transform=FALSE,
         samples=colnames(dds.s), anno=anno,
         path=tmp, folder="glimma", launch=FALSE)

Error in [.data.frame(res.df, , cols) : undefined columns selected

Many thanks in advance Nathan

rna-seq
• 5,751 views • link
written by n.tear • n.tear

Hello Nathan. What is the output of:

str(res.sLFC)
colnames(res.sLFC)
colnames(dds.s)

?

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin thank you for your reply

str(res.sLFC)
    Formal class 'DESeqResults' [package "DESeq2"] with 7 slots
      ..@ priorInfo      :List of 4
      .. ..$ type         : chr "apeglm"
      .. ..$ package      : chr "apeglm"
      .. ..$ version      :Classes 'package_version', 'numeric_version'  hidden list of 1
      .. .. ..$ : int [1:3] 1 10 0
      .. ..$ prior.control:List of 7
      .. .. ..$ no.shrink            : int [1:2] 1 2
      .. .. ..$ prior.mean           : num 0
      .. .. ..$ prior.scale          : num 0.865
      .. .. ..$ prior.df             : num 1
      .. .. ..$ prior.no.shrink.mean : num 0
      .. .. ..$ prior.no.shrink.scale: num 15
      .. .. ..$ prior.var            : num 0.749
      ..@ rownames       : chr [1:31172] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
      ..@ nrows          : int 31172
      ..@ listData       :List of 4
      .. ..$ baseMean      : num [1:31172] 30899.28 2.93 2457.81 1310.71 417.27 ...
      .. ..$ log2FoldChange: Named num [1:31172] -2.98 -3.053 0.289 -0.931 0.92 ...
      .. .. ..- attr(*, "names")= chr [1:31172] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
      .. ..$ lfcSE         : Named num [1:31172] 0.314 1.084 0.165 0.183 0.334 ...
      .. .. ..- attr(*, "names")= chr [1:31172] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
      .. ..$ svalue        : num [1:31172] 0.00011 0.03615 0.88092 0.86355 0.81761 ...
      ..@ elementType    : chr "ANY"
      ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
      .. .. ..@ rownames       : NULL
      .. .. ..@ nrows          : int 4
      .. .. ..@ listData       :List of 2
      .. .. .. ..$ type       : chr [1:4] "intermediate" "results" "results" "results"
      .. .. .. ..$ description: chr [1:4] "
      .. .. ..@ elementType    : chr "ANY"
      .. .. ..@ elementMetadata: NULL
      .. .. ..@ metadata       : list()
      ..@ metadata       :List of 6
      .. ..$ filterThreshold: Named num 0.986
      .. .. ..- attr(*, "names")= chr "11.63265%"
      .. ..$ filterTheta    : num 0.116
      .. ..$ filterNumRej   :'data.frame':  50 obs. of  2 variables:
      .. .. ..$ theta : num [1:50] 0 0.0194 0.0388 0.0582 0.0776 ...
      .. .. ..$ numRej: num [1:50] 8324 8372 8443 8492 8551 ...
      .. ..$ lo.fit         :List of 2
      .. .. ..$ x: num [1:50] 0 0.0194 0.0388 0.0582 0.0776 ...
      .. .. ..$ y: num [1:50] 8322 8379 8436 8493 8551 ...
      .. ..$ alpha          : num 0.1
      .. ..$ lfcThreshold   : num 2
    > colnames(res.sLFC)
    [1] "baseMean"       "log2FoldChange" "lfcSE"          "svalue"        
    > colnames(dds.s)
     [1] "C18"     "C54"     "C56"     "C59"     "C23"     "C27"     "C38"
• 0 views • link
written by n.tear • n.tear

Should res.sLFC not just be a data-frame? - and you will have to specify the x and y values via xval and yval, respectively.

• 0 views • link
written by Kevin Blighe Kevin Blighe

I generated the res.sLFC using:

res.sLFC <- lfcShrink(dds.s, coef="co", type="apeglm", lfcThreshold=2)

As a data frame:

> res.sLFCdf <- as.data.frame(res.sLFC)
> res.sLFCdf
                    baseMean log2FoldChange      lfcSE       svalue
ENSG00000000003 3.089928e+04   -2.980109835 0.31391008 1.098978e-04
ENSG00000000005 2.927032e+00   -3.052505162 1.08427276 3.614998e-02
ENSG00000000419 2.457810e+03    0.288835065 0.16471434 8.809206e-01
ENSG00000000457 1.310707e+03   -0.930850511 0.18271970 8.635484e-01
ENSG00000000460 4.172665e+02    0.919899394 0.33440589 8.176108e-01

Now I have changed toa data frame the error is: 

>glMDPlot(res.sLFCdf, status=status, counts=counts(dds.s,normalized=TRUE),
+          groups=dds.s$comparison, transform=FALSE,
+          samples=colnames(dds.s), anno=anno,
+          path=tmp, folder="glimma", launch=FALSE)
Error in if (isDefaultStates(x) && isDefaultCols(cols)) { : 
  missing value where TRUE/FALSE needed
• 0 views • link
written by n.tear • n.tear

I see. You will have to also set the following:

glMDPlot(..., xval <- 'baseMean', yval <- 'log2FoldChange', ...)

In fact, you may even have to create a new column for the logged base mean, and then use that:

res.sLFCdf$logmean <- log(res.sLFCdf$baseMean + 1)
glMDPlot(..., xval <- 'logmean', yval <- 'log2FoldChange', ...)
• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin, many thanks for your continued support.

Unfortunately I am still having no luck

> res.sLFCdf$logmean <- log(res.sLFCdf$baseMean + 1)
> glMDPlot(res.sLFCdf, xval <- 'logmean', yval <- 'log2FoldChange', 
+          status=status, counts=counts(dds.s,normalized=TRUE),
+          groups=dds.s$comparison, transform=FALSE,
+          samples=colnames(dds.s), anno=anno,
+          path=tmp, folder="glimma", launch=TRUE)
Error in if (isDefaultStates(x) && isDefaultCols(cols)) { : 
  missing value where TRUE/FALSE needed
• 0 views • link
written by n.tear • n.tear

Seems that there are missing or infinite values in your logmean variable. Can check the output of this:

table( is.na(res.sLFCdf$logmean) )

table( res.sLFCdf$logmean == 'Inf' )
• 0 views • link
written by Kevin Blighe Kevin Blighe

output suggests not 

> table( is.na(res.sLFCdf$logmean) )

FALSE 
31172 
> table( res.sLFCdf$logmean == 'Inf' )

FALSE 
31172
• 0 views • link
written by n.tear • n.tear

Can you literally check your input data, in both formatting (string, numeric, factor, etc) and also for other missing values? Like, literally, look through it or output it to disk and open it in Excel.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Yes looked in exel and no missing values and checked the class of each column in res.sLFCdf and they are correct

could it be an issue that not all the rownames of res.s are in org.Hs.eg.db when I make the anno?

> library(org.Hs.eg.db)
> res.s$symbol <- mapIds(org.Hs.eg.db,keys = rownames(res.s),column = "SYMBOL",keytype="ENSEMBL")
'select()' returned 1:many mapping between keys and
columns
> table(rownames(res.s) %in% keys(org.Hs.eg.db, "ENSEMBL"))

FALSE  TRUE 
10649 20523

I also have some NAs in my status object. could this be a problem?

status <- as.numeric(res.s$padj < .1)
status
   [1]  1  1  0  1  1  0  1  0  1  0  0  1  0  1  1  1  0
  [18]  0  0  0  0  0  0  1  1  0  1  0  0  1  0  0  0  1
  [35]  1  1  0  1  0  1  1  0  0  1  1  0  1  1  0  0  0
  [52]  0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  1  0
  [69]  0 NA  1  0  1  0  1  0  0 NA  0  0  0  0  0  0  0
• 0 views • link
written by n.tear • n.tear

Did you manage to get any further with this? I am limited in what I can do from a distance. Best idea would be to find the last known working example, and then work from there

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin,

Thanks for your help. Ive done exacly as you have said above and managed to work it out in the end.

Cheers

Nathan

• 0 views • link
written by n.tear • n.tear

0 answers

No answers yet.

Log in to answer this question.

More posts like this