TIL clusters: multiple assignments

Setup

library(ithi.utils)
load_base_libs()

library(methods)
library(factoextra)
library(NbClust)
library(ComplexHeatmap)
library(cluster)
library(entropy)
library(mclust)

library(ithi.meta)
library(ithi.figures)
library(ithi.utils)

ihc_table_path <- snakemake@input$ihc_table
molecular_subtype_file <- snakemake@input$molsubtypes
ipynb_path <- snakemake@input$ipynb_path

## For python
til_clusters_output <- snakemake@params$til_clusters_output
ihc_features_output <- snakemake@params$ihc_features_output

tils_for_cluster <- snakemake@params$tils_for_cluster
db_path <- snakemake@params$db

annotation_colours <- ithi.figures::get_annotation_colours()
ihc_table <- fread(ihc_table_path)
molsubtypes <- fread(molecular_subtype_file)

ipynb_base_path <- basename(ipynb_path)

til_clusters <- ithi.figures:::get_til_clusters(ihc_table, molsubtypes, tils_for_cluster = tils_for_cluster, 
    nclusts = 3)

Analysis

This reviewer comment wants us to address the issue that – like HGSC molecular subtype assignments – TIL subtype assignments may overlap in some cases. In other words, how likely is it that a sample would be assigned to more than one TIL subtype?

In the case of transcriptional subtypes, this type of analysis was done in the cited papers (e.g. Verhaak et al.) with ssGSEA. This makes sense for expression data, where gene signatures can be defined for each molecular subtype. In the case of these TIL subtypes, however, this isn’t nearly as relevant – we aren’t looking at thousands of features (rather, we’re only looking at 8) – so defining signatures is tough.

There are a few approaches I can think of to address this:

• Post-hoc cluster assessment, e.g. with the silhouette method
• Unsupervised clustering
• Supervised clustering

Post-hoc assessment of clusters: silhouette method

tilheat <- ithi.figures::plot_til_density_heatmap(ihc_table, tils_for_cluster, 
    3, molsubtypes, annotation_colours)
res.hc <- as.hclust(ComplexHeatmap::column_dend(tilheat))

fviz_dend(res.hc, palette = "jco", rect = TRUE, show_labels = FALSE)

ihc_data <- subset(ihc_table, select = c("condensed_id", tils_for_cluster))
valid_rows <- apply(ihc_data, 1, function(x) all(!is.na(x)))
ihc_data <- as.data.frame(ihc_data[valid_rows, ])

ihc_data <- ihc_data %>% plyr::rename(c(E_CD8_density = "Epithelial CD8+", E_CD4_density = "Epithelial CD4+", 
    E_CD20_density = "Epithelial CD20+", E_Plasma_density = "Epithelial Plasma", 
    S_CD8_density = "Stromal CD8+", S_CD4_density = "Stromal CD4+", S_CD20_density = "Stromal CD20+", 
    S_Plasma_density = "Stromal Plasma"))
ihc_mat <- ihc_data %>% tibble::column_to_rownames(var = "condensed_id")
ihc_mat_scaled <- scale(ihc_mat)

ihc_mat_scaled <- ithi.utils::clip_values(ihc_mat_scaled, hi = 2, lo = -2)

sil.hc <- cluster::silhouette(cutree(res.hc, 3), dist(ihc_mat_scaled, method = "euclidean"))

til_clusters_table <- table(til_clusters$cluster)

til_clusters_table

 N-TIL  S-TIL ES-TIL 
    47     30     34 

plot(sil.hc)

So all N-TILs are confidently clustered, whereas some S-TIL and ES-TIL are not as confidently clustered according to the silhouette.

sil.hc.mat <- sil.hc
class(sil.hc.mat) <- "matrix"
sil_confusion_mat <- sil.hc.mat %>% as.data.frame %>% subset(sil_width < 0)

sil_confusion_table <- table(sil_confusion_mat$cluster)

So 11 out of 30 are negative silhouette, and 7 out of NA are negative silhouette. The neighbour cluster in all but one case is N-TIL (in one case, it’s S-TIL).

In total, the fraction of negative silhouette samples is:

sum(sil_confusion_table)/sum(til_clusters_table)

[1] 0.1621622

which is pretty good. But this rate is quite high in S-TIL.

