Question: How do I get the weights of genes used in PCA in DESeq2?
1
cristian240 wrote:

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.

modified 3.1 years ago by igor9.9k • written 3.1 years ago by cristian240
2
igor9.9k wrote:

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)
``````

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)

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?

Got it!

assay() is in the package SummarizedExperiment.

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

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)

}