Analyze C. glabrata time course expression

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2     ✓ purrr   0.3.4
✓ tibble  3.0.1     ✓ dplyr   1.0.0
✓ tidyr   1.1.0     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
── Conflicts ────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

Load data

tab <- read_csv("../data/Ex009_experiment_set_up_20171019.csv")

Parsed with column specification:
cols(
  ALT_ID = col_character(),
  Sample = col_character(),
  Strain = col_character(),
  Genotype = col_character(),
  Replicate = col_double(),
  Timepoint = col_character()
)

dat <- read_tsv("../data/Ex009_normalized_log2_read_counts.zip")

Parsed with column specification:
cols(
  .default = col_double(),
  gene = col_character()
)
See spec(...) for full column specifications.

Prepare data

In this analysis we will focus on the “wt” samples,

# choose the samples to use and convert the time points to numeric variable
time2num <- c("pre" = 0, "20m" = 20, "30m" = 30, "45m" = 45, "60m" = 60,
               "90m" = 90, "120m" = 120, "150m" = 150, "180m" = 180, 
               "240m" = 240)
samples.use <- tab %>% 
  filter(Genotype == "wt", Timepoint != "del80") %>% 
  mutate(Time = time2num[Timepoint]) %>% 
  select(Sample, Time)

# extract the list of sample names
wt.list <- samples.use$Sample

# subset the data
dat1 <- dat %>% select(gene, all_of(wt.list))

# there are technical duplicates for each time point. for our purpose
# we just need one value per time point. It is sensible to compute the 
# mean value for each gene at each time point. to do this we first convert
# the data table into a long format, which allows for aggregation functions
dat1.long <- dat1 %>% 
  pivot_longer(starts_with("S"), names_to = "sample") %>% 
  left_join(samples.use, by = c("sample" = "Sample")) %>% 
  select(gene, time = Time, exn = value)

# calculate mean value for each timepoint within each gene
dat1.aggr <- dat1.long %>% 
  group_by(gene, time) %>% 
  summarize(avg.exn = mean(exn), sd.exn = sd(exn), .groups = "drop_last") %>%  
  # subtract the value of the first timepoint from the 
  # rest to form relative expression level (baseline)
  mutate(rel.exn = avg.exn - first(avg.exn))

# convert the long format back to the wide format for extracting matrix
dat1.wide <- dat1.aggr %>% 
  pivot_wider(id_cols = gene, names_from = time, values_from = avg.exn)

dat2.wide <- dat1.aggr %>% 
  pivot_wider(id_cols = gene, names_from = time, values_from = rel.exn)

# convert the result into a matrix for downstream analysis
m.dat1 <- as.matrix(dat1.wide[,-1])
rownames(m.dat1) <- dat1.wide$gene

m.dat2 <- as.matrix(dat2.wide[,-1])
rownames(m.dat2) <- dat2.wide$gene

Calculate distance (dissimilarity) between genes

Our goal is to group genes based on their temporal profile, that is, genes that share the same temporal pattern, such as increase in expression over time, should be grouped together.

For example, here are some genes known to be induced after starvation:

# plot examples
exp.list <- paste0("CAGL0",c("B02475g", "F02145g", "K10868g", "J04202g"))
exp.label <- paste(exp.list, c("PHO84","PHM2","CTA1","HSP12"), sep = " ")
names(exp.label) <- exp.list
dat1.aggr %>% 
  filter(gene %in% exp.list) %>% 
  ggplot(aes(x = time, y = rel.exn)) + geom_line() + geom_point(shape = 0) +
  geom_errorbar(aes(ymin = rel.exn - 1.96*sd.exn, 
                    ymax = rel.exn + 1.96*sd.exn)) +
  facet_wrap(~gene, labeller = labeller(gene = exp.label)) + 
  xlab("time (min)") + ylab("log2 normalized mRNA counts - time_0") +
  labs(caption = "points are average of at least two biological replicates\nerror bars are 95% confidence intervals") +
  theme(plot.caption = element_text(hjust = 0))

And here are genes that are not expected to respond to phosphate starvation:

# plot examples
exp.list <- paste0("CAGL0",c("K12694g", "K05005g", "D06138g", "D05170g"))
exp.label <- paste(exp.list, c("ACT1","ALG9","HEM2","PHO4"), sep = " ")
names(exp.label) <- exp.list
dat1.aggr %>% 
  filter(gene %in% exp.list) %>% 
  ggplot(aes(x = time, y = rel.exn)) + geom_line() + geom_point(shape = 0) +
  geom_errorbar(aes(ymin = rel.exn - 1.96*sd.exn, 
                    ymax = rel.exn + 1.96*sd.exn)) +
  facet_wrap(~gene, labeller = labeller(gene = exp.label)) + ylim(-5,5) +
  xlab("time (min)") + ylab("log2 normalized mRNA counts - time_0") +
  labs(caption = "points are average of at least two biological replicates\nerror bars are 95% confidence intervals") +
  theme(plot.caption = element_text(hjust = 0))

