【数据可视化】可视化层次聚类树

https://github.com/xieguigang/sciBASIC

层次聚类通过计算不同类别数据点间的相似度来创建一棵有层次的嵌套聚类树。基于层次聚类分析，我们可以初步可视化我们的一些原始数据：

• 例如对样本的层次聚类分类，可以让我们了解到样本在分组之间以及分组内的异质性。
• 对生物序列进行基于相似度的层次聚类分析，我们可以了解到序列之间的相似性程度或者进化关系

如果我们初步乍一看层次聚类树的外观以及其数据结构，我们就会发现，层次聚类树和二叉树非常的相似。但是层次聚类树和二叉树是不一样的两样东西哦。层次聚类树和二叉树二者之间首先在构建的算法原理上(都是基于对象相似性或者等价性)虽然有些许的相似，但是也相差很远。层次聚类树的建立过程是依据距离的极值进行不断的合并或者分裂，最后所得到的树的形状是一致的。但是对于二叉树而言，一般情况下，最终所生成的二叉树的形状与最开始的root数据点有很大的关系。

在之前的博客文章《二叉树聚类》中所展示的可视化一棵基于二叉树的聚类结果。

层次聚类的VisualBasic代码实现

对于我们所进行研究的对象实体数据点，能被看成是一类的，所以要么它们距离近，要么它们或多或少有共同的特征。在最常使用的Kmeans聚类之中，我们就是基于欧几里得距离来对数据点进行自动分类。在层次聚类之中，也是可以通过欧几里得距离来进行对象的分类的划分。

无监督机器学习中的层次聚类算法。该算法基于嵌套簇的拆分和合并。根据计算的方向，层次聚类树的建立过程可以分为：自底向上（凝聚方法，agglomerative）或者自上而下（分裂方法，divisive）这两种方法过程：

自下而上

1. 初始化，每个样本当做一个簇
2. 计算任意两簇距离，找出 距离最近的两个簇，合并这两簇
3. 重复步骤【2】直到，最远两簇距离超过阈值，或者簇的个数达到指定值，终止算法

自上而下

1. 初始化，所有样本集中归为一个簇
2. 在同一个簇中，计算任意两个样本之间的距离，找到距离最远的两个样本点a,b，将a,b作为两个簇的中心;
3. 计算原来簇中剩余样本点距离a，b的距离，距离哪个中心近，分配到哪个簇中
4. 重复步骤【2】、【3】直到，最远两簇距离不足阈值，或者簇的个数达到指定值，终止算法

层次聚类

在层次聚类树中，有一个比较重要的距离概念。样本节点之间的距离就是度量合并集群的链接标准，存在有如下所示的一些距离计算方法：

' SingleLinkageStrategy
Public Function CalculateDistance(distances As ICollection(Of Distance)) As Distance
    Dim min As Double = Double.NaN

    For Each dist As Distance In distances
        If Double.IsNaN(min) OrElse dist.Distance < min Then 
            min = dist.Distance
        End If
    Next

    Return New Distance(min)
End Function

single: 最近邻, 把类与类间距离最近的作为类间距

' WeightedLinkageStrategy
Public Function CalculateDistance(distances As ICollection(Of Distance)) As Distance
    Dim sum As Double = 0
    Dim weightTotal As Double = 0

    For Each distance As Distance In distances
        weightTotal += distance.Weight
        sum += distance.Distance * distance.Weight
    Next

    Return New Distance(sum / weightTotal, weightTotal)
End Function

weighted: WPGMA（加权分组平均）法。

' CompleteLinkageStrategy
Public Function CalculateDistance(distances As ICollection(Of Distance)) As Distance
    Dim max As Double = Double.NaN

    For Each dist As Distance In distances
        If Double.IsNaN(max) OrElse dist.Distance > max Then 
            max = dist.Distance
        End If
    Next

    Return New Distance(max)
End Function

complete: 最远邻, 把类与类间距离最远的作为类间距

' AverageLinkageStrategy
Public Function CalculateDistance(distances As ICollection(Of Distance)) As Distance
    Dim sum As Double = 0
    Dim result As Double

    For Each dist As Distance In distances
        sum += dist.Distance
    Next

    If distances.Count > 0 Then
        result = sum / distances.Count
    Else
        result = 0.0
    End If

    Return New Distance(result)
End Function

average: 平均距离, 类与类间所有pairs距离的平均。UPGMA算法（非加权组平均）法

在拥有了距离计算方法之后，我们就可以进行层次聚类分析了。在这里我们来了解自下而上的层次聚类过程。假设，我们有如下数据结构的层次聚类的Cluster对象：

