Hi, I am trying to help a labmate who is using DiffBind to do a PCA analysis on some ATAC-seq data that she generated. I have experience with PCA but not with this specific software. How can we retrieve the PCA loadings/eigenvalues for each individual for each PC? I am used to using other PCA tools that allow you to retried the spreadsheet of loadings/eigenvalues in addition to plotting PC's 1-10 against each other. As far as I can tell DiffBind allows you to plot PC1, PC2, and PC3 against each other but does not provide the actual numbers, nor does it provide access to any further PC's. Any help with this would be greatly appreciated.

dba.PCA documentation: https://www.rdocumentation.org/packages/DiffBind/versions/2.0.2/topics/dba.plotPCA

Thanks!

DiffBind does allow any components to be plotted again any other components by specifying which components to plot using components= when calling dba.plotPCA() (the default is components=1:3 but this may be changed by the user).

There is no way to directly retrieve the computed PCA data. You can retrieve the underlying read count data (normalized by default) using dba.peakset() with bRetrieve=TRUE and run a PCA on that.

What version are you using? That information is somewhat dated.

Since version 2.14, the default logic used is as follows:

• Normalized read counts for the samples being plotted are retrieved
• Any peaks with no non-zero values for these samples are removed
• The log2() values for the read counts are computed after setting values < 1 to 1
• prcomp()(not princomp()) is used to compute a principal components analysis on the transpose of these values
• The percentage variance for each component is computed by squaring the $sdev values, dividing by the total, and multiplying by 100.
• The $x values are plotted. These are the loading values that have been centered and multiplied by the rotation matrix.

Here's some code to show the essence of what is happening:

data(tamoxifen_counts)
dba.plotPCA(tamoxifen, DBA_FACTOR)

peaks <- dba.peakset(tamoxifen,bRetrieve=TRUE)
vals <- as.matrix(mcols(peaks))
vals[vals<1] <- 1
vals <- log2(vals)
pc <- prcomp(t(vals))
percent <- pc$sdev[1:2]^2 / sum(pc$sdev^2) * 100
percent <- as.character(round(percent))

lattice::xyplot(
  PC2 ~ PC1,
  data = as.data.frame(pc$x[,1:2]),
  xlab = percent[1], ylab = percent[2],
  pch = 16, cex = 1.75, aspect = 1)