---
title: "PCAtools: everything Principal Component Analysis"
author: "Kevin Blighe, Aaron Lun"
date: "`r Sys.Date()`"
package: "`r packageVersion('PCAtools')`"
output:
  html_document:
    toc: true
    toc_depth: 3
    number_sections: true
    theme: united
    highlight: tango
fig_width: 7
bibliography: library.bib
vignette: >
    %\VignetteIndexEntry{PCAtools: everything Principal Component Analysis}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
    %\usepackage[utf8]{inputenc}
---

# Introduction

Principal Component Analysis (PCA) is a very powerful technique that has wide applicability in data science, bioinformatics, and further afield. It was initially developed to analyse large volumes of data in order to tease out the differences/relationships between the logical entities being analysed. It extracts the fundamental structure of the data without the need to build any model to represent it. This 'summary' of the data is arrived at through a process of reduction that can transform the large number of variables into a lesser number that are uncorrelated (i.e. the ‘principal components'), whilst at the same time being capable of easy interpretation on the original data [@PCAtools] [@BligheK].

*PCAtools* provides functions for data exploration via PCA, and allows the user to generate publication-ready figures. PCA is performed via *BiocSingular* [@Lun] - users can also identify optimal number of principal components via different metrics, such as elbow method and Horn's parallel analysis [@Horn] [@Buja], which has relevance for data reduction in single-cell RNA-seq (scRNA-seq) and high dimensional mass cytometry data.

# Installation

## 1. Download the package from Bioconductor

```{r getPackage, eval=FALSE}

    if (!requireNamespace('BiocManager', quietly = TRUE))
        install.packages('BiocManager')
        BiocManager::install('PCAtools')

```

Note: to install development version:

```{r getPackageDevel, eval=FALSE}

    devtools::install_github('kevinblighe/PCAtools')

```

## 2. Load the package into R session

```{r Load, message=FALSE}

    library(PCAtools)

```

# Quick start

For this vignette, we will load breast cancer gene expression data with
recurrence free survival (RFS) from
[Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade To Improve Prognosis](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE2990).

First, let's read in and prepare the data:


```{r}

  library(Biobase)
  library(GEOquery)

  # load series and platform data from GEO
  gset <- getGEO('GSE2990', GSEMatrix = TRUE, getGPL = FALSE)

  x <- exprs(gset[[1]])

  # remove Affymetrix control probes
  x <- x[-grep('^AFFX', rownames(x)),]

  # extract information of interest from the phenotype data (pdata)
  idx <- which(colnames(pData(gset[[1]])) %in%
    c('age:ch1', 'distant rfs:ch1', 'er:ch1',
      'ggi:ch1', 'grade:ch1', 'size:ch1',
      'time rfs:ch1'))

  metadata <- data.frame(pData(gset[[1]])[,idx],
    row.names = rownames(pData(gset[[1]])))

  # tidy column names
  colnames(metadata) <- c('Age', 'Distant.RFS', 'ER', 'GGI', 'Grade',
    'Size', 'Time.RFS')

  # prepare certain phenotypes
  metadata$Age <- as.numeric(gsub('^KJ', NA, metadata$Age))
  metadata$Distant.RFS <- factor(metadata$Distant.RFS, levels=c(0,1))
  metadata$ER <- factor(gsub('\\?', NA, metadata$ER), levels=c(0,1))
  metadata$ER <- factor(ifelse(metadata$ER == 1, 'ER+', 'ER-'), levels = c('ER-', 'ER+'))
  metadata$GGI <- as.numeric(metadata$GGI)
  metadata$Grade <- factor(gsub('\\?', NA, metadata$Grade), levels=c(1,2,3))
  metadata$Grade <- gsub(1, 'Grade 1', gsub(2, 'Grade 2', gsub(3, 'Grade 3', metadata$Grade)))
  metadata$Grade <- factor(metadata$Grade, levels = c('Grade 1', 'Grade 2', 'Grade 3'))
  metadata$Size <- as.numeric(metadata$Size)
  metadata$Time.RFS <- as.numeric(gsub('^KJX|^KJ', NA, metadata$Time.RFS))

  # remove samples from the pdata that have any NA value
  discard <- apply(metadata, 1, function(x) any(is.na(x)))

  metadata <- metadata[!discard,]

  # filter the expression data to match the samples in our pdata
  x <- x[,which(colnames(x) %in% rownames(metadata))]

  # check that sample names match exactly between pdata and expression data 
  all(colnames(x) == rownames(metadata))

```

Conduct principal component analysis (PCA)

```{r}

  p <- pca(x, metadata = metadata, removeVar = 0.1)

```

## A scree plot

```{r ex1, fig.height = 7, fig.width = 16, fig.cap = 'Figure 1: A scree plot to show the proportion of explained variance by PC'}

  screeplot(p)

```

## A bi-plot

```{r ex2, fig.height = 7, fig.width = 8, fig.cap = 'Figure 2: A bi-plot of PC1 versus PC2'}

  biplot(p)

```

## A pairs plot

```{r ex3, fig.height = 10, fig.width = 10, fig.cap = 'Figure 3: A pairs plot, comparing PC1 - PC5 on a pairwise basis'}

  pairsplot(p)

```

## A loadings plot

```{r ex4, fig.height = 6, fig.width = 8, fig.cap = 'Figure 4: Plot the component loadings and label genes most responsible for variation'}

  plotloadings(p)

```

## An eigencor plot

```{r ex5, fig.height = 4, fig.width = 8, fig.cap = 'Figure 5: Correlate PCs to metadata variables'}

  eigencorplot(p,
    metavars = c('Age','Distant.RFS','ER','GGI','Grade','Size','Time.RFS'))

```

# Advanced features

All plots in PCAtools are highly configurable and should cover virtually all
general usage requirements. The following sections take a look at some of these
advanced features, and form a somewhat practical example of how one can use
PCAtools to make a clinical interpretation of data.

## Determine optimum number of PCs to retain

A scree plot on its own just shows the accumulative proportion of explained variation,
but how can we determine the optimum number of PCs to retain?  *PCAtools* provides two
metrics for this purpose: elbow method and Horn's parallel analysis [@Horn] [@Buja].

Let's perform Horn's parallel analysis first:

```{r}

  horn <- parallelPCA(x)
  horn$n

```

Now the elbow method:

```{r}

  elbow <- findElbowPoint(p$variance)
  elbow

```

In most cases, the identified values will disagree. This is because finding the
correct number of PCs is a difficult task and is akin to finding the 'correct'
number of clusters in a dataset - there is no correct answer.

Taking the value from Horn's parallel analyis, we can produce a new scree plot:

```{r ex6, fig.height = 7, fig.width = 9, fig.cap = 'Figure 6: Advanced scree plot illustrating optimum number of PCs'}

  library(ggplot2)

  screeplot(p,
    components = getComponents(p, 1:20),
    vline = c(horn$n, elbow)) +
    geom_text(aes(horn$n + 1, 50, label = "Horn's", vjust = -1)) +
    geom_text(aes(elbow + 1, 50, label = "Elbow", vjust = -1))

```

If all else fails, one can simply take the number of PCs that contributes to
a pre-selected total of explained variation, e.g., in this case, 27 PCs account
for >80% explained variation.

## Modify bi-plots

The bi-plot comparing PC1 versus PC2 is the most characteristic plot of PCA.
However, PCA is much more than the bi-plot and much more than PC1 and PC2. This
said, PC1 and PC2, by the very nature of PCA, are indeed usually the most
important parts of PCA.

In a bi-plot, we can shade the points by different groups and add many more features.

### Colour by a factor from the metadata, use a custom label, add lines through center, and add legend

```{r ex7, fig.height = 6, fig.width = 8, fig.cap = 'Figure 7: adding lines and a legend to a bi-plot'}

  biplot(p,
    lab = paste0(p$metadata$Age, 'yo'),
    colby = 'ER',
    hline = 0, vline = 0,
    legendPosition = 'right')

```

### Supply custom colours, add more lines, and increase legend size

```{r ex8, fig.height = 7, fig.width = 7, fig.cap = 'Figure 8: supplying custom colours to a bi-plot'}

  biplot(p,
    colby = 'ER', colkey = c('ER+'='forestgreen', 'ER-'='purple'),
    hline = 0, vline = c(-25, 0, 25),
    legendPosition = 'top', legendLabSize = 16, legendIconSize = 8.0)

```

### Change shape based on tumour grade, remove connectors, and add titles

```{r ex9, fig.height = 7, fig.width = 7, fig.cap = 'Figure 9: supplying custom colours and shapes to a bi-plot, removing connectors, and modifying titles'}

  biplot(p,
    colby = 'ER', colkey = c('ER+'='forestgreen', 'ER-'='purple'),
    hline = 0, vline = c(-25, 0, 25),
    legendPosition = 'top', legendLabSize = 16, legendIconSize = 8.0,
    shape = 'Grade', shapekey = c('Grade 1'=15, 'Grade 2'=17, 'Grade 3'=8),
    drawConnectors = FALSE,
    title = 'PCA bi-plot',
    subtitle = 'PC1 versus PC2',
    caption = '27 PCs == 80%')

```

### Remove labels, modify line types, remove gridlines, and increase point size

```{r ex10, fig.height = 6, fig.width = 8, fig.cap = 'Figure 10: removing labels, modifying line types, removing gridlines, and increasing point size in a bi-plot'}

  biplot(p,
    lab = NULL,
    colby = 'ER', colkey = c('ER+'='royalblue', 'ER-'='red3'),
    hline = 0, vline = c(-25, 0, 25),
    vlineType = c('dotdash', 'solid', 'dashed'),
    gridlines.major = FALSE, gridlines.minor = FALSE,
    pointSize = 5,
    legendPosition = 'left', legendLabSize = 16, legendIconSize = 8.0,
    shape = 'Grade', shapekey = c('Grade 1'=15, 'Grade 2'=17, 'Grade 3'=8),
    drawConnectors = FALSE,
    title = 'PCA bi-plot',
    subtitle = 'PC1 versus PC2',
    caption = '27 PCs == 80%')

```

### Colour by a continuous variable (colour controlled by ggplot2 engine); plot other PCs

```{r ex11, fig.height = 6, fig.width = 8, fig.cap = 'Figure 11: colouring by a continuous variable and plotting other PCs in a bi-plot'}

  biplot(p, x = 'PC10', y = 'PC50',
    lab = NULL,
    colby = 'Age',
    hline = 0, vline = 0,
    hlineWidth = 1.0, vlineWidth = 1.0,
    gridlines.major = FALSE, gridlines.minor = TRUE,
    pointSize = 5,
    legendPosition = 'left', legendLabSize = 16, legendIconSize = 8.0,
    shape = 'Grade', shapekey = c('Grade 1'=15, 'Grade 2'=17, 'Grade 3'=8),
    drawConnectors = FALSE,
    title = 'PCA bi-plot',
    subtitle = 'PC10 versus PC50',
    caption = '27 PCs == 80%')

```

## Quickly explore potentially informative PCs via a pairs plot

The pairs plot in PCA unfortunately suffers from a lack of use; however, for
those who love exploring data and squeezing every last ounce of information out
of data, a pairs plot provides for a relatively quick way to explore useful
leads for other downstream analyses.

As the number of pairwise plots increases, however, space becomes limited. We
can shut off titles and axis labeling to save space. Reducing point size and
colouring by a variable of interest can additionally help us to rapidly skim
over the data.

```{r ex12, fig.height = 8, fig.width = 7, fig.cap = 'Figure 12: maximising available space in a pairs plot'}

  pairsplot(p,
    components = getComponents(p, c(1:10)),
    triangle = TRUE, trianglelabSize = 12,
    hline = 0, vline = 0,
    pointSize = 0.4,
    gridlines.major = FALSE, gridlines.minor = FALSE,
    colby = 'Grade',
    title = 'Pairs plot', plotaxes = FALSE,
    margingaps = unit(c(-0.01, -0.01, -0.01, -0.01), 'cm'))

```

We can arrange these in a way that makes better use of the screen space by
setting 'triangle = FALSE'. In this case, we can further control the layout
with the 'ncol' and 'nrow' parameters, although, the function will
automatically determine these based on your input data.

```{r ex13, fig.height = 6, fig.width = 9, fig.cap = 'Figure 13: arranging a pairs plot horizontally'}

  pairsplot(p,
    components = getComponents(p, c(4,33,11,1)),
    triangle = FALSE,
    hline = 0, vline = 0,
    pointSize = 0.8,
    gridlines.major = FALSE, gridlines.minor = FALSE,
    colby = 'ER',
    title = 'Pairs plot', plotaxes = TRUE,
    margingaps = unit(c(0.1, 0.1, 0.1, 0.1), 'cm'))

```

## Determine the variables that drive variation among each PC

If, on the bi-plot or pairs plot, we encounter evidence that 1 or more PCs
are segregating a factor of interest, we can explore further the genes that
are driving these differences along each PC.

For each PC of interest, 'plotloadings' determines the variables falling within
the top/bottom 5% of the loadings range, and then creates a final consensus list of
these. These variables are then plotted.

The loadings plot, like all others, is highly configurable. To modify the
cut-off for inclusion / exclusion of variables, we use 'rangeRetain', where
0.01 equates to the top/bottom 1% of the loadings range per PC. We can also
add a title, subtitle, and caption, and alter the shape and colour scheme.

```{r ex14, fig.height = 6, fig.width = 9, fig.cap = 'Figure 14: modifying cut-off for labeling in a loadings plot'}

  plotloadings(p,
    rangeRetain = 0.01,
    labSize = 3.0,
    title = 'Loadings plot',
    subtitle = 'PC1, PC2, PC3, PC4, PC5',
    caption = 'Top 1% variables',
    shape = 24,
    col = c('limegreen', 'black', 'red3'),
    drawConnectors = TRUE)

```

We can check the genes to which these relate by using biomaRt:

```{r eval=TRUE}

  library(biomaRt)

  mart <- useMart('ENSEMBL_MART_ENSEMBL')
  mart <- useDataset('hsapiens_gene_ensembl', mart)

  getBM(mart = mart,
    attributes = c('affy_hg_u133a', 'ensembl_gene_id',
      'gene_biotype', 'external_gene_name'),
    filter = 'affy_hg_u133a',
    values = c('215281_x_at', '214464_at', '211122_s_at', '205225_at',
      '202037_s_at', '204540_at', '215176_x_at', '205044_at', '208650_s_at',
      '205380_at'),
    uniqueRows = TRUE)

```

At least one interesting finding is 205225_at (ESR1), which is by far the gene
most responsible for variation along PC2. The previous bi-plots showed that
this PC also segregated ER+ from ER- patients. The other results could be
explored.

With the loadings plot, in addition, we can instead plot absolute values and
modify the point sizes to be proportional to the loadings. We can also switch
off the line connectors and plot the loadings for any PCs

```{r ex15, fig.height = 9, fig.width = 11, fig.cap = 'Figure 15: plotting absolute component loadings'}

  plotloadings(p,
    components = getComponents(p, c(4,33,11,1)),
    rangeRetain = 0.1,
    labSize = 3.0,
    absolute = FALSE,
    title = 'Loadings plot',
    subtitle = 'Misc PCs',
    caption = 'Top 10% variables',
    shape = 23, shapeSizeRange = c(1, 16),
    col = c('white', 'pink'),
    drawConnectors = FALSE)

```

## Correlate the principal components back to the clinical data

Further exploration of the PCs can come through correlations with clinical data.
This is also a mostly untapped resource in the era of 'big data' and can help
to guide an analysis down a particular path (or not!).

We may wish, for example, to correlate all PCs that account for 80% variation in
our dataset and then explore further the PCs that have statistically significant
correlations.

'eigencorplot' is built upon another function by the *PCAtools* developers, namely
[CorLevelPlot](https://github.com/kevinblighe/CorLevelPlot). Further examples
can be found there.

```{r ex16, fig.height = 4, fig.width = 12, fig.cap = 'Figure 16: correlating PCs that account for at least 80% variation to clinical variables'}

  eigencorplot(p,
    components = getComponents(p, 1:27),
    metavars = c('Age','Distant.RFS','ER','GGI','Grade','Size','Time.RFS'),
    col = c('darkblue', 'blue2', 'black', 'red2', 'darkred'),
    cexCorval = 0.7,
    colCorval = 'white',
    fontCorval = 2,
    posLab = 'bottomleft',
    rotLabX = 45,
    posColKey = 'top',
    cexLabColKey = 1.5,
    scale = TRUE,
    main = 'PC1-27 clinical correlations',
    colFrame = 'white',
    plotRsquared = FALSE)

```

We can also supply different cut-offs for statistical significance, plot
R-squared values, and specify correlation method:

```{r ex17, fig.height = 5, fig.width = 12, fig.cap = 'Figure 17: modifying cut-offs and symbols for statistical significance in eigencorplot'}

  eigencorplot(p,
    components = getComponents(p, 1:horn$n),
    metavars = c('Age','Distant.RFS','ER','GGI','Grade','Size','Time.RFS'),
    col = c('white', 'cornsilk1', 'gold', 'forestgreen', 'darkgreen'),
    cexCorval = 1.2,
    fontCorval = 2,
    posLab = 'all',
    rotLabX = 45,
    scale = TRUE,
    main = bquote(Principal ~ component ~ Pearson ~ r^2 ~ clinical ~ correlates),
    plotRsquared = TRUE,
    corFUN = 'pearson',
    corUSE = 'pairwise.complete.obs',
    signifSymbols = c('****', '***', '**', '*', ''),
    signifCutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1))

```

Clearly, PC2 is coming across as the most interesting PC in this experiment,
with highly statistically significant correlation (p<0.0001) to ER status,
tumour grade, and GGI (genomic Grade Index), an indicator of response.
It comes as no surprise that the gene driving most variationn along PC2 is
*ESR1*, identified from our loadings plot.

This information is, of course, not new, but shows how PCA is much more than
just a bi-plot used to identify outliers!

## Plot the entire project on a single panel

```{r ex18, fig.height = 10, fig.width = 15, fig.cap = 'Figure 18: a merged panel of all PCAtools plots'}

  pscree <- screeplot(p, components = getComponents(p, 1:30),
    hline = 80, vline = 27, axisLabSize = 10, returnPlot = FALSE) +
    geom_text(aes(20, 80, label = '80% explained variation', vjust = -1))

  ppairs <- pairsplot(p, components = getComponents(p, c(1:3)),
    triangle = TRUE, trianglelabSize = 12,
    hline = 0, vline = 0,
    pointSize = 0.8, gridlines.major = FALSE, gridlines.minor = FALSE,
    colby = 'Grade',
    title = '', titleLabSize = 16, plotaxes = FALSE,
    margingaps = unit(c(0.01, 0.01, 0.01, 0.01), 'cm'),
    returnPlot = FALSE)

  pbiplot <- biplot(p, lab = NULL,
    colby = 'ER', colkey = c('ER+'='royalblue', 'ER-'='red3'),
    hline = 0, vline = c(-25, 0, 25), vlineType = c('dotdash', 'solid', 'dashed'),
    gridlines.major = FALSE, gridlines.minor = FALSE,
    pointSize = 2, axisLabSize = 12,
    legendPosition = 'left', legendLabSize = 10, legendIconSize = 3.0,
    shape = 'Grade', shapekey = c('Grade 1'=15, 'Grade 2'=17, 'Grade 3'=8),
    drawConnectors = FALSE,
    title = 'PCA bi-plot', subtitle = 'PC1 versus PC2',
      caption = '27 PCs == 80%',
    returnPlot = FALSE)

  ploadings <- plotloadings(p, rangeRetain = 0.01, labSize = 2.5,
    title = 'Loadings plot', axisLabSize = 12,
    subtitle = 'PC1, PC2, PC3, PC4, PC5',
    caption = 'Top 1% variables',
    shape = 24, shapeSizeRange = c(4, 4),
    col = c('limegreen', 'black', 'red3'),
    legendPosition = 'none',
    drawConnectors = FALSE,
    returnPlot = FALSE)

  peigencor <- eigencorplot(p,
    components = getComponents(p, 1:10),
    metavars = c('Age','Distant.RFS','ER','GGI','Grade','Size','Time.RFS'),
    #col = c('royalblue', '', 'gold', 'forestgreen', 'darkgreen'),
    cexCorval = 0.6,
    fontCorval = 2,
    posLab = 'all', 
    rotLabX = 45,
    scale = TRUE,
    main = "PC clinical correlates",
    cexMain = 1.5,
    plotRsquared = FALSE,
    corFUN = 'pearson',
    corUSE = 'pairwise.complete.obs',
    signifSymbols = c('****', '***', '**', '*', ''),
    signifCutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1),
    returnPlot = FALSE)

    library(cowplot)
    library(ggplotify)

    top_row <- plot_grid(pscree, ppairs, pbiplot,
      ncol = 3,
      labels = c('A', 'B  Pairs plot', 'C'),
      label_fontfamily = 'serif',
      label_fontface = 'bold',
      label_size = 22,
      align = 'h',
      rel_widths = c(1.05, 0.9, 1.05))

    bottom_row <- plot_grid(ploadings,
      as.grob(peigencor),
      ncol = 2,
      labels = c('D', 'E'),
      label_fontfamily = 'serif',
      label_fontface = 'bold',
      label_size = 22,
      align = 'h',
      rel_widths = c(1.5, 1.5))

    plot_grid(top_row, bottom_row, ncol = 1, rel_heights = c(1.0, 1.0))

```

# Acknowledgments

The development of *PCAtools* has benefited from contributions
and suggestions from:

Krushna Chandra Murmu

Jinsheng

# Session info

```{r}

sessionInfo()

```

# References