hierarchical clustering

So far we have been looking at how to plot raw data, wrangle and summarize data, and work with edgelist data. This chapter continues the search for relationships between the samples in a data set, but now by sorting those samples into a branching tree of nested groups, and then deciding where to cut it.

linkage dendrograms

Consider the Alaska lakes dataset: which lake is most similar, chemically speaking, to Lake Narvakrak? Answering this requires calculating numeric distances between samples based on their chemical properties, just as we did for networks, and then building a tree from those distances. We can do it all in one step by using analysis = "hclust":

alaska_lake_data_wide <- tidyr::pivot_wider(dplyr::select(alaska_lake_data, -element_type), names_from = 'element', values_from = 'mg_per_L')

AK_lakes_clustered <- runMatrixAnalysis(
    data = alaska_lake_data_wide,
    analysis = c("hclust"),
    scale_variance = TRUE,
    agglomeration_method = "ward.D2",
    columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:15],
    columns_w_sample_ID_info = c("lake", "park")
)
AK_lakes_clustered
## # A tibble: 39 × 26
##    sample_unique_ID   lake  park  parent  node branch.length
##    <chr>              <chr> <chr>  <int> <int>         <dbl>
##  1 Devil_Mountain_La… Devi… BELA      31     1         0.973
##  2 Imuruk_Lake_BELA   Imur… BELA      39     2         0.820
##  3 Kuzitrin_Lake_BELA Kuzi… BELA      38     3         0.703
##  4 Lava_Lake_BELA     Lava… BELA      33     4         0.743
##  5 North_Killeak_Lak… Nort… BELA      22     5         4.09 
##  6 White_Fish_Lake_B… Whit… BELA      22     6         4.09 
##  7 Iniakuk_Lake_GAAR  Inia… GAAR      29     7         1.27 
##  8 Kurupa_Lake_GAAR   Kuru… GAAR      36     8         0.954
##  9 Lake_Matcharak_GA… Lake… GAAR      36     9         0.954
## 10 Lake_Selby_GAAR    Lake… GAAR      37    10         1.12 
## # ℹ 29 more rows
## # ℹ 20 more variables: label <chr>, isTip <lgl>, x <dbl>,
## #   y <dbl>, branch <dbl>, angle <dbl>, bootstrap <dbl>,
## #   water_temp <dbl>, pH <dbl>, C <dbl>, N <dbl>, P <dbl>,
## #   Cl <dbl>, S <dbl>, F <dbl>, Br <dbl>, Na <dbl>,
## #   K <dbl>, Ca <dbl>, Mg <dbl>

It works! Now we can plot our cluster diagram with a ggplot add-on called ggtree. We’ve seen that ggplot takes a “data” argument (i.e. ggplot(data = <some_data>) + geom_*() etc.). In contrast, ggtree takes an argument called tr, though with the output of the runMatrixAnalysis() function, these two (data and tr) can be treated the same, so, use: ggtree(tr = <output_from_runMatrixAnalysis>) + geom_*() etc.

Note that ggtree also comes with several great new geoms: geom_tiplab() and geom_tippoint(). Let’s try those out:

library(ggtree)
AK_lakes_clustered %>%
ggtree() +
  geom_tiplab() +
  geom_tippoint() +
  theme_classic() +
  scale_x_continuous(limits = c(0,10))

A first, unstyled dendrogram of the Alaskan lakes. A tree diagram in which each tip is one lake and each join is the merge of two groups, with the horizontal position of a join recording how far apart those two groups were when they merged; joins near the root therefore mark very different groups. Tip labels are the combined lake and park identifiers, and every tip is drawn at the same horizontal position because a hierarchical clustering dendrogram is ultrametric. Data are the 'alaska_lake_data' dataset used in UMD CHEM5725, clustered on thirteen scaled analytes with Ward linkage.

Figure 2.9: A first, unstyled dendrogram of the Alaskan lakes. A tree diagram in which each tip is one lake and each join is the merge of two groups, with the horizontal position of a join recording how far apart those two groups were when they merged; joins near the root therefore mark very different groups. Tip labels are the combined lake and park identifiers, and every tip is drawn at the same horizontal position because a hierarchical clustering dendrogram is ultrametric. Data are the ‘alaska_lake_data’ dataset used in UMD CHEM5725, clustered on thirteen scaled analytes with Ward linkage.

Cool! Though that plot could use some tweaking… let’s try:

AK_lakes_clustered %>%
ggtree() +
    geom_tiplab(aes(label = lake), offset = 0.2, align = TRUE) +
    geom_tippoint(shape = 21, aes(fill = park), size = 4) +
    scale_x_continuous(limits = c(0,12)) +
    scale_fill_brewer(palette = "Set1") +
    # theme_classic() +
    theme(
      legend.position = c(0.2,0.8)
    )

The same dendrogram, with lake names as tip labels and park identity shown on the tips. A tree diagram in which each tip is one of twenty Alaskan lakes, the tip label gives the lake name, and the fill color of the tip point gives the national park the lake sits in (BELA, GAAR or NOAT). Each join marks the merge of two groups and is positioned at the distance at which they merged. North Killeak and White Fish join only at the far right, which marks them as the two most distinctive lakes in the set. Data are the 'alaska_lake_data' dataset used in UMD CHEM5725, clustered on thirteen scaled analytes with Ward linkage.

