布隆过滤器在宏基因组测序reads数据分类处理中的作用

宏基因组测序直接从环境样本获取所有生物的遗传物质，产生的海量短读序列（reads）需要被快速准确地分类到不同物种或功能类别。然而，宏基因组数据具有复杂性高、物种多样且未知序列多等特点，这给分类算法带来了巨大挑战。传统的序列比对方法虽然准确，但在面对庞大的参考数据库时计算开销巨大，难以满足实时分析的需求。因此，研究者开发了多种基于k-mer（长度为k的子序列）的快速分类方法，其中布隆过滤器（Bloom Filter）作为一种高效的概率数据结构，在针对测序reads做物种上的快速分类这项工作中起到了一些关键作用。

布隆过滤器最初由Burton Bloom于1970年提出（Space/time trade-offs in hash coding with allowable errors，Communications of the ACM, Volume 13, Issue 7；DOI: 10.1145/362686.362692），用于以极小的空间开销快速判断元素是否属于一个集合。它通过牺牲一定的准确性（允许一定概率的假阳性）来换取内存占用的大幅降低和查询速度的显著提升。在宏基因组reads分类中，布隆过滤器常被用于高效存储和查询海量的k-mer集合，从而加速序列相似性判断和物种注释过程。

布隆过滤器原理

布隆过滤器由一个位数组（bit array）和一组哈希函数构成。初始化时，所有位都置为0。插入一个元素时，将该元素经过多个哈希函数映射到位数组的若干位置，并将这些位置置为1。查询时，同样用哈希函数计算待查元素的位置，若所有对应位均为1，则“可能”存在该元素（因为可能被其他元素映射到相同位置而置1）；若存在任一位置为0，则“一定”不存在该元素。这种机制保证了布隆过滤器不会出现假阴性（即不存在的东西不会被误判为存在），但可能产生假阳性（存在的东西可能被误判为不存在）。

布隆过滤器的突出优势在于其极高的空间效率。相比直接存储所有元素（如使用哈希表或数组），布隆过滤器仅用少量位就能表示海量元素集合。例如，在宏基因组应用中，一个包含数百万k-mers的集合，用布隆过滤器存储可能只需几MB内存，而传统方法可能需要数百MB。同时，查询操作的时间复杂度为O(k)（k为哈希函数个数），与元素数量无关，因此查询速度极快。这使得布隆过滤器非常适合需要频繁判断“是否存在”的场景，例如检查一个k-mer是否属于某个参考基因组集合。

宏基因组分类常涉及判断reads中的k-mers是否与已知物种的k-mers匹配。直接存储所有参考k-mers会占用巨大内存，而布隆过滤器可以紧凑地表示这些k-mers集合。通过将参考序列的所有k-mers插入布隆过滤器，可以快速过滤掉大量不相关的reads，仅保留可能与参考匹配的候选序列，从而大幅加速后续分析。此外，布隆过滤器对插入顺序不敏感，支持动态添加元素，这在构建可扩展的分类索引时非常有用。

从头实现一个布隆过滤器

布隆过滤器算法的精髓就是通过哈希算法将我们的数据输入映射为序列中的某一个位置，然后通过检查所映射到的对应的位置上是否存在对应的数据标记即可了解我们的目标对象是否存在构建布隆过滤器的背景数据之中了。在VisualBasic.NET中，可以直接使用BitArray来表示这样子的一个可以记录对应的位置上是否存在数据标记的数据结构：

