企业官网建设流程全解析

做微生物组和宏基因组分析的人，几乎没有谁能绕开“物种注释”这一步。拿到测序数据之后，大家最关心的事情其实就两个：这一段序列来自哪个物种？这个样本里不同物种的相对丰度是多少？Kraken2+Bracken就是目前社区里解决这两个问题最常用的组合之一。Kraken2用k-mer精确匹配做超快速分类，速度比传统的序列比对工具快一个数量级；Bracken在这个基础上做丰度重估，把“这条reads归到哪个节点”的离散判断，换算成“这个物种占了多大比例”的丰度估计。这套组合的安装本身不复杂，但真要从零开始跑顺、跑出可靠结果，还是有不少细节值得掰开揉碎讲一讲。这篇文章就按我的实操路径，从工具选型、安装建库、物种分类、丰度估计到各种报错排查，一步步给你过一遍。

1. 物种注释选型：为什么是Kraken2+Bracken这套班子

1.1 Kraken2的快，建立在k-mer精确匹配上

先聊一个很多人问过我的问题：同样是做宏基因组物种注释，为什么不用BLAST或者DIAMOND，非得用Kraken2？

核心差别在算法路径上。BLAST这类比对工具，是拿每一条reads去和参考数据库做局部联配，算相似性得分，然后再通过得分阈值判定分类。这个过程准确是准确，但计算量非常大，一条双端150bp的reads拆成两段，每段都要和几万个基因组做比对，几百万条reads跑下来，计算集群也得等上大半天。

Kraken2没有走这条老路。它把参考基因组先切成固定长度的k-mer（默认是35bp），存入一个哈希表，然后对每条reads也做同样的切分，逐个k-mer去哈希表里查命中。每个k-mer命中后会指向一个分类学节点，算法会把整条reads上所有k-mer命中的节点汇总，再取这些节点的“最低共同祖先”（LCA，Lowest Common Ancestor）作为这条reads的分类结果。整个过程没有联配计算，全是哈希查询，所以速度能比传统比对快出几个量级。

但这里藏着一个先天的短板。LCA策略本身是为了保守，避免把reads错误地分到某个具体的种。它一旦发现这句话的k-mer在多个物种里都出现了，就会往上一级、两级甚至更多级回溯，最后给出一个“属”“科”甚至“目”级别的注释。你看着结果里一大片“Escherichia属”“Bacteroides属”，心里会犯嘀咕：这到底是个什么丰度？

1.2 Bracken补上的，正是丰度估计这一环

Kraken2给出的分类结果，本质上是“哪条reads大致来自哪个分类节点”，不是“某个物种在样本里占多少比例”。如果你直接拿Kraken2报告里的reads数去算丰度，那些被LCA挂到属或科一级的reads就全丢了，物种层面的占比会被明显低估。

Bracken的出现就是为了解决这件事。它的思路是：同一套数据库里，每个物种的基因组大小不同、k-mer组成不同，那在固定读长下，理论上每个物种能被“命中”的reads数也是有规律的。Bracken会读入Kraken2的分类报告，再结合数据库中每个物种的k-mer分布信息，用一个类似贝叶斯重分配的方式，把那些被归类到高层级节点的reads，按比例重新“摊还”到下属的各个物种头上，最终输出物种级别的相对丰度估计。

如果你还是觉得抽象，可以这样理解：Kraken2像个快递分拣员，看到地址不完整只能给你送到“市”这一级；Bracken则拿着各“小区”的住户名单和建筑面积，重新估算每个小区实际住着多少人。

所以这套组合的定位非常明确：Kraken2负责“快、准地把reads归到某个分类层级”，Bracken负责“在分类基础上补全丰度信息”。两者不是二选一的关系，而是前后串联的一条流水线。

2. 安装前的功课：环境依赖与三条可行路径

2.1 最省心的方式：conda/mamba一条命令

如果你不是有特殊限制，我强烈建议直接用conda或者mamba装，别在这上面浪费时间。创建一个独立环境，把两个工具都放进去，避免污染基础环境，也方便以后整体删除重建：

conda create -n tax -y -c bioconda -c conda-forge kraken2 bracken conda activate tax

conda解析依赖比较慢是出了名的，等得人心烦。推荐先装一个mamba，然后：

mamba create -n tax -y -c bioconda -c conda-forge kraken2 bracken mamba activate tax

装完先确认一下工具是否真的可用：

kraken2 --version bracken -h

如果conda源里有了对应的二进制包，这两条命令基本不会出错。唯一要留意的是，Kraken2和Bracken的版本要匹配，conda默认帮你处理好了，不太需要操心。

2.2 源码编译：适合离线服务器和追新版的情况

