Hi there,

I have run DeSeq2 in R with my RNA-expression data (no problem at all with this) and I have create a PC plot as one of the results of these analyses. The data looks good and they seem to be grouped quite clearly. But now I wonder whether there is a way os statistically test that the PC component (axis values) are statistically different.

I know that with "normal" PCA one can run "constrained rda analyses" by using the package vegan but I am not sure whether there is something similar for PC plot creating with DeSeq2.

Here is the code:

library(DESeq2) 

## Load file cts
## Load file coldata

## Run the analyses

dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition")

##  Transformation function
vsd <- varianceStabilizingTransformation(dds, blind = TRUE)

## Principal component plot of the samples
plotPCA(vsd, intgroup=c("condition","moment"))  ## visualizing the overall effect of experimental covariates and batch effects

## Customize the PCA plot using the ggplot function.
pcaData <- plotPCA(vsd, intgroup=c("condition","moment"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar.exp"))

## Do the plot
ggplot(pcaData, aes(PC1, PC2, color=condition, shape=moment)) +
  geom_point(size=5) +
  xlab(paste0("PC1: 59% variance")) +
  ylab(paste0("PC2: 13% variance")) 

Thank you so much in advance