RPointCloud: CLL Clinical Data

Introduction

We want to illustrate the RPointCloud package (Version 0.8.0) with a clinical data set. Not surprisingly, we start by loading the package.

library(RPointCloud)

We also load several other useful packages (some of which may eventually get incorporated into the package requirements).

suppressMessages( library(Mercator) )
library(ClassDiscovery)
library(Polychrome)
data(Dark24)
data(Light24)
suppressMessages( library(igraph) )
suppressMessages( library("ape") )
suppressPackageStartupMessages( library(circlize) )

Now we fetch the sample clinical data set that is included with the package.

data(CLL)
ls()
## [1] "Dark24"    "Light24"   "clinical"  "daisydist" "oopt"      "ripdiag"
dim(clinical)
## [1] 266  29
colnames(clinical)
##  [1] "AgeAtDx"                    "Sex"                        "Race"                      
##  [4] "mutation.status"            "Light.chain.subtype"        "Zap70Protein"              
##  [7] "Rai.Stage"                  "CatRAI"                     "Serum.beta.2.microglobulin"
## [10] "LogB2M"                     "CatB2M"                     "Serum.LDH"                 
## [13] "LogLDH"                     "White.blood.count"          "LogWBC"                    
## [16] "CatWBC"                     "CatCD38"                    "Hypogammaglobulinemia"     
## [19] "Massive.Splenomegaly"       "Matutes"                    "Hemoglobin"                
## [22] "Platelets"                  "Prolymphocytes"             "stat13"                    
## [25] "stat11"                     "stat12"                     "stat17"                    
## [28] "Dohner"                     "Purity"

The clinical object is a numeric matrix containing the clinical data (binary values have been converted to 0-1; categorical values to integers). The daisydist is a distance matrix (stored as a dist object). Coombes and colleagues [J Biomed Inform. 2021 Jun;118:103788] showed that this is the best way to measure distances between mixed-type clinical data. The ripdiag object is a “Rips diagram” produced by applying the TDA algorithm to the daisy distances.

TDA Built-in Visualizations of the Rips Diagram

Here are some plots of the TDA results using tools from the original package. (I am not sure what any of these are really good for.)

diag <- ripdiag[["diagram"]]
opar <- par(mfrow = c(1,2))
plot(diag, barcode = TRUE, main = "Barcode")
plot(diag, main = "Rips Diagram")
Figure 1 : The Rips barcode diagram from TDA.

Figure 1 : The Rips barcode diagram from TDA. 

par(opar)
rm(opar)

Mercator Visualizations of the Underlying Data and Distance Matrix

Now we use our Mercator package to view the underlying data.

mercury <- Mercator(daisydist, metric = "daisy", method = "hclust", K = 8)
mercury <- addVisualization(mercury, "mds")
mercury <- addVisualization(mercury, "tsne")
mercury <- addVisualization(mercury, "umap")
mercury <- addVisualization(mercury, "som")
opar <- par(mfrow = c(3,2), cex = 1.1)
plot(mercury, view = "hclust")
plot(mercury, view = "mds", main = "Mult-Dimensional Scaling")
plot(mercury, view = "tsne", main = "t-SNE")
plot(mercury, view = "umap", main = "UMAP")
barplot(mercury, main = "Silhouette Width")
plot(mercury, view = "som", main = "Self-Organizing Maps")
Figure 2 : Mercator Visualizations of the distance matrix.

Figure 2 : Mercator Visualizations of the distance matrix. 

par(opar)
rm(opar)

Dimension Zero

Here is a picture of the “zero-cycle” data, which can also be used ultimately to cluster the points (where each point is a patient). The connected lines are similar to a single-linkage clustering structure, showing when individual points are merged together as the TDA parameter increases.

nzero <- sum(diag[, "dimension"] == 0)
cycles <- ripdiag[["cycleLocation"]][2:nzero]
W <- mercury@view[["umap"]]$layout
plot(W, main = "Connected Zero Cycles")
for (cyc in cycles) {
  points(W[cyc[1], , drop = FALSE], pch = 16,col = "red")
  X <- c(W[cyc[1], 1], W[cyc[2],1])
  Y <- c(W[cyc[1], 2], W[cyc[2],2])
  lines(X, Y)
}
Figure 3 : Hierarchical connections between zero cycles.

Figure 3 : Hierarchical connections between zero cycles.

Using iGraph

We can convert the 0-dimensional cycle structure into a dendrogram, by first passing them through the igraph package. We start by putting all the zero-cycle data together, which can be viewed as an “edge-list” from the igraph perspective.

edges <- t(do.call(cbind, cycles)) # this creates an "edgelist"
G <- graph_from_edgelist(edges)
G <- set_vertex_attr(G, "label", value = attr(daisydist, "Labels"))

Note that we attached the sample names to the graph, obtaining them from the daisy distance matrix. Now we use two different algorithms to decide how to layout the graph.