有些内网集群不能访问公网conda源，或者你想用最新开发版，那就得走源码编译。Kraken2的编译过程很轻量，主要依赖g++、make、rsync和zlib，系统一般都会带：

git clone https://github.com/DerrickWood/kraken2.git cd kraken2 ./install_kraken2.sh /path/to/kraken2-bin

编译完成后会在指定目录生成kraken2、kraken2-build、kraken2-inspect三个可执行文件。记得把目录加进PATH：

export PATH=/path/to/kraken2-bin:$PATH

Bracken的源码编译同样简单：

git clone https://github.com/jenniferlu717/Bracken.git cd Bracken bash install_bracken.sh

这个脚本会编译几个perl的C扩展模块，用来加速k-mer分布统计。如果报错，十有八九是系统缺少g++或perl的开发头文件，装上对应包再跑一遍就行。

这里有个非常容易踩的坑：bracken-build脚本在运行时会调用kraken2和kraken2-inspect，这两个程序必须在PATH里能找到。很多人单独编译了kraken2却没有写进PATH，结果Bracken建库时一直报“Kraken2 executable not found”，排查半天才发现是环境变量的问题。

2.3 数据库选型：这一步直接决定你结果的上限

安装好工具只是热身，真正影响结果的是数据库。Kraken2官方提供多种预建数据库，选哪个完全取决于你的数据来源和分析目标。

最常用的标准库（Standard）包含RefSeq里面细菌、古菌、病毒、质粒、人类基因组以及UniVec_Core等数据，建好后数据库文件在8GB左右，分类时峰值内存大约需要30GB到50GB，适合绝大多数宏基因组（WGS）样本。如果你想同时覆盖真菌和原生动物，比如做环境样本或某些特殊临床样本，可以下载PlusPF库，这个库在标准库基础上加入了真菌、原生动物数据，体积和内存需求都会翻好几倍。

我的建议是：普通肠道菌群、土壤宏基因组，先跑标准库；水体、植物根际或者你明确知道样本里有大量真菌的，再考虑PlusPF。库选大了，不仅建库慢、占内存，分类时很多随机噪声也会跟着进来，反而增加误注释率。

自建库则适合有明确目标的场景，比如只关心某个物种属内的鉴定，或者有参考基因组不在NCBI库里。自建库的流程后面单独讲。

3. Kraken2实操：建库、分类与结果文件解读

3.1 标准库的分步构建

官方一条命令就能建标准库：

kraken2-build --standard --threads 32 --db ~/db/kraken2_std

这条命令会先下载NCBI的分类学信息（taxonomy），再下载RefSeq里各个库的序列，最后切k-mer建索引。整个过程对网络要求很高，NCBI官方推荐至少有50GB可用磁盘空间，我实际跑下来，包括临时文件和数据库本体，预留30GB以上比较稳妥。构建期间的内存波动很大，我自己在一台48GB内存的服务器上跑过，高峰期差点被OOM杀掉。

更可控的做法是分步执行。网络经常会在下载大文件时断掉，一条龙命令一旦中断就要从头再来，非常折磨人。我习惯这样拆成三步：

kraken2-build --download-taxonomy --db ~/db/kraken2_std kraken2-build --download-library bacteria --db ~/db/kraken2_std kraken2-build --download-library archaea --db ~/db/kraken2_std kraken2-build --download-library viral --db ~/db/kraken2_std kraken2-build --download-library plasmid --db ~/db/kraken2_std kraken2-build --download-library human --db ~/db/kraken2_std kraken2-build --build --db ~/db/kraken2_std --threads 32

注意：--download-library可以反复执行，它会断点续传；已经下好的库不会重复下载。建库完成后，数据库目录里会出现hash.k2d、opts.k2d、taxo.k2d三个文件，后续分类全靠这三个文件。如果缺了任何一个，分类时Kraken2会直接报错，所以平时要么不挪动这些文件，要么三个一起拷贝。

3.2 单端、双端与批量样本的分类命令

数据库就绪后，分类就是一条命令的事。单端reads这样跑：

kraken2 --db ~/db/kraken2_std --threads 16 \ --output sample.kraken --report sample.report \ sample.fq.gz

双端reads加一个--paired参数：

kraken2 --db ~/db/kraken2_std --threads 16 --paired \ --output sample.kraken --report sample.report \ sample_R1.fq.gz sample_R2.fq.gz

如果你有一批样本要跑，别一个个手动敲命令，直接用for循环处理：

for sample in $(cat sample_list.txt); do kraken2 --db ~/db/kraken2_std --threads 16 --paired \ --output ${sample}.kraken --report ${sample}.report \ ${sample}_R1.fq.gz ${sample}_R2.fq.gz done

这里我强烈建议加上一个参数：

--confidence 0.2