Public Class Cluster : Implements INamedValue

    Public Property Distance As Distance

    ''' <summary>
    ''' value of <see cref="Distance"/>
    ''' </summary>
    ''' <returns></returns>
    Public ReadOnly Property DistanceValue As Double
        Get
            Return Distance.Distance
        End Get
    End Property

    Public Property Parent As Cluster

    ''' <summary>
    ''' 名称是唯一的？
    ''' </summary>
    ''' <returns></returns>
    Public Property Name As String Implements INamedValue.Key
    Public ReadOnly Property Children As IList(Of Cluster)
    Public ReadOnly Property LeafNames As List(Of String)

    ''' <summary>
    ''' 是否是一个叶节点？
    ''' </summary>
    ''' <returns></returns>
    Public ReadOnly Property isLeaf As Boolean
        Get
            Return Children.Count = 0
        End Get
    End Property

    Public ReadOnly Property TotalDistance As Double
        Get
            Dim dist As Double = If(Distance Is Nothing, 0, Distance.Distance)
            If Children.Count > 0 Then
                dist += Children(0).TotalDistance
            End If
            Return dist
        End Get
    End Property

End Class

则可以有，在最开始的时候，我们是将每一个数据点都当作为一个聚类簇的，所以我们可以在最开始对所有的数据点产生有聚类簇：

Private Shared Function createClusters(clusterNames As String()) As IList(Of Cluster)
    Return clusterNames _
        .Select(Function(clusterName) New Cluster(clusterName)) _
        .AsList
End Function

接着呢，我们对初始的所有聚类簇进行两两比较计算出任意两个cluster之间的距离：

Dim linkages As New DistanceMap

For col As Integer = 0 To clusters.Count - 1
    For row As Integer = col + 1 To clusters.Count - 1
        Dim link As New HierarchyTreeNode
        Dim lCluster As Cluster = clusters(col)
        Dim rCluster As Cluster = clusters(row)

        link.LinkageDistance = distances(col)(row)
        link.Left = (lCluster)
        link.Right = (rCluster)

        linkages.Add(link)
    Next
Next

接下来，我们就需要找出距离最近的两个簇，合并为一个簇：

Dim minDistLink As HierarchyTreeNode = Distances.RemoveFirst()

Clusters.Remove(minDistLink.Right())
Clusters.Remove(minDistLink.Left())

Dim oldClusterL As Cluster = minDistLink.Left()
Dim oldClusterR As Cluster = minDistLink.Right()
Dim newCluster As Cluster = minDistLink.Agglomerate(Nothing)
Dim distanceValues As New List(Of Distance)

For Each iClust As Cluster In Clusters
    Dim link1 As HierarchyTreeNode = findByClusters(iClust, oldClusterL)
    Dim link2 As HierarchyTreeNode = findByClusters(iClust, oldClusterR)

    If link1 IsNot Nothing Then
        Dim distVal As Double = link1.LinkageDistance
        Dim weightVal As Double = link1.GetOtherCluster(iClust).WeightValue
        distanceValues.Add(New Distance(distVal, weightVal))
        Distances.Remove(link1)
    End If

    If link2 IsNot Nothing Then
        Dim distVal As Double = link2.LinkageDistance
        Dim weightVal As Double = link2.GetOtherCluster(iClust).WeightValue
        distanceValues.Add(New Distance(distVal, weightVal))
        Distances.Remove(link2)
    End If

    Dim newLinkage As New HierarchyTreeNode With {
        .Left = iClust,
        .Right = newCluster,
        .LinkageDistance = linkageStrategy _
            .CalculateDistance(distanceValues) _
            .Distance
    }

    Call distanceValues.Clear()
    Call Distances.Add(newLinkage, direct:=True)
Next

Call Distances.Sort()
Call Clusters.Add(newCluster)

在上面的代码之中，Distances对象是一个排序表，即每调用一次Distances.Add方法添加一条连接进去，就会自动按照链接两端的节点的距离排序。所以在Distances对象中，排序表中的第一个元素总是距离最近的节点链接。

在上面的代码之中，我们首先将Distances对象中距离最小的链接拿出来，将链接两端的节点合并到一个Cluster之中。接着将这个Cluster中的左右树枝节点从原来的聚类簇列表中删除（因为被合并了，取而代之的是一个新的Cluster）。重新计算整个Cluster列表中的链接距离。

至此，我们完成了自下而上的层次聚类算法中的步骤2。

最后，我们通过While循环不断的重复上面的算法过程，直到我们所建立的层次聚类树满足条件：

Do While Not builder.TreeComplete
    Call builder.Agglomerate(linkageStrategy)
Loop

可视化层次聚类树

在层次聚类树的绘制过程中，有两点是需要我们进行额外关注的：

1. 就是聚类簇节点之间的距离应该要体现在可视化作图结果之中。即Cluster节点在可视化图中的坐标位置是与节点之间的距离有关的。
2. 因为层次聚类结果是一个树形结构，所以我们会需要使用递归方法进行绘制，方便编写代码。

