摘要: karioCaS 是一个基于 R 语言开发的宏基因组数据分析工具包，旨在提升 Kraken2（及其他基于 k-mer 分类器）分类结果的可靠性。其核心目标是通过系统性地优化置信度评分（Confidence Score，CS）参数，在统计噪声过滤与生物信号保留之间取得平衡，从而有效抑制假阳性分类结果的过度膨胀。

该工具包的方法论建立在两个核心原则之上：其一为领域特异性分析，即针对古菌、细菌、真核生物和病毒四大生物域分别独立执行分析与可视化，摒弃"一刀切"的传统做法；其二为置信度评分探索，通过数据驱动的可视化手段，帮助用户识别每个生物域所对应的最优置信度阈值，精确定位区分统计噪声与真实生物信号（即"生物暗物质"）的数学拐点。

在实现层面，karioCaS 要求用户采用严格的目录结构和文件命名规范，支持多种置信度评分表示格式（百分比、遗留简写及显式小数），并采用多策略模型框架进行综合评估。该包目前作为 R 开发版本提供，并正准备提交至 Bioconductor 平台，适用于需要高可靠性微生物群落组成推断的科研项目。

karioCaS：基于置信度分数的宏基因组域特异性微生物群分类优化工具¶

概述¶

宏基因组学研究中，基于 k-mer 的分类工具（如 Kraken2）在高灵敏度设置下极易产生大量假阳性分类结果，造成统计噪声与真实生物信号的混淆。如何在保留真实生物信号（即"生物暗物质"，Biological Dark Matter）的同时有效过滤统计噪声，是宏基因组数据分析中的核心挑战。

karioCaS 是一个专为解决上述问题而设计的 R 包。其全称为 Kraken Confidence Scores for Reliable Domain-Specific Microbiota Inference and Discovery，核心思路是充分利用 Kraken2 中一个常被忽视的关键参数——置信度分数（Confidence Score，CS）——通过跨多个严格度水平（stringency levels）评估分类结果的稳定性，为每个生物域独立确定最优 CS 阈值。

该工具的两大设计原则为：域特异性分析（Domain-Specific Analysis）与置信度分数探索（Confidence Score Exploration）。前者承认不同生物域（古菌 Archaea、细菌 Bacteria、真核生物 Eukaryota、病毒 Viruses）具有各异的分类特性，拒绝"一刀切"策略；后者则通过数据驱动的可视化和多策略模型，精确定位噪声与信号之间的数学拐点。

karioCaS 输出高置信度的精炼 .mpa 文件，可直接用于下游分析，显著提升宏基因组研究的可重复性与生物学可靠性。

核心原理与功能¶

置信度分数（Confidence Score）的作用机制¶

Kraken2 的置信度分数（CS）决定了将一条 reads 分配给某分类单元所需的最低 k-mer 匹配比例：

• CS = 0.0：只要有单个 k-mer 匹配即可完成分类（最高灵敏度，噪声最大）
• CS = 0.1：至少 10% 的 k-mer 必须与某基因组完全匹配才能完成分配
• CS = 1.0：所有 k-mer 必须一致（最高严格度，噪声最小）

karioCaS 的核心假设是：随着 CS 逐步提升，真实的生物信号（稳定存在的分类单元）会在某个阈值之后趋于稳定，而统计噪声（假阳性分类）则会在低 CS 段快速消失。通过捕捉这一"稳定拐点"，可以客观确定每个域的最优 CS。

域特异性分析框架¶

所有分析和可视化均独立在四个生物域上执行：

• Archaea（古菌）
• Bacteria（细菌）
• Eukaryota（真核生物）
• Viruses（病毒）

这一设计避免了跨域 CS 阈值混用导致的偏差，体现了"one size does not fit all"（不同域需要不同策略）的原则。

工作流程与核心功能模块¶

第 000 步：数据协调（Data Harmonization）¶

import_karioCaS()：将来自不同 CS 严格度设置下的原始 Metaphlan 风格输出文件标准化，整合为统一的 Bioconductor TreeSummarizedExperiment（TSE）对象，作为所有后续分析的数据基础。

第 001–005 步：可视化探索与客观阈值确定（Visual Exploration & Objective Thresholding）¶

taxa_retention()：量化随严格度增加各域保留分类单元的百分比，同时计算每个域的最优 CS（稳定性指数，Stability Index）。其组图（group plot）以虚线标记每个域的中位最优 CS，并生成供第 1000 步使用的 SI_Audit_<rank> 表格。

支持以下四种阈值策略（通过 method= 参数指定）：

taxa_resolution()：评估分类深度，即亲代到子代（Parent-to-Child）的分类解析度。默认分析 retrieve_selected_taxa() 的最终拼图输出，每样本生成一张图；传入 CS= 参数则可分析导入数据在单一 CS 下的情况。