这个confidence参数的范围是0到1，代表了一条reads上至少有多少比例的k-mer命中某个分类节点，才会保留该分类结果。默认值是0，意思是只要有任意k-mer命中就直接分类，这样很容易被基因组间的共有序列带偏。宏基因组数据我自己习惯用0.2，如果你觉得unclassified比例太高，损失了太多信息，可以降到0.1；如果是纯菌株鉴定或者对特异性要求高，用0.3甚至0.4都可以。

另一个值得了解的参数是--minimum-hit-groups，默认值2。它要求一条reads上至少要有两组相邻的k-mer都命中同一个分类，才算有效分类，用来抑制随机偶发命中造成的假阳性。大部分情况下保持默认即可，不需要动。

3.3 Kraken2输出文件和report报告怎么读

Kraken2的--output文件默认是制表符分隔，每条reads一行，共5列：

1. 分类标记：C表示已分类，U表示未分类
2. 序列ID
3. NCBI分类号，未分类时为0
4. 序列长度
5. LCA匹配信息，形如0:35 279010:18 279010:17，表示这条reads上每个分类节点命中了多少个k-mer

这个LCA匹配信息是调试的好帮手。如果一条reads被标注为U，你可以看第5列是空的还是有命中但低于confidence阈值；如果被标记为C但挂到了很上级的分类，说明这条reads的k-mer大量分散在不同物种里，这是数据库本身的局限，不是代码的问题。

--report报告文件则是另一个视角，它按分类层级汇总了所有reads，每行6列：

看这份报告时，我最常做的一件事是先扫一遍第2列占比最高的分类。宏基因组样本里出现一个注释为unclassified的比例特别高，通常是数据库覆盖不够，而不是样本本身没有已知物种。

4. Bracken实操：从Kraken2报告到物种丰度表

4.1 为什么不能直接拿Kraken2报告当丰度

这个问题我几乎每次培训都会被问。你看，Kraken2的报告第二列写着某个物种有5000条reads，那它的丰度不就是5000除以总reads数吗？

问题出在那部分被LCA挂到更高层级的reads上。比如一条reads明明来自大肠杆菌，但它的k-mer在志贺氏菌、沙门氏菌里也能命中，Kraken2为了不犯错，会把它标记为肠杆菌科（Enterobacteriaceae）而不是大肠杆菌。这样你在种水平统计reads数时，这部分reads就被漏掉了。样本里亲缘关系越近的物种越多，漏掉的比例就越大。

Bracken干的事情，就是把这些被“保守处理”的reads，按照数据库里各物种的k-mer出现概率，重新分配到具体物种头上，最后给出一个校正后的丰度估计。

4.2 先给数据库生成Bracken需要的索引

Bracken不是直接读Kraken2数据库就能跑的，它需要先基于Kraken2数据库生成一份k-mer分布统计文件。这一步必须做，而且要和Kraken2建库时的k-mer长度保持一致：

bracken-build -d ~/db/kraken2_std -t 32 -k 35 -l S -x /path/to/kraken2-bin

参数说明：

• -d：Kraken2数据库目录路径
• -t：线程数
• -k：Kraken2建库时的k-mer长度，默认就是35
• -l：你想把reads重分配到哪个分类层级。S代表种，G代表属，F代表科；也可以按需用O、C、P
• -x：kraken2的安装目录，脚本需要调用kraken2和kraken2-inspect

这一步会遍历整个数据库里所有物种的k-mer信息，所以比较耗时。我建标准库的Bracken索引，在32线程的机器上大约要跑1到2个小时。如果之后你换了k-mer长度重建Kraken2数据库，Bracken索引也必须用相同的-k参数重新生成，否则后面会报一大串“file not found”或者干脆给出完全错误的结果。

4.3 运行Bracken：记好-r这个关键参数

索引生成好之后，命令非常简洁：

bracken -d ~/db/kraken2_std \ -i sample.report \ -o sample.bracken \ -r 150 -l S -t 32

-i输入Kraken2的report文件，-o输出结果文件，-d指向同一个数据库目录。这里重点说两个参数。

-r是reads的平均长度。为什么它重要？因为Kraken2命中k-mer的数量和reads长度直接相关，150bp的reads比100bp的reads平均命中更多的k-mer，Bracken在做丰度重分配时，需要知道这个长度来校正概率。如果你测序读长是PE150，就写150；如果是PE100，就写100。写了不对的读长，丰度估计会有系统偏差。

-l指定你要输出到哪个分类层级。如果你在bracken-build阶段先建了-l S的索引，这里就可以填S。两者不要求完全一致，bracken运行时会自动调用对应的索引文件，但前提是你提前建好了这个层级的索引。稳妥的做法是建库时就把S和G都建一遍，分析时灵活切换，不用重新跑。

