生物序列图嵌入算法

假若现在有两条Fasta序列放在你面前，现在需要你进行这两条Fasta序列的相似度计算分析。如果对于我而言，大学刚毕业刚入门生物信息学的时候，可能只能够想到通过blast比对的方式进行序列相似性计算分析。基于blast比对方式可以找到生物学意义上的序列相似性结果，但是计算的效率会比较低。假设现在让你使用这些序列进行机器学习建模分析，或者基于传统数学意义上的基于相似度的无监督聚类分析的时候，面对这些长度上长短不一的生物序列数据，可能会比较蒙圈，因为传统的数学分析方法都要求我们分析的目标至少应该是等长的向量数据。

使用图理论描述序列

在接触过图分析理论，对序列数据有一个比较深刻的认识之后，实现上面的分析目标，实际上我们可以将目标生物序列转换为图对象，然后基于图对象进行一维嵌入表示的方式将长度上长短不一的序列转换为等长的数值向量。最后基于得到的向量数据，就可以进行基于常规的数学方法做机器学习建模或者无监督的数据挖掘分析了。

生物序列图嵌入方法就是这样一种基于图论的方法将变长的序列数据转换为等长的向量的方法，其主要是通过图理论将不同长度的两条生物序列嵌入表示为两个等长的向量，当任意的生物序列都转换为等长的两个向量之后，我们就可以通过传统的聚类算法，机器学习算法来进行生物序列数据的分析，例如无监督聚类，机器学习分类。

那这个方法的实现原理是怎样的？

序列的图对象表示

其实这个序列图嵌入方法与自然语言处理分析之中的TextRank算法中对自然语言文本的图嵌入过程有一些类似，只不过TextRank处理的自然语言文本是没有太多重复单词的文本数据，而对于序列图嵌入方法而言，其所进行图嵌入的目标则是存在有大量重复单词出现的序列数据。序列图嵌入和TextRank自然语言文本处理方法，在最开始对所输入的目标文本的图对象表示，都可以基于下面的方法描述来进行：

1. 对于任意一个给定的文本数据，我们会需要首先进行文本的分词处理，对于TextRank算法所处理的自然语言文本，如果是英文文本的话，会非常容易，直接按照空格做拆分。对于中文还需要一些额外的算法做分词。对于生物序列数据，则是非常简单的吧序列文本里面的每一个字符拿出来就是一个单词对象了
2. 接着我们就可以在分词结果的基础上进行图对象的构建了：在图对象里面，都是存在有节点和节点之间的边连接关系这两种类型的数据组成的。在我们前面所分词处理后的文本数据中，所分出来的每一个单词实际上都相当于图对象之中的节点对象；那么，这个时候，源文本中的前一个单词和后一个单词二者之间就可以添加一条边连接。通过单词和单词之间的前后链接关系我们就可以将任意一个文本数据转换为基于单词节点的图对象了

在完成了上面所描述的图对象的构建之后，我们接下来就可以将这个图对象进行向量化表示就完成了将我们所获得的序列数据描述为等长的向量数据这一目标了。那么如何将我们的图对象进行向量化描述呢？

> example_sequence
ATGCCCGTCCCGTAAA

我们可以考虑使用kmer计数的方法来做简单的将我们的需要做分析的生物序列转化为向量化的描述，例如在下面表格中所展示的以k=2的kmer片段的获取然后做计数：

通过序列中，残基元素之间的边连接关系得到上面的相邻碱基计数表格后，实际上就完成了对一条任意长度的Fasta序列嵌入为了等长元素向量；在完成对原始Count值进行序列长度归一化之后，就可以与其他的任意Fasta序列的嵌入向量结果之间进行比较计算

上面的kmer计数做法非常简单直接，但是仅仅考虑了元素的出现次数，对于具有生物学意义的生物序列上的可能存在的长距离上的关联却无法体现出来。虽然直接使用上面的信息做为嵌入结果，效果可能会不怎么样，但是假若我们观察上面的kmer信息就可以发现，如果我们将前一个元素作为网络边链接的源节点，后一个元素作为网络边链接中的目标节点，那实际上我们就是针对我们的序列做了一次网络图转换。现在好歹我们有了一个比较容易理解的学习开端了。