A toy example

In order to cluster genes based on this criteria, we need to first quantify the dissimilarity between genes. Here we are less interested in the “magnitude” of gene expression than the “profile” or temporal dynamics of gene expression. In the above example, notice that the absolute maximum expression levels of the four genes differ a lot. However, they all show the same “pattern”, i.e. an early rise followed by sustained expression. Let’s look at a simple example consisting of four genes and four time points:

toy <- matrix(c(2,4,4,4,6,6,3,2,20,40,40,40,45,45,35,20), 
              byrow = T, nrow = 4, 
              dimnames = list(gene = paste("gene", 1:4), 
                              time = c(10,20,30,40)))
toy

        time
gene     10 20 30 40
  gene 1  2  4  4  4
  gene 2  6  6  3  2
  gene 3 20 40 40 40
  gene 4 45 45 35 20

plot(x = 1, type = "n", xlim = c(10,50), ylim = c(0, 50))
t = c(10,20,30,40)
for(i in 1:4){
  points(t, toy[i,], pch = i)
  lines(t, toy[i,])
}

legend("topright", legend = rownames(toy), pch = 1:4)

Scaling to the rescue

Although the scales make it a bit difficult to see, but one can tell that genes 1 and 3 “rise” over time while genes 2 and 4 “fall” over time. The question is, how can we turn that intuition into a number that can be calculated from the data? One idea is to “scale” the data by subtracting the mean expression level across time of each gene from each time point, and divide the result by the standard deviation of the four time points. This will “center” and standardize the genes.

toy.scaled <- scale(t(toy)) # note that most r functions operate on columns
# since we want to operate on the time series, we transpose the matrix so
# that each gene is in a column
matplot(t, toy.scaled, type = "b", pch = 1:4, lty = 1, col = 1)
legend("topright", colnames(toy.scaled), pch = 1:4)

dist() function

Now the (dis)similarity patterns are more clear. But we still don’t have a number. We can now calculate the “Euclidian distance” between each pair of genes. In a 2D or 3D space, this measure would be the familiar “distance”. In higher dimensions such as here (each gene is a data point in a 4-D space spanned by their expression levels at each of the four time points), the idea would be the same, where \(D := \sqrt{\sum_{i=1}^{4}(x_i-y_i)^2}\). Using the dist() function, we found the result to be

dist(t(toy.scaled)) # dist() computes the distance between rows

          gene 1    gene 2    gene 3
gene 2 3.0652078                    
gene 3 0.0000000 3.0652078          
gene 4 2.9937354 0.5955187 2.9937354

Compare this to the result had we applied the function on the unscaled raw data

dist(toy)

         gene 1   gene 2   gene 3
gene 2  5.00000                  
gene 3 64.89992 64.53681         
gene 4 68.89848 66.25708 32.78719

Pearson’s Correlation Coefficient as another way to measure similarity

There is another way to quantify the dissimilarity between the expression patterns of the genes: calculate the Pearson’s correlation coefficient between the time series and subtract it from 1. The Pearson’s correlation coefficient describes the colinearity between two variables. Often written as “r”, it is closely related to linear regression. In fact, you may well have seen “r” or “\(R^2\)” written on plots with linear regression lines. Intuitively, two time series are perfectly positively correlated if they rise or fall together, and perfectly anti-correlated if they have the opposite pattern. In the example above, gene 1 and 3 would have a Pearson’s correlation coefficient of +1. Now let’s compute the results for all pairs:

Note

cor() computes the pairwise correlation between COLUMNS of a matrix and therefore we need to transpose the original matrix

cor(t(toy))

           gene 1     gene 2     gene 3     gene 4
gene 1  1.0000000 -0.5659165  1.0000000 -0.4937419
gene 2 -0.5659165  1.0000000 -0.5659165  0.9408929
gene 3  1.0000000 -0.5659165  1.0000000 -0.4937419
gene 4 -0.4937419  0.9408929 -0.4937419  1.0000000

Compare this to the result using the untransposed matrix, which shows the correlations between time points (this may be of interest as well).

cor(toy)

          10        20        30        40
10 1.0000000 0.8948817 0.7860212 0.5019804
20 0.8948817 1.0000000 0.9788613 0.8351485
30 0.7860212 0.9788613 1.0000000 0.9288416
40 0.5019804 0.8351485 0.9288416 1.0000000

Note that we don’t have to scale the original data because the calculation of correlation includes the scaling (see wiki)

Clustering genes

Armed with two different ways to measure the distance or similarity between genes, we can now move on to cluster the genes and visualize them with the aheatmap() function from the NMF package.

Use Euclidean distance

