Question: Principal components in DESeq2
1
6 months 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?

deseq2 • 210 views
modified 6 months ago • written 6 months ago by Mozart10
Answer: C: Principal components in DESeq2
2
6 months ago by
swbarnes2340
swbarnes2340 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 months ago by swbarnes2340 So, you can confirm that apart from the 'intgroup.df' everything is fine? ADD REPLYlink written 6 months 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 months ago by Kevin Blighe380 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!

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