文章阅读目录大纲
宏基因组测序直接从环境样本获取所有生物的遗传物质，产生的海量短读序列（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位置的计算。
那，现在基于上面的哈希计算函数，我们就可以实现下面所示的两个布隆过滤器的工作逻辑：
- 向布隆过滤器中添加背景元素：将输入的kmer字符串数据，通过上面的哈希函数计算出对应的bitmap位置，然后将对应的bitmap位置的值设置为True，即可将目标背景中的特定kmer在布隆过滤器中标记上
- 检测目标元素是否存在于布隆过滤器中：现在，我们将某一个需要做检查的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。这两个关键条件是：
- 预期要存储的元素数量 (n)：你估计最多会有多少个元素被加入到布隆过滤器中。
- 期望的假阳性率：你愿意接受的“误判”概率。即，一个不在集合中的元素，被布隆过滤器判断为“可能存在”的概率。
因此，bitarray 长度 m 的计算，本质上是一个在给定 n 和 p 的前提下，求解最优 m 的数学问题。这个数学问题的原理可以用下面的几个步骤来描述：
第一步：计算一个比特位在插入 n 个元素后仍然为 0 的概率
- 初始状态：一个空的布隆过滤器，其
m个比特位全部为 0。 - 插入一个元素：当我们插入一个元素时，会使用
k个哈希函数计算出k个位置，并将这k个位置的比特位设置为 1。 - 单个比特位的命运：对于
bitarray中的任意一个特定比特位，它被这一个元素的哈希函数击中（设置为1）的概率是 k/m。那么，它没有被这一个元素击中的概率就是 1 - k/m。 - 插入
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 时。
- 我们对
x使用k个哈希函数，得到k个位置。 - 如果这
k个位置上的比特位恰好都为 1，布隆过滤器就会错误地判断x"可能存在"，这就是一次假阳性。 - 由于哈希函数的独立性，这
k个比特位都为 1 的概率就是它们各自为 1 的概率的乘积。 -
因此，假阳性率
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/eln(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)
Dim kmers As String() = KSeq.KmerSpans(read.GetSequenceData, k).ToArray
' scan every genome its bloom filter model
' and get kmer hits counter result
For Each genome As KmerBloomFilter In genomes
Dim numHits As Integer = genome.KmerHitNumber(kmers)
Dim unmap As Integer = kmers.Length - numHits
If hits.ContainsKey(genome.ncbi_taxid) Then
hits(genome.ncbi_taxid) += numHits
Else
hits(genome.ncbi_taxid) = numHits
End If
If hits.ContainsKey(0L) Then
hits(0L) += unmap
Else
hits(0L) = unmap
End If
Next
' filter of the genome hits via reads sequence coverage filter
Dim numCov As Integer = kmers.Length * coverage
Dim tax_id = From tax As KeyValuePair(Of Integer, Integer)
In hits
Where tax.Key > 0 AndAlso tax.Value >= numCov
Select tax.Key
' measure taxonomy classification via LCA method
Dim lcaNode As LcaResult = NcbiTaxonomyTree.GetLCAForMetagenomics(tax_id, min_supports)
If lcaNode.LCATaxid = 0 Then
Dim topHit = hits.Where(Function(a) a.Key > 0).OrderByDescending(Function(a) a.Value).FirstOrDefault
If topHit.Key > 0 AndAlso topHit.Value > 0 Then
If topHit.Value / kmers.Length > 0.975 Then
lcaNode = New LcaResult With {
.lcaNode = NcbiTaxonomyTree.GetNode(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(hits.Keys.Any(Function(t) t > 0), "C", "U"),
.TaxID = lcaNode.LCATaxid
}
End Function在上面所展示的测序reads分类函数中，展示了基于布隆过滤器做分类的三个主要的过程：
- 我们将测序reads分解为kmer序列片段集合，然后使用这个集合快速的扫描每一个基因组的布隆过滤器模型，得到reads在每一个基因组模型上的匹配结果情况
- 基于匹配的数量，按照一个reads片段覆盖度的阈值做过滤，过滤出能够基本上覆盖完全整条reads的结果
- 经过覆盖度的过滤操作后，可能会留下几个基因组的匹配结果，则这个时候我们就可以通过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");
- 布隆过滤器在宏基因组测序reads数据分类处理中的作用 - 2025年12月18日
- limma程序包在RNA-seq差异表达分析中的数学算法原理与实现详解 - 2025年12月16日
- 零分布原理及其在生物信息学中的应用 - 2025年12月15日
No responses yet