假设我们知道了每一个节点的父节点的坐标位置，则我们可以根据距离计算出当前所需要进行绘制的节点的坐标位置，计算代码如下：

Dim x = plotRegion.Left + plotRegion.Right - scaleX(partition.DistanceValue)

对于进行树形结构绘制的递归过程，我们可以有：

Dim n As Integer = 0

parentPt = New PointF(x, y)

For Each part As Cluster In orders
    DendrogramPlot(part, unitWidth, g, plotRegion, i + n, scaleX, parentPt, labelPadding, charWidth)
    n += part.Leafs
Next

在是上面所示的代码之中，parentPt是父节点的坐标位置。则每一次递归绘制都会递归的以parentPt为参考，根据距离计算出新的parentPt坐标点位置用于子节点的布局计算。现在通过上面的递归调用，树形结构布局已经可以被正常的计算出来了。接着我们来看一下节点的绘制：

If Not parentPt.IsEmpty Then
    ' 绘制连接线
    Call g.DrawLine(linkColor, parentPt, New PointF(parentPt.X, y))
    Call g.DrawLine(linkColor, New PointF(x, y), New PointF(parentPt.X, y))
End If

If (partition.isLeaf OrElse showAllNodes) AndAlso theme.PointSize > 0 Then
    Call g.DrawCircle(New PointF(x, y), theme.PointSize, pointColor)
End If

If showLeafLabels AndAlso (partition.isLeaf OrElse showAllLabels) Then
    Dim lsize As SizeF = g.MeasureString(partition.Name, labelFont)
    Dim lpos As New PointF(x + labelPadding, y - lsize.Height / 2)

    Call g.DrawString(partition.Name, labelFont, Brushes.Black, lpos)
End If

在递归调用的DendrogramPlot过程之中，我们在代码中会根据计算出来的节点布局位置，得到三个点（父节点，子节点，连接线转折点）。则使用这三个点可以通过DrawLine绘制出两个节点之间的连接线。对于节点的可视化，我们直接在计算出来的（x，y）位置上画一个圆点即可。当然我们也可以用饼图，条形图之类的替代圆点，进行更加高级的可视化图形的绘制。

上面所描述的完整的层次聚类树的可视化代码，大家可以阅读源代码文件：【DendrogramPanelV2.vb】。

R#脚本中可视化层次聚类树

在R#脚本之中，可以通过hclust函数进行层次聚类分析。对于hclust函数而言，其接受一个距离矩阵，然后依据这个距离矩阵进行层次聚类计算得到层次聚类树。在得到层次聚类树之后，我们可以直接将得到的树通过plot函数可视化出来：

imports "clustering" from "MLkit";
imports "visualPlots" from "graphics";

setwd(dirname(@script));

const sampleInfo       = read.csv("../../sampleInfo.csv");
const labels as string = sampleInfo[, "ID"];
const class as string  = `#${as.character(sampleInfo[, "color"])}`;
# by samples
const d = "../../metabolome.txt"
|> read.csv(
    row.names   = 1,
    tsv         = TRUE,
    check.names = FALSE
)
|> t
|> dist
|> hclust
;

print(d);

bitmap(file = "./hclust_samples.png", size = [2100, 3300]) {
    plot(d,
        class       = lapply(1:length(labels), i -> class[i], names = labels),
        padding     = "padding: 200px 400px 200px 200px;",
        axis.format = "G3",
        links       = "stroke: darkblue; stroke-width: 8px; stroke-dash: dash;",
        pt.color    = "gray",
        label       = "font-style: normal; font-size: 13; font-family: Bookman Old Style;",
        ticks       = "font-style: normal; font-size: 10; font-family: Bookman Old Style;"
    );
}

通过上面的层次聚类树可视化结果可以看到，QC样本可以很好的聚类在一块。这说明实验的系统误差比较好。CK分组的样本也可以聚类在一块，说明CK分组组内的差异比较小。对于A01分组而言，其在可视化结果之中被分裂为不同的分支，说明A01分组内的样本的重复一致性较低。

• Author
• Recent Posts

谢桂纲

高级数据科学家 at 苏州帕诺米克

Working on Engineered bacteria CAD design on its genome from scratch. Writing scientific computing software for Tianhe & Sunway TaihuLight supercomputer. Do scientific computing programming in R/R# language, he is also the programming language designer of the R# language on the .NET runtime.

Latest posts by 谢桂纲 (see all)

• 【MZKit】简单自动化组织分区 - 2023年11月5日
• 【MZKit教程】质谱成像原始数据文件查看 - 2023年6月29日
• 生物序列图嵌入算法 - 2023年6月29日