reads_per_taxa()：在对数 reads 轴上进行饱和度分析（saturation analysis），通过与最优 CS 相同的算法引擎计算每个域的最优最低 reads 数（饱和曲线肘点），在组图上标记每个域的中位值，并生成 Reads_Audit_<rank> 表格，为排除低丰度背景/假阳性分类单元提供定量阈值。

upset_kariocas()：识别"瞬时"（transient）与"持续"（persistent）分类单元，揭示哪些分类单元在不同 CS 水平下始终存在（真实信号），哪些仅在低 CS 时出现（噪声）。

heatmaps_karioCaS()：绘制详细的丰度热图，展示随 CS 提升而发生的分类单元消亡（extinction）模式。

分组叠加机制（Group Overlays）¶

taxa_retention() 和 reads_per_taxa() 默认按生物学分组（biological group）绘图，而非逐样本输出：每个样本以淡线显示，组均值（± SD）高亮展示，并按域分面（facet）。分组规则为去除样本名末尾数字，例如：

• SAMPLE33、SAMPLE34 → 分组 SAMPLE
• CONTROL01、TREATED01 → 分组 CONTROL、TREATED

通过 detail_samples= 参数控制详细样本输出：

• NULL（默认）：仅输出分组图
• "all"：为每个样本生成详细面板（输出至 per_sample/ 子文件夹）
• "SAMPLE33, SAMPLE45"：仅为指定样本生成详细面板

第 1000 步：最终生物拼图（The Ultimate Biological Mosaic）¶

retrieve_selected_taxa()：综合使用两个数据驱动阈值——来自 taxa_retention() 的最优 CS（SI_Audit）和来自 reads_per_taxa() 的最优最低 reads 数（Reads_Audit）——提取每个域的存活分类单元，生成经过统计噪声清洗的高置信度 .mpa 文件，可直接用于下游分析。

CS_* 和 reads_min_* 参数各接受以下值：

• "auto"：自动从审计表中读取主策略阈值
• "secondary"：读取次策略阈值
• 手动数值：用户自定义覆盖

安装与使用¶

安装¶

karioCaS 可直接从 GitHub 安装，当前正在准备提交至 Bioconductor。

# 安装 devtools（如已安装可跳过）
install.packages("devtools")

# 从 GitHub 安装 karioCaS 开发版
devtools::install_github("thiagoparentefiocruz/karioCaS")

目录结构要求（严格遵守）¶

karioCaS 要求项目目录符合特定结构。原始 Metaphlan 风格输出文件必须放置于基础目录下的 000_mpa_original 文件夹中：

YOUR_PROJECT_DIR/
 └── 000_mpa_original/
      ├── SAMPLE01_CS00.mpa       # 样本01，CS = 0.0（无严格度过滤）
      ├── SAMPLE01_CS05.mpa       # 样本01，CS = 0.5（中等严格度）
      └── SAMPLE01_CS10.mpa       # 样本01，CS = 1.0（最高严格度）

文件命名规则¶

文件名格式为 {SAMPLE}_CS{XX}.mpa，各部分说明如下：

1. SAMPLE：样本名称，可任意命名，但样本名本身不得包含下划线 _
2. _CS：固定字段，必须存在
3. XX：置信度分数，支持三种等价写法：

函数参数中的 CS 值格式¶

confidence_score= 和 CS_B= 等参数接受以下两种格式：

• Kraken 分数形式，范围 (0, 1]，例如 1.0 表示最高严格度
• 百分比形式，值 > 1，例如 40 等同于 CS = 0.4

实战示例¶

以下为完整的 karioCaS 工作流示例，涵盖从数据导入到最终拼图输出的全部步骤：

library(karioCaS)  # 加载 karioCaS 包

# 设置项目根目录路径
proj_dir <- "path/to/your_project"

# ==============================
# 第 1 步：数据协调（Step 000）
# ==============================
# 将不同 CS 下的原始 .mpa 文件整合为统一的 TSE 对象
import_karioCaS(project_dir = proj_dir)

# ==============================
# 第 2 步：可视化探索与阈值确定（Steps 001-005）
# ==============================

# 绘制分组叠加的分类单元保留率曲线，
# 同时计算每个域的稳定性指数（Stability Index），
# 默认使用 Kneedle 肘点检测方法
taxa_retention(project_dir = proj_dir)

# 针对特定样本生成详细 PDF（输出至 per_sample/ 子文件夹）
taxa_retention(project_dir = proj_dir, detail_samples = "SAMPLE33, SAMPLE45")

# 评估分类单元在 CS 升高过程中的"持续性"（UpSet 图）
upset_kariocas(project_dir = proj_dir)

# 分析 NGS Reads 保留率（分组叠加图；
# 使用 detail_samples = "all" 可为每个样本生成详细面板）
reads_per_taxa(project_dir = proj_dir)

