How do I get the weights of genes used in PCA in DESeq2?

Question

[biomagician]

Dear all,

I have a matrix of reads counts for each gene for each strain of an animal. I am using DESeq2 to analyse the clustering of the strains. I have performed a Variance Stabilising Transform on my DESeq2 object and can do a PCA:

p <- DESeq2::plotPCA(vsd, intgroup=c('Condition'))

p <- p + labs(title = 'PCA on gene expression levels')

p

Now, I would like to know the coefficients/weights of each gene for each principal component. Does anyone know how to do this?

Thanks.

C.

Answer 1

[igor]

The plotPCA function performs PCA like this:

plotPCA.DESeqTransform = function(object, intgroup="condition", ntop=500, returnData=FALSE)
{
  # calculate the variance for each gene
  rv <- rowVars(assay(object))

  # select the ntop genes by variance
  select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]

  # perform a PCA on the data in assay(x) for the selected genes
  pca <- prcomp(t(assay(object)[select,]))

  ...

So object is whatever object you feed it, such as vsd. It doesn't actually store the full PCA results anywhere as far as I can tell, so you would have to repeat the PCA manually (using those three commands).

To get the loadings with prcomp:

pca <- prcomp(data)
loadings <- as.data.frame(pca$rotation)

More detailed instructions: http://stackoverflow.com/questions/12760108/principal-components-analysis-how-to-get-the-contribution-of-each-paramete

[biomagician]

Hi, Thanks so much. Do you know how I can get the matrix out of vsd and feed it to prcom()?

That's the code I used to get vsd:

dds <- DESeq2::DESeqDataSetFromHTSeqCount(sampleTable=sampleData, directory = paste0('output/counts/htseq/tophat/', sampleType, '/'), design=~Condition)

dds2 <- DESeq2::DESeq(dds)

vsd <- DESeq2::varianceStabilizingTransformation(dds2)

[biomagician]

I see that it should be the assay() function but when I use it, R does not find it:

assay(vsd)

Error: could not find function "assay"

It is not part of DESeq2 either:

DESeq2::assay(vsd) Error: 'assay' is not an exported object from 'namespace:DESeq2'

Weird, how can the plotPCA() function work then?

[biomagician]

Got it!

assay() is in the package SummarizedExperiment.

[igor]

Are you not loading DESeq2 with library(DESeq2)? That should set up the namespace properly.

[biomagician]

I do something different actually. I put all the names of the packages I need in a character vector and then supply that vector to this function:

ipak <- function(pkg){

new.pkg <- pkg[!(pkg %in% installed.packages()[, 'Package'])]

if (length(new.pkg))

install.packages(new.pkg, dependencies = TRUE)

sapply(pkg, require, character.only = TRUE)

}