说实话第一次接触宏基因组数据的时候看着几千万条reads完全不知道从哪下手。想知道一个样本里有哪些微生物、它们的丰度是多少光靠比对NCBInt或者BLAST这种传统方式一跑就是好几天而且blast的结果处理起来也麻烦。后来同事给我推荐了Kraken2说这个软件速度快得离谱比对完还有Bracken做丰度校正两条命令就能出物种组成结果。我当时半信半疑装上之后跑了一个肠道菌群样本几百万条reads十分钟不到就分完了确实被惊到了。这个组合现在是宏基因组物种注释领域的标配工具也是大多数生信流程的首选方案。Kraken2负责把测序reads快速分类到物种分类学节点上Bracken则在其基础上估计每个物种的真实丰度。整套流程对新手非常友好安装不复杂使用命令也不多最关键的是运行速度快普通服务器就能跑。这篇文章就围绕安装与使用这条主线把我从零开始踩过的坑、调过的参数、看过的文档系统整理一遍给正准备入坑的朋友一份可以直接照着做的参考。1. 为什么是Kraken2Bracken这对组合1.1 物种注释到底做什么物种注释简单说就是拿到测序数据之后判断每一条reads是来自哪个物种。听起来很简单但实际操作里有两个核心难点。第一个难点是数量。一个宏基因组样本动辄几千万条reads每条reads要做一次序列比对如果用传统BLAST或者Bowtie2时间和计算资源都会爆炸。第二个难点是精度。很多微生物之间序列相似度很高尤其在不同菌株之间区分起来并不容易。Kraken2的核心思路不是把所有序列全部两两比对而是先把参考基因组切碎成固定长度的短序列片段为每个片段建立索引计算出一个k-mer集合。然后每来一条reads就把它也切成一堆k-mer去索引里查这些k-mer分别落在哪些物种上再基于所有k-mer的综合投票来确定这条reads的分类归属。这种方式的优势在于查询过程不需要逐步扫描数据库只需要一次哈希查找所以速度极快。Kraken2的分类结果还会同时给出该分类在分类学层级上的位置从界、门、纲、目、科、属到种。1.2 Kraken2的分类精度问题但Kraken2有个明显的短板它输出的只是reads被分到了哪个taxon节点上而不是这个物种在样本里的真实丰度。举个实际例子如果给Kraken2输入100条reads其中80条分给了物种A20条分给了物种B那么Kraken2报告里写的就是A占80%、B占20%。这看上去挺合理的但问题在于不同物种的基因组大小差异很大。有的细菌基因组只有2Mb有的真菌基因组有30Mb也就是说在丰度相同的情况下基因组大的物种自然会被分配到更多reads。另外Kraken2分类结果里还有一类特殊标记比如一条reads同时匹配到多个物种时它会被分配到它们共同的分类学祖先节点上送到属甚至科这个级别这类reads如果直接算进物种丰度会出现虚高的问题。这就是Bracken存在的意义。Bracken全称是Bayesian Reestimation of Abundance with KrakEN它能利用参考数据库里各物种的k-mer分布信息把被分配到属、科等更高级别节点上的reads重新分布到具体物种上同时校正基因组大小带来的偏差最终输出更接近真实情况的物种丰度表格。1.3 在什么场景下选这对组合Kraken2Bracken最适合的场景是宏基因组样本的物种组成分析典型用途有肠道菌群研究、环境微生物多样性调查、临床样本病原体快速筛查。如果目标是分析16S扩增子数据Kraken2也可以跑但像QIIME2加Greengenes这类专门流程可能更适合如果目标是做菌株级别的精细分型那Kraken2就有些吃力了需要用到StrainPhlAn这类工具。不过做绝大多数属种水平的组成分析Kraken2Bracken几乎是最省心、最高效的选项。2. 安装之前先想清楚这三件事2.1 硬件条件评估安装和使用Kraken2之前先评估一下自己的机器配置。Kraken2本身是C写的对CPU的利用率很高支持多线程一般8核16线程的服务器就够了。内存方面需要特别注意因为k-mer索引是直接加载到内存里的数据库大小直接决定内存需求量。Kraken2官方提供了几种预构建数据库数据库类型磁盘占用内存需求说明MiniKraken (8GB)约8GB约8GB小规模测试或快速预览Standard (标准库)约150GB约75GB常规宏基因组分析覆盖细菌、古菌、病毒等PlusPF约200GB以上约100GB以上标准库加真菌、原生生物使用范围更广我自己第一次用的时候只准备了32GB内存的机器跑标准库加载到一半直接OOM内存溢出后来换了能扩容的机器才解决。如果你用笔记本电脑跑小样本可以先从MiniKraken库开始如果跑正式项目建议至少64GB内存和几百GB剩余磁盘空间。2.2 数据库选型思路数据库的选择直接关系到分析结果的质量这块建议做正式分析时不要图省事。标准数据库包含细菌、古菌、病毒以及人类基因组的参考序列覆盖宏基因组研究中最常见的微生物类群。如果样本可能含有真菌使用PlusPF数据库更合适它额外加入了真菌和原生生物的数据。MiniKraken数据库虽然下载快、占用小但覆盖物种有限不建议用于正式科研数据。另一个容易被忽略的点是数据库版本。Kraken2官方会定期更新数据库内容不同版本的数据库之间分类结果可能存在差异。同一项目中的所有样本务必使用同一版本的数据库否则不同批次之间的比较会有系统误差。我自己的建议是建库时间和版本记录下来写论文方法学部分也得写清楚。2.3 安装方式取舍Kraken2和Bracken的安装方式主要有两种conda包管理器安装和源码编译安装。conda方式最省事几分钟装完我推荐大多数人用。源码编译也简单主要价值在于自定义编译参数优化性能不过对普通用户来说收益不大。Bracken需要注意版本匹配问题Bracken必须和Kraken2的主版本兼容。比如Kraken2是2.x版本Bracken也要用对应的2.x版本。如果直接用conda同时安装Kraken2和Brackenconda会自动解决版本依赖比较省心。3. Kraken2与Bracken完整安装步骤记录3.1 用conda安装Kraken2如果你的分析环境还没有配置conda需要先装一个Miniconda或Anaconda。安装过程不展开直接说装好conda之后的操作。先创建一个独立环境避免不同软件依赖互相干扰conda create -n kraken2_env python3.8 conda activate kraken2_env然后安装Kraken2conda install -c bioconda kraken2装完验证一下kraken2 --version正常情况下会输出类似Kraken version 2.1.3这样的信息。如果提示找不到命令检查一下环境是否激活或者用which kraken2看看安装路径。3.2 源码编译方式安装源码方式适合没有conda或需要自定义优化的场景。先下载源码git clone https://github.com/DerrickWood/kraken2.git cd kraken2 ./install_kraken2.sh /path/to/install/directory脚本执行完后Kraken2会安装到指定目录下的bin文件夹里。为了使用方便把bin目录加入环境变量echo export PATH/path/to/install/directory/bin:$PATH ~/.bashrc source ~/.bashrc源码编译的好处是可以指定编译参数比如启用某些CPU指令集优化但对绝大多数流程影响不大我平时还是用conda版本。3.3 安装BrackenBracken的官方GitHub仓库是jeniferjs/bracken。conda安装最简单conda install -c bioconda bracken装完检查bracken -v输出版本号说明安装成功。源码安装也简单git clone https://github.com/jeniferjs/bracken.git cd bracken chmod x install_bracken.sh ./install_bracken.shBracken安装完成后bin目录下会生成几个脚本bracken、bracken-build、est_abundance.py等。bracken-build这个脚本是用来为Kraken2数据库构建丰度分布文件的后面建库时要用到所以路径一定要记清楚。3.4 快速验证安装是否可用安装完成后最稳妥的方式是跑一个最小测试。先下载MiniKraken数据库试跑一下或者用自带的测试数据。Kraken2源码目录里通常有一个example文件夹里面放了示例数据。kraken2 --db minikraken_db_v2 --threads 4 example.fa能正常输出分类结果文件说明整个安装链路没问题。我在第一次安装时跳过这个测试直接跑正式数据结果建库路径写错了排查了半天才发现问题。装完先跑最小测试这个习惯能省很多时间。4. 数据库下载与构建这里最耗时4.1 使用官方预构建数据库Kraken2官方提供预构建数据库直接下载解压就能用省去自己构建的漫长时间。下载地址在Kraken2的GitHub仓库“Database”部分有提供或者用下面的命令从官方AWS存储拉取# 下载标准库大小约150GB wget https://genome-idx.s3.amazonaws.com/kraken/k2_standard_2023.tar.gz # 解压到指定目录 mkdir -p /path/to/kraken2_db tar -xvzf k2_standard_2023.tar.gz -C /path/to/kraken2_db下载标准库需要较长时间网络条件差的话建议使用支持断点续传的下载工具或者分时段下载。解压后的目录里包含hash.k2d、opts.k2d、taxo.k2d三个核心文件还有一个seqid2taxid.map映射文件。使用预构建数据库时需要确保Bracken的kmer分布文件也一起下载。最新的标准库压缩包中通常已经包含了Bracken需要用的database150mers.kmer_distrib等文件。如果没有还需要额外运行Bracken的构建步骤。有一个容易踩的坑下载的数据库版本和Bracken的期望版本不匹配运行bracken的时候会报错中断所以下载时留意一下数据库自带文件的完整性。4.2 手动构建自定义数据库如果研究物种比较特殊或者数据里有大量参考数据库没覆盖的物种就需要自己构建数据库。手动构建分三步走。第一步下载分类学信息kraken2-build --download-taxonomy --db $DBNAME这里的$DBNAME是自定义的数据库目录名要先创建好。这个命令会从NCBI下载taxdump文件构建分类学树。第二步下载参考序列按分类学组别下载细菌库kraken2-build --download-library bacteria --db $DBNAME可选类别包括bacteria、viral、fungi、archaea、protozoa等。如果下载中断重新执行命令会继续下载未完成的部分。第三步构建k-mer数据库kraken2-build --build --db $DBNAME --threads 16这一步会把下载下来的参考序列切碎成k-mer并建立索引过程非常消耗内存和CPU。标准库构建时内存不够就会直接报错建议至少64GB。构建完成后数据库目录下会出现hash.k2d文件就说明建好了。4.3 为Bracken构建kmer分布文件Bracken要想做丰度校正必须先针对Kraken2数据库构建一个k-mer分布文件。这一步使用bracken-build脚本完成。bracken-build -d $DBNAME -t 16 -k 35 -l 150参数解释-d指定数据库目录-t指定线程数-k指定k-mer长度必须和Kraken2建库时使用的k-mer长度一致默认是35-l指定reads长度根据测序reads的实际长度来一般宏基因组是150bp构建完成后数据库目录里会出现类似kmer_distrib_150的文件。Bracken构建过程会先对数据库里的参考序列重新做一次Kraken2分类生成所有物种的k-mer分布信息这个过程同样需要较长时间建议后台运行并记录日志。nohup bracken-build -d $DBNAME -t 16 -k 35 -l 150 bracken_build.log 21 期间用tail -f bracken_build.log观察进度。构建失败时日志里会有明确报错信息最常见的是Kraken2版本不兼容报错内容类似于Unknown kraken2 database format这时候要检查Kraken2版本是否太旧Bracken要求2.0.7以上。4.4 数据库管理的几个习惯Kraken2数据库文件很大管理不好容易出问题。先说磁盘规划预构建库解压时需要的临时空间差不多是压缩包两倍大小比如150GB的压缩包解压时需要预留300GB以上空间。其次数据库路径建议用固定目录并把路径写入环境变量或者分析流程配置里避免每次运行命令时手写路径出现差错。最后不同版本的数据库不要放在同一个目录里命名时带上版本号比如k2_standard_2023这样的命名时间久了也不容易搞混。5. 实战完整跑通一个Kraken2Bracken分析5.1 准备工作与测试数据为了演示完整流程我准备了一份模拟的宏基因组双端测序数据大约100万条reads模拟成一个肠道菌群样本。实际操作时可以先从SRA下载公开宏基因组数据或者自己生成模拟数据做测试。输入数据的格式通常是FASTQ。如果数据是原始下机数据带有接头建议先用fastp或Trimmomatic做质控和去接头不要拿原始数据直接跑分类。质控干净虽然会损失一部分reads但能减少假阳性分类。我自己的经验是接头残留确实会导致一些reads错误地比对到人工序列上进而影响物种注释。准备一个样本清单文件每个样本一行格式为样本名加路径方便批量处理。5.2 Kraken2分类命令详解Kraken2分类的标准命令kraken2 --db /path/to/k2_standard_2023 \ --threads 16 \ --paired \ --output sample1.kraken \ --report sample1.kraken_report \ --use-mpa-style \ sample1_R1.fastq.gz sample1_R2.fastq.gz逐项解释参数--db指定数据库路径--threads线程数按机器配置调整一般设为核心数的一半或全部--paired输入为双端reads此模式下两条reads都匹配到同一物种才算数--output分类结果输出文件每行一条reads格式为“分类状态 taxid readID 序列 分类路径”--report汇总报告文件展示每个taxon的reads数和比例--use-mpa-style以MPA样式输出报告便于后续合并多个样本运行结束后sample1.kraken_report的文件内容类似94.32 943201 943201 U 0 unclassified 5.68 56782 45210 R 1 root 5.68 56782 0 R1 131567 cellular organisms ...第一列是比例第二列是该分类及其子分类覆盖的reads总数第三列是该分类本身直接分配的reads数第四列是分类学级别U表示未分类、R表示根节点、G表示属、S表示种。--use-mpa-style生成的报告格式有些不同更适合下游合并和热图可视化。如果只跑一个样本两种格式差别不大如果跑多个样本要合并成丰度表MPA格式会更方便。5.3 Bracken丰度校正命令详解得到Kraken2报告文件后接下来用Bracken做丰度重估。运行前先确定目标分类层级。通常分析物种组成用种水平命令如下bracken -d /path/to/k2_standard_2023 \ -i sample1.kraken_report \ -o sample1.bracken \ -r 150 \ -l S \ -t 16参数解释-d数据库目录-i输入Kraken2报告文件-o输出Bracken丰度结果-rreads长度必须与构建Bracken kmer分布时设置的-l一致-l分类层级S表示种SpeciesG表示属Genus-t线程数Bracken运行完成后会生成两个文件sample1.bracken是最终的丰度估算文件另一个是sample1.bracken_report后期可视化常用到的就是后一个。看一个输出片段示例name taxonomy_id taxonomy_lvl kraken_assigned_reads added_reads new_est_reads fraction_total_reads Escherichia coli 562 S 12000 3500 15500 0.1123 Bacteroides vulgatus 821 S 8000 1200 9200 0.0667每一行代表一个物种重点是new_est_reads和fraction_total_reads这两列前者是该物种校正后的reads数后者是该物种在样本中的相对丰度。把fraction_total_reads加总就能得到完整的物种组成谱。5.4 多样本批量处理与合并实际项目中很少只测一个样本几十个样本批量跑时需要对处理流程做点小小优化比如用循环把每个样本提交给后台任务执行。先写一个Shell循环分别对每个样本运行Kraken2和Bracken为了不让流程卡在单条命令上后台任务加日志记录for sample in sample1 sample2 sample3; do kraken2 --db /path/to/k2_standard_2023 \ --threads 8 \ --paired \ --output ${sample}.kraken \ --report ${sample}.kraken_report \ ${sample}_R1.fastq.gz ${sample}_R2.fastq.gz done全部跑完后用KrakenTools里的combine_kreports.py脚本把多个Kraken2报告合并成一个表生成以分类学单元为行、样本为列的丰度矩阵。combine_kreports.py --report sample1.kraken_report sample2.kraken_report sample3.kraken_report --output combined_report.tsvcombined_report.tsv就是下游做热图、PCoA、差异分析的标准输入表格。这里列一下KrakenTools的安装方式它是个Python工具包直接从GitHub拉下来就能用不依赖额外的包git clone https://github.com/jenniferlu717/KrakenTools.git对丰度表做后续分析时我一般先用vegan包或phyloseq包做多样性分析再用ggplot2画堆叠柱状图展示物种组成。这一步已经跨入统计可视化范畴不做展开提醒一点Bracken输出的丰度是相对丰度做差异比较时需要考虑闭合效应必要时转成绝对丰度或者用ALDEx2这类专门工具。6. 常见问题排查与经验技巧6.1 运行过程中的典型报错把这一年多被问得最多的问题和解决思路整理成一个速查表报错信息可能原因解决办法Error: Unable to open database数据库路径错误或文件损坏检查路径、重新解压数据库Out of memory内存不足数据库加载失败换大内存机器或用MiniKrakenk-mer file hash.k2d not found数据库不完整确认构建步骤完成hash.k2d必须存在Bracken: Error reading distribution fileBracken的kmer分布文件未构建或不匹配运行bracken-build生成对应reads长度的分布文件Unknown kraken2 database formatKraken2版本过旧升级Kraken2到2.0.7以上Cannot open taxo.k2d数据库目录结构不全重新下载taxo.k2d文件6.2 运行速度的调优思路Kraken2的运行速度受三个因素影响CPU线程数、数据库类型、输入文件大小。线程数上来之后速度提升明显但也不是无限制增长一般16核以内线性扩展还行再往上提升就变缓了。磁盘IO也是一个影响点机械硬盘跑大型数据库时启动阶段会明显偏慢有条件用SSD能省不少时间。输入文件建议用压缩格式FASTQ.GZ直接输入Kraken2支持自动解压不用手动解压可以省下大量临时磁盘空间和IO时间。我第一次跑的时候为了保险把数据都解压成了明文的fastq白白占了三倍磁盘空间后续发现直接喂压缩文件完全没问题。6.3 Bracken分类层级选择的细节Bracken默认支持从域到种多个层级实际分析中物种水平-l S用得最多属水平-l G次之。做物种水平分析时要注意参考数据库里有些条目标注不完整比如某些菌株在NCBI的分类学树里只定位到属这时Bracken无法把reads校正到种只能归到上一级结果里就会出现“Genus”和“Species”混着的情况。处理这种结果时可以把不同分类层级分开看或者用物种注释工具提供的分类学过滤器过滤过低支持度的分类。6.4 踩过几次坑之后的几点心得首先数据库下载别图快官方标准库最好下完整包只用一部分参考序列建库看起来省了时间实际跑出来的结果覆盖度差很多。其次跑正式数据前一定先用小规模测试数据验证流程。最后做好数据库版本记录。同一项目里用不同版本的数据库结果不一致是一个很尴尬的问题尤其投稿返修时审稿人问起物种注释版本如果拿不出记录会非常被动。另外分享一个日常小技巧写分析流程时把Kraken2和Bracken的命令参数写进配置文件里用Snakemake或Nextflow管理的话每个样本的输入输出参数自动记录这样不仅完整复现方便出了问题排查日志也快。Kraken2Bracken的高效和易用让物种注释从一个耗时的瓶颈环节变成了分析流程里最顺畅的一步。只要数据库选好、参数配好从原始数据到物种丰度表半天就能跑完一轮。实际使用中如果遇到数据库下载或内存不足这类问题先按上面表格排查大多数情况下都能解决。后面如果再深入一些可以把Kraken2的结果和金标准比对工具做交叉验证或者把注释结果接入功能预测流程这些都是已经比较成熟的方向了。