【网络可视化】基于Louvain算法的网络集群发现

https://github.com/xieguigang/sciBASIC

Louvain算法是基于模块度的网络节点集群发现算法。该算法在效率和效果上都表现较好，并且能够发现层次性的网络节点集群结构，其优化目标是最大化整个网络集群模块的模块度（Modularity）。

我们在进行一些生物学网络数据可视化的时候，网络集群划分是一个比较频繁遇到的共通性的问题。例如：

• 将基因蛋白以及代谢物都联通在一起的网络中，我们可以分出哪些代谢途径模块；
• 对单细胞测序的结果进行分型分析，得到的细胞集群数据可以划分为哪些细胞分型
• 在宏基因组的数据处理结果中，我们可以将所检测到的环境样本中的细菌都划分为哪些相关性高的集群

Figure 2 Metavirome network. A visualization of the all-verses-all BLAST network. Viral groups are assigned a unique color and were visualized in Gephi. Separation of the viral groups was accomplished using Gephi’s Force Atlas 2 plugin, a force-directed layout algorithm. Ovals highlight viral groups whose members matched against seeded viruses in a parallel network analysis. Despite the large number of seeded viruses (1812), only 26 viruses were assigned to viral groups. All 26 viruses were crenarchaeal viruses. Edge connections are lines connecting members within a viral group and between members of differing viral groups. Care should be taken to not associate distance between viral groups as indicative of sequence similarity. To enhance clarity and distinctions between viral groups, only those containing 50 contigs or more are included in the figure. Viral groups numbering: group 0 = SIRV-1,2; group 23 = ASV-1, SSV-1-2, 4-9; group 26 = ATV; group 28 = AFV-1; group 29 = STIV-1,2; group 31 = AFV2-3, 6-9, SIFV; group 32 = STST-1,2 and ARSV-1, SRV.
Viral assemblage composition in Yellowstone acidic hot springs assessed by network analysis. The ISME Journal (2015) 9, 2162–2177; doi:10.1038/ismej.2015.28

上面的问题，我们都可以通过Louvain算法来回答。

Louvain算法

Louvain算法的思想很简单：

1. 将图中的每个节点看成一个独立的集群簇，这个时候集群的数目与节点个数相同；
2. 对每个节点i，依次尝试把节点i分配到其每个邻居节点所在的集群。计算分配前与分配后的模块度变化，并记录最大的那个邻居节点，如果，则把节点i分配最大的那个邻居节点所在的集群，否则保持不变；
3. 重复【2】，直到所有节点的所属集群簇不再变化；
4. 对图进行压缩，将所有在同一个集群簇中的节点压缩成一个新节点。集群内节点之间的边的权重转化为新节点的环的权重，集群间的边权重转化为新节点间的边权重；
5. 重复【1】直到整个图的模块度不再发生变化。

下面我们来一次看看上面每一个步骤是怎么实现为代码的吧：

首先，我们需要将网络中的每一个节点都看作为一个独立的集群簇，这一步实现起来非常的简单：

Sub setCluster0()
    cluster = New Integer(n - 1) {}

    For i As Integer = 0 To n - 1
        ' 一个节点一个簇
        cluster(i) = i
    Next
End Sub

接着呢，我们会尝试将第i个节点分配到其邻居节点所属的集群，并计算分配前与分配后的模块度变化：

Dim bestCluster = -1  ' 最优的簇号下标(先默认是自己)
Dim maxx_deltaQ = 0.0 ' 增量的最大值
Dim vis = New Boolean(n - 1) {}
Dim cur_deltaQ As Double

cluster_weight(cluster(i)) -= node_weight(i)
j = head(i)

While j <> -1
    ' l代表領接点的簇号
    Dim l = cluster(edge(j).v)

    If vis(l) Then
        ' 一个領接簇只判断一次
        Exit While
    Else
        vis(l) = True
    End If

    cur_deltaQ = edgeWeightPerCluster(l)
    cur_deltaQ -= node_weight(i) * cluster_weight(l) * resolution

    If cur_deltaQ > maxx_deltaQ Then
        bestCluster = l
        maxx_deltaQ = cur_deltaQ
    End If

    edgeWeightPerCluster(l) = 0
    j = edge(j).next
End While

在上面的代码之中，我们通过一个While循环遍历i所指向的节点的所有邻居节点，然后分别计算i变成了邻居节点所属的集群分类之后的模块度Q值。与原先的集群分类下的Q值作比较，得到了cur_deltaQ。如果这个cur_deltaQ比最大的Q值maxx_deltaQ要高的话，说明当前的集群移动是有利的。这个时候我们就记录下这个最大化的新集群编号到bestCluster之中