''' <summary>
''' 使用 BitArray 来高效地存储位数组
''' </summary>
ReadOnly bitmap As BitArray

那现在我们就只需要添加一个哈希函数，将我们的数据输入转换为对应的位置索引号即可实现一个布隆过滤器了。在这里，我们针对kmer字符串进行哈希计算的对应的哈希函数如下所示：

''' <summary>
''' 获取元素在位数组中映射的 k 个位置。
''' 这里使用双哈希技术来生成 k 个独立的哈希值。
''' </summary>
''' <param name="item">要哈希的元素。</param>
''' <returns>包含 k 个位置的整数数组。</returns>
Private Function GetHashPositions(item As String) As Integer()
    ' 1. 将字符串转换为字节数组
    Dim data As Byte() = Encoding.UTF8.GetBytes(item)

    ' 2. 使用两个不同的种子进行双哈希
    Dim h1 As UInteger = MurmurHash.MurmurHashCode3_x86_32(data, 0) ' 种子 0
    Dim h2 As UInteger = MurmurHash.MurmurHashCode3_x86_32(data, &HFFFFFFFFUI) ' 种子 FFFFFFFF
    Dim positions(_k - 1) As Integer

    ' 3. 组合哈希以生成 k 个位置
    For i As Integer = 0 To _k - 1
        ' 组合哈希：Hash_i = (h1 + i * h2) mod m
        ' 确保 h2 不为 0，以避免所有位置都相同
        If h2 = 0 Then
            h2 = 1
        End If

        positions(i) = (h1 + i * h2) Mod _m
    Next

    Return positions
End Function

在上面的计算函数中，我们通过相应的哈希算法，将输入的kmer字符串数组计算为对应的哈希值。然后基于计算出来的哈希值针对bitmap的长度进行取余，就完成了对应的bitmap位置的计算。

那，现在基于上面的哈希计算函数，我们就可以实现下面所示的两个布隆过滤器的工作逻辑：

1. 向布隆过滤器中添加背景元素：将输入的kmer字符串数据，通过上面的哈希函数计算出对应的bitmap位置，然后将对应的bitmap位置的值设置为True，即可将目标背景中的特定kmer在布隆过滤器中标记上
2. 检测目标元素是否存在于布隆过滤器中：现在，我们将某一个需要做检查的kmer字符串数据，仍然是通过上面的哈希函数计算出其在布隆过滤器中的对应的bitmap位置，然后检查所映射的bitmap位置是否是True，即可完成检测目标元素是否存在于布隆过滤器中的检查工作了

基于上面的算法描述，下面我们就可以实现对应的两个方法函数：

''' <summary>
''' 将指定的元素添加到布隆过滤器中。
''' </summary>
''' <param name="item">要添加的元素。不能为 Nothing。</param>
Public Sub Add(item As String)
    Dim positions = GetHashPositions(item)

    For Each pos As Integer In positions
        bitmap(pos) = True
    Next
End Sub

''' <summary>
''' 检查布隆过滤器中是否可能包含指定的元素。
''' </summary>
''' <param name="item">要检查的元素。不能为 Nothing。</param>
''' <returns>如果元素可能存在，则为 true；如果元素绝对不存在，则为 false。</returns>
Public Function Contains(item As String) As Boolean
    Dim positions = GetHashPositions(item)

    For Each pos As Integer In positions
        ' 只要有一个哈希位置不为1，就说明元素绝对不存在
        If Not bitmap(pos) Then
            Return False
        End If
    Next

    ' 如果所有位置都为1，则元素可能存在（存在假阳性可能）
    Return True
End Function

布隆过滤器位映射长度的计算统计学原理

到现在这里，你可能会有一个疑问：我的bitmap的长度会是多少？首先，我们需要明确布隆过滤器初始化时的目标。我们不是随便设置一个 bitarray 长度，而是希望在满足两个关键条件的前提下，找到最优的（即最小的）bitarray 长度 m。这两个关键条件是：

1. 预期要存储的元素数量 (n)：你估计最多会有多少个元素被加入到布隆过滤器中。
2. 期望的假阳性率：你愿意接受的“误判”概率。即，一个不在集合中的元素，被布隆过滤器判断为“可能存在”的概率。

因此，bitarray 长度 m 的计算，本质上是一个在给定 n 和 p 的前提下，求解最优 m 的数学问题。这个数学问题的原理可以用下面的几个步骤来描述：

