1.1 Party trick

A fun party trick is to take vector and order it using order()

x <- runif(5)
x

## [1] 0.7809753 0.5210077 0.8932159 0.5110514 0.2235770

xo <- x[order(x)]
xo

## [1] 0.2235770 0.5110514 0.5210077 0.7809753 0.8932159

The original order can the be retrieved by applying order() twice, and using the result to order xo

xo[order(order(x))]

## [1] 0.7809753 0.5210077 0.8932159 0.5110514 0.2235770

Neat, eh? The ordered vector xo can be transformed in a way that might be difficult for the original data, for instance finding differences between ordered values

diff(c(0, xo))[order(order(x))]

## [1] 0.259967575 0.009956273 0.112240639 0.287474455 0.223576978

or transforming numeric values into letters, such that the smallest value of x is given the first letter, etc.

LETTERS[seq_along(xo)][order(order(x))]

## [1] "D" "C" "E" "B" "A"

1.2 Quantile normalization

Create a matrix with two columns representing ‘gene expression’ of 100 (actually, 10,000’s of) genes in two samples. The idea is that ‘base-line’ expression follows some kind of exponential distributition

set.seed(123)
n_genes <- 1000
base <- 2 + rexp(n_genes)

and that the expression observed in each sample is the baseline plus some (normally distributed) error. The following assumes that the size of the measurement error decreases as the gene expression increases; sample y has (e.g., because of sample processing) a consistently stronger signal than sample x.

expression <- data.frame(
    x = base + rnorm(n_genes, 0.3) / base,
    y = 1.1 * base + rnorm(n_genes, 0.4) / base
)

Visualize the expression of each gene in the two samples…

library(ggplot2)
ggplot(expression, aes(x, y)) + geom_point()

and the distribution of expression values across all genes

df <- reshape::melt(expression)

## Using  as id variables

ggplot(df, aes(x = value, colour = factor(variable))) +
    geom_density()

A commmon hypothesis in, e.g., RNA-seq, is that the expression of most genes is unchanged, and differences between distributions are due to technical artifacts. Under this assumption, we expect the distributions in the second figure, above, to be identical. So we could try to normalize samples so that the distributions of each sample look about the same. Quantile normalization is an extreme form of this normalization, and makes the distributions exactly the same.

1. Place each column in order

   x <- expression$x
   xo <- x[order(x)]
   y <- expression$y
   yo <- y[order(y)]

   For our sample of 100, x[1] is the first quantile, x[2] the second, etc.

2. Calculate the mean (or more robust, e.g., median) statistic for each quantile. I’ll do this by using cbind() to create a two-column matrix, and then apply() to iterate over each row to calculate the mean

   mo <- cbind(xo, yo)
   row_mean <- apply(mo, 1, mean)

3. Use the mean of each quantile in the original vector

   normalized <- data.frame(
       x = row_mean[order(order(x))],
       y = row_mean[order(order(y))]
   )

Visualize normalized data – same overall relationship between expression in two samples

ggplot(normalized, aes(x, y)) + geom_point()

…but now identical distributions

df <- reshape::melt(normalized)

## Using  as id variables

ggplot(df, aes(x = value, colour = factor(variable))) +
    geom_density()

2.1 Functions

We could write a script, copy/pasting each time we wanted to use this

x <- expression$x
xo <- x[order(x)]
y <- expression$y
yo <- y[order(y)]
mo <- cbind(xo, yo)
row_mean <- apply(mo, 1, mean)
normalized <- data.frame(
    x = row_mean[order(order(x))],
    y = row_mean[order(order(y))]
)

It would be better to write a function

quantile_normalize <-
    function(expression)
{
    x <- expression$x
    xo <- x[order(x)]
    y <- expression$y
    yo <- y[order(y)]

    mo <- cbind(xo, yo)
    row_mean <- apply(mo, 1, mean)

    data.frame(
        x = row_mean[order(order(x))],
        y = row_mean[order(order(y))]
    )
}

and use that on different data sets.

normalized <- quantile_normalize(expression)

2.2 Testing

We’re going to update our function to be more efficient, general, and robust, but we want to have some confidence that the function still works as expected. So we’ll make a much smaller data set that we can verify ‘by hand’ and use as a test case.

expression <- data.frame(
    x = c(1, 2, 4, 7),
    y = c(3, 2, 5, 8)
)
quantile_normalize(expression)

##     x   y
## 1 1.5 2.5
## 2 2.5 1.5
## 3 4.5 4.5
## 4 7.5 7.5

We can use a formal concept called a ‘unit test’ to check that our code is correct

library(testthat)
test_that("quantile_normalize() works", {
    ## calculate simple example 'by hand'...
    x <- c(1, 2, 4, 7)
    y <- c(3, 2, 5, 8)
    xo <- x[order(x)]
    yo <- y[order(y)]
    mo <- apply(cbind(xo, yo), 1, mean)
    expected <- data.frame(x = mo[order(order(x))], y = mo[order(order(y))])

    ## Compare to outcome of our function
    expression <- data.frame(x = x, y = y)
    observed <- quantile_normalize(expression)
    expect_equal(observed, expected)      # other expectations possible
})

3.1 General

Our function assumes that there are exactly two samples, and that the two samples are labelled x and y. Let’s try and generalize our function using apply() to order an arbitrary number of columns,

m <- apply(expression, 2, function(v) v[order(v)])
m

##      x y
## [1,] 1 2
## [2,] 2 3
## [3,] 4 5
## [4,] 7 8

