lexicmap genome search

$ lexicmap genome search -h
Search genomes against an index, with ANI and AF computed

Algorithm:
  1. Genome screening:
     Candidate genomes are screened by the total length of shared seeds (k-mers longer
     than the value of -p/--seed-min-prefix) between the query genome and the genomes
     in the index, with lexichash masking. Seeds never cross query contig boundaries.

     Flag -n/--top-n-genomes keeps score groups in descending order until their total
     genome count reaches or exceeds N, including the entire cutoff group (0 for all).
     For scores 10, 10, 10, 10, 9, 9, 9, 8, 7, 6, -n 3 keeps the four genomes scoring 10;
     -n 5 keeps those four plus the three genomes scoring 9.

     For highly redundant databases such as AllTheBacteria, which may contain hundreds of
     thousands of highly similar genomes, the screening parameters can be moderately tuned
     together to improve discrimination. For example: -p 21 -m 4096 -n 100

     -----------------------------------------------------------------------------------
      Flag   Controls               Effect of a larger value
     -----------------------------------------------------------------------------------
      -p     Minimum seed prefix    Fewer seed matches and lower memory use; may lose
             length                 shorter matches and reduce nearest-neighbor recall.
      -m     Number of masks        Better discrimination; slower screening and much
                                    higher memory use.
      -n     Candidate genomes      Higher sensitivity but slower alignment; all ties
             retained               at the cutoff score are retained.
     -----------------------------------------------------------------------------------

  2. Alignment:
     a) (default): cut the query genome into non-overlapping fragments and align them to
        the candidate genomes with an approach similar to that used in 'lexicmap search',
        with steps of seed matching, chaining, pseudoalignment, and base-level alignment.
     b) (OrthoANI mode): cut both the query genome and candidate genomes into non-overlapping
        fragments and only orthologous fragment pairs are used for calculating the ANI and
        AF values, which is similar to the algorithm of OrthoANI.

     Flag -n/--top-n-genomes keeps all matches with the top N distinct ANI percentages
     as printed with four decimal places. Both stages may retain more than N genomes.

     All chunks of each retained genome are included in alignment and genome-size
     denominators, including chunks without screening seed hits.

     Results are sorted by ANI percentages rounded to four decimal places, then by
     full-precision qAF and sAF (descending).

     Repeated sequences can produce many k-mer matches, so alignment time can differ
     substantially between genomes of similar sizes.

     Query-fragment indexes are reused across candidate genomes within
     --max-genome-cache-memory (default 4G, shared across concurrent queries;
     0 disables reuse). Increase it for large or concurrent queries when memory
     permits. Other search memory is additional.

Attention:
  1. Input should be (gzipped) FASTA records from files or stdin, with one genome per file.
     Experimental feature: Long-read FASTQ files are also accepted. For these inputs,
     --adjust-ani-by-quality can be used to estimate ANI corrected for sequencing errors:
         e = mean_i(10^(-Q_i/10)) and ANIadj = (ANI - e/3) / (1 - 4e/3),
     where Q_i are the Phred scores in aligned query regions and ANI is expressed as a
     fraction. The correction assumes independent, symmetric substitution errors.
  2. One or more input files are accepted, via positional parameters
     and/or a file list via the flag -X/--infile-list.
  3. For multiple queries, the order of queries in output might be different from the input.
  4. With an index built using --soft-masking, lowercase query bases are converted
     to A for both seed screening and final alignment. This can affect ANI and AF
     for soft-masked genomes, such as fungal genomes; use unmasked queries when possible.
     Ambiguous bases (including N) are also replaced with A/C/G/T in queries and
     stored genomes, so ANI and AF in ambiguous regions require caution.