图对象嵌入表示

因为我们创建图对象的本意就是用来描述任意两个节点之间的关联关系的，所以对于任意图对象而言，我们使用两两组合排列的方式吧图对象里面所有节点的连接关系给列举出来实际上就可以将图对象转换为向量了。当组合出来的两个节点之间存在有边连接关系的时候，就把边连接关系的权重放入到向量中对应的位置；当组合出来两个节点之间没有边连接关系的时候，就放一个零，这样子，把图对象表示为向量是不是非常简单。

对于通过自然语言文本生成的图对象而言，上面的方法可能会有些行不通，因为自然语言分词出来的单词可能是无限数量类型的，所以假若描述为向量的话，向量理论上也是无限长的。但是对于生物序列数据所产生的图对象而言，则情况简单的多了。

其实我们观察生物序列的文本特点就可以了解到，任意一条生物序列数据，里面都是由有限的基础的若干元素所构成的：DNA序列ATGC，RNA序列AUGC，蛋白序列则是20种常见氨基酸的简写字母。那么我们通过上面所描述的方法进行生物序列图对象的向量化表示，得到的就是有限长度的向量结果。

那么现在既然有了可以用来描述对应生物序列的向量数据了，那现在我们就可以基于这些向量进行常规数学方法的机器学习建模和数据挖掘分析了。

使用SGT算法进行任意序列的图嵌入

我在github上发现了一个比较有趣的代码库，可以通过上面所讲解的图嵌入的方式针对我们的生物序列做向量化嵌入：

Sequence Graph Transform (SGT) — Sequence Embedding for Clustering, Classification, and Search

Sequence Graph Transform (SGT) is a sequence embedding function. SGT extracts the short- and long-term sequence features and embeds them in a finite-dimensional feature space. The long and short term patterns embedded in SGT can be tuned without any increase in the computation.

https://github.com/cran2367/sgt/blob/25bf28097788fbbf9727abad91ec6e59873947cc/python/sgt-package/sgt/sgt.py

在github上面的代码库中，算法的描述原文为：

Compute embedding of a single or a collection of discrete item sequences. A discrete item sequence is a sequence made from a set discrete elements, also known as alphabet set. For example, suppose the alphabet set is the set of roman letters, {A, B, ..., Z}. This set is made of discrete elements. Examples of sequences from such a set are AABADDSA, UADSFJPFFFOIHOUGD, etc. Such sequence datasets are commonly found in online industry, for example, item purchase history, where the alphabet set is the set of all product items. Sequence datasets are abundant in bioinformatics as protein sequences.

Using the embeddings created here, classification and clustering models can be built for sequence datasets. Read more in https://arxiv.org/pdf/1608.03533.pdf

翻译成中文就是，SGT 把一条离散符号序列，变成一个固定维度的向量 / 矩阵，这个向量编码了“不同符号之间在不同距离下的关系强度”，并可以通过参数 kappa 控制更关注“近程关系”还是“远程关系”。在SGT算法中，并不是简单地数“某个词出现几次”，而是刻画“不同符号前后如何组合、以及组合的距离模式”。通俗的来讲，就是 SGT算法 想要做什么？：

• 对一个字母表里的每一对符号 (u, v)，比如 (A, B)、(C, D)……
• SGT 计算这一对符号在序列中“以各种距离出现”的平均强度；
• 再对这些“距离关系”做一个带指数衰减的加权，让远距离的影响可控；
• 最终得到一个“符号对关系矩阵 / 向量”，可以用来做聚类、分类、相似度搜索等。

以我们使用SGT算法来进行任意长度的生物序列的向量化嵌入工作为例，我们在这里的向量化嵌入工作中，可以给定：

• 一个有限字母表（符号集合）：Σ = {a₁, a₂, …, aₙ}
  1. 对于核酸：Σ = {A, C, G, T, U}
  2. 对于蛋白质：Σ = {20 种氨基酸单字母缩写}
