Question

Principal components in DESeq2

1

Entering edit mode

Mozart ▴ 30

@mozart-20625

Last seen 4.1 years ago

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 • 2.6k views

ADD COMMENT • link 5.5 years ago Mozart ▴ 30

score 2 · Accepted Answer · 2019-05-31

2

Entering edit mode

swbarnes2 ★ 1.4k

@swbarnes2-14086

Last seen 9 minutes ago

San Diego

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 COMMENT • link 5.5 years ago swbarnes2 ★ 1.4k

0

Entering edit mode

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

ADD REPLY • link 5.5 years ago Mozart ▴ 30

2

Entering edit mode

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 REPLY • link 5.5 years ago Kevin Blighe ★ 4.0k

0

Entering edit mode

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 REPLY • link 5.5 years ago Mozart ▴ 30

0

Entering edit mode

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

ADD REPLY • link 5.5 years ago Kevin Blighe ★ 4.0k