---
title: "karioCaS: Kraken Confidence Scores Exploration"
author:
- name: Thiago Parente
  affiliation: Instituto Oswaldo Cruz, Fiocruz, Rio de Janeiro, Brazil
  email: thiago.parente@fiocruz.br
package: karioCaS
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{karioCaS: Kraken Confidence Scores Exploration}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup_knitr, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    warning = FALSE,
    message = FALSE,
    fig.width = 10,
    fig.height = 7
)
```

# Introduction

Metagenomic classification with k-mer based tools such as Kraken2 is fast, but prone to false-positive assignments. Kraken2's `--confidence` parameter (the Confidence Score, CS) controls how much k-mer evidence a read needs before it is assigned to a taxon: higher values remove spurious assignments, but also discard true signal. The usual practice of applying a single, global CS to a whole sample ignores that biological domains (Bacteria, Archaea, Eukaryota and Viruses) differ widely in database representation, genome size and k-mer uniqueness, so **one threshold rarely fits all**.

**karioCaS** (Kraken Confidence Scores for Reliable Domain-Specific Microbiota Inference and Discovery) takes the same sample classified at several CS values and evaluates how taxa and reads are lost as stringency increases, separately for each domain. From these decay curves it derives objective, domain-specific thresholds for the CS and for the minimum number of reads, and assembles a final "mosaic" profile in which each domain is filtered at its own optimal stringency.

## karioCaS and Bioconductor

karioCaS works upstream of most microbiome analysis. The reports produced at the different CS values are imported into a single `r BiocStyle::Biocpkg("TreeSummarizedExperiment")` object (taxa as rows, one column per sample and CS), so the data can be inspected or passed on with standard `r BiocStyle::Biocpkg("SummarizedExperiment")` accessors and with the wider ecosystem of packages built on these classes, such as `r BiocStyle::Biocpkg("mia")`. Existing Bioconductor packages for microbiome data start from an already chosen taxonomic profile; to our knowledge, none of them addresses how the Kraken2 confidence score used to produce that profile affects the result. karioCaS fills this gap and hands a curated, high-confidence profile to those downstream tools.

# Installation

karioCaS can be installed from Bioconductor with:

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("karioCaS")
```

```{r setup}
library(karioCaS)
```

# Input data and project folder