• 一条序列：S = (s₁, s₂, …, s_L)，其中 s_t ∈ Σ

然后我们基于上面的计算的目标字符集，定义了针对目标生物序列的计算目标：在工作代码中定义出一个映射：S → F(S) ∈ ℝ^{n×n} 或 ℝ^{n²} 使得 F(S) 能编码 S 中的短程和长程模式，且维度只依赖于 |Σ|，不直接依赖于 L。由于我们是经常以稀疏矩阵的形式来表示SGT算法中所使用到的图对象，所以自然而然的，SGT 的输出就是是一个 n×n 矩阵，我们需要在输出之前再做 flatten 成一个 n² 维的向量给使用者。

SGT算法的核心思想：用“有向完全图”表示符号对关系

在SGT算法之中，SGT 把字母表 Σ 看作一个有向完全图的节点集合：节点为我们所限定的目标字符集，每个符号 a ∈ Σ，然后通过这个字符集中的字符在我们的目标序列上的先后顺序就可以建立起有向边，在这里所构建的有向图中任意有序对 (a_i, a_j)，都有一条边 a_i → a_j。在进行序列嵌入的时候，基于所建立的有向图，SGT 会为每一条边 (a_i, a_j) 计算一个“强度值”，这个强度值综合了：

• a_i 后面出现 a_j 的“频率”
• a_i 与 a_j 在序列中的“距离”分布
• 通过参数 κ 控制对“距离”的敏感度

最终得到的 n×n 稀疏矩阵（网络图结构对象） M[i, j] 就是边 a_i → a_j 的强度值。在vb.net代码中，我们可以使用下面的LINQ查询代码来从我们的生物序列数据中建立起上面所描述的有向图中的边连接集合信息：

From ai As Integer In U
From bj As Integer In V
Where bj > ai
Select (i:=ai, j:=bj)

额外的，我们也可以定义出其他的LINQ查询，实现出具有不同效果的邻接图用于作为进行SGT算法的图嵌入的数据来源，例如：

' just find for pattern AB in current tuple graph
From ai As Integer In U
From bj As Integer In V
Where bj = ai + 1
Select (i:=ai, j:=bj)

数学计算原理

为了讲清楚SGT算法的数学公式原理，我们先定义几组变量。

• 位置集合 - 设：P(a) = { t ∈ {1,…,L} | s_t = a }即符号 a 在序列 S 中所有出现的位置集合。
• 出现对集合 - 对两个符号 u, v ∈ Σ，定义在 S 中形成的所有“出现对”为网络图中的边(i,j), 如上面的LINQ代码所示，每一对 (i, j) 需要满足下面的属性：i ∈ P(u)（位置 i 上是 u）；j ∈ P(v)（位置 j 上是 v）；并且 j > i（v 在 u 之后）。然后，我们就可以记这个有序对集合为C(u, v) = { (i, j) | i ∈ P(u), j ∈ P(v), j > i }
• 两个核心矩阵 W₀ 和 Wₖ - 在SGT算法中，矩阵W₀[u, v]可以被简单理解为“u → v”的出现频率 / 规模，而Wₖ[u, v]则是对这些出现对做一个“带距离衰减的加权总和”。这样子SGT 最终的嵌入值是从 Wₖ 和 W₀ 组合出来的。

现在，我们有了上面定义的三组变量之后，就可以开始学习SGT的具体的实现公式了。下面给出 SGT 的具体公式（在VB.NET代码中和Python版本的代码相比较，不同实现会有微小差异，但核心结构我保证是一致的）。

Dim size = alphabets.Length
Dim l = 0
Dim W0 As NumericMatrix = NumericMatrix.Zero(size, size)
Dim Wk As NumericMatrix = NumericMatrix.Zero(size, size)
Dim positions = get_positions(sequence, alphabets)
Dim alphabets_in_sequence = sequence.Distinct.ToArray
Dim cu, cv As Vector
Dim c As (i As Integer, j As Integer)()
Dim V2 As Integer()

