---
title: 'Multivariate Analysis: PCA and batch removal'
author: "Susan Holmes"
date: "Pune, 17/12/2019"
output: 
  BiocStyle::html_document:
    number_sections: no
    toc: yes
    toc_depth: 4
    toc_float: true
editor_options: 
  chunk_output_type: console

---


## Goal

In this lab we will learn the basics of Multivariate Analysis and PCA using 
a few simple examples.

Work through this lab by running all the R code to your computer and making sure 
that you understand the input and the output. Make alterations where you seem
fit. We encourage you to work through  this lab with a partner. 

This Lab is part of two book chapters from Modern Statistics for Modern Biology by
Susan Holmes and Wolfgang Huber @HolmesHuber2019 which are available freely [online Chapter 7](http://web.stanford.edu/class/bios221/book/Chap-MultivaHetero.html)
[online Chapter 9](http://web.stanford.edu/class/bios221/book/Chap-MultivaHetero.html).

```{r setup, warning=FALSE, message=FALSE, results='hide'}
knitr::opts_chunk$set(echo = TRUE)
pkgs_needed = c("tidyverse","GGally", "factoextra", "ade4")
setdiff(pkgs_needed, installed.packages())
#If this is not empty you made need to execute the following line:
#BiocManager::install(setdiff(pkgs_needed, installed.packages()))
library("tidyverse")
```


## Obtain Data 

In this lab we will be working with the following datasets:

```{r}
#turtles = read.table(url("https://web.stanford.edu/class/bios221/data/PaintedTurtles.txt"),
#                 header=TRUE)
turtles<-read.table("../../data/PaintedTurtles.txt",header=TRUE)
head(turtles)
```

```{r}
#download.file(url = "https://web.stanford.edu/class/bios221/data/athletes.RData",
#              destfile = "athletes.RData",mode = "wb")
load("athletes.RData")
athletes[1:3,]
```

## Low dimensional data summaries and preparation

It is instructive to first consider 2-dimensional summaries of the datasets:

```{r}
library("GGally")
ggpairs(turtles[,-1], axisLabels="none")
```

Can you do the same for the athletes data?

```{r}
#TODO

```

Correlations can be displayed on a color scale by a simple call to the `pheatmap` function:

```{r heatmapathletes}
library("pheatmap")
pheatmap(cor(athletes),cell.width=10,cell.height=10)
```

## Preprocessing the data

Our first task in data analysis is to transform the data: standardizing the data 
to a common standard deviation. This rescaling is done using the `scale` 
function which makes every column have a variance of 1 (and also mean 0).

```{r turtlesDim12}
scaledTurtles=data.frame(scale(turtles[,-1]),sex=turtles[,1])
ggplot(scaledTurtles,aes(x=width,y=height, group =sex)) +
  geom_point(aes(color=sex))
```

Can you compute the standard deviation and mean of each column in the `turtles`
data frame? Can you do the same on the scaled dataset, i.e. on `scaledturtles`?

**Question**: What was the mean of turtles' heights before standardizing?

**Answer:**

```{r}
mean(turtles$length)
```


**Question**: What was the standard deviation of turtles' widths before standardizing?

**Answer:**
```{r}
sd(turtles$width)
```

**Question**: What was the standard deviation of turtles' widths after standardizing?

**Answer:**
```{r}
sd(scaledTurtles$width)
```


## Dimension reduction

```{r}
library("ggplot2")
athletes = scale(athletes)
n = nrow(athletes)
athletes = data.frame(athletes)
```


```{r SimpleScatter}
p = ggplot(athletes, aes(x = weight,y=  disc)) +
  geom_point(size = 2, shape=21)
p + geom_point(aes(y = rep(0, n)), colour="red") +
  geom_segment(aes(xend = weight, yend = rep(0,n)), linetype = "dashed")
```

Now try to do the following:

**Question**: Calculate the variance of the red points in the above
figure. (points are projected onto the weight-axis).

**Answer:**
```{r}
var(athletes$weight)
```


## Summarize 2D-data by a line

We regress `disc` on `weight` with the `lm` function (linear model) to find the regression line; its slope (a) is given by the second coefficient in the output of `lm` and its intercept (b) is the first:

```{r Reg1}
ath_gg = ggplot(athletes, aes(x = weight, y = disc)) +
  geom_point(size = 2, shape = 21)
reg1 = lm(disc ~ weight, data = athletes)
a1 = reg1$coefficients[1] # intercept
b1 = reg1$coefficients[2] # slope
pline1 = ath_gg + geom_abline(intercept = a1, slope = b1,
    col = "blue", lwd = 1.5)
pline1 + geom_segment(aes(xend = weight, yend = reg1$fitted),
    colour = "red", arrow = arrow(length = unit(0.15, "cm")))
```

Can you regress `weight` on `discs` and generate a similar plot? 

```{r}
reg2 = lm(weight ~ disc, data = athletes)
a2 = reg2$coefficients[1] # intercept
b2 = reg2$coefficients[2] # slope
pline2 = ath_gg + geom_abline(intercept = -a2/b2, slope = 1/b2,
    col = "darkgreen", lwd = 1.5)
pline2 + geom_segment(aes(xend=reg2$fitted, yend=disc),
    colour = "orange", arrow = arrow(length = unit(0.15, "cm")))
```

Can you create a plot that shows all points, as well as both regression lines, i.e., a plot that show both the line you get from `lm(disc ~ weight)` and `lm(weight ~ disc)`?

```{r}
pline1 + geom_segment(aes(xend = weight, yend = reg1$fitted), colour = "blue", alpha = 0.35) +
  geom_abline(intercept = -a2/b2, slope = 1/b2, col = "darkgreen", lwd = 1.5, alpha = 0.8) +
  geom_segment(aes(xend = reg2$fitted, yend = disc), colour = "orange", alpha = 0.35) +
   coord_fixed()
```

### A line that minimizes distances in both directions

Now we will plot the line chosen to minimize the sum of squares of the orthogonal (perpendicular) projections of data points onto it; we call this the principal component line.

```{r PCAmin}
X = cbind(athletes$disc, athletes$weight)
svda = svd(X)
pc = X %*% svda$v[, 1] %*% t(svda$v[, 1])
bp = svda$v[2, 1] / svda$v[1, 1]
ap = mean(pc[, 2]) - bp * mean(pc[, 1])

p + geom_segment(xend = pc[,1], yend = pc[,2]) + 
  geom_abline(intercept = ap, slope = bp, col = "purple", lwd = 1.5) + 
  coord_fixed()
```

Can you create a plot that includes both the line from the plot above, plus the two regression lines `lm(disc ~ weight)` and `lm(weight ~ disc)`?

```{r}
pline1 + geom_segment(aes(xend = weight, yend = reg1$fitted), colour = "blue", alpha = 0.35) +
  geom_abline(intercept = -a2/b2, slope = 1/b2, col = "darkgreen", lwd = 1.5, alpha = 0.8) +
  geom_segment(aes(xend = reg2$fitted, yend = disc), colour = "orange", alpha = 0.35) +
  geom_abline(intercept = ap, slope = bp, col = "purple", lwd = 1.5, alpha = 0.8) +
  geom_segment(xend = pc[, 1], yend = pc[, 2], colour = "purple", alpha = 0.35) + coord_fixed()
```


If we rotate the `(discus, weight)` plane with this change of coordinates making the purple line the horizontal $x$ axis, we obtain what is know as the first principal plane:

```{r CompareSDs}
ppdf = data.frame(PC1n = -svda$u[,1]*svda$d[1], PC2n=svda$u[,2] * svda$d[2])

ggplot(ppdf,aes(x = PC1n, y=PC2n)) + geom_point() + ylab("PC2 ") +
  geom_hline(yintercept = 0, color = "purple", lwd = 1.5, alpha = 0.5) +
  geom_point(aes(x = PC1n, y = 0),color = "red") + xlab("PC1 ")+
  xlim(-3.5, 2.7)+ylim(-2, 2) + coord_fixed() +
  geom_segment(aes(xend = PC1n,yend = 0), color = "red")
```

**Question**: What are the sums of squares of the red segments equal to?

**Answer:**
```{r}
sum(ppdf$PC2n^2)
```

**Question**: What is the variance of this new set of red points?

**Answer:**
```{r}
var(ppdf$PC1n)
```

**Question**: What is the sum of the variances of `ppdf$PC1n` and `ppdf$PC2n`?

**Answer:**
```{r}
var(ppdf$PC1n) + var(ppdf$PC2n)
```
We could have gotten the same results using the `princomp` command as follows:

```{r}
pca_athletes = princomp(X)
```

Now compare (note that e.g. loadings are not unique up to sign, but the lines they define are the same):

```{r}
svda$v
pca_athletes$loadings
```

```{r}
head(pca_athletes$scores)
```

```{r}
head(svda$u %*% diag(svda$d))
```

**Question**: Which field in `pca_athletes` contains approximately the same object as `c(sd(ppdf$PC1n), sd(ppdf$PC2n))`? 

**Answer:**
```{r}
pca_athletes$sdev
```

**Question**: Unfortunately the results stored in the above field do not perfectly match  `c(sd(ppdf$PC1n), sd(ppdf$PC2n))`. If you multiply by which correction factor will you get the results to match `c(sd(ppdf$PC1n), sd(ppdf$PC2n))`?

**Answer:**
The factor is $\sqrt(33/32)$.

```{r}
sqrt(33/32)*pca_athletes$sdev
```
Not that the above is the same as:
```{r}
c(sd(ppdf$PC1n), sd(ppdf$PC2n))
```
The difference is that ``princomp`` returns unbiased estimates of sample standard deviations.

## Turtle PCA

Now let's continue inspecting the turtles data.

```{r PCAturtlesunscaled}
turtles3var = turtles[, -1]
apply(turtles3var, 2, mean)
```

We start by looking at the variances of the three components in the **un**standardized case:

```{r simplecomp}
apply(turtles3var, 2, var)
```

Next we see that basically all 3 variables are very strongly correlated:

```{r PCAturtles}
turtlesc = scale(turtles3var)
cor(turtlesc)
```

Because of the strong correlations, we would expect that the data matrix can be well approximated by a rank 1 matrix. Let's do the PCA:

```{r}
library("factoextra")
pca1 = princomp(turtlesc)
# or alternatively:
#pca1 = ade4::dudi.pca(turtlesc, scannf = FALSE)
pca1
fviz_eig(pca1, geom = "bar", width = 0.4)
```

The screeplot showing the eigenvalues for the standardized data: one very large component in this case and two very small ones. In this case the data are (almost) one dimensional.


**Question**: What is the percentage of variance
explained by the first PC?

**Answer:**

```{r}
100*pca1$sdev[1]^2/sum(pca1$sdev^2)
```


```{r turtlesbiplot}
fviz_pca_biplot(pca1, label = "var", col.ind = turtles[,1]) 
```

Add ellipses for female and male groups to the plot above.

```{r}
fviz_pca_biplot(pca1, label = "var", col.ind = turtles[,1], addEllipses=TRUE) 
```

## Back to the athletes

Now let us try to interpret another scree plot with more dimensions.

```{r}
library("ade4")
pca.ath = dudi.pca(athletes, scannf = FALSE)
pca.ath$eig
```

```{r}
fviz_eig(pca.ath, geom = "bar", bar_width = 0.3) + ggtitle("")
```

The screeplot make a clear drop after the second eigenvalue. This indicates a good approximation will be obtained at rank 2. Let’s look at an interpretation of the first two axes by projecting the loadings of the original (old) variables as they project onto the two new ones.

```{r}
fviz_pca_var(pca.ath, col.circle = "black") + ggtitle("")
```

**Note**

It can seem paradoxical that the m variables are opposed to the others.

**Question: Why does this occur?**

We can make the variables align and give the left direction on PCA 1 to
be an axis of athletic ability by changing the signs:

```{r}
athletes[, c(1, 5, 6, 10)] = -athletes[, c(1, 5, 6, 10)]
cor(athletes) %>% round(1)
```

```{r}
pcan.ath = dudi.pca(athletes, nf = 2, scannf = FALSE)
pcan.ath$eig
```

```{r}
fviz_pca_var(pcan.ath, col.circle="black") + ggtitle("")
```

```{r}
fviz_pca_ind(pcan.ath) + ggtitle("") + ylim(c(-2.5,5.7))
```

What do you notice about the numbers?

## Contiguous or supplementary information

Many programs and workflows for biological sequence analysis or assays separate the environmental and contextual information, which they call **metadata**, from the assay data or sequence reads. We discourage such practice as the exact connections between the samples and covariates are important. A lost connection between the assays and covariates makes later analyses impossible. Covariates such as clinical history, time, batch or location are important and should be considered components of the data.

Extra information about sample batches, dates of measurement, different protocols are often named *metadata*; this can be misnomer if it is implied that metadata are somehow less important. Such information is *real data* that need to be integrated into the analyses. We typically store it in a *data.frame* or a similar R class and tightly link it to the primary assay data.

### Known batches in data

Here we show an example of an analysis that was done by @holmes2011psb on bacterial abundance data from *Phylochip* microarrays. The experiment was designed to detect differences between a group of healthy rats and a group who had Irritable Bowel Disease. 
This example shows a case where the nuisance batch effects become apparent in the analysis of experimental data. It is an illustration of the fact that best practices in data analyses are sequential and that it is better to analyse data as they are collected to adjust for severe problems in the experimental design *as they occur*, instead of having to deal with them *post mortem*.

When data collection started on this project, days 1 and 2 were delivered and we made the plot that appears in the Figure below. This showed a definite *day* effect. When investigating the source of this effect, we found that both the protocol and the array were different in days 1 and 2. This leads to uncertainty in the source of variation, we call this **confounding** of effects.

This section will use a new package from Bioconductor called "sva", you need to install it first (take away the comments in the following code, also please load the data from the website).

```{r vsn28Exprs}
#IBDchip = readRDS(url("https://web.stanford.edu/class/bios221/Pune/data/vsn28Exprd.rds"))
IBDchip = readRDS("../../data/vsn28Exprd.rds")
#BiocManager::install("sva")
library("sva")
```

**Question**
What class is the `IBDchip` ? Look at the last row of the matrix, 
what do you notice? 

**Answer**

```{r IBDchip}
class(IBDchip)
dim(IBDchip)
tail(IBDchip[,1:3])
summary(IBDchip[nrow(IBDchip),])
```


The data are normalized abundance measurements of 8634 taxa measured on 28 samples. We use a rank-threshold transformation, giving the top 3000 most abundant tax scores from 3000 to 1, and letting the 5634 least abundant all have a score of 1. We separate out the assay data from the day variable which should be considered a factor:

```{r assayIBD}
assayIBD = IBDchip[-nrow(IBDchip), ]
day = factor(IBDchip[nrow(IBDchip), ])
```

```{r screepc12, fig.keep = 'high', fig.cap = "The screeplot shows us that the phenomenon can be usefully represented in two dimensions.", echo = FALSE, fig.width = 8, fig.height = 6}
rankthreshPCA = function(x, threshold = 3000) {
  ranksM = apply(x, 2, rank)
  ranksM[ranksM < threshold] = threshold
  ranksM = threshold - ranksM
  dudi.pca(t(ranksM), scannf = FALSE, nf = 2)
}
pcaDay12 = rankthreshPCA(assayIBD[,day!=3])
fviz_eig(pcaDay12, bar_width = 0.6) + ggtitle("")
```

Instead of using the continuous, normalized data, we use a robust analysis replacing the values by their ranks. The lower values are considered ties encoded as a threshold chosen to reflect the number of expected taxa thought to be present:

```{r rankthreshPCA-code, warning = FALSE}
rankthreshPCA = function(x, threshold = 3000) {
  ranksM = apply(x, 2, rank)
  ranksM[ranksM < threshold] = threshold
  ranksM = threshold - ranksM
  dudi.pca(t(ranksM), scannf = FALSE, nf = 2)
}
pcaDay12 = rankthreshPCA(assayIBD[,day!=3])
day12 = day[day!=3]
fviz(pcaDay12, element="ind", axes=c(1,2), geom=c("point","text"),
  habillage = day12, repel = TRUE, palette = "Dark2",
  addEllipses = TRUE, ellipse.type = "convex") + ggtitle("") +
  coord_fixed()
```

```{r rankthreshPCA, fig.keep = 'high', fig.cap = "We have used colors to identify the different days and have kept the sample labels as well. We have also added convex hulls for each day. The group mean is identified as the point with the larger symbol (circle, triangle or square).", warning = FALSE, echo = FALSE}
rtPCA1 <- fviz(pcaDay12, element="ind", axes=c(1,2), geom=c("point","text"),
  habillage = day12, repel = TRUE, palette = "Dark2",
  addEllipses = TRUE, ellipse.type = "convex") + ggtitle("") +
  coord_fixed()
rtPCA1
```

**Question**
Why do we use a threshold for the ranks? 

**Answer**
Low abundances, at noise level occur for species that are not really present, of which there are more than half. A large jump in rank for these observations could easily occur without any meaningful reason. Thus we create a large number of ties for low abundance. 

**Question** 
What can you say about the PCA plot?

**Answer**
The figure shows that the sample arrange themselves naturally into two different groups according to the day of the samples. After discovering this effect, we delved into the differences that could explain these distinct clusters. There were two different protocols used (protocol 1 on day 1, protocol 2 on day 2) *and* unfortunately two different provenances for the arrays used on those two days (array 1 on day 1, array 2 on day 2).

A third set of data of four samples had to be collected to deconvolve the confounding effect. Array 2 was used with protocol 2 on Day 3, 

The new figure below shows the new PCA plot with all the samples created by the following:

```{r Threesetspca123-code, warning=FALSE}
pcaDay123 = rankthreshPCA(assayIBD)
fviz(pcaDay123, element="ind", axes=c(1,2), geom=c("point","text"),
  habillage = day, repel=TRUE, palette = "Dark2",
  addEllipses = TRUE, ellipse.type = "convex") + ggtitle("") +
  coord_fixed()
```

```{r Threesetspca123, fig.keep = 'high', fig.show = "hold", fig.cap = "When comparing the three day analysis to that of the first two days, we notice the inversion of signs in the coordinates on the second axis: this has no biological relevance. The important finding is that group 3 overlaps heavily with group 1 indicating that it was the protocol change on Day 2 which created the variability.", fig.margin = FALSE, warning=FALSE,echo = FALSE}
rtPCA1
fviz(pcaDay123, element="ind", axes=c(1,2), geom=c("point","text"),
  habillage = day, repel=TRUE, palette = "Dark2",
  addEllipses = TRUE, ellipse.type = "convex") + ggtitle("") +
  coord_fixed()
```



In which situation would it be preferable to make confidence ellipses around the group means using the following code?

```{r exCodeEllipse}
fviz_pca_ind(pcaDay123, habillage = day, labelsize = 3,
  palette = "Dark2", addEllipses = TRUE, ellipse.level = 0.69)
```


```{r screepc123, fig.keep = 'high', fig.cap = "The eigenvalue screeplot the case of 3 groups is extremely similar to that with two groups shown in Figure \\@ref(fig:screepc12)."}
fviz_eig(pcaDay123,bar_width=0.6) + ggtitle("")
```

Through this visualization we were able to uncover a flaw in the original experimental design. The first two batches shown in green and brown were both balanced with regards to IBS and healthy rats. They do show very different levels of variability and overall multivariate coordinates. In fact, there are two **confounded** effects. Both the arrays and protocols were different on those two days. We had to run a third batch of experiments on day 3, represented in purple, this used protocol from day 1 and the arrays from day 2. The third group faithfully overlaps with batch 1, telling us that the change in protocol was responsible for the variability.

### Removing batch effects

Through the combination of the continuous measurements from `assayIBD` and the **supplementary** batch number as a factor, the PCA map has provided an invaluable investigation tool. 

This is a good example of the use of **supplementary points**. 

The mean-barycenter points are created by using the group-means of points in each of the three groups and serve as extra markers on the plot.

We can decide to re-align the three groups by subtracting the group means so that all the batches are centered on the origin. A slightly more effective way is to use the `ComBat` function available in the **[sva](https://bioconductor.org/packages/sva/)** package. This function uses a similar, but slightly more sophisticated method (Empirical Bayes mixture approach [@Leek:2010:batch]). We can see its effect on the data by redoing our robust PCA (see the result in this Figure:

```{r CombatIBD, fig.keep = 'high', fig.cap = "The modified data with the batch effects removed now show three batch-groups heavily overlapping and centered almost at the origin.",message = FALSE, warning = FALSE}
model0 = model.matrix(~1, day)
combatIBD = ComBat(dat = assayIBD, batch = day, mod = model0)
pcaDayBatRM = rankthreshPCA(combatIBD)
fviz(pcaDayBatRM, element = "ind", geom = c("point", "text"),
  habillage = day, repel=TRUE, palette = "Dark2", addEllipses = TRUE,
  ellipse.type = "convex", axes =c(1,2)) + coord_fixed() + ggtitle("")
```

Also see a quick overview of the important steps:

[10 tips for dimension reduction](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006907)


## References


