50-normalisation.Rmd

# Data normalisation: centring, scaling, quantile normalisation {#sec-norm}

**Learning objectives**

- What is data normalisation and why do we need it?
- What are centring and scaling?
- What is quantile normalisation?
- How to assess if normalisation worked.
- Learn about handling of missing values.

## Introduction

In a typical high throughput experiment, we assay thousands of
features (gene transcripts, proteins, metabolites, ...) in a certain
number of biologically diverse samples (from about 6 to hundreds or
thousands).

- But what do you think would happen is we took the same sample, such
  as a suspension of U2OS cells (a cell line stemming from bone
  cancer), split it in two, and measured it twice using the exact same
  protocol? Will the measurements be exactly identical? Why?

  These two samples would be **technical replicates** and the observed
  differences would be solely the result of **technical variability**.

- If we now take two independent different cell cultures of the same
  cell line, do you expect there to be larger, similar or smaller
  difference than in the example above? Why?

  One would refer to these two samples as **biological replicates**,
  and the observed differences would be due to **technical** and
  **biological variability**.

- In biomedical research, experiments aim at measuring biological
  variability by comparing two or more biological conditions in a
  controlled setting. If we now take two independent cell cultures and
  treat one of them using a drug, and then measure the two samples, do
  you expect there to be larger, similar or smaller difference than in
  the examples above? Why?

  The observed differences between the two samples would comprise
  **technical variability**, **biological variability** and additional
  **biological variability due to the treatment**. To be able to
  measure any biologically signal in the data, the biological
  variability of interest, i.e. the one produced by the treatment,
  must be larger than the technical variability.

However, before analysing any data, it is often (always) necessary to
make the samples as comparable as possible by removing the unwanted
technical variability (that should be shared among all samples)
without removing biological variability (that will differentiate the
samples biologically).

## Data transformation

The very first step in preparing a dataset is to visualise the
distribution of the values. Below we see the density plots for 6
proteomics samples and about 4000 peptides from the CPTAC study that
compares the reproducibility of quantitative mass spectrometry using
human proteins spiked into a yeast proteome background.

```{r logtans, echo = FALSE, message = FALSE, fig.cap = "Raw quantiation values from the proteomics CPTAC data."}
library("rWSBIM1322")
library("tidyverse")
library("SummarizedExperiment")

data(cptac_se)
as_tibble(assay(cptac_se)) %>%
    pivot_longer(names_to = "sample", values_to = "expression",
                 everything()) %>%
    ggplot(aes(x = expression, colour = sample)) +
    geom_density()
```

We see that the vast majority of the data are at very low values with
some very high values:

```{r, echo = FALSE}
summary(assay(cptac_se))
```

Such data are difficult to visualise and to analyse, because of the skewness of their
distribution. The first step is thus generally (but
not always[^negbinom]) to log-transform the data, as shown below.