For Each char_i As SeqValue(Of Char) In alphabets_in_sequence.SeqIterator
    Dim i As Integer = char_i.i
    Dim u As Char = char_i.value
    Dim Upos As Integer() = positions(u)

    For Each char_j As SeqValue(Of Char) In alphabets_in_sequence.SeqIterator
        Dim j As Integer = char_j.i
        Dim v As Char = char_j.value
        Dim pos_i = _alphabets.IndexOf(u)
        Dim pos_j = _alphabets.IndexOf(v)

        If positions.ContainsKey(v) Then
            V2 = positions(v)
        Else
            Call $"KeyNotFound: The given key '{v}' was not present in the dictionary.".Warning
            V2 = {}
        End If

        c = combine(Upos, V2).ToArray
        cu = c.Select(Function(ic) ic.i).ToArray
        cv = c.Select(Function(ic) ic.j).ToArray

        W0(pos_i, pos_j) = c.Length
        Wk(pos_i, pos_j) = (-kappa * (cu - cv).Abs).Exp().Sum
    Next

    l += Upos.Length
Next

If lengthsensitive Then
    W0 /= l
End If

' avoid divide by 0
W0(W0 = 0.0) = 1.0E-18

Dim sgt = (Wk / W0) ^ (1 / kappa)
Dim sgtv As Double() = sgt.ArrayPack.IteratesALL.ToArray

Return sgtv

在上面所实现的SGT计算代码中，代码W0(pos_i, pos_j) = c.Length针对W0矩阵进行边连接出现次数的计数。在这里，最基本的计数是“有多少个 (i, j) 满足条件”：记：N{u,v} = |C(u, v)| = |{ (i, j) | i ∈ P(u), j ∈ P(v), j > i }|，则 W0 定义为（未归一化时）：W₀[u, v] = N{u,v}。在这里实现的SGT算法中，会根据 lengthsensitive 参数配置做一个针对输入的目标序列其自身长度的“归一化”操作，若 在代码中我们启用了 lengthsensitive = True（长度敏感模式），则通常会除以序列长度 L 或总出现数，使得不同长度的序列在同一尺度上：W₀[u, v] = N{u,v} / L，反之若 lengthsensitive = False（默认），则 W₀ 就是原始计数 N{u,v}。

接着，我们在得到了W0计数后，就可以计算出Wk矩阵，完成距离加权总和。在 SGT图嵌入算法 其关键之处在于“距离加权”，即我们对每个 (i, j) ∈ C(u, v)，用“距离”的指数衰减来加权。从上面的构建连接图的LINQ查询代码中，我们可以了解到，对于任意两个节点在序列上的距离，可以直接基于d = j − i（从位置 i 到位置 j）计算出来，那对每个 (i, j)我们就可以定义出一个权重 w(i, j) = exp( −κ · |i − j| ) = exp( −κ·d )，其中 κ > 0 是超参数 kappa。然后把所有 (i, j) 的权重求和，得到：Wₖ[u, v] = Σ_{(i, j) ∈ C(u, v)} exp( −κ·|i − j| )。

在这里有一个关键的参数kappa（κ > 0，常用取值：1, 5, 10 等），其作用就是为我们控制“长程依赖”被捕捉的程度。我们直观的对这个参数做理解就是：

• κ 越大，上面的计算代码中，exp(−κ·d) 在 d 稍微变大时就会变成极小值，这就意味着几乎只有相邻 / 非常近的 u, v 对才会对 Wₖ 有贡献，即我们所进行SGT图嵌入将会更偏向“局部模式”（n-gram 风格）
• 而当κ参数越小的时候，对应的节点间相互作用也就会衰减慢，远距离的 (i, j) 也有较明显的权重。这个时候嵌入会更偏向“全局结构 / 长程依赖”

最后，我们就可以把上面得到的两个通过输入的序列建立的图对象矩阵，进行组合，构建出SGT的嵌入结果M[u, v]。在SGT 的原始论文中，把 Wₖ 和 W₀ 组合成一个“平均衰减强度的 κ 次根”：M[u, v] = ( Wₖ[u, v] / W₀[u, v] )^{1/κ}。对应的计算代码就是很简单的一行矩阵运行代码：

