plotPCA in DESeq2
1
0
Entering edit mode
yueli7 ▴ 10
@yueli7-8401
Last seen 3 months ago
China

Hello, I used plotPCA function in DESeq2.

Thanks in advance for great help!

Best,

Yue

> plotPCA(vsd, intgroup=c("condition"))
Error in (function (classes, fdef, mtable)  : 
unable to find an inherited method for function 'exprs' for signature '"DESeqTransform"'

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
[9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] affycoretools_1.58.4        RColorBrewer_1.1-2          pheatmap_1.0.12            
[4] vsn_3.54.0                  ggplot2_3.3.1               DESeq_1.38.0               
[7] lattice_0.20-41             locfit_1.5-9.4              data.table_1.12.8          
[10] IHW_1.14.0                  airway_1.6.0                pasilla_1.14.0             
[13] tximeta_1.4.5               DESeq2_1.26.0               SummarizedExperiment_1.16.1
[16] DelayedArray_0.12.3         BiocParallel_1.20.1         matrixStats_0.56.0         
[19] Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.1        
[22] IRanges_2.20.2              S4Vectors_0.24.4            BiocGenerics_0.32.0        
[25] tximportData_1.14.0         readr_1.3.1                 tximport_1.14.2            
[28] limma_3.42.2               

loaded via a namespace (and not attached):
[1] R.utils_2.9.2            tidyselect_1.1.0         RSQLite_2.2.0            AnnotationDbi_1.48.0    
[5] htmlwidgets_1.5.1        grid_3.6.1               munsell_0.5.0            codetools_0.2-16        
[9] preprocessCore_1.48.0    withr_2.2.0              colorspace_1.4-1         Category_2.52.1         
[13] OrganismDbi_1.28.0       knitr_1.28               rstudioapi_0.11          labeling_0.3            
[17] slam_0.1-47              bbmle_1.0.23.1           GenomeInfoDbData_1.2.2   lpsymphony_1.14.0       
[21] mixsqp_0.3-43            hwriter_1.3.2            bit64_0.9-7              farver_2.0.3            
[25] coda_0.19-3              vctrs_0.3.1              generics_0.0.2           xfun_0.14               
[29] biovizBase_1.34.1        BiocFileCache_1.10.2     R6_2.4.1                 Glimma_1.14.0           
[33] apeglm_1.8.0             invgamma_1.1             AnnotationFilter_1.10.0  bitops_1.0-6            
[37] reshape_0.8.8            assertthat_0.2.1         scales_1.1.1             nnet_7.3-14             
[41] gtable_0.3.0             affy_1.64.0              ggbio_1.34.0             ensembldb_2.10.2        
[45] rlang_0.4.6              genefilter_1.68.0        splines_3.6.1            rtracklayer_1.46.0      
[49] lazyeval_0.2.2           acepack_1.4.1            dichromat_2.0-0          hexbin_1.28.1           
[53] checkmate_2.0.0          BiocManager_1.30.10      reshape2_1.4.4           GenomicFeatures_1.38.2  
[57] backports_1.1.7          Hmisc_4.4-0              RBGL_1.62.1              tools_3.6.1             
[61] affyio_1.56.0            ellipsis_0.3.1           gplots_3.0.3             ff_2.2-14.2             
[65] Rcpp_1.0.4.6             plyr_1.8.6               base64enc_0.1-3          progress_1.2.2          
[69] zlibbioc_1.32.0          purrr_0.3.4              RCurl_1.98-1.2           prettyunits_1.1.1       
[73] rpart_4.1-15             openssl_1.4.1            ashr_2.2-47              cluster_2.1.0           
[77] magrittr_1.5             truncnorm_1.0-8          mvtnorm_1.1-1            SQUAREM_2020.3          
[81] ProtGenerics_1.18.0      hms_0.5.3                xtable_1.8-4             XML_3.99-0.3            
[85] emdbook_1.3.12           jpeg_0.1-8.1             gcrma_2.58.0             gridExtra_2.3           
[89] compiler_3.6.1           biomaRt_2.42.1           bdsmatrix_1.3-4          tibble_3.0.1            
[93] KernSmooth_2.23-17       crayon_1.3.4             ReportingTools_2.26.0    R.oo_1.23.0             
[97] htmltools_0.5.0          GOstats_2.52.0           Formula_1.2-3            geneplotter_1.64.0      
[101] DBI_1.1.0                dbplyr_1.4.4             MASS_7.3-51.6            rappdirs_0.3.1          
[105] Matrix_1.2-17            R.methodsS3_1.8.0        gdata_2.18.0             pkgconfig_2.0.3         
[109] GenomicAlignments_1.22.1 numDeriv_2016.8-1.1      foreign_0.8-76           foreach_1.5.0           
[113] annotate_1.64.0          XVector_0.26.0           AnnotationForge_1.28.0   stringr_1.4.0           
[117] VariantAnnotation_1.32.0 digest_0.6.25            graph_1.64.0             Biostrings_2.54.0       
[121] htmlTable_1.13.3         edgeR_3.28.1             GSEABase_1.48.0          curl_4.3                
[125] Rsamtools_2.2.3          gtools_3.8.2             lifecycle_0.2.0          jsonlite_1.6.1          
[129] PFAM.db_3.10.0           askpass_1.1              BSgenome_1.54.0          pillar_1.4.4            
[133] GGally_2.0.0             httr_1.4.1               survival_3.2-3           GO.db_3.10.0            
[137] glue_1.4.1               fdrtool_1.2.15           iterators_1.0.12         png_0.1-7               
[141] bit_1.1-15.2             Rgraphviz_2.30.0         stringi_1.4.6            blob_1.2.1              
[145] oligoClasses_1.48.0      latticeExtra_0.6-29      caTools_1.18.0           memoise_1.1.0           
[149] dplyr_1.0.0              irlba_2.3.3          


> sessionInfo
function (package = NULL) 
{
z <- list()
z$R.version <- R.Version()
z$platform <- z$R.version$platform
if (nzchar(.Platform$r_arch)) 
    z$platform <- paste(z$platform, .Platform$r_arch, sep = "/")
z$platform <- paste0(z$platform, " (", 8 * .Machine$sizeof.pointer, 
    "-bit)")
z$locale <- Sys.getlocale()
z$running <- osVersion
z$RNGkind <- RNGkind()
if (is.null(package)) {
    package <- grep("^package:", search(), value = TRUE)
    keep <- vapply(package, function(x) x == "package:base" || 
        !is.null(attr(as.environment(x), "path")), NA)
    package <- .rmpkg(package[keep])
}
pkgDesc <- lapply(package, packageDescription, encoding = NA)
if (length(package) == 0) 
    stop("no valid packages were specified")
basePkgs <- sapply(pkgDesc, function(x) !is.null(x$Priority) && 
    x$Priority == "base")
z$basePkgs <- package[basePkgs]
if (any(!basePkgs)) {
    z$otherPkgs <- pkgDesc[!basePkgs]
    names(z$otherPkgs) <- package[!basePkgs]
}
loadedOnly <- loadedNamespaces()
loadedOnly <- loadedOnly[!(loadedOnly %in% package)]
if (length(loadedOnly)) {
    names(loadedOnly) <- loadedOnly
    pkgDesc <- c(pkgDesc, lapply(loadedOnly, packageDescription))
    z$loadedOnly <- pkgDesc[loadedOnly]
}
z$matprod <- as.character(options("matprod"))
es <- extSoftVersion()
z$BLAS <- as.character(es["BLAS"])
z$LAPACK <- La_library()
class(z) <- "sessionInfo"
z
}
<bytecode: 0x55cc2a8c1198>
<environment: namespace:utils>   

> plotPCA
nonstandardGenericFunction for "plotPCA" defined from package "BiocGenerics"

function (object, ...) 
{
standardGeneric("plotPCA")
}
<bytecode: 0x55cbfa16ddc0>
<environment: 0x55cbfa166f00>
Methods may be defined for arguments: object
Use  showMethods("plotPCA")  for currently available ones.
software error deseq2 • 517 views
ADD COMMENT
0
Entering edit mode

vsd is supposed to be a DESeqTransform object. From the manual:

a DESeqTransform object, with data in assay(x), produced for example by either rlog or varianceStabilizingTransformation.

So, something like:

vsd <- varianceStabilizingTransformation(dds)
plotPCA(vsd, intgroup=c("condition"))
ADD REPLY
0
Entering edit mode

I guess the function is masked by any package for microarray analysis because it tries to call exprs() on vsd.

ADD REPLY
0
Entering edit mode

Could be - moved back to comment. plotPCA is overloaded between DESeq2 and affycoretools

ADD REPLY
0
Entering edit mode

One note (more for developers than users) is that if a method masks plotPCA or other function names that are listed as generics here, they can avoid this by exporting a method for those generics.

ADD REPLY
0
Entering edit mode

Thank you all of yours' great help! Best,

Yue

> rld<-rlogTransformation(dds2)
rlog() may take a few minutes with 30 or more samples,
vst() is a much faster transformation
> plotPCA(rld, intgroup=c("condition"))
ADD REPLY
0
Entering edit mode

Try to avoid loading DESeq. This version is official deprecated.

ADD REPLY
0
Entering edit mode

Hello Kevin Blighe, Thank you so much for your great help! Best, Yue

> vsd <- varianceStabilizingTransformation(dds)
Error in getVarianceStabilizedData(cds) : 
is(cds, "CountDataSet") is not TRUE
ADD REPLY
0
Entering edit mode

Hello Kevin Blighe, Thank you so much for your great help! Best, Yue

ADD REPLY
0
Entering edit mode
ATpoint ▴ 850
@atpoint-13662
Last seen 11 hours ago
Germany

Try DESeq2::plotPCA(vsd, intgroup=c("condition"))

ADD COMMENT
0
Entering edit mode

It's not necessary to specify an S4 function using the double colon, so long as there isn't an object dispatch clash (which there isn't in the case of DESeq2 and affycoretools).

> library(DESeq2)
> library(affycoretools)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

> showMethods(plotPCA)
Function: plotPCA (package BiocGenerics)
object="DESeqTransform"
object="ExpressionSet"
object="matrix"

> example(varianceStabilizingTransformation)

vrncST> dds <- makeExampleDESeqDataSet(m=6)

vrncST> vsd <- varianceStabilizingTransformation(dds)

vrncST> dists <- dist(t(assay(vsd)))

vrncST> # plot(hclust(dists))
vrncST> 
vrncST> 
vrncST> 
vrncST> 
> plotPCA(vsd)
> 

ADD REPLY
0
Entering edit mode

I did not say it was that particular package, but there is a conflict here since plotPCA (from whatever package it comes) tries to call exprs() on the input object, right?

ADD REPLY
0
Entering edit mode

True. You only said it was a microarray analysis package. Kevin Blighe said it was affycoretools that overloaded DESeq2, which isn't a thing, since both packages use the method defined by BiocGenerics and export correctly.

It's either another package that has a bare or S3 function, or perhaps the OP wrote their own? Anyway, without the output from sessionInfo it's not possible to say.

ADD REPLY
0
Entering edit mode

Was the effect of responding in haste while trying to do all of my 1 million jobs. plotPCA() is actually used in many different packages, it seems. Won't know much more until user gets back (currently middle of the night in China).

ADD REPLY
1
Entering edit mode

Yes, and me going all 'Grrr, talking smack about my package, are ya!' probably without warrant. ;-D

ADD REPLY
0
Entering edit mode

Thank you all of yours' great help!

> rld<-rlogTransformation(dds2)
rlog() may take a few minutes with 30 or more samples,
vst() is a much faster transformation
> plotPCA(rld, intgroup=c("condition"))
ADD REPLY
0
Entering edit mode

You need to supply the results from

sessionInfo()

What you have supplied is the function body, which you get by doing


sessionInfo

Which is a different thing.

ADD REPLY
0
Entering edit mode

You could also do

plotPCA


Which will give the function body, and where it comes from, which would be pretty helpful as well.

ADD REPLY
0
Entering edit mode

Hello James W. MacDonald, Thank you so much for your great help! I have already added these information.

Thank you again! Best, Yue

ADD REPLY
0
Entering edit mode

You seem to have too many packages open in the same workspace. Do you really need all of these in the same workflow:

affycoretools              DESeq       
tximeta     DESeq2
tximportData  tximport
limma

You should close your workspace when not using it, and only load packages that you need. For example, there is no major reason why you should have DESeq and DESeq2 loaded together (DESeq is also being deprecated, as per Michael's comment above).

Be selective in the packages that you load, and close your workspace when not using it. One can run into all sorts of problems by leaving packages loaded and old variables in the same workspace.

ADD REPLY
1
Entering edit mode

Hello, Kevin Blighe, Thank you so much for your great help!

Best, Yue

ADD REPLY

Login before adding your answer.

Traffic: 242 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