稀疏矩阵表示：处理稀疏kmer计数矩阵
做宏基因组分析的人，几乎都跟kmer计数矩阵打过交道。这东西有多大？我见过一个土壤宏基因组样本，建出来的kmer计数矩阵有800G。说实在的，一开始我也想过「硬吃」，结果内存直接爆了。
后来我想明白了——kmer矩阵最大的特点是「稀疏」。一个40G的样本里，出现过的kmer种类可能只有总数的5%~10%。剩下90%以上的格子全是0。你花大价钱买的内存，存的全是0，这不科学。
那怎么办？稀疏矩阵表示就是用来解决这个问题的。今天我就带你把主流方案过一遍。
为什么kmer矩阵极度稀疏
先看个实际例子。假设你有一个宏基因组样本，测序后得到10亿条reads，每条150bp。你建k=31的kmer集。
理论上，4^31 ≈ 4.6×10¹⁸种可能的kmer。但你的样本里，实际出现的kmer数量通常只有10⁷~10⁹的数量级。覆盖率是多少？
10⁹ / 4.6×10¹⁸ ≈ 2×10⁻¹⁰
说白了，99.99999998%的kmer根本没有出现。这就是稀疏的根源。我以前跟刚入行的同事解释这事，他们说「不可能这么夸张」。但现实就是如此，你的样本再复杂，也拼不过4^31这个天文数字。
技术路线对比
项目里我用到过三种主流的稀疏表示方式。这里给你列个表，你可以根据自己的场景直接选。
| 表示方式 | 存储结构 | 查询效率 | 适合场景 |
|---|---|---|---|
| 字典（Dict of kmer→count） | 哈希表 | O(1) | 小样本，快速原型 |
| 压缩稀疏行（CSR） | 数组+偏移 | O(log n) | 大规模矩阵运算 |
| 布隆过滤+计数器 | 位向量+哈希 | O(k) | 存在性判断，内存极度受限 |
我个人习惯：如果样本量在100G以下，用Dict。超过100G而且要做矩阵运算，必须上CSR。布隆过滤我一般只用在预处理阶段，比如先过滤掉罕见kmer。
先看最简单的：基于字典的稀疏表示
Python里直接建一个字典，key是kmer字符串，value是计数。这是最直觉的做法。
from collections import defaultdict
def build_kmer_dict(reads, k=31):
kmer_count = defaultdict(int)
for read in reads:
for i in range(len(read) - k + 1):
kmer = read[i:i+k]
kmer_count[kmer] += 1
return kmer_count这段代码看着简单，但坑不少。我遇到过最典型的：kmer字符串本身占内存太狠了。一个31bp的kmer，Python字符串开销接近50字节，加上哈希表本身的负载因子，实际内存可能是理论值的5倍。
那怎么办？做整数映射。
def kmer_to_int(kmer):
mapping = {'A': 0b00, 'C': 0b01, 'G': 0b10, 'T': 0b11}
val = 0
for base in kmer:
val = (val << 2) | mapping[base]
return val
def build_kmer_dict_int(reads, k=31):
kmer_count = defaultdict(int)
for read in reads:
for i in range(len(read) - k + 1):
kmer = read[i:i+k]
key = kmer_to_int(kmer)
kmer_count[key] += 1
return kmer_count把kmer编码成整数，每个31mer只需要62位，Python里用一个int就够了。内存直接降低到原来的1/5左右。这是我优化内存的第一步，百试百灵。
再来看看CSR格式，它才是做大规模运算的利器
CSR全称Compressed Sparse Row。不扯太多理论，直接上核心思路：
- data数组：存放所有非零元素的值
- indices数组：存放每个非零元素所在的列索引
- indptr数组：存放每行在data中的起始和结束位置
比如一个3行×5列的矩阵，只有4个非零元素。用CSR存，内存占用大概是原来的1/100。
Python里用scipy.sparse的csr_matrix可以直接搞定。
from scipy.sparse import csr_matrix
import numpy as np
# 假设我们有一个密集矩阵
dense_matrix = np.array([
[0, 0, 5, 0, 0],
[0, 2, 0, 0, 0],
[0, 0, 0, 0, 9]
])
# 转成CSR
sparse_matrix = csr_matrix(dense_matrix)
print("data:", sparse_matrix.data)
print("indices:", sparse_matrix.indices)
print("indptr:", sparse_matrix.indptr)
print("shape:", sparse_matrix.shape)data: [5 2 9]
indices: [2 1 4]
indptr: [0 1 2 3]
shape: (3, 5)你看，原来需要15个格子，现在只存3个非零值加一些索引信息。空间省了多少？自己算算。
在宏基因组里，我们通常把样本当行，kmer当列。比如100个样本，每个样本有1000万个独特kmer的计数。用CSR存，数据量可能从1000GB直接掉到10GB以内。
我的避坑指南
这些东西我第一次用的时候，踩过的坑比你们想象的要多。这里列几个典型的：
我踩过的坑之一：CSR做行切片时，indptr会被复制
我曾经在一个200G的CSR矩阵上做行切片操作，以为只是「取子集」，结果内存暴涨了3倍。后来查文档才发现，CSR的切片会生成一个新的indptr数组。如果行数特别多，这个数组本身可能占好几G。
解决方案：用 sparse_matrix.tocsr()[rows, :] 之前，先确认有多少行。如果只需要少数行，不如直接读原始数据重新构建。
我踩过的坑之二：dict的碰撞导致性能崩溃
有一回跑kmer计数，字典里塞了5亿个key。Python的哈希表开始疯狂rehash，CPU占用冲到99%，但进度条一动不动。
解决方案：提前估算kmer种类数量，初始化时设置 initial_size。比如预估5亿种，上来就 dict(initial_size=5_000_000_000)，能省掉一半的运行时rehash开销。
我踩过的坑之三：布隆过滤的假阳性在计数场景下会放大错误
布隆过滤只能告诉你「可能存在」或「肯定不存在」。如果你用它来做计数缩减，假阳性会把原本不存在的kmer「误判」为存在，导致后续的计数矩阵里出现「幽灵kmer」。
解决方案：只在预处理筛掉「肯定不存在」的低频kmer，不要用布隆过滤做精确计数。
三种方案的决策流程图
我整理了一张图，把选择逻辑画清楚了。你拿到项目，照着这个图走，基本不会错。
实际项目中的应用案例
给你说个真实的项目。当时我们在分析一个海洋微生物宏基因组，100个样本，总数据量约1.2T。kmer矩阵如果用密集存储，算下来至少要600G内存。服务器内存是512G的，装不下。
我们的做法是：
- 先用整数编码把kmer映射到64位整数
- 对每个样本构建字典，但内存限制在2G以内，超过就写入磁盘
- 所有样本的字典合并后，构建CSR矩阵
- 用CSR矩阵做PCA和聚类分析
最终内存峰值控制在120G以内，跑通了整个流程。我去跟课题组汇报的时候，他们说「原来我们之前超内存是因为没做稀疏化」。就这么简单的一个转变，效果天壤之别。
关键知识点总结：
- kmer矩阵天然稀疏，密度通常低于0.0001%
- 小样本（<100G）用字典+整数编码
- 大样本+矩阵运算用CSR
- 布隆过滤只做预处理筛选，不做精确计数
- Python里优先用scipy.sparse，少自己造轮子
对付超大的kmer计数矩阵，别想着一口气全读到内存里。稀疏表示就是让你用更少的内存，装更多的信息。下次遇到内存爆表，先问问自己——「我是不是存了太多0？」。答案往往是肯定的。