Hi, I am trying out some packages for PCA in combination with the reads-to-genes edgeR-workflow for my RNA-seq analysis. I used the directions given from this post: Trying to do a PCA using edgeR with expression data

Group is equivalent to the Biomarkerquartiles_Q variable.

When I use the plotMDS function with the following code:

keep <- filterByExpr(y, group = y$samples$group)

y <- y[keep, , 
       keep.lib.sizes = F]

y <- calcNormFactors(y)

pch <- c(0,1,2,15,16,17)
colors <- rep(c("darkgreen", "red", "blue"), 2)
plotMDS(y, col=colors[factor(targets$Biomarkerquartiles_Q)],
        pch=pch[factor(targets$Biomarkerquartiles_Q)],
        gene.selection="common")

legend("topleft", legend=levels(factor(targets$Biomarkerquartiles)),
       pch=pch, col=colors, ncol=2)

I get this plot:

And although there is no clear clustering of the quartiles, it is pretty clear that there is some clustering going on.

When I use the PCAtools package with the following code:

p <- pca(cpm(y, log = T),
         metadata = y$samples,
         removeVar = 0.1)
PCA.plot <- biplot(p,
    lab = paste0(p$metadata$group),
    colby = "group"
)

I get the following plot:

I know both packages use different approaches but is it expected that the result can differ so much? Intriguingly, the x-axis and y-axis descriptions of both approaches seem to report comparable amounts of variations explained by the PCs.

Additionally, where could I begin to learn how to interpret the results from the plotMDS output?

The reason why the two plots are different is that the first plot is a PCA plot while the second is a biplot. Quite different plots.

The PCA plot from edgeR plotMDS is showing very clearly that your data has a massive batch effect, which is something that you need to investigate. The biplot doesn't seem to be showing any useful information.

Any batch effect should be included in the design matrix.

The results from plotMDS seem pretty clear. You can see examples of MDS plots and their interpretation in the limma and edgeR case studies. As the documentation explains, distances on the MDS plot correspond to log2 leading fold-changes.