Taxonomic operations:
  1. Taxonomy data, including NCBI-format taxdump files (-T/--taxdump) and a genome-ID-to-TaxId
     mapping file (-G/--genome2taxid), are needed for filtering genomes by TaxId(s) and/or
     showing taxonomic names.

     Taxdump files can be createted from any taxonomy data with TaxonKit, see
     https://bioinf.shenwei.me/taxonkit/usage/#create-taxdump

     If -T and -G are not provided, it will try to find them in the index directory,
     in the subdirectory "taxdump" and the file "taxdump/taxid.map", respectively.

  2. Users can limit search by TaxId(s) via -t/--taxids or --taxid-file.
     Only genomes with descendant TaxIds of the specific ones or themselves are searched,
     in a similar way with BLAST+ 2.15.0 or later versions.
     Negative values are allowed as a black list.

     For example, searching non-Escherichia (561) genera of Enterobacteriaceae (543) family with
     -t 543,-561.

  3. Flag --show-genome-name and --show-species-name can be used to show the taxonomic/species
     name of the subject genome in the output.

Output format:
  Tab-delimited format with 10 columns. --adjust-ani-by-quality appends
  ANIadj as the last column.

    1.  query,    Query genome ID.
    2.  subject,  Subject genome ID.
    3.  ANI,      Average nucleotide identity.
    4.  qAF,      Align fraction of the query genome (sum_aligned_length / query_genome_size).
    5.  sAF,      Align fraction of the subject genome (sum_aligned_length / subject_genome_size).
                  Terminal pieces shorter than half a fragment join the preceding fragment
                  of the same contig. Contigs shorter than the alignment seed length are
                  not aligned but remain included in genome sizes for AF.
    6.  qctgs,    Number of contigs in the query genome.
    7.  qsize,    Size of the query genome.
    8.  sctgs,    Number of contigs in the subject genome.
    9.  ssize,    Size of the subject genome.
    10. sname,    Taxonomic name of the subject name.
    11. ANIadj,   Quality-adjusted ANI (optional, experimental feature).

Usage:
  lexicmap genome search [flags] -d <index path> [query.fasta[.gz] ...] [-o result.tsv[.gz]]

