4、功能基因注释与丰度矩阵归纳总结
功能基因注释，说白了就是把测序得到的基因片段“翻译”成生物学功能。你拿到的原始数据是一堆ATCG，但你想知道的是——这些微生物在干嘛？产甲烷？降解污染物？还是产生抗生素？
我做宏基因组分析头几年，最怕的就是这步。注释跑完了，结果一看，全是“hypothetical protein”。嗯，那是数据库该更新了。今天咱们把这条路彻底走通。
4.1 功能注释的核心逻辑
序列比对是基础。将预测出的基因（ORF）蛋白序列，怼到已知功能数据库里，找同源——说白了就是“看它像谁”。像KEGG就告诉你代谢通路，COG告诉你直系同源簇，Pfam告诉你蛋白结构域。
我给自己定了个习惯：同时跑至少两个数据库，再用投票机制综合。单个数据库总有盲区，交互相容才靠谱。
4.2 常用数据库与工具选择
下面这张表是我反复对比后的推荐配置。你直接拿来用就行。
| 数据库 | 特点 | 适用场景 | 工具推荐 |
|---|---|---|---|
| KEGG | 通路层级清晰，K编号富集 | 代谢通路分析 | KOBAS, ghostKOALA |
| COG/eggNOG | 直系同源簇分类，速度快 | 整体功能谱 | eggNOG-mapper |
| Pfam | 结构域级，深度大 | 蛋白家族注释 | hmmscan (HMMER) |
| CAZy | 碳水化合物活性酶专精 | 糖类降解/合成 | dbCAN2 |
| ResFinder/CARD | 抗性基因专库 | 耐药基因筛查 | ABRicate, RGI |
我的经验：别只用一个数据库。KEGG注释不到的那些基因，往往在Pfam或CAZy里能找到线索。我曾经被一个未知基因卡了三天，最后是Pfam的“glycosyl transferase family 2”给了解答。
4.3 丰度矩阵的构建
注释完了，你得到的是每个基因的功能标签。但问题来了：这个功能在样本里有多少？ 丰度矩阵就是回答这个问题的。
最直接的算法——基于比对覆盖度。比如你用Bowtie2把reads map回基因，然后用featureCounts计数，每个基因得到一条counts。然后按功能聚合：同属K00123的基因，counts加起来，就是K00123在样本里的丰度。
这里有个关键：要不要做长度校正？我建议必须做。长基因天然有更多read覆盖，不做校正，你的丰度其实反映的是基因长度，不是真实的丰度。常用方法是RPKM（Reads Per Kilobase per Million），或者FPKM、TPM。我个人偏爱TPM，跨样本可比性最好。
实际跑起来，我一般用这一套命令（假设你已经有了gene_count矩阵和基因长度文件）：
# 计算RPKM，然后按功能聚合
awk 'NR==1{print; next} { len=$3; for(i=4;i<=NF;i++) { rpkm=$i*1e9/(len*sum) }; print }' gene_counts.tsv > gene_rpkm.tsv
# 用Python或R批量按KEGG_ko映射聚合
python3 aggregate_by_ko.py -i gene_rpkm.tsv -m ko_map.txt -o ko_matrix.tsv注意：聚合之前，一定检查多对多映射。一个基因可能对应多个KO——这种情况我通常取主要KO。曾经有个样本因为用了“全保留”，把泛素化通路和碳代谢搅在一起，下游分析完全走偏。
4.4 规范化与归一化——避坑指南
你拿到的原始丰度矩阵，直接扔进随机森林？千万别。
我曾经犯过这个错误。第一版模型准确率87%，我高兴坏了。后来发现，全是文库大小差异导致的假关联。丰度随文库大小线性变化，模型学到的其实是“这个样本测了多少数据”。
正确的做法分两步：
- 文库归一化： 用总counts或中位数ratio将样本拉到同一尺度。常用TMM（来自edgeR）、CSS（来自metagenomeSeq）或简单的CPS（counts per sample）。
- 功能级标准化： 每个功能的丰度除以其在所有样本的中位数或均值。这一步消除功能间量级差异，让模型更关注相对变化。
我个人的标准流程：
- 先用脱盐包（DESeq2的varianceStabilizingTransformation）做方差稳定变换。它能同时处理文库大小和低计数功能。
- 再对每个功能做
scale(center=T, scale=T)，转为Z-score。这样随机森林的特征重要性才真正可比。
一句话总结：归一化不是可选项，是必选项。不做归一化的模型，就是在过拟合文库深度。
4.5 实践流程：从注释到矩阵
下面这张SVG图展示了从原始序列到功能丰度矩阵的完整流程——我自己的项目里就贴这个图，方便跟团队对齐。
4.6 矩阵质量评估
矩阵建好之后，别急着建模。先做一轮质量控制：
- 稀疏度： 超过70%的零值？那说明你的功能列表可能太细，建议聚合到更高级别（比如从KO到Pathway）。
- 离群样本： 做一次PCA快速扫描。如果某个样本单独挂在一个角落，很可能是文库问题或污染。
- 功能覆盖率： 看看有多少功能在所有样本中的表达量都低于1。这些低信噪比功能，我建议直接剔除。不然它们会稀释重要特征的权重。
一个小技巧：把矩阵输出前，加一列“功能名+注释来源”，类似“K00123|KEGG|alcohol dehydrogenase”。这样后续做功能富集或通路映射时，你不用来回查表。
好了，到现在你手里应该有一张干净、归一化、高覆盖的功能丰度矩阵。下一站就是随机森林了。不过我觉得——先把这步走稳，模型的事我们回头再聊。