Dim sgt = (Wk / W0) ^ (1 / kappa)

看完上面的代码后，我们在这里再继续深入解释以下这个嵌入的实现形式：

• 对于代码中的Wₖ[u, v] / W₀[u, v]矩阵除法运算部分，可以看作“所有 u→v 出现对（i,j）的 exp(−κ·d) 的平均值”，即Wₖ 是加权和，W₀ 是计数（相当于权重都为 1 的和），那总和除以总数就是得到平均值了
• 接着，我们对这个“平均衰减”再开 κ 次方：当 κ = 1 时，就是直接取平均衰减；当 κ 增大时，相当于把衰减强度“放大”，使得远程模式的影响力相对更小。这是因为在这里当k越大的时候exp(−κ·d) 在较大 d 时迅速接近 0，意味着平均衰减值变小，我们再针对其开方之后数值还会进一步被“压缩”；当k越小的时候，两个节点再序列上间隔的距离影响就会被弱化，此时，算法更像是“只关心出现对是否存在，不太在乎多远”

简单的例子来理解SGT算法

假设我们有一个有限的字母表 Σ = {A, B, C}，可以存在有序列S = ABACAB（长度L=6），那按照上面所展示的代码逻辑我们可以先算各符号位置，有：

• P(A) = {1, 3, 5}
• P(B) = {2, 6}
• P(c) = {4}

我们以 (A, B) 两个节点构建出来的网络边连接为例来进行理解：对于(A,B)，其满足 i ∈ P(A), j ∈ P(B), j > i 的 (i, j)位置的关系对包含有3个(1,2)，(3,6)，(5,6)。因此，N_{A,B} = 3 → W₀[A, B] = 3（或 /L 做 lengthsensitive）。这个时候我们实际上就得到了上面所提到的W0和Wk矩阵中的元素值，即：

• 假设 κ = 1，则：d 分别是：1, 3, 1Wₖ[A, B] = exp(−1) + exp(−3) + exp(−1) ≈ 0.3679 + 0.0498 + 0.3679 ≈ 0.7856
• 则嵌入值：M[A, B] = (Wₖ[A, B] / W₀[A, B])^{1/κ} = (0.7856 / 3)^1 ≈ 0.2619

类似地，我们通过编写代码对图中的每对 (u, v) 都算一次，就得到 3×3 的 SGT 矩阵，然后拉平成 9 维向量。这个9个维度的向量就是我们针对上面所示的序列S进行的向量化嵌入结果。

使用GCModeller进行Fasta序列的图嵌入

基于上面我们所熟知后的的SGT算法，现在我们在GCModeller之中定义下一个帮助类，用来进行任意Fasta序列的图嵌入操作。这个帮助类的代码如下：