karioCaS works on a project folder. Place one Kraken2 report per sample and CS, converted to MPA format (for example with KrakenTools' `kreport2mpa.py`), inside a sub-folder named `000_mpa_original`. File names must follow `SAMPLENAME_CSXX.mpa`, where `XX` encodes the CS: `SAMPLE01_CS09.mpa` is sample `SAMPLE01` classified with `--confidence 0.9` (`CS00` to `CS10` are read as tenths; `CS40`, `CS100` or `CS0.05` are also accepted). Include a report at CS 0 for every sample: it is the reference for the detailed per-sample retention panels.

Each analysis step writes its results to its own numbered sub-folder of the project. This vignette uses the example data shipped with the package, a single sample classified at six CS values (0, 0.2, 0.4, 0.6, 0.8 and 0.9). How these files were produced is described in `inst/script/` of the package source. We copy them to a temporary folder, because the installed package folder must not be written to.

```{r prep_data}
proj_dir <- file.path(tempdir(), "karioCaS_Tutorial")
dir.create(file.path(proj_dir, "000_mpa_original"), recursive = TRUE)
mock_files <- list.files(
    system.file("extdata", "your_project_name", "000_mpa_original",
        package = "karioCaS"
    ),
    full.names = TRUE
)
file.copy(mock_files, file.path(proj_dir, "000_mpa_original"))
basename(mock_files)
```

# Importing the reports

`import_karioCaS()` reads every report in `000_mpa_original` and assembles them into a `TreeSummarizedExperiment`:

* each report becomes one column, identified as `SAMPLE_CSxx`, with the sample name and the CS (as an integer percentage, e.g. `CS09` becomes 90) stored in `colData()`;
* each distinct taxonomic path becomes one row; a taxon absent from a report gets a count of 0 in that column, and identical paths within one report are summed;
* the lineage string is split into eight rank columns (Domain, Kingdom, Phylum, Class, Order, Family, Genus, Species) in `rowData()`, plus `Rank`, the lowest rank present in the path.

No normalisation or transformation is applied: the assay holds the read counts exactly as in the reports. Note that MPA counts are *cumulative*: the count of a genus row includes the reads assigned to its species.

```{r import_data}
tse_object <- import_karioCaS(project_dir = proj_dir)
tse_object
head(SummarizedExperiment::colData(tse_object))
```

The object is also saved, together with a tab-separated copy of the matrix, in `001_imported_matrix/`, where the following steps read it from.

# Analysis results: the `kariocas_result` object

Every analysis step returns a `kariocas_result` object:

* `$data`: the computed table;
* `$plots`: the figures, as `ggplot`/`patchwork` objects (or UpSetR plots) that can be printed or modified;
* `$paths`: the files written to the project folder.

By default (`export = TRUE`) results are also written to the step's sub-folder as PDF, TSV and RDS files, with a log in `logs/`. With `export = FALSE` nothing is written, which is convenient for interactive exploration, for example when comparing `method` or `tax_level` values. The final mosaic (see *The final mosaic* below) reads the threshold tables *exported* by the previous steps, so run those with `export = TRUE` before building it.

# Objective thresholds

## Confidence Score: the Stability Index

`taxa_retention()` counts, for each sample and domain, how many distinct taxa at `tax_level` are detected (at least one read) at each CS, and expresses it as a percentage of the maximum across CS values, usually the count at CS 0. This decay curve drops steeply while spurious assignments are removed and flattens once only well-supported taxa remain. The **Stability Index (SI)** engine locates that transition and reports it as the *Primary SI*, the optimal CS for that domain. Several strategies are available through `method`:

* **Kneedle** (default): draws a straight line between the first and last points of the curve and picks the CS where the curve lies furthest below it, i.e. the elbow. Parameter-free.
* **Post-cliff**: the first CS after the steepest drop whose step loss falls below a tolerance estimated from the tail of the curve (CS >= 50); more conservative. With `kneedle`, this value is also reported as a *Secondary SI* when it is stricter than the elbow.
* **Segmented**: fits two straight lines (broken-stick regression) and uses the CS at the best breakpoint.
* **Dynamic** / **Manual**: the first CS whose step loss is within the tail-derived tolerance, or within a loss (%) given by the user.

```{r optimize_cs}
retention <- taxa_retention(project_dir = proj_dir, tax_level = "Species")
retention

audit_df <- retention$data
audit_df[
    audit_df$SI_Type %in% "Primary_SI",
    c("Sample", "Domain", "CS", "Taxa_Count", "Pct_Retained", "Step_Loss_Pct")
]
```

`Taxa_Count` is the number of taxa detected at that CS, `Pct_Retained` its percentage of the maximum, and `Step_Loss_Pct` the percentage points lost since the previous CS. The figure overlays the curves of all samples of a biological group (groups are inferred from sample names by removing trailing digits, e.g. `SAMPLE33` and `SAMPLE34` form group `SAMPLE`); the bold line is the group mean with a +/- SD band, and the dashed line is the median Primary SI of each domain.

```{r retention_plot}
retention$plots[[1]]
```

Detailed per-sample panels (all ranks, taxa versus reads) are produced on request with `detail_samples = "all"` or a list of sample names.

A domain needs a decay curve with at least three CS values to be assessed. In the example data, viral species are detected only at CS 0 (16 species) and CS 0.2 (one species), so no Stability Index is reported for Viruses. If `CS_V = "auto"` were requested in the next step, CS 0 would be used for Viruses and a warning would say so; the example below therefore sets the viral thresholds manually.

## Minimum reads: saturation analysis

Low-count taxa are a second source of false positives. `reads_per_taxa()` raises a minimum-read cutoff step by step (1, 2, 3, 4, 5, 7, 10, then log-spaced values up to the largest count) and counts how many taxa at `analysis_level` still have at least that many reads, at each CS. The elbow of this saturation curve, located with the same engine on a log10 read axis, is the optimal minimum number of reads for the domain at that CS.

```{r reads_audit}
reads <- reads_per_taxa(project_dir = proj_dir, analysis_level = "Species")
head(reads$data[reads$data$SI_Type %in% "Primary_SI", ])
```

# The final mosaic

`retrieve_selected_taxa()` builds, for each sample, a profile in which each domain is taken from the report at its own CS and filtered by its own minimum reads. Each threshold can be `"auto"` (the Primary SI or optimal reads from the exported audits), `"secondary"` (the more conservative value, when available) or a number. The mosaic keeps all taxonomic ranks, like the original reports.

```{r retrieve_taxa}
mosaic <- retrieve_selected_taxa(
    project_dir = proj_dir,
    CS_B = "auto", reads_min_B = "auto", # optimal CS + optimal reads
    CS_A = "secondary", reads_min_A = "auto", # conservative CS, optimal reads
    CS_E = 40, reads_min_E = 10, # manual values for Eukaryota
    CS_V = 0, reads_min_V = 0 # no filtering for Viruses
)
mosaic

# Thresholds applied to each domain and where they came from
unique(mosaic$data[, c(
    "Domain", "CS", "CS_Source", "Min_Reads", "Min_Reads_Source"
)])
```

The mosaic is written to `004_final_mosaic/` in MPA (`mpa/`) and tab-separated (`tsv/`) formats, ready for downstream analysis.

# Visual exploration

The following functions help to inspect how individual taxa behave across CS values. Here they are run with `export = FALSE`, so the figures are only returned.

## Taxa persistence across CS

`upset_kariocas()` shows, for each sample and domain, which taxa are detected at which CS values. Taxa detected at every CS form the stable core, and taxa present only at low CS are the ones stringency removes.

```{r upset}
ups <- upset_kariocas(proj_dir, tax_level = "Genus", export = FALSE)
ups$plots[["SAMPLE01_Bacteria_Genus"]]
```

## Relative abundance heatmaps

`heatmaps_karioCaS()` shows the relative abundance (%) of the most abundant taxa surviving at a target CS, across all CS values up to that target. The remaining taxa are aggregated into rows according to the CS at which they were last detected.

```{r heatmap, fig.height = 10}
hm <- heatmaps_karioCaS(
    proj_dir,
    analysis_rank = "Species", confidence_score = 40, top_n = 10,
    export = FALSE
)
hm$plots[[1]]
```

## Taxonomic resolution

`taxa_resolution()` compares, for the most abundant parent taxa (e.g. genera), the reads resolved to the child rank (e.g. species) with those that stay at the parent level. By default it analyses the final mosaic; here it is run on the imported data at CS 40.

```{r resolution}
res <- taxa_resolution(proj_dir, CS = 40, export = FALSE)
head(res$data)
res$plots[[1]]
```

## Core and unique taxa across samples

When a study has several samples per biological group, `group_upset()` compares the taxa detected across the samples of each group and classifies them as *Core* (present in every sample), *Shared* or *Unique* (present in one sample only), a useful pattern to spot sample-specific signals and false positives. It needs at least two samples per group, which the example data does not have, so the call is shown without being run:

```{r group_upset, eval = FALSE}
group_upset(project_dir = proj_dir)
```

```{r cleanup, include = FALSE}
unlink(proj_dir, recursive = TRUE)
```

# Acknowledgements

The scientific concept, the method design and its validation are the author's. The package code was written with the assistance of large language models, which also contributed to the formulation of some of the calculations. All code was reviewed and tested by the author.

# Session information

```{r sessionInfo, echo = FALSE}
sessionInfo()
```