Trouble with Custom PC Plot
2
0
Entering edit mode
newsomew13 • 0
@newsomew13-9225
Last seen 7.9 years ago
United States

Hello,

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?

Thanks!
 

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 <- as.data.frame(colData(rld)[, intgroup, drop=FALSE])
  group <- if (length(intgroup) > 1) {
    factor(apply( intgroup.df, 1, paste, collapse=" : "))
  } else {
    colData(rld)[[intgroup]]
  }
  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]
    return(d)
  }
 
  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)

locale:
[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.5k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 8 hours ago
United States

This line looks malformed

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

What do you want to do here?

ADD COMMENT
0
Entering edit mode
newsomew13 • 0
@newsomew13-9225
Last seen 7.9 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.

ADD COMMENT
0
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]

ADD REPLY

Login before adding your answer.

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