Figure 2.10: The same dendrogram, with lake names as tip labels and park identity shown on the tips. A tree diagram in which each tip is one of twenty Alaskan lakes, the tip label gives the lake name, and the fill color of the tip point gives the national park the lake sits in (BELA, GAAR or NOAT). Each join marks the merge of two groups and is positioned at the distance at which they merged. North Killeak and White Fish join only at the far right, which marks them as the two most distinctive lakes in the set. Data are the ‘alaska_lake_data’ dataset used in UMD CHEM5725, clustered on thirteen scaled analytes with Ward linkage.

Very nice! Since North Killeak and White Fish are so different from the others, we could re-analyze the data with those two removed:

alaska_lake_data_wide %>%
  filter(!lake %in% c("North_Killeak_Lake","White_Fish_Lake")) -> alaska_lake_data_wide_filtered

runMatrixAnalysis(
    data = alaska_lake_data_wide_filtered,
    analysis = c("hclust"),
    scale_variance = TRUE,
    agglomeration_method = "ward.D2",
    columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide_filtered)[3:15],
    columns_w_sample_ID_info = c("lake", "park")
) %>%
ggtree() +
    geom_tiplab(aes(label = lake), offset = 0.2, align = TRUE) +
    geom_tippoint(shape = 21, aes(fill = park), size = 4) +
    scale_x_continuous(limits = c(0,9)) +
    scale_fill_brewer(palette = "Set1") +
    # theme_classic() +
    theme(
      legend.position = c(0.1,0.85)
    )
## Replacing NAs in your data with mean

The Alaskan lakes re-clustered with the two most distinctive lakes removed. A tree diagram of the same form as the previous figure but with eighteen tips: North Killeak and White Fish were dropped before the distances were recomputed. The tip label gives the lake name and the tip fill color gives the national park. Removing the two outliers rescales the horizontal axis, which makes the structure among the remaining lakes readable. Data are the 'alaska_lake_data' dataset used in UMD CHEM5725, clustered on thirteen scaled analytes with Ward linkage.

Figure 2.11: The Alaskan lakes re-clustered with the two most distinctive lakes removed. A tree diagram of the same form as the previous figure but with eighteen tips: North Killeak and White Fish were dropped before the distances were recomputed. The tip label gives the lake name and the tip fill color gives the national park. Removing the two outliers rescales the horizontal axis, which makes the structure among the remaining lakes readable. Data are the ‘alaska_lake_data’ dataset used in UMD CHEM5725, clustered on thirteen scaled analytes with Ward linkage.

concept check

Fill in the blank with the analysis that builds a tree by step-by-step merging. Press Run. This one takes a few seconds, and the object it makes is used again later in the chapter.

AK_tree <- runMatrixAnalysis(
  data = alaska_lake_data_wide,
  analysis = "____",
  scale_variance = TRUE,
  agglomeration_method = "ward.D2",
  columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:15],
  columns_w_sample_ID_info = c("lake", "park")
)
nrow(AK_tree)
sum(AK_tree$isTip)

Self-check: a colleague re-runs this clustering with agglomeration_method = “single” instead of “ward.D2” and gets a visibly different tree, from the same data and the same distances. How is that possible?

It is not; a different tree means the distances must have been recomputed differently The two rules measure the distance to a group differently, so they merge the same distance matrix in a different order Single linkage discards samples that sit too far from any group, which changes the shape of the tree

Fill in the blank with the column to sort by, so that the closest pair of lakes comes to the top. Press Run.

runMatrixAnalysis(
  data = alaska_lake_data_wide,
  analysis = "dist",
  output_format = "long",
  scale_variance = TRUE,
  columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:15],
  columns_w_sample_ID_info = c("lake", "park")
) %>%
  arrange(____) %>%
  head(2)

Self-check: the pair at the top of that table is the first merge the algorithm makes, and it is the same first merge under single, complete, average and Ward linkage. Why do all four agree here?

Because Ward’s method runs first and the other three copy its opening move Because the two lakes sit in the same national park Because every group still holds a single sample, so distance-to-a-group is just distance-between-samples and the linkage rule has nothing to choose between yet

agglomeration methods

Hierarchical clustering builds its tree from the leaves to the root using a very simple procedure:

Start with every sample in a group of its own.
Find the two groups that are closest to each other, and merge them into one group.
Repeat step 2 until everything is in a single group.

Each merge becomes a join in the tree, and the position of the join (a node in the tree) records how far apart the two groups were when they merged. Branches that join close to the tips were very similar. Branches that only join near the root were very different. This merging procedure is what runMatrixAnalysis() does by default, and it is the only tree-building rule we use in this chapter. (For completeness: the function also takes a tree_method argument, which defaults to "linkage_dendrogram" — the procedure just described. Its one alternative, "neighbor_joining", is a tree-building rule borrowed from phylogenetics that does not merge groups a pair at a time, and so ignores agglomeration_method entirely. We meet it properly in the phylogenies chapter.)

Step 2 hides a question, though. At the start, every group is a single sample, so the distance between two groups is simply the distance between two samples, read straight from the distance matrix. But after the first merge, how far is a group of two samples from a third sample? Or from a group of five? There is no single right answer. The rule chosen is called the agglomeration method, or linkage, and different rules build different trees from exactly the same distance matrix. In practice, the choice of linkage often changes the tree more than the choice of distance measure does. (We will meet other ways of measuring distance when we get to embeddings). Either way, the two extremes are the easiest to understand:

Single linkage: the distance between two groups is the distance between their closest members. Two groups merge as soon as any one sample in the first is near any one sample in the second.
Complete linkage: the distance between two groups is the distance between their furthest members. Two groups only merge once every sample in the first is reasonably close to every sample in the second.

Let’s see what difference this makes, using the solvents dataset. We will cluster the solvents on six physical properties. (We leave out vapor_pressure, because it is missing for five of the solvents.) Everything is identical between the two runs except the agglomeration method:

solvent_properties <- c(
  "boiling_point", "melting_point", "density",
  "relative_polarity", "formula_weight", "refractive_index"
)

solvents_single <- runMatrixAnalysis(
    data = solvents,
    analysis = c("hclust"),
    scale_variance = TRUE,
    agglomeration_method = "single",
    columns_w_values_for_single_analyte = solvent_properties,
    columns_w_sample_ID_info = c("solvent", "category")
)

solvents_complete <- runMatrixAnalysis(
    data = solvents,
    analysis = c("hclust"),
    scale_variance = TRUE, 
    agglomeration_method = "complete",
    columns_w_values_for_single_analyte = solvent_properties,
    columns_w_sample_ID_info = c("solvent", "category")
)

Now let’s plot the two trees side by side:

single_plot <- ggtree(solvents_single) +
    geom_tiplab(aes(label = solvent), offset = 0.05, size = 3) +
    geom_tippoint(shape = 21, aes(fill = category), size = 3) +
    scale_x_continuous(limits = c(0, 2.3)) +
    scale_fill_brewer(palette = "Dark2") +
    theme(legend.position = "none")

complete_plot <- ggtree(solvents_complete) +
    geom_tiplab(aes(label = solvent), offset = 0.1, size = 3) +
    geom_tippoint(shape = 21, aes(fill = category), size = 3) +
    scale_x_continuous(limits = c(0, 6)) +
    scale_fill_brewer(palette = "Dark2") +
    theme(legend.position = "right")

plot_grid(single_plot, complete_plot, nrow = 1, rel_widths = c(1, 1.3), labels = c("A", "B"))

Single and complete linkage build very different trees from identical distances. Two dendrograms of the same thirty-two solvents, in which each tip is one solvent, the tip fill color gives the solvent category recorded in the dataset, and each join is positioned at the distance at which the two groups it joins merged. A) Single linkage, which merges on the closest pair of members and produces a lopsided, staircase-shaped tree in which water, acetic acid and carbon disulfide are left to join at the very end. B) Complete linkage, which merges only once all members are close and produces four compact, chemically recognizable groups. Data are the 'solvents' dataset used in UMD CHEM5725, clustered on six scaled physical properties.

Figure 2.12: Single and complete linkage build very different trees from identical distances. Two dendrograms of the same thirty-two solvents, in which each tip is one solvent, the tip fill color gives the solvent category recorded in the dataset, and each join is positioned at the distance at which the two groups it joins merged. A) Single linkage, which merges on the closest pair of members and produces a lopsided, staircase-shaped tree in which water, acetic acid and carbon disulfide are left to join at the very end. B) Complete linkage, which merges only once all members are close and produces four compact, chemically recognizable groups. Data are the ‘solvents’ dataset used in UMD CHEM5725, clustered on six scaled physical properties.

The two trees were built from identical distances, but they tell very different stories. Single linkage produces a lopsided, staircase-like tree. There are a few small clusters (the small alcohols, the chlorinated solvents), but they hang one after another off a single long backbone, and water, acetic acid, and carbon disulfide are left to join on their own at the very end. Cutting this tree into four groups would give one giant group and three lone solvents. Complete linkage gives compact groups that a chemist would recognize: the small polar protic solvents (water, methanol, ethanol, the propanols, acetic acid, and acetonitrile); the volatile, low-boiling solvents (ether, acetone, ethyl acetate, THF, pentane and hexane); the heavier, higher-boiling solvents (the aromatics, DMF, pyridine, octanol and so on); and the dense chlorinated solvents, which carbon disulfide joins.

Why do we observe those differences? The reason is that solvent properties don’t fall into neat, separate boxes. They vary gradually from one solvent to the next. Single linkage follows chains of near neighbors: if A is close to B, and B is close to C, then A and C end up in the same group even if A and C are very different. This is called chaining, and on data that form a continuum it strings nearly everything into one long group. Complete linkage refuses to merge two groups until all of their members are close, so it cuts the same continuum into compact chunks.

Notice also the colors, which show the dataset’s category column. Neither tree reproduces those categories, and that is not a failure. The categories record functional groups, but we clustered on physical properties, for example, chlorobenzene behaves like the other aromatics, not like dichloromethane, so that is where it lands.

Single and complete linkage give the most different trees when:

the samples form a continuum rather than distinct groups. Single linkage chains along the gradient, giving a lopsided tree. Complete linkage chops it into compact pieces, giving a more balanced one. This is what happened with the solvents.
a few intermediate samples bridge two real groups. Single linkage uses the bridge to join the two groups early. Complete linkage keeps them apart until late, because their furthest members are still far apart.
one large, spread-out group sits next to tight ones. Complete linkage tends to split the big group into pieces the size of the tight ones. Single linkage follows its actual shape.
there are outliers. Under complete linkage, one extreme sample drags a group’s furthest-member distance upwards and delays that group’s merges. Single linkage simply leaves the outlier on a long branch of its own that joins at the end.