[^negbinom]: When analysing count data, such as in the case of RNA-Seq
    or spectral counting in proteomics, it is better not to transform
    the data and use dedicated distribution for count-based data. In
    the cases above, the [negative binomial
    distribution](https://en.wikipedia.org/wiki/Negative_binomial_distribution)
    has been shown to accurately model technical and biological
    variability of such experimental data.

```{r logtrans2, echo = FALSE, fig.cap = "CPTAC data after log transformation"}
as_tibble(assay(cptac_se)) %>%
    pivot_longer(names_to = "sample", values_to = "expression",
                 everything()) %>%
    ggplot(aes(x = log2(expression), colour = sample)) +
    geom_density()
```

`r msmbstyle::question_begin()`

Load the CPTAC data from the `rWSBIM1322` package (version 0.1.3 or
later) and reproduce the figures above. See `?cptac` for details .

`r msmbstyle::question_end()`

`r msmbstyle::solution_begin()`
```{r, eval = FALSE}
library("rWSBIM1322")
library("SummarizedExperiment")
data(cptac_se)

library("tidyverse")
as_tibble(assay(cptac_se)) %>%
    pivot_longer(names_to = "sample", values_to = "expression",
                 everything()) %>%
    ggplot(aes(x = expression, colour = sample)) +
    geom_density()

as_tibble(assay(cptac_se)) %>%
    pivot_longer(names_to = "sample", values_to = "expression",
                 everything()) %>%
    ggplot(aes(x = log2(expression), colour = sample)) +
    geom_density()
```
`r msmbstyle::solution_end()`

This is particularly important for fold-changes, as illustrated
below. Typical threshold for down- and up-regulation are 0.5 and 2,
with the absence of any change around 1. In log2 space, these
threshold become -1 and 1, and symmetrically centred around 0.


```{r lfc, echo = FALSE, fig.cap = "Non-symetric fold-changes (left) and symetric log fold-changes (right) distributions.."}
lfc <- rnorm(1000)
fc <- 2^lfc
par(mfrow = c(1, 2))
plot(density(fc), xlab = "Fold-changes", main = "")
abline(v = 1)
abline(v = c(0.5, 2), col = "red")
plot(density(lfc), xlab = "log2 fold-changes", main = "")
abline(v = 0)
abline(v = c(-1, 1), col = "red")
```

## Normalisation

In addition to possible transformation of the data, it is necessary to
further process the data to remove as much as technical variability as
possible while keeping biological variability. This step is called
**normalisation**. One of the important requirements of most
normalisation techniques is that most proteins aren't expected to
change among biological conditions. In other words, normalisation
expects only a minority of biological features to be differentially
expressed in the conditions of interest.

### Centring and scaling

`r msmbstyle::question_begin()`
We are going to start by generating a data set to precisely illustrate
the effect of the methods.

1. Use the `rnorm()` function to generate a distribution of 1000
   values centred around 0 and with a standard deviation
   of 2. Visualise these data.

2. Generate four such distribution with parameters *N(6, 2)*,
   *N(4,2)*, *N(4, 1)*, *N(7, 3)* and create a matrix or dataframe
   with rownames `gene1` to `gene1000` and colnames `sample1` to
   `sample4`. Visualise these data and discuss whether these samples
   could be compared against each other. To assure replication of this
   simulation, set the random number generation seed to 123.
`r msmbstyle::question_end()`


`r msmbstyle::solution_begin()`
```{r sln1, fig.cap = "Simulated data for 4 samples and 1000 genes."}
n <- 1000
set.seed(123)
x <- data.frame(sample1 = rnorm(n, 6, 2),
                sample2 = rnorm(n, 4, 2),
                sample3 = rnorm(n, 4, 1),
                sample4 = rnorm(n, 7, 3))
head(x)
boxplot(x)
```
`r msmbstyle::solution_end()`


**Centring** refers to the operation of modifying the mean value of a set of
values by subtracting a fixed value from each individual value. On the
figure above, this equates to shifting the values up to down. A
typical value is the mean of all the data to be centred.

`r msmbstyle::question_begin()`
For each column, calculate its mean value and subtract it from the
values. Visualise and interpret the centred data.
`r msmbstyle::question_end()`

`r msmbstyle::solution_begin()`
```{r}
cmns <- colMeans(x)
for (i in 1:ncol(x))
    x[, i] <- x[, i] - cmns[i]
boxplot(x)
```
`r msmbstyle::solution_end()`

**Scaling** refers to the operation of rescaling a set of values to
scale in the range of 0 and 1 (or -1 and 1). On the figure above, this
equates to changing the boxes so as to all have similar heights. A
typical scaling method is to dividing the values by their standard
deviations.

`r msmbstyle::question_begin()`
Calculate the standard deviation of each column and divide the values
by it. Visualise and interpret the centred data.
`r msmbstyle::question_end()`


`r msmbstyle::solution_begin()`

```{r}
csdvs <- apply(x, 2, sd)
for (i in 1:ncol(x))
    x[, i] <- x[, i] / csdvs[i]
boxplot(x)
```

`r msmbstyle::solution_end()`


`r msmbstyle::question_begin()`
The above oberations can also be performed with R's `scale`
function. Familiarise yourself with it by reading the documentation,
then regenerate the data above and repeat the scaling/centring
operations using `scale`.
`r msmbstyle::question_end()`


### Quantile normalisation

Quantile normalisation is a method that will make different data
distributions identical. An example is shown below using a small
dataset with quantitative data for three samples (S1, S2, and S3) and
4 genes (A to D) (example taken from the Wikipedia page).

```{r quantex1}
x <- cbind(S1 = c(5, 2, 3, 4),
           S2 = c(4, 1, 4, 2),
           S3 = c(3, 4, 6, 8))
rownames(x) <- LETTERS[1:4]
x
```

The first step is to rank (from lowest to largest) each value in each
sample (column). For sample S1, gene B has the lowest value, hence
rank 1, then gene C gets rank 2, gene D gets rank 3, then gene A, with
the highest value, gets rank 4. We store these in a new matrix `rnk`.

```{r quantex2}
(rnk <- apply(x, 2, rank, ties.method = "min"))
```

We now arrange the values according to their rank (i.e. sorting) and
calculate row-wise means: the mean of all lowest values, ... up to the
mean of all highest values.

```{r quantex3}
sorted_x <- apply(x, 2, sort)
ranked_means <- rowMeans(sorted_x)
cbind(sorted_x, ranked_means)
```

The final step is to replace the ranks in `rnk` by the respective
 `ranked_means`: the gene with the lowest expression in each
samples gets the lowest ranked mean, ..., the gene with the highest
expression in each sample get the highest ranked mean.

```{r quantex4}
x_norm <- matrix(ranked_means[rnk], ncol = 3)
dimnames(x_norm) <- dimnames(x)
round(x_norm, 2)
```

`r msmbstyle::question_begin()`

Using the `cptac_se` data (`SummarizedExperiment` object),
log-transform and normalise it using quantile normalisation, then
visualise the data before and after. To log-transform and normalise
the data, you can use the `logTransform()` and `normalize(method =
"quantiles")` functions from the `QFeatures` package.

`r msmbstyle::question_end()`

`r msmbstyle::solution_begin()`

```{r quantnormex, message = FALSE}
library(QFeatures)
data(cptac_se)
cptac <- logTransform(cptac_se, base = 2)
cptac_quant <- normalize(cptac, method = "quantiles")

as_tibble(assay(cptac)) %>%
    pivot_longer(names_to = "sample", values_to = "expression",
                 everything()) %>%
    ggplot(aes(x = expression, colour = sample)) +
    geom_density()

as_tibble(assay(cptac_quant)) %>%
    pivot_longer(names_to = "sample", values_to = "expression",
                 everything()) %>%
    ggplot(aes(x = expression, colour = sample)) +
    geom_density()
```
`r msmbstyle::solution_end()`

Note that the `scaleTransform()` method can also be used to scale
and/or centre `SummarizedExperiment` objects.

Whenever new data is presented, one of the first steps is to assess
the need for normalisation by verifying the data distributions. Large
variations thereof should indicate that extra care should be taken and
point to possible more serious issues in the data
acquisition. Visualisation of these distributions after normalisation
will necessarily comply with the normalisation strategy employed
(scaling, centering, quantile normalisation).

In the later chapters, we will see dimensionality reduction and
clustering methods, that can be used to assess the grouping of the
samples in an experiment. These methods are useful to check the effect
of the normalisation procedure, and for example verify it the samples
tend to cluster more according to the biological groups defined by the
experiment.

## Missing data

In some omics experiments, a substantial faction of the data can be
missing. Missing values are encoded as `NA` in R, and should be
represented as such in omics data.

Missing data can be very tricky to handle, as they can arise for
different reasons. The first one is related to the absence of the
feature, or an abundance that is below the detection limit of the
device. These data would typically be missing for biological reasons
and thus aren't appearing at random. Such missing values are often
denoted *missing not at random* (MNAR).

However, missing values can also appear randomly, due to uneven or
absent cDNA amplification in transcriptomics from low amounts of RNA
(for instance in in single cell RNA sequencing) or due the
semi-stochastic nature of mass spectrometry in proteomics. Such
missing values are denoted *missing (completely) at random* (M(C)AR).

There are two approaches to deal with missing values:

1. Filtering: remove features that have missing values, or a certain
   proportion of missing values.

2. Imputation: replacing missing values by sensible values. Some
   imputation methods can only be applied to data missing at random,
   while others can only be applied to data missing not at random.

`r msmbstyle::question_begin()`

Load the `se_na2` `SummarizedExperiment` object available in the
`QFeatures` package.

Calculate the number of missing values along the rows and sample, then
visualise these data and represent them as tables showing the number
of occurences of 1, 2, 3, ... missing values.

`r msmbstyle::question_end()`


`r msmbstyle::solution_begin()`
```{r}
library("QFeatures")
data(se_na2)

col_na <- colSums(is.na(assay(se_na2)))
col_na
barplot(col_na)

row_na <- rowSums(is.na(assay(se_na2)))
table(row_na)
hist(row_na)
```
`r msmbstyle::solution_end()`

All the concepts and methods seen in this chapter will be paramount in
the next ones, once we test, cluster or classify our data.


`r msmbstyle::question_begin()`

Impute the missing data by replacing all NAs with a zero. You can do
this with the `QFeatures::impute()` method and setting `method =
"zero"`. See the manual page and references therein (notably
`?MsCoreUtils::impute_matrix`) for help and details.

(Note that imputing by zero isn't an approach that is generally
recommended.)

`r msmbstyle::question_end()`

`r msmbstyle::solution_begin()`

```{r}
library("QFeatures")
data(se_na2)

se_zero <- impute(se_na2, method = "zero")
anyNA(assay(se_zero))
```
`r msmbstyle::solution_end()`


## Additional exercises

`r msmbstyle::question_begin()`
You are provided with a dataset shown below (2 biological groups A
(red) and B (green) with 3 replicates in each) and asked to analyse
it. What are your thoughts about these data and what are the
implications for the normalisation step?

```{r, echo = FALSE}
n <- 1000
x <- data.frame(A1 = rnorm(n, 5, 2),
                A2 = rnorm(n, 5, 2),
                A3 = rnorm(n, 5, 2),
                B1 = rnorm(n, 10, 2),
                B2 = rnorm(n, 10, 2),
                B3 = rnorm(n, 10, 2))
boxplot(x, col = rep(2:3, each = 3))
```
`r msmbstyle::question_end()`


`r msmbstyle::question_begin()`

Generate a dataset composed of 5000 rows and 6 columns, containing
values sampled from a normal distribution *N(0, 1)*. Visualise the
samples distributions. Scale and center the data, visualising the
samples distributions at each step. Interpret what you see.

`r msmbstyle::question_end()`