and apply the same sort of logic to reconstruct the quantiles

quantile_normalize <-
    function(x)
{
    m <- apply(x, 2, function(v) v[order(v)])
    row_mean <- apply(m, 1, mean)
    apply(x, 2, function(v, row_mean) row_mean[order(order(v))], row_mean)
}

Our unit test will fail, because the original function returned a data.frame, whereas the current function returns a matrix. Let’s update the unit test (we could also coerce the matrix to a data.frame before returning from quantile_normalize()).

test_that("quantile_normalize() works", {
    x <- c(1, 2, 4, 7)
    y <- c(3, 2, 5, 8)
    xo <- x[order(x)]
    yo <- y[order(y)]
    mo <- apply(cbind(xo, yo), 1, mean)
    expected <- cbind(x = mo[order(order(x))], y = mo[order(order(y))])

    expression <- data.frame(x = x, y = y)
    observed <- quantile_normalize(expression)
    expect_equal(observed, expected)
})

3.2 Efficient

apply() is an iteration – the function is called once for each element of m. We could instead use rowMeans(), which is vectorized and makes only one function call for the entire data. This is more efficient in R, so let’s try to modify our function

quantile_normalize <-
    function(x)
{
    m <- apply(x, 2, function(v) v[order(v)])
    row_mean <- rowMeans(m)
    apply(x, 2, function(v, row_mean) row_mean[order(order(v))], row_mean)
}

Our unit test still works

test_that("quantile_normalize() works", {
    x <- c(1, 2, 4, 7)
    y <- c(3, 2, 5, 8)
    xo <- x[order(x)]
    yo <- y[order(y)]
    mo <- apply(cbind(xo, yo), 1, mean)
    expected <- cbind(x = mo[order(order(x))], y = mo[order(order(y))])

    expression <- data.frame(x = x, y = y)
    observed <- quantile_normalize(expression)
    expect_equal(observed, expected)
})

so we have some confidence that we’ve changed our implementation without changing the outcome.

3.3 Robust

Good functions provide something like a contract – if certain types of arguments are provided, then the function will return a certain type of value. In our case we might think of the contract as “if you provide me with a data.frame with all numeric columns, or a numeric matrix, I’ll provide you with a matrix of quantile-normalized columns”. There could be additional constraints on the contract, for instance “the input cannot have NA values” or “dimnames of the input are the same as dimnames of the output”. Checking that the inputs satisfy the contract greatly simplifies our function, and typically provides the user with helpful indications before things go wrong in a cryptic way. So it seems like we should validate the inputs of the function as soon as possible.

quantile_normalize <-
    function(x)
{
    ## so long as the input can be coerced to a matrix...
    x <- as.matrix(x)

    ## and can be validated to conform to our contract...
    stopifnot(
        is.numeric(x),
        !anyNA(x)
    )

    ## ...we have confidence that we will satisfy the return value
    m <- apply(x, 2, function(v) v[order(v)])
    row_mean <- rowMeans(m)
    apply(x, 2, function(v, mo) mo[order(order(v))], mo)
}

We can check that our original unit test is satisfied, but also add aditional tests…

test_that("'quantile_normalize()' validates inputs", {
    m <- cbind(letters, LETTERS)
    expect_error(quantile_normalize(m))

    df <- data.frame(x=rnorm(26), y = letters)
    expect_error(quantile_normalize(df))

    m <- matrix(rnorm(10), nrow = 5)
    m[1,1] <- NA
    expect_error(quantile_normalize(m))
})

It seems like our function should preserve dimnames on the original object. The function requires modification, and we can write a new unit test

quantile_normalize <-
    function(x)
{
    ## validate inputs
    x <- as.matrix(x)
    stopifnot(
        is.numeric(x),
        !anyNA(x)
    )

    ## quantile normalization
    m <- apply(x, 2, function(v) v[order(v)])
    row_mean <- rowMeans(m)
    result <-
        apply(x, 2, function(v, row_mean) row_mean[order(order(v))], row_mean)

    ## propagate dimnames
    dimnames(result) <- dimnames(x)
    result
}

test_that("'quantile_normalize()' propagates dimnames", {
    m <- matrix(rnorm(10), 5, dimnames=list(LETTERS[1:5], letters[1:2]))
    observed <- quantile_normalize(m)
    expect_identical(dimnames(observed), dimnames(m))
})

We can also check that our function is robust to less common cases, like inputs with 0 rows and / or 0 columns. These tests lead to additional modifications – apply() does not always return a matrix!

quantile_normalize <-
    function(x)
{
    ## validate inputs
    x <- as.matrix(x)
    stopifnot(
        is.numeric(x),
        !anyNA(x)
    )

    ## quantile normalize
    m <- apply(x, 2, function(v) v[order(v)])
    dim(m) <- dim(x)    # apply() doesn't always return a matrix!

    row_mean <- rowMeans(m)

    result <- apply(
        x, 2, function(v, row_mean) row_mean[order(order(v))], row_mean
    )
    dim(result) <- dim(x)

    ## propagate dimnames
    dimnames(result) <- dimnames(x)
    result
}

test_that("'quantile_normalize()' works with edge cases", {
    m <- matrix(rnorm(5), nrow = 5)
    expect_identical(m, quantile_normalize(m))

    m <- matrix(rnorm(5), ncol = 5)
    expect_identical(matrix(mean(m), ncol = 5), quantile_normalize(m))

    m <- matrix(0, 0, 0)
    expect_identical(m, quantile_normalize(m))
})