Plot PCA for RNASeq at any principles as scatterplot. PC1-PC2, PC3-PC4 any as you wish.

plotPCA_anyPrinciples.R

# The function is the basically the same as https://github.com/mikelove/DESeq2/blob/master/R/plots.R.
# Inspired by https://www.biostars.org/p/243695/. 
# I added two variables, pp1 and pp2 to let user chose which principle to plot. pp1 and pp2 only take integer values.

# Last modified: 2021-12-03

plotPCA_jh = function(pp1=1, pp2=2, 
                      object, intgroup="condition", 
                      ntop=1000, returnData=FALSE) {
  
  
  # calculate the variance for each gene
  rv <- rowVars(assay(object))
  
  # select the ntop genes by variance
  select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]
  
  # perform a PCA on the data in assay(x) for the selected genes
  pca <- prcomp(t(assay(object)[select,]))
  
  # the contribution to the total variance for each component
  percentVar <- pca$sdev^2 / sum( pca$sdev^2 )
  
  if (!all(intgroup %in% names(colData(object)))) {
    stop("the argument 'intgroup' should specify columns of colData(dds)")
  }
  
  intgroup.df <- as.data.frame(colData(object)[, intgroup, drop=FALSE])
  
  # add the intgroup factors together to create a new grouping factor
  group <- if (length(intgroup) > 1) {
    factor(apply(intgroup.df, 1, paste, collapse=":"))
  } else {
    colData(object)[[intgroup]]
  }
  
  # assembly the data for the plot
  d <- data.frame(PC1=pca$x[,pp1], PC2=pca$x[,pp2], group=group, intgroup.df, 
                  name=colnames(object))
  
  if (returnData) {
    attr(d, "percentVar") <- percentVar[pp1:pp2]
    return(d)
  }
  
  ggplot(data=d, aes_string(x="PC1", y="PC2", color="group", label = "name")) + 
    geom_point(size=3) + 
    xlab(paste0("PC", pp1, ": ",round(percentVar[pp1] * 100),"% variance")) +
    ylab(paste0("PC", pp2, ": ",round(percentVar[pp2] * 100),"% variance")) +
    #scale_color_brewer(type = "qual", palette = "Dark2", direction = 1)+
    colorspace::scale_color_discrete_qualitative(palette = "Dark 3", rev = TRUE)+
    cowplot::theme_cowplot()+
    coord_fixed()
}