Trouble with Custom PC Plot
Entering edit mode
newsomew13 • 0
Last seen 6.4 years ago
United States


I'm attempting to produce a PC plot which includes PC3, PC4, etc. I have modified the existing plotPCA function but I keep getting the following errors:

Error in plotPCA(rld, intgroup = c("sex")) :
  no method for coercing this S4 class to a vector
In addition: Warning message:
In all(intgroup %in% "sex", (colData(rld))) :
  coercing argument of type 'S4' to logical

I find this error message confusing, as plotPCA works perfectly fine when I run plotPCA unmodified. Would anyone mind checking to see if there are any obvious errors in my code?


plotPCA <- function (rld, intgroup = "sex", ntop = 500, returnData = FALSE)
  rv <- rowVars(assay(rld))
  select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]
  pca <- prcomp(t(assay(rld)[select,]))
  percentVar <- pca$sdev^2 / sum( pca$sdev^2 )
  if (!all(intgroup %in% 'sex', (colData(rld)))) {
    stop("the argument 'intgroup' should specify columns of colData(dds)")
  intgroup.df <-[, intgroup, drop=FALSE])
  group <- if (length(intgroup) > 1) {
    factor(apply( intgroup.df, 1, paste, collapse=" : "))
  } else {
  d <- data.frame(PC3 = pca$x[, 3], PC4 = pca$x[, 4], group=group, intgroup.df, name=colnames(rld))

  if (returnData) {
    attr(d, "percentVar") <- percentVar[1:2]
  ggplot(data=d, aes_string(x="PC3", y="PC4", color="group")) + geom_point(size=3) +
    xlab(paste0("PC3: ",round(percentVar[1] * 100),"% variance")) +
    ylab(paste0("PC4: ",round(percentVar[2] * 100),"% variance"))

print(plotPCA(rld, intgroup=c("sex")))

Session info:

R version 3.3.0 (2016-05-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 10586)

[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

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

other attached packages:
 [1] matrixStats_0.50.2         ggplot2_2.1.0             
 [3] DESeq2_1.12.2              SummarizedExperiment_1.2.2
 [5] Biobase_2.32.0             GenomicRanges_1.24.0      
 [7] GenomeInfoDb_1.8.2         IRanges_2.6.0             
 [9] S4Vectors_0.10.0           BiocGenerics_0.18.0       
[11] BiocInstaller_1.22.2      

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.5          RColorBrewer_1.1-2   plyr_1.8.3          
 [4] XVector_0.12.0       tools_3.3.0          zlibbioc_1.18.0     
 [7] digest_0.6.9         rpart_4.1-10         annotate_1.50.0     
[10] RSQLite_1.0.0        gtable_0.2.0         lattice_0.20-33     
[13] Matrix_1.2-6         DBI_0.4-1            gridExtra_2.2.1     
[16] genefilter_1.54.2    cluster_2.0.4        locfit_1.5-9.1      
[19] grid_3.3.0           nnet_7.3-12          data.table_1.9.6    
[22] AnnotationDbi_1.34.2 XML_3.98-1.4         survival_2.39-4     
[25] BiocParallel_1.6.2   foreign_0.8-66       latticeExtra_0.6-28
[28] Formula_1.2-1        geneplotter_1.50.0   Hmisc_3.17-4        
[31] scales_0.4.0         splines_3.3.0        xtable_1.8-2        
[34] colorspace_1.2-6     labeling_0.3         acepack_1.3-3.3     
[37] munsell_0.4.3        chron_2.3-47  

deseq2 plotPCA PCA rna-seq • 1.0k views
Entering edit mode
Last seen 1 day ago
United States

This line looks malformed

 if (!all(intgroup %in% 'sex', (colData(rld)))) {

What do you want to do here?

Entering edit mode
newsomew13 • 0
Last seen 6.4 years ago
United States

Honestly, I'm not really sure. The source file shows this as a template:

if (!all(intgroup %in% names(colData(object)))) {

which I didn't really know what to do with.

I figured I should try "intgroup %in% 'sex'" or just "sex", but either way it produces the same error.

Entering edit mode

If you're new to R, and customizing a function for your own purposes, I'd only change the necessary parts.

Note that the % of variance annotation on your x and y axes will be incorrect because you are using percentVar[1] and percentVar[2]


Login before adding your answer.

Traffic: 234 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6