# we first use the Euclidean distance on the scaled matrix
# note that dist() operates on the rows
t(toy.scaled)

        time
gene             10        20         30        40
  gene 1 -1.5000000 0.5000000  0.5000000  0.500000
  gene 2  0.8488747 0.8488747 -0.6063391 -1.091410
  gene 3 -1.5000000 0.5000000  0.5000000  0.500000
  gene 4  0.7406129 0.7406129 -0.1058018 -1.375424
attr(,"scaled:center")
gene 1 gene 2 gene 3 gene 4 
  3.50   4.25  35.00  36.25 
attr(,"scaled:scale")
   gene 1    gene 2    gene 3    gene 4 
 1.000000  2.061553 10.000000 11.814539 

toy.dist1 <- dist(t(toy.scaled), method = "euclidean")
toy.dist1

          gene 1    gene 2    gene 3
gene 2 3.0652078                    
gene 3 0.0000000 3.0652078          
gene 4 2.9937354 0.5955187 2.9937354

toy.clust1 <- hclust(toy.dist1)
plot(toy.clust1)

Use correlation matrix

# remember that cor() operates on the columns
t(toy)

    gene
time gene 1 gene 2 gene 3 gene 4
  10      2      6     20     45
  20      4      6     40     45
  30      4      3     40     35
  40      4      2     40     20

toy.dist2 <- 1-cor(t(toy))
toy.dist2

         gene 1     gene 2   gene 3     gene 4
gene 1 0.000000 1.56591646 0.000000 1.49374193
gene 2 1.565916 0.00000000 1.565916 0.05910708
gene 3 0.000000 1.56591646 0.000000 1.49374193
gene 4 1.493742 0.05910708 1.493742 0.00000000

toy.clust2 <- hclust(as.dist(toy.dist2))
plot(toy.clust2)

Notice that the topologies are the same but the tree height units are different between the two. But that doesn’t concern us.

Before we apply what we learned to our actual gene expression dataset, let’s look at one more thing, that is, how heatmap() function works.

Dissecting the heatmap() function

Brief explanation

heatmap() actually involves a series of steps. First the function applies, by default, the hierarchical clustering algorithm on both the rows and the columns based on the input data (non-scaled). It then reorders the rows and the columns to match the order after the clustering. Note, however, the clustering is done on the unscaled data. This is especially confusing because the base R’s heatmap() function has an argument called scaled. This, however, only affects how the data are plotted, not how they are clustered. Note that this clustering step can be skipped by setting Rowv = NA and Colv = NA, each of which controls whether clustering and reordering is done for the rows and the columns respectively.

Once the rows and columns are reordered, unless disabled by the arguments shown above, the function then checks the value of scale. The default for this argument depends on the value of another argument, symm, whose default is false. If “symm” is false, scale defaults to “row”, i.e. scaling by row. Otherwise, scaling is off.

After reordering and scaling, if applicable, heatmap() then transposes the resulting matrix and plots it using the base function image(). We will see the reason for transposin the matrix later. Let’s first understand how image() works. This function simply draws rectangles and colors them by the value given. Note there is something nonintuitive about the mapping from the matrix to the image – the image() function draws a cartesian coordinate with (0,0) at the lower left corner. It then colors the square corresponding to (i,j) to the value of X[i,j] in the matrix. If you follows me so far, you realize that the while in the matrix the first index is the row and the second the column, in the image, the rows now become the x-axis and the columns the y-axis. Essentially, the matrix has been rotated 90 degrees counterclockwise.

Using the image() function

Now let’s look at the toy example to understand these intuitively. First, if we directly apply image() to the toy data, this is what we get

image(x = 1:4, y = c(10,20,30,40), z = toy, xlab = "gene", ylab = "time (min)", xaxt = "n", yaxt = "n")
axis(side = 1, at = 1:4, labels = rownames(toy))

axis(side = 2, at = c(10,20,30,40), labels = colnames(toy))

Next we use the heatmap() function without clustering or scaling

heatmap(toy, Rowv = NA, Colv = NA, scale = "none")

Once we imagine that the row labels be moved to the left side, we realize that the heatmap() function effectively swapped the rows and the columns, i.e. transposed the matrix. Presumably, this makes the rows and columns appear nearly the same as the input matrix. I say “nearly” because – notice that – the row orders were reversed, again because in a plot, (0,0) is in the lower left corner.

Adding the clustering based reordering

Now let’s add the clustering – note that we only want to cluster the genes, and would like to keep the time points in their original order, for obvious reasons. To achieve this, we realize that time series are in the columns. So we just need to set Colv = NA, while leaving Rowv = NULL, which means it uses the default clustering order.

# with row clustering
heatmap(toy, Rowv = NULL, Colv = NA, scale = "none")

Adding scaling doesn’t change the reordering

