1. Introduction

In this practical, you will apply unsupervised learning techniques, such as clustering and principal component analysis (PCA) to investigate the structure of the data. Try to understand all the code. During the hands-on, you will find a series of short questions (Q) and tasks (T) that you will have to complete. Remember that they can appear as possible questions in the weekly questionaires.

Required libraries

Before starting, install and load the required libraries:

library(ggplot2)
library(ggrepel)
library(pheatmap)
library(factoextra)
library(cowplot)
library(scatterplot3d)

The factoextra R package makes easy to extract and visualize the output of exploratory multivariate data analyses, including principal component and cluster analyses, but also many others, such as correspondence analysis or factor analysis. Although here we will just use it to plot some of our results, it is a valuable resource worth to have a look at. You can learn more about it here.

2. Datasets

Gene expression in different cell types

Holmes et al. (2005) studied gene expression profiles of sorted T-cell populations: naïve (NAI) effector (EFFE) and memory (MEM) T cells; from 10 different subjects: 5 healthy (HEA) and 5 melanoma (MEL) patients. The columns are a subset of gene (actually the associated microarray probes) expression measurements, they correspond to 156 genes that show differential expression between cell types.

3. Study cell types based on their expression patterns

First, load the dataset and inspect its content:

load("tcells.RData")

dim(tcells) # Obtain the number of rows (n) and columns (p)
## [1]  30 156
tcells[1:5, 1:5] # Visualize the first rows and columns
##                   X3968     X14831      X13492      X5108     X16348
## HEA26_EFFE_1 -2.6108361 -1.1857923 -0.05612926 -0.1470299  0.5230758
## HEA26_MEM_1  -2.2592865 -0.4718562  0.27730991  0.5438521 -0.3659554
## HEA26_NAI_1  -0.2719081  0.8170308  0.81269167  0.7198425 -0.9049818
## MEL36_EFFE_1 -2.2431022 -1.0806009 -0.23557156 -0.1760595  0.6406782
## MEL36_MEM_1  -2.6798807 -0.1463671  0.24909402  0.9519808 -0.2040837

Note that with this large number of variables, you have it very difficult to visualize the dataset. In this scenario, it can be useful to apply principal component analysis (PCA):

tcells <- scale(tcells, center = T, scale = T)  # Center and scale our variables prior to PCA
pca <- prcomp(tcells)                           # Compute PCA

var_explained <- pca$sdev^2/sum(pca$sdev^2)*100 # Percentage of variance explained
head(var_explained)
## [1] 63.814332  9.622202  3.973310  3.618176  2.395551  2.130956
PCs <- as.data.frame(pca$x[, 1:2])              # First two PCs
head(PCs)
##                     PC1       PC2
## HEA26_EFFE_1 -9.5861230 -2.654871
## HEA26_MEM_1  -2.3930805  3.012799
## HEA26_NAI_1   9.4027407 -1.553217
## MEL36_EFFE_1 -9.1913487 -2.803891
## MEL36_MEM_1  -0.7298769  6.123905
## MEL36_NAI_1  10.6834050 -1.861107
ggplot(data = PCs, aes(PC1, PC2, label = rownames(PCs))) + 
  geom_point() +
  geom_text_repel() +
  theme_bw()

Q1: Do you identify some patterns in the PCA plot? Which ones?

Q2: Which is the percentage of variance explained by PC1? and PC2?

Let’s add the information of the cell type, the disease status and the individual id:

## Preprocess the rownames to get the cell type, the disease status and the individual id (don't worry about it)
l <- lapply(strsplit(rownames(PCs), "_"), function (x){c(substr(x[1], 1,3), x[2], as.numeric(substring(x[3],1,1)))})
info <- do.call(rbind, l)
colnames(info) <- c("status", "celltype", "id")
rownames(info) <- rownames(tcells)
##

PCs <- cbind(info, PCs)
head(PCs)
##              status celltype id        PC1       PC2
## HEA26_EFFE_1    HEA     EFFE  1 -9.5861230 -2.654871
## HEA26_MEM_1     HEA      MEM  1 -2.3930805  3.012799
## HEA26_NAI_1     HEA      NAI  1  9.4027407 -1.553217
## MEL36_EFFE_1    MEL     EFFE  1 -9.1913487 -2.803891
## MEL36_MEM_1     MEL      MEM  1 -0.7298769  6.123905
## MEL36_NAI_1     MEL      NAI  1 10.6834050 -1.861107
ggplot(data = PCs, aes(PC1, PC2, col = celltype, shape = status, label = id)) + 
  geom_point() +
  geom_text_repel() +
  theme_bw()

Q3: How many samples are not where they should? What makes you think so? Do you see any subgroup structure? (e.g. within the same group, the samples from the same individual are similar, the samples from the healthy individuals are similar, etc.)

By looking at the PCA, it seems pretty clear that there are three major groups. It is important that you remember that PCA is an unsupervised method, and therefore these groups have been found without any prior information on the celltype, the individual ID or the disease status.

Can we arrive to the same conclusion using K-means? Let’s try it out:

set.seed(123)

# Empty vector to store the results
res <- c()

# Run K-means for different values of k and store J (WSS)
for (k in 1:10){
  res <- rbind(res, cbind(k, J = kmeans(tcells, centers = k, nstart = 10)$tot.withinss))
}
head(res)
##      k        J
## [1,] 1 4524.000
## [2,] 2 2090.360
## [3,] 3 1431.513
## [4,] 4 1311.958
## [5,] 5 1197.656
## [6,] 6 1108.381
ggplot(as.data.frame(res), aes(k, J)) + 
  geom_path() + 
  geom_point() +
  scale_x_continuous(breaks = 1:10) +
  theme_bw()