上面整个迭代迭代计算的过程可以通过下面的While循环来实现。我们一直计算整个移动过程，直到新的模块度Q值没有任何变化为止：

Do
    Dim i As Integer = order(point)
    point = (point + 1) Mod n

    If TryMoveNode(i) Then
        enum_time = 0
        update_flag = True
    Else
        enum_time += 1
    End If
Loop While enum_time < n

在完成了上面的一次循环计算过程之后，我们再对图进行构建，完成图压缩即可。上面的算法过程完整的实现代码大家可以阅读这个源代码文件：【Louvain.vb】

网络模块度（Modularity）的计算

我们在计算完网络节点集群之后，下面就需要计算一下我们的集群模块划分的模块度来评估我们的计算结果。模块度是评估一个网络集群模块划分好坏的度量方法，它的物理含义是集群模块内节点的连边数与随机情况下的边数的差，它的取值范围是 [−1/2,1)。我们也可以理解其是网络集群模块内部边的权重减去所有与模块节点相连的边的权重和。对于无向图我们更加容易理解：即模块内部边的度数减去模块内节点的总度数。模块度的定义如下：

我们将上面的公式转换为代码就是：

''' <summary>
''' compute the modularity of the network comminity
''' </summary>
''' <param name="g"></param>
''' <returns></returns>
Public Shared Function Modularity(g As NetworkGraph) As Double
    Dim m As Double = g.ComputeNodeDegrees.Values.Sum / 2
    Dim q As Double = 0
    Dim communitySet = GetCommunitySet(g)

    For Each entry In communitySet
        Dim Cset As Node() = entry.Value

        For i As Integer = 0 To Cset.Length - 1
            Dim u As Node = Cset(i)

            For j As Integer = 0 To Cset.Length - 1
                Dim v As Node = Cset(j)
                Dim auv As Double = 0

                If u.adjacencies.hasNeighbor(v) Then
                    auv = 1
                End If

                Dim ku = u.degree.In + u.degree.Out
                Dim kv = v.degree.In + u.degree.Out
                Dim subQ As Double = auv - ((ku * kv) / (2 * m))

                q += subQ
            Next
        Next
    Next

    Dim mov As Double = (1.0 / (2.0 * m)) * q
    Return mov
End Function

VisualBasic编程实例

我们进行网络节点集群聚类需要经过三个步骤，分别为：创建网络模型，执行计算，最后导出网络结果。在下面我给出的实例代码中，首先的For循环是我们进行网络模型创建的代码过程；之后呢，就可以通过一个公开的Communities.Analysis接口函数进行我们所构建的网络模型中的节点的聚类计算分析了；最后将网络模型转换为Tabular表格数据保存至文件即可导出分析结果。

Imports Microsoft.VisualBasic.Data.visualize.Network
Imports Microsoft.VisualBasic.Data.visualize.Network.Analysis
Imports Microsoft.VisualBasic.Data.visualize.Network.FileStream
Imports Microsoft.VisualBasic.Data.visualize.Network.Graph
Imports Microsoft.VisualBasic.Serialization.JSON

Dim links As String()() = ' load links data...
Dim g As New NetworkGraph

' build network via links
For Each line As String() In links
    If g.GetElementByID(line(0)) Is Nothing Then
        g.CreateNode(line(0))
    End If
    If g.GetElementByID(line(1)) Is Nothing Then
        g.CreateNode(line(1))
    End If

    g.CreateEdge(g.GetElementByID(line(0)), g.GetElementByID(line(1)), 1)
Next

' the original network with communities labeled
Dim clusters As NetworkGraph = Communities.Analysis(g)

Call Console.WriteLine(Communities.Modularity(clusters))
Call Console.WriteLine(Communities.Community(g).GetJson(indent:=True))

Call clusters _
    .Tabular _
    .Save("/path/to/save/folder")

在上面的代码实例之中，我们实际上所需要使用到的计算分析代码就只有一行：

' the original network with communities labeled
Dim clusters As NetworkGraph = Communities.Analysis(g)

整个的计算过程非常的简单简洁。我们在完成了网络节点集群的计算之后，将网络模型导出为文件，送到Cytoscape软件之中进行可视化验证我们的计算结果。可以看见，我们的代码所实现的网络节点集群聚类的效果是非常的好的：

从Cytoscape软件渲染出来的网络图我们可以看到，通过Louvain算法对我们的示例网络模型进行集群的划分结果是非常好的。每一个节点密集分布的簇中，相关的节点基本上都被划分到同一个集群之中了。通过对网络节点集群的Modularity值的计算结果为0.814316166040941。

• 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日