We found that the clustering is not what we have expected, i.e. gene 1/3 and gene 2/4. This is because by default heatmap() clusters the rows based on the unscaled values. So what if we set scale = "row"? Let’s see:

#heatmap with row scaling and clustering
heatmap(toy, Rowv = NULL, Colv = NA, scale = "row")

Implement the reordering step outside the heatmap() function

So now the plotted colors match the scaled values, but the clustering order is still the original! This shows that the scaling step is only for plotting, not for reordering. We can actually implement the ordering step outside of the heatmap() function

# below is how reordering is done inside the heatmap function
# recall that we have set toy.dist1 = dist(toy)
toy.dist1

          gene 1    gene 2    gene 3
gene 2 3.0652078                    
gene 3 0.0000000 3.0652078          
gene 4 2.9937354 0.5955187 2.9937354

# we can then cluster the genes using this distance matrix, which we stored as toy.clust1
plot(toy.clust1)

# we can then use this order in the heatmap
heatmap(toy, Rowv = as.dendrogram(toy.clust1), Colv = NA, scale = "row")

Now this is what we would have expected. This tells us that if we should think of heatmap and clustering as two steps, even though the heatmap() function in base R and many of its variations combine the two. This helps use clearly define

1. which variable do we want to cluster, the genes or the time series (or any other dimension)?
2. how do we want to measure the dissimilarity (distance) – by Euclidean disance or Pearson’s Correlation Coefficient?
3. how do we want to visualize the data, e.g. whether to reorder the rows or columns, whether to plot scaled or unscaled data etc. Note that scaling would remove the difference between genes in their absolute expression level, leaving only the temporal dynamics visible. Sometimes that may be exactly what we want, while in other cases that may not be desirable.

Apply to our data

plot the raw data, scaled or not

Since we are interested in genes that change their expression over the time course, it is natural to first look at the distribution of maximum change in expression for each gene over all time points, and then choose a cutoff and only select those genes with maximum change in expression greater than that cutoff.

# check distribution of changes
dat1.aggr %>% 
  group_by(gene) %>% 
  summarize(max.exn = max(abs(rel.exn))) %>% 
  ggplot(aes(max.exn)) + stat_ecdf(geom = "step") +
  scale_x_continuous(breaks = -1:15) +
  xlab("log2 maximum fold change in expression relative to time_0") +
  ylab("cumulative # of genes") +
  ggtitle("Cumulative distribution of maximum fold change in expression")

`summarise()` ungrouping output (override with `.groups` argument)

Let’s select the genes with a maximum fold change of 4 or above (log2 fold change of 2 or above)

# subset the data
morethan4fold <- apply(m.dat2, 1, function(x) any(abs(x) >= 2))
m.dat3 <- m.dat2[morethan4fold,]

Now we can plot the log2 fold changes relative to time 0 for these genes:

# we will use the pheatmap function from the pheatmap package, which is 
# similar but with better defaults than the base R's heatmap()
suppressPackageStartupMessages(library(pheatmap))
# nmf.options(grid.patch=TRUE)
colour <- colorRampPalette( c("blue", "black", "yellow"), space="rgb")(64) # for plotting
pheatmap(m.dat3, color = colour, breaks = seq(-8,8,length.out = 65), cluster_rows = F, cluster_cols = F, scale = "none", labels_row = "")

Now let’s try to cluster the genes using Pearson’s Correlation Coefficient. First, we will compute the correlation matrix and visualize it using the image() function, with reordering based on hierarchical clustering.

library(corrplot) # plot correlation matrix in heatmap format
dat3.cor <- cor(t(m.dat3), method = "pearson")
dat3.hc <- hclust(as.dist(1-dat3.cor), method = "complete")
od <- dat3.hc$order
colour <- colorRampPalette( c("blue", "black", "yellow"), space="rgb")(64)
corrplot(dat3.cor, col = colour, method = "color", order = "hclust", hclust.method = "ward.D", addrect = 4, tl.pos = "n", title = "Pairwise Pearson's Correlation Coefficient for 1519 genes with >4 fold change")

Next we apply the clustering order to the rows in the heatmap for gene expression.

pheatmap(m.dat3, color = colour, breaks = seq(-8,8,length.out = 65), cluster_rows = dat3.hc, cluster_cols = F, scale = "none", labels_row = "")

Another question is how similar are the time points with each other. We use similar approaches above, but this time applied to the time points as columns.

dat1.cor <- cor(m.dat1, method = "pearson")
dat1.hc <- hclust(as.dist(1-dat1.cor), method = "complete")
od1 <- dat1.hc$order
colour <- colorRampPalette( c("blue", "black", "yellow"), space="rgb")(64)
pheatmap(dat1.cor, col = colour, cluster_rows = F, cluster_cols = F, scale = "none")

The result is self-explanatory.