Flags:
  -O, --OrthoANI                         ► Compute OrthoANI using reciprocal best hit of fragment
                                         pairs. Type 'lexicmap genome search --help' for details.
      --adjust-ani-by-quality            ► (Experimental feature) Estimate ANI corrected for Phred+33
                                         sequencing errors in the aligned query regions and add the
                                         ANIAdjusted column. FASTQ input is required. The raw ANI is
                                         still used for filtering and sorting.
      --align-band int                   ► Band size in backtracking the score matrix (pseudo
                                         alignment phase). (default 100)
      --align-max-gap int                ► Maximum gap in a HSP segment. (default 100)
  -l, --align-min-match-len int          ► Minimum aligned length in a HSP segment. (default 30)
  -i, --align-min-match-pident float     ► Minimum base identity (percentage) in a HSP segment.
                                         (default 70)
      --debug                            ► Print debug information, including a progress bar.
                                         (recommended when searching with one query).
      --extra                            ► Show extra columns for -S/--only-genome-screening.
      --frag-size int                    ► The size of non-overlap fragments cut for ANI computation.
                                         (default 1020)
      --gc-interval int                  ► Force garbage collection every N queries (0 for disable).
                                         The value can't be too small. (default 4)
  -G, --genome2taxid string              ► Two-column tabular file for mapping genome ID to TaxId,
                                         needed for filtering results with TaxIds. Genome IDs in the
                                         index can be exported via "lexicmap utils genomes -d db.lmi/ |
                                         csvtk cut -t -f 1 | csvtk uniq -Ut". (default: <taxdump
                                         path>/taxid.map)
  -h, --help                             help for search
  -d, --index string                     ► Index directory created by "lexicmap index".
  -k, --keep-genomes-without-taxid       ► Keep genome hits without TaxId, i.e., those without TaxId
                                         in the --genome2taxid file.
      --kmer-scale int                   ► Using 1/scale of k-mers for seeding (default mode) or
                                         fragment comparison (OrthoANI mode). Available values: 2, 4, 8.
                                         (default 4)
  -w, --load-whole-seeds                 ► Load the whole seed data into memory for faster seed
                                         matching. It will consume a lot of RAM.
  -m, --masks int                        ► Only use seed data of N masks for genome
                                         prefiltering/screening. It should be 0 (for all masks in the
                                         index) or power of 4 (needs to be >= 256, e.g., 256, 1024,
                                         4096, 16384). (default 1024)
  -e, --max-evalue float                 ► Maximum evalue of a HSP segment. (default 1e-15)
      --max-genome-cache-memory string   ► Memory budget for reusing query-fragment indexes across
                                         candidate genomes, shared by concurrent queries (supports
                                         K/M/G; 0 disables reuse). Other search memory is additional.
                                         (default "4G")
      --max-open-files int               ► File budget for index and query readers, including 9
                                         reserved descriptors. Minimum: 9 + seed chunks + genome batches
                                         (9 + genome batches with -w). Every seed chunk and genome batch
                                         retains at least one reader; insufficient budgets are rejected.
                                         Query-file and index-reader concurrency are reduced to fit. Set
                                         "ulimit -n" at least this high. (default 1024)
  -J, --max-query-conc int               ► Maximum number of concurrent queries. (default 8)
      --max-subject-genome-size int      ► Maximum size of subject genomes to be considered (in MB).
                                         (default 20)
  -F, --min-af float                     ► Only output results where one genome has aligned fraction >
                                         than this value (percentage). (default 15)
  -I, --min-ani float                    ► Only output results where one genome has ANI > than this
                                         value (percentage). (default 70)
  -q, --min-qcov-per-hsp float           ► Minimum query coverage (percentage) per HSP. (default 30)
  -S, --only-genome-screening            ► Only perform genome screening, no ANI computation.
  -o, --out-file string                  ► Out file, supports a ".gz" suffix ("-" for stdout).
                                         (default "-")
      --ref-name-regexp string           ► Regular expression with at least one capturing group (the
                                         first is used) for extracting the reference name from the input
                                         filename. Attention: use double quotation marks for patterns
                                         containing commas, e.g., -p '"A{2,}"'. (default
                                         "(?i)(.+)\\.(f[aq](st[aq])?|fna)(\\.gz|\\.xz|\\.zst|\\.bz2)?$")
  -p, --seed-min-prefix int              ► Minimum prefix length of matched seeds in the genome
                                         filtering phase. (default 21)
  -g, --show-genome-name                 ► Show the taxonomic name of subject genome in field 'sname'.
                                         Flags -T/--taxdump and -G/--genome2taxid are needed.
  -s, --show-species-name                ► Show the species name of subject genome in field 'sname'.
                                         Flags -T/--taxdump and -G/--genome2taxid are needed.
  -T, --taxdump string                   ► Directory containing taxdump files (nodes.dmp, names.dmp,
                                         etc.), needed for filtering results with TaxIds. For other
                                         non-NCBI taxonomy data, please use 'taxonkit create-taxdump' to
                                         create taxdump files. (default: <index path>/taxdump)
      --taxid-file string                ► TaxIds from a file for filtering results, where the taxids
                                         are equal to or are the children of the given taxids. Negative
                                         values are allowed as a black list.
  -t, --taxids strings                   ► TaxIds(s) for filtering results, where the taxids are equal
                                         to or are the children of the given taxids. Negative values are
                                         allowed as a black list.
  -N, --top-n-chains int                 ► Keep the top N chains in a genome for the query (0 for all)
                                         in the chaining phase. Value 1 is not recommended as the best
                                         chaining result does not always bring the best alignment.
                                         (default 5)
  -n, --top-n-genomes int                ► Keep the top N screened genomes, including all ties at the
                                         cutoff score, then all matches with the top N distinct ANI
                                         percentages as printed with four decimal places (0 for all).
                                         May return more than N genomes. (default 10)
      --windows int                      ► The number of windows in lexichash masking, for genome
                                         screening (1 is sufficient). (default 1)

Global Flags:
      --force-gzip-output    ► Gzip-compress result output regardless of the filename, including stdout.
  -X, --infile-list string   ► File of input file list (one file per line). If given, they are
                             appended to files from CLI arguments.
      --log string           ► Log file.
      --quiet                ► Do not print any verbose information. But you can write them to a file
                             with --log.
  -j, --threads int          ► Number of CPU cores to use. By default, it uses all available cores.
                             (default 16)