第一步：计算一个比特位在插入 n 个元素后仍然为 0 的概率

1. 初始状态：一个空的布隆过滤器，其 m 个比特位全部为 0。
2. 插入一个元素：当我们插入一个元素时，会使用 k 个哈希函数计算出 k 个位置，并将这 k 个位置的比特位设置为 1。
3. 单个比特位的命运：对于 bitarray 中的任意一个特定比特位，它被这一个元素的哈希函数击中（设置为1）的概率是 k/m。那么，它没有被这一个元素击中的概率就是 1 - k/m。
4. 插入 n 个元素：现在我们插入了 n 个元素。由于哈希函数的独立性，这个特定比特位在 n 次插入过程中都未被击中的概率是：

   # P_bit_is_0_after_n_insertions
   (1 - 1/m)^(k*n)

   注意：这里用 1/m 而不是 k/m 是更精确的推导。一个哈希函数将其设置为1的概率是 1/m，k 个哈希函数都没设置它的概率是 (1 - 1/m)^k。n 个元素都这样操作，所以是 ((1 - 1/m)^k)^n = (1 - 1/m)^(kn)。

第二步：利用极限公式简化

当 m 很大时，我们可以使用一个重要的极限公式来简化上面的表达式：

# limit_as_x_to_infinity 
# 约等于 1/e
(1 - 1/x)^x  

其中 e 是自然常数，约等于 2.718。我们将这个近似应用到我们的概率计算中：

# approximation
# 约等于 (1/e)^(k*n/m) = e^(-k*n/m)
((1 - 1/m)^m)^(k*n/m)

所以，一个比特位在插入 n 个元素后仍然为 0 的概率约等于：

# P_bit_is_0
exp(-k*n/m)

那么，这个比特位为 1 的概率就是：

# P_bit_is_1
1 - exp(-k*n/m)

第三步：推导假阳性率 p

假阳性发生在什么情况下？当我们查询一个不在集合中的元素 x 时。

1. 我们对 x 使用 k 个哈希函数，得到 k 个位置。
2. 如果这 k 个位置上的比特位恰好都为 1，布隆过滤器就会错误地判断 x "可能存在"，这就是一次假阳性。
3. 由于哈希函数的独立性，这 k 个比特位都为 1 的概率就是它们各自为 1 的概率的乘积。

4. 因此，假阳性率 p 的近似计算公式为：

   p = (1 - exp(-k*n/m))^k

   这个公式是布隆过滤器理论的核心，它揭示了 m, k, n, p 四个变量之间的近似关系。

第四步：求解最优的 k 和 m

我们的目标是在给定 n 和 p 的情况下，找到最小的 m。要做到这一点，我们首先需要找到在 m 和 n 固定的情况下，能使假阳性率 p 最小的最优哈希函数数量 k。通过对 p 关于 k 求导并令导数为 0，可以解出最优的 k：

k_opt = (m/n) * log(2)

这个结论非常直观：位数组越大（m 越大），或者元素越少（n 越小），就可以使用越多的哈希函数来降低误判率。现在，我们将这个最优的 k 代回到假阳性率 p 的公式中，来求解 m：

p_approximation = (1 - exp(-((m/n) * log(2)) * n/m))^((m/n) * log(2))

简化指数部分：

p_simplified = (1 - exp(-log(2)))^((m/n) * log(2))

因为 e^-ln(2) = 1/e^ln(2) = 1/2，所以：

# 等于 2^(-(m/n) * log(2))
p_final = (1/2)^((m/n) * log(2))

对两边取自然对数 ln：

ln_p = -(m/n) * (log(2))^2

最后，整理得到 m 的计算公式：

m = -n * log(p) / (log(2))^2

现在，我们不要太注重研究这个公式推导的过程，直接拿出对应的计算公式的代码答案，答案如下：