Q4: Which number of clusters would you select by looking at the elbow plot? Why? What is the nstart parameter doing? Hint: check out ?kmeans.

Now let’s plot again the two first PCs, coloring by cell type, and setting a shape that corresponds to each cluster:

set.seed(123)
cluster <- kmeans(tcells, centers = 3, nstart = 10)$cluster
cluster <- paste0("cluster", cluster)
head(cluster)
## [1] "cluster2" "cluster3" "cluster1" "cluster2" "cluster3" "cluster1"
PCs <- cbind(cluster, PCs)

ggplot(PCs, aes(PC1, PC2, shape = cluster, col = celltype)) + 
  geom_point() +
  theme_bw()

Q5: Do the 3 clusters obtained by K-means match the grouping by cell type?. How many samples are not where they should? What makes you think so?

Now, let’s try to obtain a finer cluster structure so that we can see clearer whether the samples from the same individual or the same disease group look like each other. We can use hierarchical clustering for that:

D <- dist(tcells, method = "euclidean")    # n x n    
# Other distances: manhattan, canberra ...

hc <- hclust(D, method = "complete")        
# Other clustering methods: average, ward.D2 ...

fviz_dend(hc, cex = 0.4, k = 3, palette = "jco")

Here we are telling hierarchical clustering to cut the tree at the heigth it leaves three clusters, that we have colored. Again, note that this clusters have been obtained in an unsupervised manner, although they match approximately the three major cell types. Regarding the cluster substructure, note that samples from the same individual tend to cluster together.

Q6: Is there any grouping by disease status?

Note also that if we transpose our gene expression matrix: t(tcells), we can cluster also the genes:

tcells2 <- t(tcells)                       # (transposed) p x n
D2 <- dist(tcells2, method = "euclidean")  # p x p
hc2 <- hclust(D2, method = "complete")     

fviz_dend(hc2, cex = 0.5)

In this case we are not cutting the tree at any level.

T1: Use K-means to obtain the corresponding number of gene clusters. Hint: Simply repeat the code for computing the elbow plot, but now with t(tcells) instead of tcells as input. Which value of K would you select?

Often it is very useful to visualize the actual expression values (note that as we have centered and scaled, small values will be negative) using a heatmap, that incorporates both clusterings (samples and genes), and includes a color legend with additional information (e.g. the cell type). We can do that with pheatmap:

pheatmap(tcells, 
         annotation_row = as.data.frame(info), 
         show_colnames = F,
         cutree_rows = 3)

T2: Add the option cutree_cols = x, where \(x\) is the number of gene clusters that you selected in task T1. How is the expression of the genes in each gene cluster for each of the sample clusters? Fill the following table with ‘high’, ‘medium’ and ‘low’ (e.g. if you decide that there should be 3 gene clusters) according to the heatmap:

gene cluster 1 gene cluster 2 gene cluster 3
sample cluster 1 high/low high/low high/low
sample cluster 2 high/low high/low high/low
sample cluster 3 high/low high/low high/low

Interestingly, besides clustering algorithms, there are methods such as multidimensional scaling (aka Principal Coordinate Analysis, PCoA), which allow to represent graphically in a PCA-like manner a distance matrix. When the distance matrix is built using Euclidean distances, the result is identical to PCA. Check it out:

pcoa <- cmdscale(D, k = 2)
colnames(pcoa) <- c("PCo1", "PCo2")

head(pcoa)
##                    PCo1      PCo2
## HEA26_EFFE_1 -9.5861230 -2.654871
## HEA26_MEM_1  -2.3930805  3.012799
## HEA26_NAI_1   9.4027407 -1.553217
## MEL36_EFFE_1 -9.1913487 -2.803891
## MEL36_MEM_1  -0.7298769  6.123905
## MEL36_NAI_1  10.6834050 -1.861107
p1 <- ggplot(as.data.frame(pcoa), aes(PCo1, PCo2)) + 
  geom_point() +
  theme_bw() +
   ggtitle("PCoA (euclidean)")

D3 <- dist(tcells, method = "canberra")
pcoa <- cmdscale(D3, k = 2)
colnames(pcoa) <- c("PCo1", "PCo2")

p2 <- ggplot(as.data.frame(pcoa), aes(PCo1, PCo2)) + 
  geom_point() +
  theme_bw() +
  ggtitle("PCoA (canberra)")

p3 <- ggplot(data = PCs, aes(PC1, PC2)) + 
  geom_point() +
  theme_bw() +
  ggtitle("PCA")

plot_grid(p3, p1, p2, ncol = 3, labels = c("a", "b", "c"))

Back to PCA, for visualization purposes we often reduce dimensionality from \(p\) to 2 principal components, but remember that you could also visualize 3 principal components:

## Obtain the celltype colors (do not worry about this)
cl <-  c("#F8766D", "#00BA38", "#619CFF")
names(cl) <- c("EFFE", "MEM", "NAI")
cl <- cl[match(info[,"celltype"], names(cl)) ]
## 

scatterplot3d(pca$x[, 1:3], color = cl, pch = 19, angle = 60)

Finally, generate a biplot and try to interpret the results:

rownames(pca$rotation) <- paste0("gene", 1:nrow(pca$rotation)) # Change array probe names by gene1, gene2 ...
fviz_pca_biplot(pca, 
                geom = "point", 
                habillage = info[,"celltype"],    # Color and shape defined by celltype
                col.var = "violet",               # Variables (genes) shown as violet arrows
                title = NULL,                     # No title
                repel = T,                        # Improve visualization of overlapping labels
                select.var = list(contrib = 5))   # Only added to plot 5 genes

Q8: Answer the following:

Remember to answer the short questions (Q) and complete the tasks (T). They can appear as possible questions in the weekly questionaires.