kamino_cli/

lib.rs

//! # kamino
//!
//! kamino builds an amino-acid alignment in a reference-free, alignment-free manner from
//! a set of proteomes. It is not “better” than traditional marker-based pipelines, but it is
//! simpler and faster to use.
//!
//! Typical usage ranges from between-species to within-phylum phylogenetic analyses (bacteria,
//! archaea, and eukaryotes).
//!
//! ## Input modes
//! kamino accepts proteome files as input in one of two modes:
//! - **Directory mode** (`--input-directory`): a directory containing FASTA proteomes
//!   (plain text or `.gz` compressed). Each file represents one isolate. Filenames minus the
//!   extension become sequence names in the final amino-acid alignment.
//! - **Table mode** (`--input-file`): a tab-delimited file mapping a species/sample name
//!   to a proteome path (one name + path pair per line). This is useful when file names
//!   do not encode the sample name or when proteomes are located in multiple directories.
//!
//! In the directory mode, files are recognized by their extension (.fas, .fasta, .faa, .fa, .fna; gzipped ot not).
//!
//! For **bacterial** isolates, the phylogenomic alignment can also be generated directly from genome assemblies 
//! by selecting the option `--genomes` (using either `-i` or `-I`). In this case, an ultra-fast but approximate 
//! protein prediction is performed, and the predicted proteomes are written to a temporary directory.
//!
//! ## Arguments
//! - `-i`, --input-directory <INPUT>: input directory with FASTA proteomes (plain or .gz)
//! - `-I`, --input-file <INPUT_FILE>: tab-delimited file mapping species name to proteome path
//! - `--genomes`: treat input files as bacterial genomes and predict proteomes before analysis
//! - `-k`, `--k`: k-mer length (default: 14; must be within the valid range for the
//!   selected recoding scheme).
//! - `-f`, `--min-freq`: minimum fraction of samples with an amino-acid per position
//!   (default: 0.85; must be ≥ 0.6).
//! - `-d`, `--depth`: maximum traversal depth from each start node (default: 12).
//! - `-o`, `--output`: output prefix for generated files (default: `kamino`).
//! - `-c`, `--constant`: number of constant positions retained from in-bubble k-mers
//!   (default: 3; must be ≤ k-1).
//! - `-l`, `--length-middle`: maximum number of middle positions per variant group
//!   (default: 35; must be ≥ 1).
//! - `-m`, `--mask`: mask middle segments with long mismatch runs (default: 5).
//! - `-t`, `--threads`: number of threads used for graph construction and analysis
//!   (default: 1).
//! - `-r`, `--recode`: amino-acid recoding scheme (default: `sr6`).
//! - `--nj`: generate a NJ tree from kamino alignment [nj=false]
//! - `-v`, `--version`: print version information and exit.
//!
//!
//! ## Important things to optimize
//! The main parameters governing the number of phylogenetic positions in the final alignment are
//! the k-mer size (-k), the depth of the recursive graph traversal (-d), and the minimum sample
//! frequency (-f).
//!
//! The default k-mer size has already been chosen to maximise the final alignment length, and
//! increasing it usually does not substantially increase the number of variant groups. It may,
//! however, be useful to decrease the k-mer size from 14 to 13 if memory consumption is too high.
//!
//! Increasing the depth of the recursive graph traversal (e.g. from 12 to 16) generally increases
//! the size of the final alignment, as kamino detects more variant groups during graph traversal.
//! This is typically the most effective approach if the alignment is deemed too short, but results 
//! in increased runtimes as the program spends more time traversing the graph.
//!
//! Finally, larger alignments can also be produced by decreasing the minimum fraction of samples
//! required to carry an amino acid (e.g. from 0.85 to 0.8), at the cost of increased missing data
//! in the final alignment. Missing data are represented by '-' (missing amino acid) and 'X'
//! (ambiguous or masked amino acid).
//!
//!
//! ## Less important parameters
//! Besides testing/benchmarking, I would not recommend modifying these parameter values.
//!
//! The number of constant positions in the final alignment can be adjusted with the --constant parameter. These are
//! taken from the left flank of the end k-mer in each variant group, next to the middle positions. Because these
//! positions are recoded, some may become polymorphic once converted back to amino acids. Using the default c=3,
//! constant positions represent 50 to 60% of the alignment.
//!
//! The `--mask parameter` controls the amino-acid masking performed by kamino to prevent long runs of polymorphism from being
//! retained in the final alignment. These correspond to genuine but unwanted polymorphisms (e.g., micro-inversions) or,
//! less frequently, errors made by kamino (“misaligned” paths due to the presence of two consecutive indels). The minimum length
//! of polymorphism runs to be masked can be decreased using this parameter to be more stringent.
//!
//! The `--length-middle` parameter is used to filter out long variant groups. Increase this parameter to allow longer
//! variant groups to be retained in the final alignment.
//!
//! Finally, the 6-letter recoding scheme can be modified using the --recode parameter, although the default sr6 recoding
//! scheme performed the best in most of my tests (sr6 ≥ dayhoff6 ≫ kgb6).
//!
//!
//! ## Output files
//! The names of the output files are controlled by a prefix (-o; default: `kamino`). The prefix
//! may include a directory path (e.g. `-o my_analyses/taxon1`). Note that the output directory is not
//! created by kamino and must already exist.
//!
//! The three output files are:
//! - `<prefix>_alignment.fas`: FASTA amino acid alignment of all samples.
//! - `<prefix>_missing.tsv`: Tab-delimited per-sample missingness percentages.
//! - `<prefix>_partitions.tsv`: Tab-delimited variant group coordinates (0-based) in the FASTA
//!   alignment, along with consensus protein names when the input proteomes are annotated.
//!
//! Additionally, a Neighbor-Joining (NJ) tree can be produced from the amino acid alignment
//! when the `--nj` argument is specified. Pairwise distances are computed using an F81
//! correction with LG stationary amino-acid frequencies. The resulting tree provides an
//! overview of isolate relationships and is not intended for detailed phylogenetic inference.
//!
use anyhow::Context;
use clap::{ArgAction, Parser};