Single and complete linkage are two extremes, and there are actually two widely used methods sit between those extremes: Average linkage ("average", also called UPGMA) uses the average of all the distances between the members of the two groups. Ward’s method ("ward.D2"), which merges whichever pair of groups adds the least to the spread within groups. It behaves like a stricter version of complete linkage, with a strong preference for compact groups of similar size.

Last, and importantly: when the groups in the data are compact and well separated, every method gives essentially the same tree. So running two methods is a check on a result: if the trees from two different methods agree, the grouping is robust; if they disagree, the data are more of a continuum than a set of clusters, and the “groups” in any one tree should be treated with caution. No matter what, when presenting a figure that contains a tree, it is essential to state which agglomeration method was used, and, ideally, the results from multiple agglomeration methods and whether they agree.

concept check

Fill in the blank with the function that returns the smallest value, so that each solvent is paired with the distance to its nearest neighbor. Press Run.

runMatrixAnalysis(
  data = solvents,
  analysis = "dist",
  output_format = "long",
  scale_variance = TRUE,
  columns_w_values_for_single_analyte = solvent_properties,
  columns_w_sample_ID_info = c("solvent", "category")
) %>%
  group_by(solvent_sample_1) %>%
  summarize(nearest_neighbor = ____(distance)) %>%
  arrange(desc(nearest_neighbor)) %>%
  head(4)

Self-check: carbon disulfide, acetic acid, and water top that table, and under single linkage they are also the last three solvents to join. What does single linkage do with the other twenty-nine?

Chains them into one large group, because each of them has some close neighbor and single linkage merges on the closest pair of members Splits them into four compact groups, one for each solvent category in the dataset Leaves them unassigned until a cut height is chosen

dendrogram interpretation

Read across, not down. (assuming root->leaves goes left to right or right to left) All of the information is in the horizontal direction. Each node/join sits at the distance at which its two groups merged, so the horizontal distance from the tips back to a node/join indicates how far apart those groups were. Joins close to the tips are groups that merged early because they were alike; joins far out toward the root merged late, because whatever those nodes brought together was not so alike.

Compare branch lengths, rather that interpreting absolute branch lengths. These plots are drawn without a numeric axis, and that is deliberate rather than an omission. The horizontal spacing is a faithful rescaling of the distances, so every comparison read off the picture is exactly right: which of two joins is deeper, and by what factor. One group merging at twice the distance of another really does look twice as far out from the tips.

The vertical direction means nothing at all. (again, assuming root->leaves goes left to right or right to left) The tips are spread evenly down the page so the labels fit in the allocated space and can be read. There is no quantity on the y-axis. Worse, the order itself is arbitrary: at every join, either branch may be drawn above the other, and both drawings are the same tree. A dendrogram can be flipped at any of its joins and still mean exactly what it meant before.

Two tips sitting next to each other vertically are not necessarily similar. In the Alaska tree above, Lake Kangilipak and Okoklik Lake are neighbors on the page, and they are indeed the closest pair in the whole data set, 0.78 apart. But White Fish Lake and Wild Lake are also neighbors on the page too, and they sit 7.57 apart: nearly ten times as far, and in the most distant fifth of all the pairs in the data. Their lines do not meet until the very last join in the tree. So to judge how related two tips are, ignore how close together they are printed. Trace left from each of them until the two paths meet. That meeting point is their join, and how far out it sits is the answer.

One way to deploy these interpretation methods is to use a dendrogram to make a statement about an unknown sample. For example, suppose that we had a measurements from an unknown lake:

unknown_lake <- data.frame(
  lake = "unknown_lake", park = "unknown_park",
  water_temp = 8.00, pH = 7.55, C = 2.2, N = 0.000, P = 0.001, Cl = 0.65, S = 0.35, F = 0.00,
  Br = 0.00, Na = 0.92, K = 0.14, Ca = 1.65, Mg = 0.23
)
unknown_lake
##           lake         park water_temp   pH   C N     P
## 1 unknown_lake unknown_park          8 7.55 2.2 0 0.001
##     Cl    S F Br   Na    K   Ca   Mg
## 1 0.65 0.35 0  0 0.92 0.14 1.65 0.23

… and maybe we want to know “which lake in our Alaska Lakes data set is most similar to the unknown lake?”. One method of course would be a distance matrix, but we can also visualize the answer with a dendrogram. For this, we do need an important function: rbind(), which binds data together, row-wise. It only works if the items being bound all have the same set of columns. Like this:

rbind(unknown_lake, alaska_lake_data_wide)
##                   lake         park water_temp   pH    C
## 1         unknown_lake unknown_park       8.00 7.55  2.2
## 2  Devil_Mountain_Lake         BELA       6.46 7.69  3.4
## 3          Imuruk_Lake         BELA      17.38 6.44  4.7
## 4        Kuzitrin_Lake         BELA       8.06 7.45  2.0
## 5            Lava_Lake         BELA      20.18 7.42  8.3
## 6   North_Killeak_Lake         BELA      11.34 8.04  4.3
## 7      White_Fish_Lake         BELA      12.05 7.82 12.3
## 8         Iniakuk_Lake         GAAR       9.10 7.01  3.3
## 9          Kurupa_Lake         GAAR       9.30 7.03  2.1
## 10      Lake_Matcharak         GAAR      10.20 6.95  5.1
## 11          Lake_Selby         GAAR      15.10 7.15  4.2
## 12      Nutavukti_Lake         GAAR      17.60 6.88  4.5
## 13         Summit_Lake         GAAR      11.90 6.45  2.4
## 14       Takahula_Lake         GAAR       9.90 6.88  2.7
## 15         Walker_Lake         GAAR      15.30 7.22  1.3
## 16           Wild_Lake         GAAR       5.50 6.98  6.5
## 17    Desperation_Lake         NOAT       2.95 6.34  2.1
## 18         Feniak_Lake         NOAT       4.51 7.24  1.8
## 19     Lake_Kangilipak         NOAT       5.36 6.56  8.5
## 20      Lake_Narvakrak         NOAT      18.30 7.31  5.8
## 21        Okoklik_Lake         NOAT       6.46 6.87  7.8
##        N     P     Cl     S    F   Br     Na    K    Ca
## 1  0.000 0.001   0.65  0.35 0.00 0.00   0.92 0.14  1.65
## 2  0.028 0.000  10.35  0.62 0.04 0.02   8.92 1.20  5.73
## 3  0.013 0.000   1.18  0.20 0.02 0.00   1.36 0.32  0.96
## 4  0.000 0.000   0.67  0.29 0.01 0.00   0.94 0.19  1.70
## 5  0.017 0.001   2.53  0.59 0.04 0.01   2.93 0.57 11.77
## 6  0.037 0.001 337.23  0.04 0.11 1.08 161.89 8.32 29.75
## 7  0.034 0.006 105.30  0.18 0.05 0.33  54.16 3.48 18.62
## 8  0.141 0.000   0.22 12.09 0.00 0.00   0.57 0.44 40.78
## 9  0.043 0.000   0.13 12.43 0.00 0.00   4.51 0.31 14.68
## 10 0.000 0.000   1.25 13.30 0.00 0.00   5.13 0.74 36.82
## 11 0.107 0.000   0.11  7.92 0.00 0.00   0.46 0.10 12.91
## 12 0.000 0.001   0.18  2.72 0.00 0.00   0.84 0.21 11.78
## 13 0.000 0.001   0.08  3.21   NA   NA   1.07 0.17  4.69
## 14 0.014 0.000   0.23  5.53 0.00 0.00   0.62 0.93 57.64
## 15 0.190 0.001   0.19  5.77 0.00 0.00   0.57 1.18 27.68
## 16 0.130 0.001   0.31 28.68 0.00 0.00   1.40 0.40 48.17
## 17 0.005 0.000   0.20  2.73 0.01 0.00   1.11 0.16  5.87
## 18 0.000 0.000   0.21  4.93 0.01 0.00   0.72 0.17  6.81
## 19 0.005 0.000   0.55  0.55 0.02 0.00   2.20 0.35  6.51
## 20 0.000 0.000   0.76  1.38 0.02 0.00   0.82 0.56 11.52
## 21 0.000 0.000   0.76  1.01 0.02 0.01   1.97 0.48  8.74
##       Mg
## 1   0.23
## 2   4.04
## 3   0.86
## 4   0.26
## 5   2.77
## 6  37.68
## 7  16.16
## 8  11.44
## 9   8.74
## 10  8.56
## 11  2.81
## 12  2.08
## 13  1.86
## 14  7.49
## 15  3.50
## 16 16.42
## 17  3.03
## 18  7.66
## 19  2.63
## 20  1.84
## 21  2.57

Then we can run the hierarchical clustering analysis (below). Nice! Looks like the unknown lake is most simlar to Kuzitrin Lake.

runMatrixAnalysis(
    data = rbind(unknown_lake, alaska_lake_data_wide),
    analysis = c("hclust"),
    scale_variance = TRUE,
    agglomeration_method = "ward.D2",
    columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:15],
    columns_w_sample_ID_info = c("lake", "park")
) %>%
ggtree() +
    geom_tiplab(aes(label = lake), offset = 0.2, align = TRUE) +
    geom_tippoint(shape = 21, aes(fill = park), size = 4) +
    scale_x_continuous(limits = c(0,9)) +
    scale_fill_brewer(palette = "Set1") +
    theme(
      legend.position = c(0.1,0.85)
    )
## Replacing NAs in your data with mean

concept check

Fill in the blank with the lake that sits next to White Fish Lake on the tree. Press Run, and compare the two distances that come back — both pairs are neighbors on the page.

runMatrixAnalysis(
  data = alaska_lake_data_wide,
  analysis = "dist",
  output_format = "long",
  scale_variance = TRUE,
  columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:15],
  columns_w_sample_ID_info = c("lake", "park")
) %>%
  filter((lake_sample_1 == "Lake_Kangilipak" & lake_sample_2 == "Okoklik_Lake") |
         (lake_sample_1 == "White_Fish_Lake" & lake_sample_2 == "____")) %>%
  select(lake_sample_1, lake_sample_2, distance)

