karioCaS: Kraken Confidence Scores Exploration
14 October 2026
Package
karioCaS 0.99.16
1 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.
1.1 karioCaS and Bioconductor
karioCaS works upstream of most microbiome analysis. The reports produced at the different CS values are imported into a single TreeSummarizedExperiment object (taxa as rows, one column per sample and CS), so the data can be inspected or passed on with standard SummarizedExperiment accessors and with the wider ecosystem of packages built on these classes, such as 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.
2 Installation
karioCaS can be installed from Bioconductor with:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("karioCaS")library(karioCaS)3 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.
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"))
#> [1] TRUE TRUE TRUE TRUE TRUE TRUE
basename(mock_files)
#> [1] "SAMPLE01_CS00.mpa" "SAMPLE01_CS02.mpa" "SAMPLE01_CS04.mpa"
#> [4] "SAMPLE01_CS06.mpa" "SAMPLE01_CS08.mpa" "SAMPLE01_CS09.mpa"4 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.CS09becomes 90) stored incolData(); - 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(), plusRank, 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.
tse_object <- import_karioCaS(project_dir = proj_dir)
tse_object
#> class: TreeSummarizedExperiment
#> dim: 2030 6
#> metadata(0):
#> assays(1): counts
#> rownames(2030): d__Archaea d__Archaea|k__Methanobacteriati ...
#> d__Bacteria|p__Pseudomonadota|c__Betaproteobacteria|o__Burkholderiales|f__Comamonadaceae|g__Hydrogenophaga|s__Hydrogenophaga
#> crassostreae
#> d__Bacteria|p__Pseudomonadota|c__Alphaproteobacteria|o__Sphingomonadales|f__Sphingomonadaceae|g__Novosphingobium|s__Novosphingobium
#> sp. KCTC 2891
#> rowData names(10): Taxonomy_Full Domain ... Species Rank
#> colnames(6): SAMPLE01_CS00 SAMPLE01_CS20 ... SAMPLE01_CS80
#> SAMPLE01_CS90
#> colData names(2): Sample_ID Confidence_Score
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> rowLinks: NULL
#> rowTree: NULL
#> colLinks: NULL
#> colTree: NULL
head(SummarizedExperiment::colData(tse_object))
#> DataFrame with 6 rows and 2 columns
#> Sample_ID Confidence_Score
#> <character> <integer>
#> SAMPLE01_CS00 SAMPLE01 0
#> SAMPLE01_CS20 SAMPLE01 20
#> SAMPLE01_CS40 SAMPLE01 40
#> SAMPLE01_CS60 SAMPLE01 60
#> SAMPLE01_CS80 SAMPLE01 80
#> SAMPLE01_CS90 SAMPLE01 90The object is also saved, together with a tab-separated copy of the matrix, in 001_imported_matrix/, where the following steps read it from.
5 Analysis results: the kariocas_result object
Every analysis step returns a kariocas_result object:
$data: the computed table;$plots: the figures, asggplot/patchworkobjects (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.
6 Objective thresholds
6.1 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.
retention <- taxa_retention(project_dir = proj_dir, tax_level = "Species")
retention
#> <kariocas_result> step: 002_taxa_retention
#> data : 15 row(s) x 7 column(s)
#> plots : 1 [SAMPLE_Group_Retention]
#> export: 3 file(s) in /tmp/RtmpvXuoop/karioCaS_Tutorial/002_taxa_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")
]
#> # A tibble: 3 × 6
#> Sample Domain CS Taxa_Count Pct_Retained Step_Loss_Pct
#> <chr> <chr> <int> <int> <dbl> <dbl>
#> 1 SAMPLE01 Bacteria 60 48 70.6 22.1
#> 2 SAMPLE01 Archaea 20 3 5 95
#> 3 SAMPLE01 Eukaryota 40 7 31.8 22.7Taxa_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.
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.
6.2 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.
reads <- reads_per_taxa(project_dir = proj_dir, analysis_level = "Species")
head(reads$data[reads$data$SI_Type %in% "Primary_SI", ])
#> Sample CS Domain Cutoff Taxa_Count Pct_Retained Step_Loss_Pct
#> 25 SAMPLE01 0 Bacteria 50000 6 8.823529 4.411765
#> 40 SAMPLE01 0 Archaea 70 2 3.333333 11.666667
#> 48 SAMPLE01 0 Viruses 4 8 50.000000 12.500000
#> 65 SAMPLE01 0 Eukaryota 100 6 27.272727 13.636364
#> 96 SAMPLE01 20 Bacteria 3000 10 14.705882 5.882353
#> 105 SAMPLE01 20 Archaea 2 2 66.666667 33.333333
#> SI_Type
#> 25 Primary_SI
#> 40 Primary_SI
#> 48 Primary_SI
#> 65 Primary_SI
#> 96 Primary_SI
#> 105 Primary_SI7 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.
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
#> <kariocas_result> step: 004_final_mosaic
#> data : 238 row(s) x 8 column(s)
#> plots : 0
#> export: 2 file(s) in /tmp/RtmpvXuoop/karioCaS_Tutorial/004_final_mosaic
# Thresholds applied to each domain and where they came from
unique(mosaic$data[, c(
"Domain", "CS", "CS_Source", "Min_Reads", "Min_Reads_Source"
)])
#> Domain CS CS_Source Min_Reads Min_Reads_Source
#> d__Archaea Archaea 40 [SI: Secondary] 4 [Reads:auto]
#> d__Bacteria Bacteria 60 [SI: Primary] 2000 [Reads:auto]
#> d__Eukaryota Eukaryota 40 [Manual] 10 [Manual]
#> d__Viruses Viruses 0 [Manual] 0 [Manual]The mosaic is written to 004_final_mosaic/ in MPA (mpa/) and tab-separated (tsv/) formats, ready for downstream analysis.
8 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.
8.1 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.
ups <- upset_kariocas(proj_dir, tax_level = "Genus", export = FALSE)
ups$plots[["SAMPLE01_Bacteria_Genus"]]
8.2 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.
hm <- heatmaps_karioCaS(
proj_dir,
analysis_rank = "Species", confidence_score = 40, top_n = 10,
export = FALSE
)
hm$plots[[1]]
8.3 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.
res <- taxa_resolution(proj_dir, CS = 40, export = FALSE)
head(res$data)
#> # A tibble: 6 × 7
#> sample Source Parent_Name Parent_Cumulative_Total Child_Sum_Resolved
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 SAMPLE01 CS40 Novosphingobium 99135 15271
#> 2 SAMPLE01 CS40 Aquirufa 54003 11768
#> 3 SAMPLE01 CS40 Ideonella 46194 44620
#> 4 SAMPLE01 CS40 Rhodoferax 14924 4035
#> 5 SAMPLE01 CS40 Hydrogenophaga 13678 2873
#> 6 SAMPLE01 CS40 Croceicoccus 168 167
#> # ℹ 2 more variables: Parent_Exclusive <dbl>, Domain <chr>
res$plots[[1]]
8.4 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:
group_upset(project_dir = proj_dir)9 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.
10 Session information
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] karioCaS_0.99.16 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 dplyr_1.2.1
#> [3] farver_2.1.2 Biostrings_2.81.9
#> [5] S7_0.2.2 fastmap_1.2.0
#> [7] SingleCellExperiment_1.35.2 lazyeval_0.2.3
#> [9] digest_0.6.39 lifecycle_1.0.5
#> [11] tidytree_0.4.8 magrittr_2.0.5
#> [13] compiler_4.6.1 rlang_1.3.0
#> [15] sass_0.4.10 tools_4.6.1
#> [17] utf8_1.2.6 yaml_2.3.12
#> [19] knitr_1.52 labeling_0.4.3
#> [21] S4Arrays_1.13.2 bit_4.6.0
#> [23] DelayedArray_0.39.8 plyr_1.8.9
#> [25] xml2_1.6.0 RColorBrewer_1.1-3
#> [27] TreeSummarizedExperiment_2.21.0 abind_1.4-8
#> [29] BiocParallel_1.47.0 withr_3.0.3
#> [31] purrr_1.2.2 BiocGenerics_0.59.12
#> [33] grid_4.6.1 stats4_4.6.1
#> [35] ggplot2_4.0.3 scales_1.4.0
#> [37] tinytex_0.61 dichromat_2.0-1
#> [39] SummarizedExperiment_1.43.0 cli_3.6.6
#> [41] UpSetR_1.4.1 rmarkdown_2.32
#> [43] crayon_1.5.3 ragg_1.5.2
#> [45] treeio_1.37.1 generics_0.1.4
#> [47] otel_0.2.0 tzdb_0.5.0
#> [49] commonmark_2.0.0 ape_5.8-1
#> [51] cachem_1.1.0 stringr_1.6.0
#> [53] parallel_4.6.1 BiocManager_1.30.27
#> [55] XVector_0.53.0 matrixStats_1.5.0
#> [57] yulab.utils_0.2.5 vctrs_0.7.3
#> [59] Matrix_1.7-6 jsonlite_2.0.0
#> [61] litedown_0.11 bookdown_0.48
#> [63] IRanges_2.47.5 hms_1.1.4
#> [65] patchwork_1.3.2 S4Vectors_0.51.10
#> [67] bit64_4.8.6 magick_2.9.1
#> [69] systemfonts_1.3.2 tidyr_1.3.2
#> [71] jquerylib_0.1.4 glue_1.8.1
#> [73] codetools_0.2-20 ggtext_0.2.0
#> [75] stringi_1.8.9 gtable_0.3.6
#> [77] GenomicRanges_1.65.4 tibble_3.3.1
#> [79] pillar_1.11.1 rappdirs_0.3.4
#> [81] htmltools_0.5.9 Seqinfo_1.3.2
#> [83] R6_2.6.1 textshaping_1.0.5
#> [85] vroom_1.7.1 evaluate_1.0.5
#> [87] lattice_0.23-1 Biobase_2.73.2
#> [89] markdown_2.0 readr_2.2.0
#> [91] gridtext_0.1.6 bslib_0.12.0
#> [93] Rcpp_1.1.2 gridExtra_2.3.1
#> [95] SparseArray_1.13.4 nlme_3.1-171
#> [97] xfun_0.61 fs_2.1.0
#> [99] MatrixGenerics_1.25.0 forcats_1.0.1
#> [101] pkgconfig_2.0.3