mod bubbles;
mod filter_groups;
mod graph;
mod io;
mod output;
mod phylo;
mod proba_filter;
mod protein_prediction;
mod recode;
mod revert_aminoacid;
mod traverse;
//mod util;

pub use recode::RecodeScheme;

struct TraversalArtifacts {
    groups: traverse::VariantGroups,
    species_names: Vec<String>,
    n_species: usize,
    recode_scheme: RecodeScheme,
}

#[allow(clippy::too_many_arguments)]
fn run_traverse(
    species_inputs: &[io::SpeciesInput],
    k: usize,
    min_freq: f32,
    depth: usize,
    length_middle: usize,
    num_threads: usize,
    recode_scheme: RecodeScheme,
) -> anyhow::Result<TraversalArtifacts> {
    // Build global graph
    let mut g = graph::Graph::new(k, recode_scheme);
    io::build_graph_from_inputs(
        species_inputs,
        k,
        min_freq,
        length_middle,
        &mut g,
        num_threads,
    )?;

    // Basic stats
    io::print_graph_size(&g);
    //io::print_outdegree_histogram(&g);

    // Bubble endpoints
    let (start_kmers, end_kmers) = bubbles::find_bubble_endpoints(&g, min_freq, num_threads);
    eprintln!(
        "bubble endpoints: start={} end={}",
        start_kmers.len(),
        end_kmers.len()
    );

    // Traverse VGs
    let groups = traverse::find_variant_groups(
        &g,
        &start_kmers,
        &end_kmers,
        depth,
        min_freq,
        length_middle,
        num_threads,
    );
    let total_paths: usize = groups.values().map(|v| v.len()).sum();
    eprintln!(
        "variant groups: groups={} paths={}",
        groups.len(),
        total_paths
    );

    // Emit amino-acid sequences per path for each variant group (call disabled, keep for debugging)
    // traverse::print_variant_group_sequences(&g, &groups);

    Ok(TraversalArtifacts {
        groups,
        species_names: g.species_names.clone(),
        n_species: g.n_species,
        recode_scheme: g.recode_scheme,
    })
}

/// Build a node-based, colored de Bruijn graph from amino-acid proteomes and analyze bubbles.
#[derive(Parser, Debug)]
#[command(name = "kamino", author, version, about, disable_version_flag = true)]
#[command(
    group = clap::ArgGroup::new("input_source")
        .required(true)
        .multiple(true)
        .args(["input", "input_file"])
)]
pub struct Args {
    /// Input directory with FASTA proteomes (plain or .gz)
    #[arg(short, long = "input-directory")]
    pub input: Option<std::path::PathBuf>,

    /// Tab-delimited file mapping species name to proteome path
    #[arg(short = 'I', long = "input-file")]
    pub input_file: Option<std::path::PathBuf>,

    /// Treat input files as bacterial genomes and predict proteomes before analysis
    #[arg(long = "genomes", default_value_t = false, hide_default_value = true)]
    pub genomes: bool,

    /// K-mer length [k=14]
    #[arg(short, long)]
    pub k: Option<usize>,

    /// Minimal fraction of samples with an amino-acid per position [m=0.85]
    #[arg(
        short = 'f',
        long = "min-freq",
        default_value_t = 0.85,
        hide_default_value = true
    )]
    pub min_freq: f32,

    /// Maximum traversal depth from each start node [d=12]
    #[arg(short, long, default_value_t = 12, hide_default_value = true)]
    pub depth: usize,

    /// Output prefix [o=kamino]
    #[arg(short, long, default_value = "kamino")]
    pub output: std::path::PathBuf,

    /// Number of constant positions to keep from the in-bubble k-mer [c=3]
    #[arg(short, long)]
    pub constant: Option<usize>,

    /// Maximum number of middle positions per variant group [l=35]
    #[arg(short = 'l', long = "length-middle", default_value_t = 35, hide_default_value = true)]
    pub length_middle: usize,

    /// Mask middle segments with long mismatch runs [m=5]
    #[arg(
        short = 'm',
        long = "mask",
        default_value_t = 5,
        hide_default_value = true
    )]
    pub mask: usize,

    /// Number of threads [t=1]
    #[arg(short = 't', long, default_value_t = 1, hide_default_value = true)]
    pub threads: usize,

    /// Recoding scheme [r=sr6]
    #[arg(short = 'r', long = "recode", value_enum, default_value_t = RecodeScheme::SR6, hide_default_value = true)]
    pub recode: RecodeScheme,

    /// Generate a NJ tree from kamino alignment [nj=false]
    #[arg(long = "nj", default_value_t = false, hide_default_value = true)]
    pub nj: bool,

    /// Display version information.
    #[arg(short = 'v', long = "version", action = ArgAction::Version)]
    pub version: (),
}

