做微生物组和宏基因组分析的人几乎没有谁能绕开“物种注释”这一步。拿到测序数据之后大家最关心的事情其实就两个这一段序列来自哪个物种这个样本里不同物种的相对丰度是多少Kraken2Bracken就是目前社区里解决这两个问题最常用的组合之一。Kraken2用k-mer精确匹配做超快速分类速度比传统的序列比对工具快一个数量级Bracken在这个基础上做丰度重估把“这条reads归到哪个节点”的离散判断换算成“这个物种占了多大比例”的丰度估计。这套组合的安装本身不复杂但真要从零开始跑顺、跑出可靠结果还是有不少细节值得掰开揉碎讲一讲。这篇文章就按我的实操路径从工具选型、安装建库、物种分类、丰度估计到各种报错排查一步步给你过一遍。1. 物种注释选型为什么是Kraken2Bracken这套班子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命中的节点汇总再取这些节点的“最低共同祖先”LCALowest 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 taxconda解析依赖比较慢是出了名的等得人心烦。推荐先装一个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三个可执行文件。记得把目录加进PATHexport PATH/path/to/kraken2-bin:$PATHBracken的源码编译同样简单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列分类标记C表示已分类U表示未分类序列IDNCBI分类号未分类时为0序列长度LCA匹配信息形如0:35 279010:18 279010:17表示这条reads上每个分类节点命中了多少个k-mer这个LCA匹配信息是调试的好帮手。如果一条reads被标注为U你可以看第5列是空的还是有命中但低于confidence阈值如果被标记为C但挂到了很上级的分类说明这条reads的k-mer大量分散在不同物种里这是数据库本身的局限不是代码的问题。--report报告文件则是另一个视角它按分类层级汇总了所有reads每行6列列含义1该分类节点及其所有子节点占全部reads的百分比2归属该节点及其子节点的reads数3仅归属该节点本身的reads数4分类层级代码S种、G属、F科、O目、C纲、P门等5NCBI分类号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参数说明-dKraken2数据库目录路径-t线程数-kKraken2建库时的k-mer长度默认就是35-l你想把reads重分配到哪个分类层级。S代表种G代表属F代表科也可以按需用O、C、P-xkraken2的安装目录脚本需要调用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-merBracken在做丰度重分配时需要知道这个长度来校正概率。如果你测序读长是PE150就写150如果是PE100就写100。写了不对的读长丰度估计会有系统偏差。-l指定你要输出到哪个分类层级。如果你在bracken-build阶段先建了-l S的索引这里就可以填S。两者不要求完全一致bracken运行时会自动调用对应的索引文件但前提是你提前建好了这个层级的索引。稳妥的做法是建库时就把S和G都建一遍分析时灵活切换不用重新跑。4.4 读懂Bracken输出列衔接下游分析Bracken的输出文件同样以制表符分隔核心列如下列名含义name物种学名taxonomy_idNCBI分类号taxonomy_lvl分类层级kraken_assigned_readsKraken2原本直接分到这个节点的reads数added_readsBracken从上级节点重新分配下来的reads数est_total_reads校正后的总reads数等于前两列之和fraction_total_reads该物种相对丰度所有物种合计为1拿到这张表之后下游的分析就自由了。只想看相对丰度直接取fraction_total_reads列排序想转为BIOM格式给QIIME2用可以借助kraken-biom这类小工具想画堆叠柱状图常见的做法是先按属水平汇总再挑出丰度最高的前20个属awk -F \t $3G {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 常见报错与排查速查表现象常见原因解决办法Error: Unable to open database数据库目录不完整或路径写错检查hash.k2d、opts.k2d、taxo.k2d三个文件是否都在A database was not found和上面一样多半是Kraken2和Bracken用了不同的数据库路径把-d参数统一避免一个指向绝对路径、一个指向相对路径taxonomy file is missing构建数据库时taxonomy没下载成功重跑kraken2-build --download-taxonomy --dbBracken报找不到kmer2read_distr文件Bracken索引没有生成或-k参数和建库时不一致重新运行bracken-build确认-k值Kraken2 executable not found编译安装的kraken2没有加入PATH确认which kraken2能找到程序再跑bracken-build并指定-x--paired报两端reads数不一致R1和R2序列数量对不上可能是QC时没过滤干净用seqkit stats检查两端reads数必要时重新跑fastp分类结果里unclassified占比异常高数据库和样本类型不匹配或序列质量太差换PlusPF或自建库先做质量控制再降低confidence阈值建库过程中途被OOM杀掉内存不够换更大的内存机器或用--fast-build占内存但速度更快以上这些报错里我遇到最多的就是第一种和第四种。尤其是Bracken索引的问题很多人以为bracken命令能跑就能出结果结果运行到一半才报缺文件浪费了不少时间。我的习惯是数据库构建完成之后立刻把hash.k2d、opts.k2d、taxo.k2d以及所有kmer2read_distr文件放到同一个目录并且坚决不随意挪动后面一劳永逸。5.4 一条我保留了很久的小习惯最后再分享一个很多老玩家不会写在文档里的细节正式跑到大批量样本之前先拿一个样本走通全流程并且把中间文件都检查一遍。我一般会先看Kraken2报告的前几行确认已知高丰度物种是不是真的出现在结果里再看Bracken输出的est_total_reads有没有明显的异常值。如果第一步就发现数据库选择有问题还能及时止损总比全批次跑完再回炉重做强。另外分类结果里如果某个物种的added_reads远大于kraken_assigned_reads你就要留意了这通常意味着这个物种和它的近缘物种共享了大量k-mer丰度估计的置信度其实并不高写文章时要慎下结论。这一套流程跑顺之后你会发现Kraken2Bracken这条流水线真是又稳又省心后续换数据、换读长也只需要在参数层面微调而已。