set.seed(2734)
Lt <- layout_as_tree(G)
L <- layout_with_fr(G)
opar <- par(mfrow = c(1,2), mai = c(0.01, 0.01, 1.02, 0.01))
plot(G, layout = Lt, main = "As Tree")
plot(G, layout = L, main = "Fruchterman-Reingold")
Figure 4 : Two igraph depictions of the zero cycle structure.

Figure 4 : Two igraph depictions of the zero cycle structure. 

par(opar)

Note that the Fruchterman-Reingold layout gives the most informative structure.

Community Structure

There are a variety of community-finding algorithms that we can apply. (Communities in grpah theory are similar to clusters in other machine learning areas of study.) “Edge-betweenness” seems to work best.

keg <- cluster_edge_betweenness(G) # 20
table(membership(keg)) 
## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
## 10 21 22 11 21 23  7  8  9  9 20  8 11 26 14 11  8 15  7  5
pal <- Dark24[membership(keg)]

The first line in the next code chunk shows that we did actually produce a tree. We explore three different ways ro visualize it

is.hierarchical(keg)
## Warning: `is.hierarchical()` was deprecated in igraph 2.0.0.
## ℹ Please use `is_hierarchical()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
## [1] TRUE
H <- as.hclust(keg)
H$labels <- attr(daisydist, "Labels")
K <-  7
colset <- Light24[cutree(H, k=K)]
G2 <- set_vertex_attr(G, "color", value = colset)
e <- 0.01
opar <- par(mai = c(e, e, e, e))
plot(G2, layout = L)
Figure 5 : Community structure, simplified.

Figure 5 : Community structure, simplified. 

par(opar)
P <- as.phylo(H)
opar <- par(mai = c(0.01, 0.01, 1.0, 0.01))
plot(P, type = "u", tip.color = colset, cex = 0.8, main = "Ape/Cladogram")
Figure 6 : Cladogram realization, from the ape package.

Figure 6 : Cladogram realization, from the ape package. 

par(opar)
rm(opar)

Visualizing Features

In any of the “scatter plot views” (e.g., MDS, UMAP, t-SNE) from Mercator, we may want to overlay different colors to represent different clinical features.

U <- mercury@view[["mds"]]
V <- mercury@view[["tsne"]]$Y
W <- mercury@view[["umap"]]$layout
featMU <- Feature(clinical[,"mutation.status"], "Mutation Status", c("cyan", "red"),
                  c("Mutated", "Unmutated"))
featRai <- Feature(clinical[,"CatRAI"], "Rai Stage", c("green", "magenta"), c("High", "Low"))
opar <- par(mfrow = c(1,2))
plot(W, main = "UMAP; Mutation Status", xlab = "U1", ylab = "U2")
points(featMU, W, pch = 16, cex = 1.4)
plot(W, main = "UMAP; Rai Stage", xlab = "U1", ylab = "U2")
points(featRai, W, pch = 16, cex = 1.4)
Figure 7 : UMAP visualizations with clinical features.

Figure 7 : UMAP visualizations with clinical features. 

par(opar)
rm(opar)

Significance

We have a statistical approach to deciding which of the detected cycles are statistically significant. Empirically, the persistence of 0-dimensional cycles looks like a gamma distribution, while the persistence of higher dimensional cycles looks like an exponential distribution. In both cases, we use an empirical Bayes approach, treating the observed distribution as a mixture of either gamma or exponential (as appropriate) with an unknown distribution contributing to heavier tails.

persistence <- diag[, "Death"] - diag[, "Birth"]

Zero-Cycles (Connected Components)

d0 <- persistence[diag[, "dimension"] == 0]
d0 <- d0[d0 < 1]
summary(d0)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.003652 0.054455 0.076682 0.084373 0.108333 0.227561
mu <- mean(d0)
nu <- median(d0)
sigma <- sd(d0)
shape <- mu^2/sigma^2
rate <- mu/sigma^2
xx <- seq(0, 0.23, length = 100)
yy <- dgamma(xx, shape = shape, rate = rate)
hist(d0, breaks = 123, freq = FALSE)
lines(xx, yy, col = "purple", lwd = 2)
Figure 8 : Histogram of the duration of zero cycles, with an overlaid gammm distibution.

Figure 8 : Histogram of the duration of zero cycles, with an overlaid gammm distibution.

One-Cycles (Loops)

Now we want to determine if there are significant “loops” in the data, and, if so, how many?

d1 <- persistence[diag[, "dimension"] == 1]
ef <- ExpoFit(d1) # should be close to log(2)/median? 
eb <- EBexpo(d1, 200)
opar <- par(mfrow = c(1,3))
plot(ef)
hist(eb)
plot(eb, prior = 0.56)
Figure 9 : Empirical Bayes detection of significant one-cycles.

Figure 9 : Empirical Bayes detection of significant one-cycles. 