''' <summary>
''' 根据给定的元素数量和假阳性率，计算最优的位数组大小。
''' </summary>
''' <param name="n">预计插入的元素数量。</param>
''' <param name="p">期望的假阳性率。</param>
''' <returns>最优的位数组大小。</returns>
''' 
<MethodImpl(MethodImplOptions.AggressiveInlining)>
Public Shared Function OptimalM(n As Integer, p As Double) As Integer
    ' 公式: m = - (n * ln(p)) / (ln(2))^2
    Dim m As Double = -n * std.Log(p) / (std.Log(2) ^ 2)
    Dim mint As Integer = CInt(m)
    Return mint
End Function

''' <summary>
''' 根据给定的位数组大小和元素数量，计算最优的哈希函数数量。
''' </summary>
''' <param name="m">位数组的大小。</param>
''' <param name="n">预计插入的元素数量。</param>
''' <returns>最优的哈希函数数量。</returns>
''' 
<MethodImpl(MethodImplOptions.AggressiveInlining)>
Public Shared Function OptimalK(m As Integer, n As Integer) As Integer
    ' 公式: k = (m/n) * ln(2)
    Dim k As Double = (m / n) * std.Log(2)
    Dim kint As Integer = CInt(k)
    Return kint
End Function

通过上面的两个函数，我们计算出来了两个参数值：m为我们的bitmap的位数的长度，k则是哈希位置的数量，即对于某一个需要做检测的kmer元素字符串数据，其需要通过检查k个bitmap中的位置都要为True才会判断目标kmer存在于我们的背景基因组中。

