Chapter 5 Data normalisation: centring, scaling, quantile normalisation

5.1 Introduction

Copy link

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).

5.2 Data transformation

Copy link

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.

Figure 5.1: Raw quantiation values from the proteomics CPTAC data.

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

##       6A_7                6A_8                6A_9          
##  Min.   :    15014   Min.   :    14940   Min.   :    19946  
##  1st Qu.:   154475   1st Qu.:   155875   1st Qu.:   169510  
##  Median :   345620   Median :   322680   Median :   362710  
##  Mean   :  1750921   Mean   :  1542153   Mean   :  1675242  
##  3rd Qu.:  1020500   3rd Qu.:   930635   3rd Qu.:  1044450  
##  Max.   :107120000   Max.   :124530000   Max.   :133480000  
##       6B_7                6B_8               6B_9          
##  Min.   :    12167   Min.   :   14447   Min.   :    18980  
##  1st Qu.:   145530   1st Qu.:  144330   1st Qu.:   172565  
##  Median :   325910   Median :  309180   Median :   366900  
##  Mean   :  1677453   Mean   : 1454185   Mean   :  1743025  
##  3rd Qu.:   948920   3rd Qu.:  881775   3rd Qu.:  1060350  
##  Max.   :112820000   Max.   :94944000   Max.   :128870000

Such data are difficult to visualise and to analyse, because of the skewness of their distribution. The first step is thus generally (but not always77 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 has been shown to accurately model technical and biological variability of such experimental data.) to log-transform the data, as shown below.

Figure 5.2: CPTAC data after log transformation

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()

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.

Figure 5.3: Non-symetric fold-changes (left) and symetric log fold-changes (right) distributions..

5.3 Normalisation

Copy link

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.

5.3.1 Centring and scaling

Copy link

► Question

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.

► Solution

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)

##    sample1   sample2  sample3   sample4
## 1 4.879049  2.008403 3.488396  6.549078
## 2 5.539645  1.920090 4.236938  6.016729
## 3 9.117417  3.964040 3.458411  2.655504
## 4 6.141017  3.735650 5.219228  4.908146
## 5 6.258575 -1.098686 4.174136 14.795471
## 6 9.430130  6.081147 3.384732  6.887755

boxplot(x)

Figure 5.4: Simulated data for 4 samples and 1000 genes.

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.

► Question

For each column, calculate its mean value and subtract it from the values. Visualise and interpret the centred data.

► Solution

cmns <- colMeans(x)
for (i in 1:ncol(x))
    x[, i] <- x[, i] - cmns[i]
boxplot(x)

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.

► Question

Calculate the standard deviation of each column and divide the values by it. Visualise and interpret the centred data.

► Solution

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

► Question

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.

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

##   S1 S2 S3
## A  5  4  3
## B  2  1  4
## C  3  4  6
## D  4  2  8

(rnk <- apply(x, 2, rank, ties.method = "min"))

##   S1 S2 S3
## A  4  3  1
## B  1  1  2
## C  2  3  3
## D  3  2  4

sorted_x <- apply(x, 2, sort)
ranked_means <- rowMeans(sorted_x)
cbind(sorted_x, ranked_means)

##      S1 S2 S3 ranked_means
## [1,]  2  1  3     2.000000
## [2,]  3  2  4     3.000000
## [3,]  4  4  6     4.666667
## [4,]  5  4  8     5.666667

x_norm <- matrix(ranked_means[rnk], ncol = 3)
dimnames(x_norm) <- dimnames(x)
round(x_norm, 2)

##     S1   S2   S3
## A 5.67 4.67 2.00
## B 2.00 2.00 3.00
## C 3.00 4.67 4.67
## D 4.67 3.00 5.67

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()

5.4 Missing data

Copy link

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.

library("QFeatures")
data(se_na2)

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

##  M1F1A  M1F4A  M1F7A M1F11A  M1F2B  M1F5B  M1F8B M1F11B  M2F1A  M2F4A  M2F7A 
##     34     45     56     39     47     52     49     61     41     42     55 
## M2F11A  M2F2B  M2F5B  M2F8B M2F11B 
##     45     51     43     57     53

barplot(col_na)

row_na <- rowSums(is.na(assay(se_na2)))
table(row_na)

## row_na
##   0   1   2   3   4   8   9  10 
## 301 247  91  13   2  23  10   2

hist(row_na)

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

library("QFeatures")
data(se_na2)

se_zero <- impute(se_na2, method = "zero")
anyNA(assay(se_zero))

## [1] FALSE

5.5 Additional exercises

Copy link