Question: How do I get the weights of genes used in PCA in DESeq2?
1
gravatar for cristian
2.5 years ago by
cristian230
cristian230 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.

loadings pca deseq2 gene list • 3.1k views
ADD COMMENTlink modified 2.5 years ago by igor8.1k • written 2.5 years ago by cristian230
2
gravatar for igor
2.5 years ago by
igor8.1k
United States
igor8.1k 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)
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

ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by igor8.1k

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)

ADD REPLYlink modified 2.5 years ago • written 2.5 years ago by cristian230

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?

ADD REPLYlink written 2.5 years ago by cristian230

Got it!

assay() is in the package SummarizedExperiment.

ADD REPLYlink written 2.5 years ago by cristian230

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

ADD REPLYlink written 2.5 years ago by igor8.1k

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)

}

ADD REPLYlink modified 2.5 years ago • written 2.5 years ago by cristian230
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1687 users visited in the last hour