Unsupervised clustering

Fuzzy clustering

We’ll do fuzzy c-means and see how well that recapitulates the hierarchical clustering results.

res.fanny <- fanny(ihc_mat_scaled, k = 3, metric = "euclidean", stand = FALSE)

fanny_clusters <- data.frame(res.fanny$clustering) %>% rownames_to_column(var = "condensed_id") %>% 
    plyr::rename(c(res.fanny.clustering = "fanny_cluster"))
fanny_clusters$fanny_cluster <- fanny_clusters$fanny_cluster %>% plyr::mapvalues(from = c(1:3), 
    to = c("N-TIL", "S-TIL", "ES-TIL"))

hc_fanny_clusters <- til_clusters %>% plyr::join(fanny_clusters)

with(hc_fanny_clusters, table(cluster, fanny_cluster))

        fanny_cluster
cluster  ES-TIL N-TIL S-TIL
  N-TIL       0    47     0
  S-TIL       9     4    17
  ES-TIL     31     0     3

with(hc_fanny_clusters, adjustedRandIndex(cluster, fanny_cluster))

[1] 0.6893836

with(hc_fanny_clusters, mi.empirical(table(cluster, fanny_cluster)))

[1] 0.6852884

So generally we’re doing pretty well on the consistency of our clustering between hierarchical and fuzzy c-means clustering.

fviz_cluster(res.fanny, ellipse.type = "norm", repel = TRUE, palette = "jco", 
    ggtheme = theme_minimal(), legend = "right")

fviz_silhouette(res.fanny, palette = "jco", ggtheme = theme_minimal())

  cluster size ave.sil.width
1       1   51          0.59
2       2   20          0.09
3       3   40         -0.02

The silhouette plot looks similar, albeit quite a bit worse for ES-TIL.

Let’s now plot the clustering probabilities.

pheatmap(res.fanny$membership, show_rownames = FALSE)

What this means is that fuzzy c-means is doing a terrible job distinguishing S-TIL from ES-TIL.

I think the silhouette stuff is probably sufficient to address this comment, but what do you guys think? The fuzzy c-means clustering is bad. Alternatively I can look into using NMF, etc, since the data is non-negative.

NMF

When I have time? Worth doing?

Supervised clustering

We can alternatively look at supervised methods – if we work from the known TIL clusters, can we distinguish N/S/ES-TIL with methods that either assign cluster membership probabilities, or with methods that can designate multiple labels (multilabel classification)?

samples_intersect <- intersect(rownames(ihc_mat_scaled), til_clusters$condensed_id)

X <- ihc_mat_scaled[samples_intersect, ]
y <- (til_clusters %>% column_to_rownames(var = "condensed_id"))[samples_intersect, 
    ]

write.table(X, file = ihc_features_output, sep = "\t", quote = FALSE, row.names = FALSE, 
    col.names = TRUE)
write.table(data.frame(cluster = as.numeric(y)), file = til_clusters_output, 
    sep = "\t", quote = FALSE, row.names = FALSE, col.names = TRUE)

Now, look in the iPython notebook for the results of this stuff.

Supervised clustering analysis

The summary of this analysis was that, with a simple multinomial-objective logistic regression model and 10-fold cross-validation, we achieve 100% validation accuracy. Furthermore, 80% of samples are classified with >80% probabillity into the correct cluster. Consistent with our other results above, S-TIL is the most difficult to cluster to call, with only 63% of samples classified with >80% probability. However, 86% of S-TIL samples can be correctly classified with >50% probability.

Summary

16% have negative silhouette scores (some S-TIL that could also be assigned to N-TIL is the most common misclassification). This is comparable to the % of negative silhouette scores for TCGA molecular subtypes (~17% or so). The conclusion from Roel’s paper was:

“The overlap in gene signature scores suggests that HGS-OvCa does not consist of mutually exclusive expression subtypes, but that each tumor sample is represented by multiple signatures at different levels of activation.”

Unsupervised fuzzy clustering methods, like fuzzy c-means, do poorly and are unable to distinguish S-TIL and ES-TIL at all. This is not necessarily surprising, as these are entirely different clustering methods that do not know the assigned true labels.

However, we show that with supervised clustering methods, more specifically logistic regression, we can get very competitive accuracy and confidence rates for our classes.