# 在指定置信度分数（CS = 0.4）下评估亲代-子代分类解析度
taxa_resolution(project_dir = proj_dir, CS = 40)

# 绘制分类单元消亡模式的丰度热图
heatmaps_karioCaS(project_dir = proj_dir)

# ==============================
# 第 3 步：提取最终生物拼图（Step 1000）
# ==============================
# CS_B / reads_min_B：细菌（Bacteria）的 CS 和最低 reads 阈值，均使用自动推断
# CS_A / reads_min_A：古菌（Archaea）的 CS 和最低 reads 阈值，均使用自动推断
# CS_E = 40：真核生物（Eukaryota）手动指定 CS = 0.4
# reads_min_E = 10：真核生物手动指定最低 reads 数 = 10
# CS_V = 0 / reads_min_V = 0：病毒（Viruses）不施加任何过滤
retrieve_selected_taxa(
  project_dir = proj_dir,
  tax_level = "Species",       # 分类层级：物种水平
  CS_B = "auto",   reads_min_B = "auto",   # 细菌：完全数据驱动
  CS_A = "auto",   reads_min_A = "auto",   # 古菌：完全数据驱动
  CS_E = 40,       reads_min_E = 10,       # 真核生物：手动覆盖
  CS_V = 0,        reads_min_V = 0         # 病毒：无过滤
)

# ==============================
# 第 4 步：评估最终拼图的分类解析度
# ==============================
# 默认分析 retrieve_selected_taxa() 的输出（最终拼图）
# 若需保留亲代分类层级，可在上一步中设置 tax_level = NULL
taxa_resolution(project_dir = proj_dir)

上述示例展示了一个典型的混合策略场景：对细菌和古菌采用完全数据驱动的自动阈值，对真核生物采用专家经验手动覆盖（CS = 0.4，最低 10 条 reads），对病毒则不施加任何过滤限制（适用于低丰度病毒研究场景）。

常见问题¶

Q1：为什么分析必须按生物域分开进行，而不能对所有分类单元统一处理？

不同生物域的基因组复杂度、序列相似性分布和数据库覆盖度差异显著。例如，病毒基因组高度多样且数据库不完整，适合较低的 CS；而细菌基因组数据库完善，较高 CS 仍能保留大量真实信号。使用统一 CS 阈值会导致某些域过度过滤，另一些域噪声残留，因此域特异性分析是更合理的策略。

Q2：文件命名中 CS05 到底代表 0.5 还是 0.05？

在旧式 0.1 步长缩写模式下，CS05 解析为 0.5，而非 0.05。若需精确表示 0.05，应使用显式小数写法 CS0.05。建议在实际使用中优先采用百分比写法（如 CS50 代表 0.5）或显式小数写法（如 CS0.50）以避免歧义。

Q3：retrieve_selected_taxa() 中 "auto" 与 "secondary" 有何区别？

• "auto"：从 SI_Audit（由 taxa_retention() 生成）或 Reads_Audit（由 reads_per_taxa() 生成）中读取主策略（默认为 Kneedle）计算出的阈值。
• "secondary"：读取次策略阈值，适用于主策略结果偏于保守或激进、需要备选参考的情况。

Q4：样本名中为什么不能使用下划线 _？

文件命名规则使用 _CS 作为固定分隔符来解析样本名和置信度分数。若样本名本身包含下划线，解析器将无法正确区分样本名与 _CS 分隔符，导致文件读取错误。应使用连字符 - 或驼峰命名法等替代方案。

Q5：taxa_retention() 的四种策略（Kneedle、Post-cliff、Segmented、Dynamic）如何选择？

• 若无特殊先验知识，优先使用默认的 Kneedle 方法，无需参数调整，跨域鲁棒性强。
• 若研究目标是最大程度过滤噪声、偏保守，选择 Post-cliff（拐点更靠后的高原区）。
• 若数据包含明显的机制转变（如生态梯度研究、暗物质探索），选择 Segmented（断棍回归）。
• 若有明确的专家经验或特定损耗容忍度要求，选择 Manual 并手动指定阈值。

总结¶

karioCaS 提供了一套完整的、数据驱动的宏基因组分类可靠性增强框架。其核心价值在于：将 Kraken2 置信度分数（CS）从一个常被忽视的参数转化为客观的噪声过滤工具；通过域特异性分析避免了跨域阈值混用的偏差；通过多策略拐点检测（Kneedle、Post-cliff、Segmented、Dynamic）灵活适应不同研究场景。

工作流从数据协调（import_karioCaS()）出发，经过多维度可视化探索（分类单元保留率、持续性、reads 饱和度、分类深度、丰度热图），最终通过 retrieve_selected_taxa() 整合双重数据驱动阈值，输出去噪后的高置信度 .mpa 文件。该文件同时支持全自动推断与手动覆盖，兼顾分析灵活性与可重复性，是宏基因组数据质量控制流程中的有力工具。