pub fn run_with_args(args: Args) -> anyhow::Result<()> {
    // k defaults and limits for SR6
    let default_k = 14usize;
    let max_k = 21usize;
    let default_constant = 3usize;

    anyhow::ensure!(
        args.min_freq >= 0.6,
        "min_freq ({}) must be ≥ 0.6.",
        args.min_freq
    );
    let min_freq = args.min_freq;

    let k = args.k.unwrap_or(default_k);
    anyhow::ensure!(
        (2..=max_k).contains(&k),
        "k={} is invalid for {}: allowed range is 2..={} (default {})",
        k,
        args.recode,
        max_k,
        default_k
    );

    // Validate constant against k-1
    let k1 = k - 1;
    let constant = args.constant.unwrap_or_else(|| default_constant.min(k1));
    anyhow::ensure!(
        constant <= k1,
        "constant ({}) must be ≤ k-1 ({}).",
        constant,
        k1
    );

    anyhow::ensure!(
        args.length_middle >= 1,
        "length_middle ({}) must be ≥ 1.",
        args.length_middle
    );

    // Decide how many threads to use (default: 1)
    anyhow::ensure!(args.threads >= 1, "threads must be ≥ 1");

    eprintln!("kamino v{}", env!("CARGO_PKG_VERSION"));
    let mut species_inputs = Vec::new();
    let mut input_labels: Vec<String> = Vec::new();

    if let Some(table) = args.input_file.as_ref() {
        input_labels.push(format!("input_file={}", table.display()));
        species_inputs.extend(io::collect_species_inputs_from_table(table)?);
    }

    if let Some(dir) = args.input.as_ref() {
        input_labels.push(format!("input={}", dir.display()));
        species_inputs.extend(io::collect_species_inputs_from_dir(dir)?);
    }

    let output_dir = args
        .output
        .parent()
        .map(std::path::Path::to_path_buf)
        .unwrap_or_else(|| std::path::PathBuf::from("."));
    let mut genomes_tmpdir = None;

    let input_label = input_labels.join(" ");

    eprintln!(
        "parameters: k={} constant={} min_freq={} depth={} length_middle={} mask={} threads={} recode={} {} output={}",
        k,
        constant,
        min_freq,
        args.depth,
        args.length_middle,
        args.mask,
        args.threads,
        args.recode,
        input_label,
        args.output.display()
    );

    if args.genomes {
        eprintln!("genome files: {}", species_inputs.len());

        let tmpdir = tempfile::Builder::new()
            .prefix("tmp_kamino_")
            .rand_bytes(6)
            .tempdir_in(&output_dir)
            .with_context(|| {
                format!(
                    "create temporary directory for predicted proteomes in {}",
                    output_dir.display()
                )
            })?;

        let predicted =
            protein_prediction::predict_proteomes(&species_inputs, tmpdir.path(), args.threads)?;
        species_inputs = predicted;
        input_labels.push("genomes=true".to_string());
        genomes_tmpdir = Some(tmpdir);
    }

    anyhow::ensure!(
        !species_inputs.is_empty(),
        "At least one input source with files must be provided."
    );

    let traversal = run_traverse(
        &species_inputs,
        k,
        min_freq,
        args.depth,
        args.length_middle,
        args.threads,
        args.recode,
    )?;

    let (alignment_len, alignment_missing_pct) = output::write_outputs_with_head(
        &species_inputs,
        &args.output,
        &traversal.species_names,
        traversal.n_species,
        traversal.recode_scheme,
        &traversal.groups,
        k,
        constant,
        min_freq,
        args.mask,
        args.threads,
        args.nj,
    )?;

    let (fas_path, tsv_path, partitions_path, tree_path) = output::output_paths(&args.output);
    eprintln!(
        "alignment: length={} missing={:.1}%",
        alignment_len, alignment_missing_pct
    );

    if args.nj {
        eprintln!(
            "output files:  {}, {}, {}, and {}",
            fas_path.display(),
            tsv_path.display(),
            partitions_path.display(),
            tree_path.display()
        );
    } else {
        eprintln!(
            "output files:  {}, {}, and {}",
            fas_path.display(),
            tsv_path.display(),
            partitions_path.display()
        );
    }

    drop(genomes_tmpdir);

    Ok(())
}