''' <summary>
''' 根据期望的元素数量和目标假阳性率，自动计算最优参数并初始化 BloomFilter 类的新实例。
''' </summary>
''' <param name="expectedItems">预计将要插入的元素数量。必须为正数。</param>
''' <param name="desiredFalsePositiveRate">期望的假阳性率，范围在 (0, 1) 之间。</param>
Public Shared Function Create(expectedItems As Integer, desiredFalsePositiveRate As Double) As BloomFilter
    If expectedItems <= 0 Then
        Throw New ArgumentOutOfRangeException(NameOf(expectedItems), "期望元素数量必须为正数。")
    End If
    If desiredFalsePositiveRate <= 0.0 OrElse desiredFalsePositiveRate >= 1.0 Then
        Throw New ArgumentOutOfRangeException(NameOf(desiredFalsePositiveRate), "假阳性率必须在 0 和 1 之间。")
    End If

    ' 根据公式计算最优的 m 和 k
    Dim optimalM = BloomFilter.OptimalM(expectedItems, desiredFalsePositiveRate)
    Dim optimalK = BloomFilter.OptimalK(optimalM, expectedItems)

    ' 调用主构造函数
    Return New BloomFilter(optimalM, optimalK)
End Function

哈希算法实现

通过阅读上面的GetHashPositions这个函数，你会发现，我们目前还缺少一个东西：就是需要一个哈希算法将目标输入输入转化为对应的哈希值用来做bitmap的位置计算。在这里我们使用了MurmurHash3_x86_32 哈希算法来实现这个哈希计算过程。

MurmurHash3_x86_32 是由 Austin Appleby 于 2008 年设计的一种非加密型哈希函数，属于 MurmurHash 家族的最新版本。该算法专为 32 位环境优化，使用 32 位算术运算，生成 32 位哈希值，适用于高性能、低碰撞率的通用哈希场景。其名称“Murmur”来源于“multiply”和“rotate”两个核心操作，结合异或（xor）运算，实现了快速且高质量的哈希值生成。

MurmurHash3_x86_32 是一种高效、分布优良、易于集成的通用哈希算法，是大数据、分布式系统和搜索引擎等领域的经典选择。在这里我们并不具体的来讨论这个哈希算法的计算原理，直接贴上对应的实现代码：

Public Module MurmurHash

    ' MurmurHash3 的常量
    Const c1 As UInteger = &HCC9E2D51UI
    Const c2 As UInteger = &H1B873593UI
    Const r1 As Integer = 15
    Const r2 As Integer = 13
    Const m As UInteger = 5UI
    Const n As UInteger = &HE6546B64UI

    ''' <summary>
    ''' 计算给定数据的32位MurmurHash3值。
    ''' </summary>
    ''' <param name="data">输入数据。</param>
    ''' <param name="seed">哈希种子。</param>
    ''' <returns>32位无符号哈希值。</returns>
    Public Function MurmurHashCode3_x86_32(data As Byte(), seed As UInteger) As UInteger
        If data Is Nothing Then Throw New ArgumentNullException(NameOf(data))

        Dim length As Integer = data.Length
        Dim h As UInteger = seed
        Dim i As Integer = 0
        Dim k As UInteger = 0

        ' --- 处理主体部分（4字节块） ---
        While length - i >= 4
            k = CUInt(data(i)) Or
                            CUInt(data(i + 1)) << 8 Or
                            CUInt(data(i + 2)) << 16 Or
                            CUInt(data(i + 3)) << 24

            k *= c1
            k = (k << r1) Or (k >> (32 - r1))
            k *= c2

            h = h Xor k
            h = (h << r2) Or (h >> (32 - r2))
            h = h * m + n
            i += 4
        End While

        ' --- 处理尾部剩余字节（修正部分） ---
        k = 0
        Select Case length - i
            Case 3
                k = k Xor CUInt(data(i + 2)) << 16
                ' 注意：这里没有 Exit Select，会继续执行 Case 2
                k = k Xor CUInt(data(i + 1)) << 8
                ' 继续执行 Case 1
                k = k Xor CUInt(data(i))
            Case 2
                k = k Xor CUInt(data(i + 1)) << 8
                ' 继续执行 Case 1
                k = k Xor CUInt(data(i))
            Case 1
                k = k Xor CUInt(data(i))
        End Select

        If length - i > 0 Then
            k *= c1
            k = (k << r1) Or (k >> (32 - r1))
            k *= c2
            h = h Xor k
        End If

        ' --- 最终混合 ---
        h = h Xor CUInt(length)
        h = h Xor (h >> 16)
        h *= &H85EBCA6BUI
        h = h Xor (h >> 13)
        h *= &HC2B2AE35UI
        h = h Xor (h >> 16)

        Return h
    End Function
End Module

至此，我们已经实现了一个布隆过滤器，下面我们就需要将上面所实现的这个布隆过滤器应用到宏基因组测序reads分类的工作之中。

布隆过滤器在宏基因组测序reads分类中的应用

将布隆过滤器应用到宏基因组的测序reads的分类工作之中，会遵循着两个步骤原则：1. 首先我们会需要有用于分类判别的背景模型数据，即将背景基因组序列转换为一个布隆过滤器；2. 在得到了背景基因组的布隆过滤器模型后，对目标待分类的测序reads做kmer查询即可大致了解目标reads在我们的背景基因组中的匹配情况。

为基因组序列创建一个布隆过滤器

在程序中完成了目标基因组的fasta序列加载工作后，我们按照设定的kmer长度对输入的基因组fasta序列生成kmer序列片段，这个生成的kmer序列片段的长度一般为fasta序列长度减去k个。然后基于这个数据量和我们所期望的最高的假阳性率，就可以为布隆过滤器提供了模型创建所需要的两个参数：

Dim pool As Fasta() = genomics.ToArray
Dim estimatedKmers As Integer = Math.Max(0, Math.Min(spanSize, pool.Max(Function(s) s.length)) - k + 1)
Dim filter As BloomFilter = BloomFilter.Create(estimatedKmers, desiredFPR)

接着，基于这个已经确定的kmer的长度参数k值，从目标基因组序列中提取出kmer序列集合，然后将每一个所提取到的kmer序列添加到布隆过滤器中即可：

For Each nt As Fasta In TqdmWrapper.Wrap(pool, bar:=bar)
    Dim ntseq As String = nt.GetSequenceData

    Call names.Add(nt.title)
    Call bar.SetLabel(nt.title)

    For i As Integer = 0 To ntseq.Length Step spanSize
        Dim len As Integer = spanSize

        If i + len > ntseq.Length Then
            len = ntseq.Length - i
        End If

        For Each kmer As String In KSeq.KmerSpans(ntseq.Substring(i, len), k)
            Call filter.Add(kmer)
        Next
    Next
Next

基于布隆过滤器的测序reads分类

那，现在通过上面的代码操作，我们在得到了一个背景基因组序列的布隆过滤器之后，就可以非常简单的做测序得到的reads上提取出来的kmer序列片段的查询了：

Public Function KmerHits(kmers As IEnumerable(Of String)) As Dictionary(Of String, Integer)
    Dim hits As New Dictionary(Of String, Integer)

    For Each kmer As String In kmers
        If bloomFilter(kmer) Then
            If Not hits.ContainsKey(kmer) Then
                hits.Add(kmer, 1)
            Else
                hits(kmer) += 1
            End If
        End If
    Next

    Return hits
End Function

通过上面的工作函数我们可以看得到，我们针对测序得到的reads中的每一个kmer片段都通过目标背景基因组序列的布隆过滤器做了一次检查，基于布隆过滤器可以让我们快速的了解到对应的kmer序列片段是否存在于目标基因组序列上。然后基于成功匹配到的结果对对应的kmer序列片段做计数，这样子就可以产生了我们的测序reads上的kmer的分布情况了。理论上基于这些kmer分布的信息就可以帮助我们借助于一些统计学方法来确认reads的分类归属。

布隆分类器的实现

上面的例子是针对某一条测序的reads针对某一个特定的背景基因组的匹配的操作。那现在假设我们针对很多个背景基因组都构建了对应了布隆过滤器模型，然后我们将测序的reads都在这些背景基因组序列所对应的布隆过滤器做一次查询，然后针对统计计数结果按照从大到小的顺序做了一次降序排序，那么得到的最高匹配数量的基因组就很有可能是当前的这条测序reads的来源基因组了。

下面的函数代码实现了上面所描述的一个测序reads的基于布隆过滤器的分类计算过程：

Public Function MakeClassify(read As IFastaProvider) As KrakenOutputRecord
    Dim hits As New Dictionary(Of Integer, Integer)
    ' split reads sequence as kmers
    Dim kmers As String() = KSeq.KmerSpans(read.GetSequenceData, k).ToArray

    ' scan bloom filter model for each background genome
    ' get kmer hits number
    For Each genome As KmerBloomFilter In genomes
        Dim numHits As Integer = genome.KmerHitNumber(kmers)

        If numHits = 0 Then
            ' current genome no hits
            ' do nothing
        Else
            If hits.ContainsKey(genome.ncbi_taxid) Then
                hits(genome.ncbi_taxid) += numHits
            Else
                hits(genome.ncbi_taxid) = numHits
            End If
        End If
    Next

    ' filter low hits genome, use high hits genome for LCA evaluation
    Dim numCov As Integer = kmers.Length * coverage
    Dim candidateTaxIds = From tax As KeyValuePair(Of Integer, Integer)
                            In hits
                            Where tax.Key > 0 AndAlso tax.Value >= numCov
                            Select tax.Key
    Dim lcaNode As LcaResult = NcbiTaxonomyTree.GetLCAForMetagenomics(candidateTaxIds, min_supports)

    ' lca node not found
    If lcaNode.LCATaxid = 0 Then
        ' filter out key(0) and then sort in desc 
        Dim descSortedHits = hits.Where(Function(a) a.Key > 0).OrderByDescending(Function(a) a.Value).ToArray
        Dim topHit = descSortedHits.FirstOrDefault

        If topHit.Key > 0 AndAlso topHit.Value > 0 Then
            ' 1 top hits must be nearly identical to the target reads
            ' 2 top hits is unique hits
            ' or top hits is significant greater than other genome hits

            ' 检查是否存在一个明确、占主导地位的命中
            ' 条件1: Read中很大一部分k-mers (例如 >99%) 都命中了这一个基因组
            Dim isDominantHit As Boolean = (topHit.Value / kmers.Length) > 0.99
            ' 条件2: 这个最佳命中是唯一的，或者显著优于第二名
            Dim isUniqueOrSignificant As Boolean = False

            If descSortedHits.Length = 1 Then
                ' 只有一个分类单元被命中
                isUniqueOrSignificant = True
            ElseIf descsortedHits.Length > 1 Then
                ' 检查最佳命中的数量是否至少是第二名的两倍
                isUniqueOrSignificant = (topHit.Value / descSortedHits(1).Value) > 2
            End If

            ' 如果同时满足两个条件，则将该最佳命中作为分类结果
            If isDominantHit AndAlso isUniqueOrSignificant Then
                lcaNode = New LcaResult With {
                    .lcaNode = NcbiTaxonomyTree.GetNode(topHit.Key),
                    .supportRatio = 1,
                    .supportedTaxids = {topHit.Key}
                }
            End If
        End If
    End If

    Return New KrakenOutputRecord With {
        .LcaMappings = hits _
            .ToDictionary(Function(a) a.Key.ToString,
                            Function(a)
                                Return a.Value
                            End Function),
        .ReadLength = read.length,
        .ReadName = read.title,
        .StatusCode = If(lcaNode.LCATaxid > 0, "C", "U"), ' C - classfied, U - unmapped, no hits
        .TaxID = lcaNode.LCATaxid,
        .Taxonomy = If(lcaNode.LCATaxid = 0, "n/a", lcaNode.lcaNode.name & $"({lcaNode.lcaNode.rank})"),
        .LCASupport = lcaNode.supportRatio,
        .LCATaxids = lcaNode.supportedTaxids
    }
End Function

在上面所展示的测序reads分类函数中，展示了基于布隆过滤器做分类的三个主要的过程：

1. 我们将测序reads分解为kmer序列片段集合，然后使用这个集合快速的扫描每一个基因组的布隆过滤器模型，得到reads在每一个基因组模型上的匹配结果情况
2. 基于匹配的数量，按照一个reads片段覆盖度的阈值做过滤，过滤出能够基本上覆盖完全整条reads的结果
3. 经过覆盖度的过滤操作后，可能会留下几个基因组的匹配结果，则这个时候我们就可以通过LCA共同最近祖先算法得到一个物种分类结果，这个物种分类结果就是目标reads的物种分类结果。

测试布隆分类器

下面的R#脚本代码展示了如何通过GCModeller所提供的宏基因组处理模块进行bloom分类器模型的加载，然后针对读取的测序reads样本数据进行分类分析，最后导出为结果表格的整个过程。

let ncbi_tax = Ncbi.taxonomy_tree("/opt/ncbi_taxonomy/");
let blooms = system.file("/opt/environment_blooms")
   |> bloom_filters(
       ncbi_tax,
       min_supports = min_supports);
let sample_reads = FastQ::read.fastq("/data.fq");
let result = blooms |> make_classify(reads = sample_reads);
let find_reads = [result]::StatusCode == "C";

result = as.data.frame(result);

message("inspect of the bllom classify result:");
print(result, max.print = 6);

write.csv(result, file = "./bloom_result.csv");

从测试的结果来看，算法代码已经可以成功的将大部分来源于Influenza A virus基因组序列的模拟测序reads都分类到了Alphainfluenzavirus influenzae(species)这个物种下：

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

• 建立KEGG的KO序列数据库 - 2026年1月4日
• 环境微生物群落GEMs建模综述 - 2025年12月29日
• 【虚拟细胞】翻译事件建模 - 2025年12月21日