par(opar)
rm(opar)
sum(d1 > cutoffSignificant(eb, 0.8, 0.56)) # posterior 80%, prior 0.56
## Warning in d1 > cutoffSignificant(eb, 0.8, 0.56): longer object length is not a multiple of
## shorter object length
## [1] 145
sum(d1 > 0.065) # post 90%
## [1] 1
which(d1 > 0.047)
## [1]  50  91  93 123 214 236 249
which(d1 > 0.065)
## [1] 236

Let’s pick out the most persistent 1-cycle.

cyc1 <- Cycle(ripdiag, 1, 236, "forestgreen")
cyc1@index
##       [,1] [,2]
##  [1,]   99  212
##  [2,]    4    7
##  [3,]    2   99
##  [4,]  110  212
##  [5,]   44   59
##  [6,]    4   75
##  [7,]    7   59
##  [8,]   75  110
##  [9,]    2   44

Each row represents an edge, by listing the IDs of the points at either end of the line segment. In this case, there are nine edges that link together to form a connected loop (or topological circle).

cyc2 <- Cycle(ripdiag, 1, 123, "red")
cyc3 <- Cycle(ripdiag, 1, 214, "purple")

opar <- par(mfrow = c(1, 3))
plot(cyc1, W, lwd = 2, pch = 16, col = "gray", xlab = "U1", ylab = "U2", main = "UMAP")
lines(cyc2, W, lwd=2)
lines(cyc3, W, lwd=2)

plot(U, pch = 16, col = "gray", main = "MDS")
lines(cyc1, U, lwd = 2)
lines(cyc2, U, lwd = 2)
lines(cyc3, U, lwd = 2)

plot(V, pch = 16, col = "gray", main = "t-SNE")
lines(cyc1, V, lwd = 2)
lines(cyc2, V, lwd = 2)
lines(cyc3, V, lwd = 2)
Figure 10 : Three views of three one-cycles.

Figure 10 : Three views of three one-cycles. 

par(opar)
rm(opar)
poof <- angleMeans(W, ripdiag, cyc1, clinical)
poof[is.na(poof)] <- 0
angle.df <- poof[, c("mutation.status", "CatB2M", "CatRAI", "CatCD38",
                     "Massive.Splenomegaly", "Hypogammaglobulinemia")]
colorScheme <- list(c(M = "green", U = "magenta"),
                    c(Hi = "cyan", Lo ="red"),
                    c(Hi = "blue", Lo = "yellow"),
                    c(Hi = "#dddddd", Lo = "#111111"),
                    c(No = "#dddddd", Yes = "brown"),
                    c(No = "#dddddd", Yes = "purple"))
annote <- LoopCircos(cyc1, angle.df, colorScheme)
image(annote)
Figure 11 : Circos plot of features varying around the most persistent cycle.

Figure 11 : Circos plot of features varying around the most persistent cycle.

Two-Cycles (Voids)

Now we want to determine if there are significant “voids” (empty interiors of spheres) in the data, and, if so, how many?

d2 <- persistence[diag[, "dimension"] == 2]
ef <- ExpoFit(d2) # should be close to log(2)/median? 
eb <- EBexpo(d2, 200)
opar <- par(mfrow = c(1, 3))
plot(ef)
hist(eb)
plot(eb, prior = 0.75)
Figure 12 : Empirical Bayes detection of significant two-cycles.

Figure 12 : Empirical Bayes detection of significant two-cycles. 

par(opar)
rm(opar)
sum(d2 > cutoff(0.8, 0.75, eb)) # posterior 80%, prior 0.56
## Warning in d2 > cutoff(0.8, 0.75, eb): longer object length is not a multiple of shorter object
## length
## [1] 87
sum(d2 > cutoff(0.95, 0.75, eb)) # posterior 90%, prior 0.56
## Warning in d2 > cutoff(0.95, 0.75, eb): longer object length is not a multiple of shorter
## object length
## [1] 82
cutoff(0.95, 0.75, eb)
## $low
## [1] 8.703499e-05
## 
## $high
## [1] 0.03263812
sum(d2 > 0.032) # post 90%
## [1] 2
which(d2 > 0.034)
## [1] 95
vd <- Cycle(ripdiag, 2, 95, "purple")
mds <- cmdscale(daisydist, k = 3)
voidPlot(vd, mds)
voidFeature(featMU, mds, radius = 0.011) # need to increase radius when you overlay one sphere on another
rgl::rglwidget()
#htmlwidgets::saveWidget(rglwidget(), "mywidget.html")
ob <- Projection(vd, mds, featMU, span = 0.2)
opar <- par(mfrow = c(1,2))
plot(ob)
image(ob, col = colorRampPalette(c("cyan", "gray", "red"))(64))
Figure 13 : Planar projection of mutation status around void.

Figure 13 : Planar projection of mutation status around void. 

par(opar)
rm(opar)
options(oopt)
#rm(list = ls())