Multi-view clustering

library(plyr)
library(dplyr)
library(data.table)
library(ggplot2)
library(pheatmap)
library(methods)
library(grid)
library(reshape2)
library(nlme)
library(stringr)
library(psych)
library(tibble)
library(cluster)

library(metautils)
library(ithi.meta)
library(clonecluster)
library(paperutils)

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db
cluster_label_file <- snakemake@input$cluster_labels
## For sample labels and plotting
data_matrix_file <- snakemake@input$data_matrix

Metadata

db <- src_sqlite(db_path, create = FALSE)
samples <- collect(tbl(db, "samples"))

Clusters

cluster_labels <- fread(cluster_label_file)
data_matrix <- fread(data_matrix_file)

mat <- as.matrix(cluster_labels)
rownames(mat) <- data_matrix$condensed_id

cluster_pct <- mat/rowSums(mat, na.rm = TRUE)
remove_ind <- apply(cluster_pct, 1, function(x) all(is.na(x)))
cluster_pct <- cluster_pct[!remove_ind, ]

Number of samples: 56.

Selecting k

cluster_counts <- 2:7
kmeans_results <- lapply(cluster_counts, function(k) {
    kmeans(cluster_pct, centers = k, iter.max = 500, nstart = 50)
})
silwidths <- sapply(kmeans_results, function(x) {
    mean(silhouette(x$cluster, dist(cluster_pct))[, "sil_width"])
})
withinss <- sapply(kmeans_results, function(x) {
    sum(x$withinss)
})

clustdf <- data.frame(nclust = cluster_counts, silhouette = silwidths, withinss = withinss)
clustdf_melted <- melt(clustdf, id.vars = "nclust", measure.vars = c("silhouette", 
    "withinss"), variable.name = "metric", value.name = "value")

ggplot(clustdf_melted, aes(x = nclust, y = value)) + geom_point() + theme_bw() + 
    theme_Publication() + facet_wrap(~metric, scales = "free")

It’s quite abundantly clear that 3 is the optimal number of clusters, by both the silhouette and elbow rule.

cluster_labels <- kmeans_results[[which(cluster_counts == 3)]]$cluster
cluster_labels_df <- rownames_to_column(data.frame(cluster_labels), "condensed_id")

Heatmap

data_matrix <- merge(data_matrix, cluster_labels_df, by = c("condensed_id"))
data_matrix$cluster_labels <- factor(data_matrix$cluster_labels)

data_matrix_subset <- subset(data_matrix, select = -c(condensed_id, cluster_labels))
datmat <- as.matrix(clip_values(scale(data_matrix_subset), 2, -2))
rownames(datmat) <- data_matrix$condensed_id

tmp <- scale(data_matrix_subset)
dists <- as.dist((1 - cor(tmp, use = "pairwise.complete"))/2)
hc <- hclust(dists, method = "ward.D")

annotation_row <- data.frame(data_matrix$cluster_labels)
rownames(annotation_row) <- data_matrix$condensed_id

annotation_colors <- list(cluster_labels = get_colour_palette(data_matrix$cluster_labels))

pheatmap(datmat[order(data_matrix$cluster_labels), hc$order], cluster_rows = FALSE, 
    cluster_cols = FALSE, annotation_row = annotation_row, annotation_colors = annotation_colors)

This doesn’t look that great …

Might have to do with the fact that the clonal analysis view isn’t accurate still (divergence will throw things off).

Plain hierarchical clustering

dists2 <- as.dist((1 - cor(t(tmp), use = "pairwise.complete"))/2)
hc2 <- hclust(dists2, method = "ward.D")
pheatmap(datmat[hc2$order, hc$order], cluster_rows = FALSE, cluster_cols = FALSE)