Principal components in DESeq2
1
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
2
Entering edit mode
swbarnes2 ★ 1.4k
@swbarnes2-14086
Last seen 20 hours 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
0
Entering edit mode

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

ADD REPLY
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
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
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

Login before adding your answer.

Traffic: 751 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6