''' <summary>
''' Create sgt embedding matrix
''' </summary>
Public Class CreateMatrix

    ReadOnly sgt As SequenceGraphTransform

    ''' <summary>
    ''' get the dimension size of the generated biological sequence 
    ''' matrix via the function <see cref="ToMatrix(FastaSeq)"/>
    ''' </summary>
    ''' <returns>The CNN input size</returns>
    Public ReadOnly Property dimension As Size
        <MethodImpl(MethodImplOptions.AggressiveInlining)>
        Get
            Return New Size(sgt.alphabets.Length, sgt.alphabets.Length)
        End Get
    End Property

    Sub New(Optional mol As SeqTypes = SeqTypes.Protein)
        Dim allChars As String() = GetChars(mol)

        sgt = New SequenceGraphTransform
        sgt.set_alphabets(allChars)
    End Sub

    Private Shared Function GetChars(mol As SeqTypes) As String()
        Select Case mol
            Case SeqTypes.Protein
                Return AminoAcidObjUtility _
                    .AminoAcidLetters _
                    .JoinIterates("-") _
                    .AsCharacter _
                    .ToArray
            Case Else
                Return mol.GetVector _
                    .JoinIterates({"-"c, "N"c}) _
                    .AsCharacter _
                    .ToArray
        End Select
    End Function

    Public Function ToMatrix(seq As FastaSeq) As Double()()
        Dim v = sgt.fit(seq.SequenceData)
        Dim m = sgt.TranslateMatrix(v)

        Return m
    End Function

    ''' <summary>
    ''' make the given fasta sequence embedding as vector
    ''' </summary>
    ''' <param name="seq"></param>
    ''' <returns></returns>
    <MethodImpl(MethodImplOptions.AggressiveInlining)>
    Public Function ToVector(seq As FastaSeq) As Double()
        Return sgt.fitVector(seq.SequenceData)
    End Function

End Class

可以看得到，这个帮助类的构造函数接受一个序列类型的参数，然后在构造函数中根据序列类别加载不同的序列字母集合到SGT算法模块之中。在创建了这个帮助类对象实例之后，我们就可以通过ToVector函数将fasta序列进行嵌入为向量了。基于这个得到的向量数据，因为其是等长的，所以就可以非常容易的应用到下游的数据建模分析之中：例如聚类，机器学习训练之类的要求结构化数据输入的工作场景。

在此帮助类的基础之上，GCModeller的序列工具包中定义了针对生物序列进行图嵌入操作的工具函数，可以在R#语言之中进行使用。

''' <summary>
''' Create algorithm for make sequence embedding
''' </summary>
''' <param name="moltype"></param>
''' <returns></returns>
<ExportAPI("seq_sgt")>
Public Function seq_sgt(Optional moltype As SeqTypes = SeqTypes.Protein) As CreateMatrix
    Return New CreateMatrix(moltype)
End Function

''' <summary>
''' embedding the given fasta sequence as vector
''' </summary>
''' <param name="sgt"></param>
''' <param name="seqs"></param>
''' <param name="env"></param>
''' <returns></returns>
<ExportAPI("seq_vector")>
Public Function seq_vector(sgt As CreateMatrix, 
                           <RRawVectorArgument> 
                           seqs As Object, 
                           Optional env As Environment = Nothing) As Object

    Dim seq_pool = GetFastaSeq(seqs, env).ToArray

    If seq_pool.Length = 1 Then
        Return sgt.ToVector(seq_pool(0))
    Else
        Dim vec As list = list.empty

        For Each seq As FastaSeq In seq_pool
            Call vec.add(seq.Title, sgt.ToVector(seq))
        Next

        Return vec
    End If
End Function

可以了解到，在GCModeller中对外导出了两个函数，一个函数seq_sgt是用于创建针对特定的序列类型的图嵌入算法对象。第二个函数就可以基于前面所创建的图嵌入算法模块，针对输入的fasta序列进行向量嵌入操作。下面是一个R#脚本的小小的针对fasta序列向量化嵌入的示例：

require(GCModeller);

imports "bioseq.fasta" from "seqtoolkit";

# get fasta sequence data
let seqs = read.fasta("./proteins.fa");
let sgt  = seq_sgt(moltype = "prot", kappa = 1,lengthsensitive= FALSE);
# get sequence embedding result
let vec  = sgt |> seq_vector(seqs, as.dataframe=TRUE);

# run data analysis on the generated embedding vectors
# embedding the raw matrix from high dimensional space
# into latent space
let latent = prcomp(vec, scal=TRUE,center=TRUE, pc=6);
# run clustering
let clusters = as.data.frame(
     kmeans(latent, centers = 6, bisecting = TRUE));

print(clusters, max.print = 6);

bitmap(file = "protein_classification.png", size = c( 1200,  800)) {
    plot(clusters$dimension_1, clusters$dimension_2,
        class = clusters$Cluster,
        colors = "paper",
        point.size = 8,
        grid.fill = "white",
        padding = "padding: 5% 10% 15% 15%;"
    );
}

在这里，我们将源代码库中的蛋白序列数据集下载下来做测试：

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

• 环境微生物群落GEMs建模综述 - 2025年12月29日
• 【虚拟细胞】翻译事件建模 - 2025年12月21日
• Motif的SequenceLogo图的绘制原理 - 2025年12月20日