Self-check: in the Alaska tree, Lake Kangilipak sits next to Okoklik Lake and White Fish Lake sits next to Wild Lake. Both pairs are neighbors on the page. What does that tell you about the two pairs?

Both pairs are similar, because neighboring tips are always each other’s closest match Nothing — vertical position carries no quantity, and these two pairs are in fact the most similar and the least similar in the data The first pair is similar and the second is not, which can be read from how far apart they are printed

cutting a dendrogram

A dendrogram shows how everything in a data set relates to everything else, but it does not, by itself, say which samples belong to which group. Getting groups out of a tree means cutting across it at a chosen position (often called a “height”, which I think is confusing since trees are often drawn with root on the left and tips on the right, but anyway…) Every join “below” the cut has already happened, so the samples under it stay together; every join above the cut is severed, so the branches it would have merged stay apart. The height is the real control, and the number of groups is what falls out of it: cut low, near the tips, and there are many small groups; cut high, near the root, and there are a few large ones. runMatrixAnalysis() cuts the tree when it is given a parameters argument — the same argument k-means takes in the next chapter — and adds a cluster column to the output, one label per sample. Asking for a height looks like this:

solvents_cut <- runMatrixAnalysis(
    data = solvents,
    analysis = c("hclust"),
    scale_variance = TRUE,
    agglomeration_method = "complete",
    parameters = c(height = 5),
    columns_w_values_for_single_analyte = solvent_properties,
    columns_w_sample_ID_info = c("solvent", "category")
)

solvents_cut %>%
  filter(isTip) %>%
  group_by(cluster) %>%
  summarize(n_solvents = n())
## # A tibble: 4 × 2
##   cluster   n_solvents
##   <chr>          <int>
## 1 cluster_1          7
## 2 cluster_2          8
## 3 cluster_3         13
## 4 cluster_4          4

Cutting at 5 leaves four groups, of 13, 8, 7 and 4 solvents. Note that cluster is filled in for the tips only: the internal nodes of a tree are joins, not samples, so they have no group to belong to and are left as NA. This is the same convention the bootstrap column follows in the other direction, where it is the tips that are NA. Coloring the tips by that column shows what the cut actually did:

ggtree(solvents_cut) +
    geom_tiplab(aes(label = solvent), offset = 0.1, size = 3) +
    geom_tippoint(shape = 21, aes(fill = cluster), size = 3) +
    geom_cut(height = 5) +
    scale_x_continuous(limits = c(0, 6)) +
    scale_fill_brewer(palette = "Set2") +
    theme(legend.position = c(0.15, 0.85))

The solvent dendrogram cut at a height of 5. A tree diagram in which each tip is one of thirty-two solvents, the tip label gives the solvent name, and the fill color of the tip point gives the group it falls into when the tree is cut at that height. The dashed vertical line is the cut itself: the three joins to its left were severed, leaving the four groups, and every join to its right had already happened and holds its group together. Data are the 'solvents' dataset used in UMD CHEM5725, clustered on six scaled physical properties with complete linkage.

Figure 2.13: The solvent dendrogram cut at a height of 5. A tree diagram in which each tip is one of thirty-two solvents, the tip label gives the solvent name, and the fill color of the tip point gives the group it falls into when the tree is cut at that height. The dashed vertical line is the cut itself: the three joins to its left were severed, leaving the four groups, and every join to its right had already happened and holds its group together. Data are the ‘solvents’ dataset used in UMD CHEM5725, clustered on six scaled physical properties with complete linkage.

The dashed line drawn by geom_cut() is the cut. Everything to its right has already merged, so each color is one unbroken piece of the tree; the three joins to its left were severed, which is what left four groups rather than one. Count them on the figure: severing three joins always leaves four pieces. The groups here are the compact families complete linkage found — the dense chlorinated solvents on their own, the polar protic solvents together, and so on.

geom_cut() is worth a word, because the obvious thing to reach for does not work. The horizontal axis ggtree draws is not the cut height: it measures distance from the root, and runs to only half the height of the tree, so geom_vline(xintercept = 5) puts the line nowhere near the cut. geom_cut() does that conversion, which is why it takes the same height that was passed to parameters — write the height once, and the picture and the cluster column are guaranteed to agree. When the tree was cut by number of groups instead, geom_cut(k = 4) draws the line in the matching place.

Moving the cut changes the answer, and the whole point of cutting by height is that the number of groups is something the data hands back, not something chosen in advance:

n_groups_at <- function(cut_at) {
  runMatrixAnalysis(
      data = solvents,
      analysis = c("hclust"),
      scale_variance = TRUE,
      agglomeration_method = "complete",
      parameters = c(height = cut_at),
      columns_w_values_for_single_analyte = solvent_properties,
      columns_w_sample_ID_info = c("solvent", "category")
  ) %>%
    filter(isTip) %>%
    pull(cluster) %>%
    unique() %>%
    length()
}

data.frame(height = seq(2,6,0.1)) %>%
  mutate(n_groups = sapply(height, n_groups_at)) %>%
  ggplot(aes(x = height, y = n_groups)) + geom_point() + geom_line()

The number of groups is not something chosen in advance -- it depends on where the tree is cut. A scatter plot with a connecting line, in which the x axis is the height at which the solvent dendrogram was cut and the y axis is the number of groups that cut leaves. Each point is one cut, taken every 0.1 units from a height of 2 to a height of 6. The curve falls in steps rather than smoothly, and the wide flat stretch from about 4.5 to 5.9 is the stability that argues for reporting four groups. Data are the 'solvents' dataset used in UMD CHEM5725, clustered on six scaled physical properties with complete linkage.

