Question: Principal components in DESeq2
1
gravatar for Mozart
6 weeks ago by
Mozart10
Mozart10 wrote:

Hello there, after reading a few posts here about this particular topic, I guess I solved this 'puzzling' issue and I just need to check that this is the correct way to do this. Essentially, I need to output the PC1 and PC2 from plotPCA function in order to let someone else to play around with an external software these data.

rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(500,length(rv)))]
pca <- prcomp(t(assay(rld)[select, ]))
intgroup.df <- as.data.frame(colData(rld)[, intgroup, drop = FALSE])
pca$x[,1]
pca$x[,2]

Is this the correct way to proceed?

Thanks in advance.

deseq2 • 113 views
ADD COMMENTlink modified 6 weeks ago • written 6 weeks ago by Mozart10
Answer: C: Principal components in DESeq2
2
gravatar for swbarnes2
6 weeks ago by
swbarnes2210
swbarnes2210 wrote:

intgroup.df isn't doing anything there. I like to throw in pca_merged <- merge(colData(dds), pca$x, by.x = 0, by.y = 0) to get all the info in one place.

ADD COMMENTlink written 6 weeks ago by swbarnes2210

So, you can confirm that apart from the 'intgroup.df' everything is fine?

ADD REPLYlink written 6 weeks ago by Mozart10
2

That is most likely what your colleagues want - yes. pca$x contains the rotated component loadings. In the code that you have pasted, though, you are also pre-selecting the genes with highest variance prior to performing PCA with prcomp(); thus, you are somewhat biasing it.

You can also get the same data with my own PCAtools package:

p <- pca(x, metadata = NULL, removeVar = NULL)
p$rotated[,1:2]

You can remove variables based on variance prior to performing PCA via the removeVar parameter

ADD REPLYlink written 6 weeks ago by Kevin Blighe170

Thanks Kevin, now I have understood why this DESeq2 plot is defined 'biased'. So, rather than doing the code above I could simply do in the following way:

rv <- rowVars(assay(rld))
p <- pca(x, metadata = NULL, removeVar = NULL)
p$rotated[,1:2]

Is this what you were referring to? Sorry, I am bit rusty!

ADD REPLYlink written 6 weeks ago by Mozart10

Yes that is the code, if you are using my PCAtools. You do not have to calculate rv, though.

ADD REPLYlink written 6 weeks ago by Kevin Blighe170
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 16.09
Traffic: 313 users visited in the last hour