4.4 读懂Bracken输出列，衔接下游分析

Bracken的输出文件同样以制表符分隔，核心列如下：

拿到这张表之后，下游的分析就自由了。只想看相对丰度，直接取fraction_total_reads列排序；想转为BIOM格式给QIIME2用，可以借助kraken-biom这类小工具；想画堆叠柱状图，常见的做法是先按属水平汇总，再挑出丰度最高的前20个属：

awk -F '\t' '$3=="G" {print $1, $NF}' sample.bracken | sort -k2 -rn | head -20

这里$NF就是最后一列的fraction_total_reads，$1是物种名，$3判断层级是否属于属。你要是愿意，也可以直接在R里读入表格，用dplyr按taxonomy_lvl分组后summarise，道理是一样的。

5. 实战进阶：自建数据库、资源调优与报错排查

5.1 自建数据库：把自己关注的参考基因组加进来

有时候NCBI的库太大，跑起来费劲；有时候你又想加入一批私有参考基因组。这时候自建一个小而精的数据库是正解。流程并不复杂：

第一步，准备FASTA文件。如果你希望这条序列被正确识别到某个物种，建议把序列头部写成Kraken2能识别的格式：

>seqid|kraken:taxid|561 ATCGATCG...

这里的561是大肠杆菌的NCBI分类号。如果你的序列没有对应的NCBI分类号，也可以单独提供一个映射文件，在构建时用--taxids选项指定。

第二步，把序列加进数据库并构建：

kraken2-build --add-to-library my_genomes.fna --db ~/db/custom_db kraken2-build --build --db ~/db/custom_db --threads 32

第三步，生成对应的Bracken索引：

bracken-build -d ~/db/custom_db -t 32 -k 35 -l S -x /path/to/kraken2-bin

自建库有一个很容易踩的坑：如果你加入的序列没有正确映射到物种层级的分类号，构建时Kraken2会把它挂到root节点，导致分类时大量reads被标记为unclassified。用kraken2-inspect可以快速检查一下建好的库：

kraken2-inspect --db ~/db/custom_db | head -20

这个命令能预览数据库中每个分类节点下有多少序列，如果发现某个物种下面序列数明显不对，多半是taxid映射出了问题。

5.2 资源和性能的实用经验值

Kraken2分类时的内存消耗主要来自数据库加载。标准库分类时大约需要30GB到50GB内存，这是哈希表常驻内存的开销，和你开多少个线程没关系。所以一台32GB内存的机器跑标准库会非常吃力，要么换库，要么用--memory-mapping参数，让数据库通过内存映射的方式从磁盘读取，速度会慢一些，但内存占用能明显降下来。

我实测下来，一台16线程、64GB内存的机器，跑一个5000万条reads的宏基因组样本，Kraken2分类大约需要20到30分钟，Bracken重估丰度只需要2到3分钟。瓶颈往往不在CPU，而在读取压缩FASTQ时的解压速度。如果你有大量样本要跑，建议先把文件从机械盘挪到SSD上，或者用zcat预解压到临时目录，IO等待能少很多。

线程数也不是越大越好。超过32线程后，Kraken2的性能提升变得很不明显，因为哈希查询本身是内存密集型的，内存带宽会先到瓶颈。与其无脑堆线程，不如同时跑两三个样本，每个分配16线程，整体吞吐会更高。

5.3 常见报错与排查速查表

以上这些报错里，我遇到最多的就是第一种和第四种。尤其是Bracken索引的问题，很多人以为bracken命令能跑就能出结果，结果运行到一半才报缺文件，浪费了不少时间。我的习惯是，数据库构建完成之后，立刻把hash.k2d、opts.k2d、taxo.k2d以及所有kmer2read_distr文件放到同一个目录，并且坚决不随意挪动，后面一劳永逸。

5.4 一条我保留了很久的小习惯

最后再分享一个很多老玩家不会写在文档里的细节：正式跑到大批量样本之前，先拿一个样本走通全流程，并且把中间文件都检查一遍。我一般会先看Kraken2报告的前几行，确认已知高丰度物种是不是真的出现在结果里；再看Bracken输出的est_total_reads有没有明显的异常值。如果第一步就发现数据库选择有问题，还能及时止损，总比全批次跑完再回炉重做强。另外，分类结果里如果某个物种的added_reads远大于kraken_assigned_reads，你就要留意了，这通常意味着这个物种和它的近缘物种共享了大量k-mer，丰度估计的置信度其实并不高，写文章时要慎下结论。这一套流程跑顺之后，你会发现Kraken2+Bracken这条流水线真是又稳又省心，后续换数据、换读长，也只需要在参数层面微调而已。