Figure 2.14: The number of groups is not something chosen in advance – it depends on where the tree is cut. A scatter plot with a connecting line, in which the x axis is the height at which the solvent dendrogram was cut and the y axis is the number of groups that cut leaves. Each point is one cut, taken every 0.1 units from a height of 2 to a height of 6. The curve falls in steps rather than smoothly, and the wide flat stretch from about 4.5 to 5.9 is the stability that argues for reporting four groups. Data are the ‘solvents’ dataset used in UMD CHEM5725, clustered on six scaled physical properties with complete linkage.

Fifteen groups at a height of 2, nine at 3, six at 4, four at 5, and only two at 6. Notice that the our number of groups at four is stable over a wide stretch: anywhere from a height of about 4.5 to 5.9 gives the same four groups. That stability is the argument for reporting four. A gap like that in the tree means the next merge happens a long way above the ones before it, so four groups is not a delicate result that a slightly different cut would overturn. A number of groups that only survives a narrow band of heights is one to be suspicious of.

When the number of groups is already known (for example, because the data are clearly telling us so, or because a downstream step needs a fixed number), that number can be asked for directly instead, and runMatrixAnalysis() finds the height that delivers it:

runMatrixAnalysis(
    data = solvents,
    analysis = c("hclust"),
    scale_variance = TRUE,
    agglomeration_method = "complete",
    parameters = c(n_clusters = 4),
    columns_w_values_for_single_analyte = solvent_properties,
    columns_w_sample_ID_info = c("solvent", "category")
) %>%
  filter(isTip) %>%
  group_by(cluster) %>%
  summarize(n_solvents = n())
## # A tibble: 4 × 2
##   cluster   n_solvents
##   <chr>          <int>
## 1 cluster_1          7
## 2 cluster_2          8
## 3 cluster_3         13
## 4 cluster_4          4

The same four groups, because a height of 5 and a request for four groups are two ways of describing one cut. The difference is which one is the assumption. Asking for four groups will always return four, even from data with no group structure at all; cutting at a height can return one group, or twenty, and that is information. The next chapter takes up the question of how to choose a number of groups when there is no tree to read a height off.

concept check

Fill in the blank with the height that leaves only two groups. Press Run.

runMatrixAnalysis(
  data = solvents,
  analysis = "hclust",
  scale_variance = TRUE,
  agglomeration_method = "complete",
  parameters = c(height = ____),
  columns_w_values_for_single_analyte = solvent_properties,
  columns_w_sample_ID_info = c("solvent", "category")
) %>%
  filter(isTip) %>%
  group_by(cluster) %>%
  summarize(n_solvents = n())

Self-check: cutting at 5 gave groups of 13, 8, 7 and 4. Cutting at 6 gives 28 and 4. What does raising the cut do?

It re-runs the clustering with different settings, so every solvent may be reassigned It splits the larger groups into smaller ones, because a higher cut crosses more joins It can only merge groups, never split them — three of the four joined into one group of 28, and the group of 4 was left untouched

annotating dendrograms

Overlaying sample traits on a ggtree-based plot is straightforward when we combine ggtree with ggplot2. We begin by running the hierarchical clustering analysis and keeping its output for later plotting.

hclust_out <- runMatrixAnalysis(
  data = chemical_blooms,
  analysis = c("hclust"),
  scale_variance = TRUE,
  agglomeration_method = "ward.D2",
  columns_w_values_for_single_analyte = colnames(chemical_blooms)[2:10],
  columns_w_sample_ID_info = "label"
)

The object returned by runMatrixAnalysis() already contains the branch coordinates, so we can pass it directly to ggtree() and add tip labels while keeping the tree readable. Note that we should deliberately control the y-axis here, that will be key for aligning the tree with other plots later.

tree_plot <- ggtree(hclust_out) +
  geom_tiplab(size = 2, align = TRUE) +
  scale_x_continuous(limits = c(0, 10)) +
  scale_y_continuous(limits = c(0, 80)) +
  theme_classic()
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
tree_plot

A dendrogram of the chemical bloom samples, drawn on a deliberately fixed y scale. A tree diagram in which each tip is one of seventy-eight samples, the tip label is the sample identifier, and each join is positioned at the distance at which the two groups it joins merged. The y axis is fixed to the range 0 to 80 so that this panel can later be aligned tip-for-tip with a heat map of the same samples. Data are the 'chemical_blooms' dataset used in UMD CHEM5725, clustered on nine scaled compound classes with Ward linkage.

Figure 2.15: A dendrogram of the chemical bloom samples, drawn on a deliberately fixed y scale. A tree diagram in which each tip is one of seventy-eight samples, the tip label is the sample identifier, and each join is positioned at the distance at which the two groups it joins merged. The y axis is fixed to the range 0 to 80 so that this panel can later be aligned tip-for-tip with a heat map of the same samples. Data are the ‘chemical_blooms’ dataset used in UMD CHEM5725, clustered on nine scaled compound classes with Ward linkage.

Next, reshape the tip-level measurements to long form so each chemical becomes its own column of tiles. Because we reuse the y coordinate supplied by ggtree, the tiles inherit the same vertical order as the tips in the tree. Note that we remove the other columns in the hclust output for simplicity - they are only needed if we want to draw the full tree. Note that we also control the y-axis here to make sure it has the same bounds (limits) as the tree we made previously.

heat_plot <- hclust_out %>%
  filter(isTip) %>%
  select(-parent, -node, -branch.length, -label, -isTip, -x, -branch, -angle, -bootstrap) %>%
  pivot_longer(cols = 3:11, names_to = "chemical", values_to = "abundance") %>%
  ggplot(aes(x = chemical, y = y, fill = abundance)) +
  geom_tile() +
  scale_y_continuous(limits = c(0, 80)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))
heat_plot

Compound abundances for the same samples, in the tree's own tip order. A heat map in which each column is one of nine compound classes, each row is one sample, and fill color gives that compound class's abundance in that sample. The y axis reuses the y coordinate supplied by ggtree rather than the sample name, so the rows appear in the order the tips appear in the dendrogram. Data are the 'chemical_blooms' dataset used in UMD CHEM5725.

Figure 2.16: Compound abundances for the same samples, in the tree’s own tip order. A heat map in which each column is one of nine compound classes, each row is one sample, and fill color gives that compound class’s abundance in that sample. The y axis reuses the y coordinate supplied by ggtree rather than the sample name, so the rows appear in the order the tips appear in the dendrogram. Data are the ‘chemical_blooms’ dataset used in UMD CHEM5725.

With matching y scales, plot_grid() can align the tree and the heat map so the tiles line up with the corresponding samples. Using align = "h" snaps them together horizontally, and axis = "tb" keeps the panel heights consistent.

plot_grid(tree_plot, heat_plot, axis = "tb", align = "h", labels = c("A", "B"))

The dendrogram and the heat map aligned. A two-panel figure combining the tree and the compound heat map on a shared y scale, so that each row of tiles sits beside the tip it belongs to and blocks of chemically similar samples can be read off directly. A) The dendrogram, each tip one sample. B) The heat map, each column one compound class and fill color giving abundance. Data are the 'chemical_blooms' dataset used in UMD CHEM5725.

Figure 2.17: The dendrogram and the heat map aligned. A two-panel figure combining the tree and the compound heat map on a shared y scale, so that each row of tiles sits beside the tip it belongs to and blocks of chemically similar samples can be read off directly. A) The dendrogram, each tip one sample. B) The heat map, each column one compound class and fill color giving abundance. Data are the ‘chemical_blooms’ dataset used in UMD CHEM5725.

Note: if we were to instead build the heat map directly from the raw chemical_blooms table, the rows fall back to their alphabetical order and the heat map no longer matches the dendrogram ordering:

chemical_blooms %>%
  pivot_longer(cols = 2:10, names_to = "chemical", values_to = "abundance") %>%
  ggplot(aes(x = chemical, y = label, fill = abundance)) +
  geom_tile() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

What happens when the heat map is built from the raw table instead. A heat map of the same nine compound classes and the same seventy-eight samples, but with the y axis mapped to the sample name, which orders the rows alphabetically. The blocks visible in the previous figure are gone, because the row order now has nothing to do with the clustering. Data are the 'chemical_blooms' dataset used in UMD CHEM5725.

Figure 2.18: What happens when the heat map is built from the raw table instead. A heat map of the same nine compound classes and the same seventy-eight samples, but with the y axis mapped to the sample name, which orders the rows alphabetically. The blocks visible in the previous figure are gone, because the row order now has nothing to do with the clustering. Data are the ‘chemical_blooms’ dataset used in UMD CHEM5725.

concept check

Using AK_tree from the first concept check, fill in the blank with the coordinate that puts the tiles in the tree’s own tip order. Press Run.

AK_tree %>%
  filter(isTip) %>%
  select(lake, y, Ca, Mg, Na, K) %>%
  pivot_longer(cols = c(Ca, Mg, Na, K), names_to = "element", values_to = "mg_per_L") %>%
  ggplot(aes(x = element, y = ____, fill = mg_per_L)) +
  geom_tile()

Self-check: mapping y to lake instead of to y also draws a perfectly good-looking heat map. Why is that version wrong?

Because lake names are text, and geom_tile needs a numeric axis Because the lake column is dropped when the clustering is run Because it orders the rows alphabetically, which has nothing to do with the order ggtree placed the tips in, so the tiles no longer line up with the tree

exercises

The Confluence

Hierarchical clustering · merge heights + cutting the tree + group summaries

You came down into the red canyon to survey its springs, and the flood season has come early. Climb the ladder to the high springs to learn which waters are kin, noting the height of every confluence as you pass, then come back down and read the shape of the land itself before the water climbs to meet you.

The Register of the Gods

Hierarchical clustering · distances + cutting the tree + group summaries

The keeper of the temple is dead, and two registers died with him: which offering each shrine takes, and which god each shrine serves. Rebuild them both from the residues left in the vessels before the sun reaches its zenith, light the ring of fires at the altar, and cross the canopy bridge to carry the word.

Integrated Bioanalytics

hierarchical clustering

linkage dendrograms

concept check

agglomeration methods

concept check

dendrogram interpretation

concept check

cutting a dendrogram

concept check

annotating dendrograms

concept check

exercises

further reading