Add multi-hash support for common Bloom filter (#88)

* Add multi-hash to ntsynt_make_common_bf.cpp

* Add hashes option to python driver and snakemake files

* Update tests for multi-hash

* Update CI

* Update README.md

* Update btllib version in README.md

* Update default hash number in ntsynt_make_common_bf.cpp

Author: lcoombe

README.md

@@ -41,7 +41,7 @@ Lauren Coombe, Parham Kazemi, Johnathan Wong, Inanc Birol, René L. Warren. 2025
 ## Usage
 
 ```
-usage: ntSynt [-h] [--fastas_list FASTAS_LIST] -d DIVERGENCE [-p PREFIX] [-k K] [-w W] [-t T] [--fpr FPR] [-b BLOCK_SIZE] [--merge MERGE]
+usage: ntSynt [-h] [--fastas_list FASTAS_LIST] -d DIVERGENCE [-p PREFIX] [-k K] [-w W] [-t T] [--fpr FPR] [--hashes HASHES] [-b BLOCK_SIZE] [--merge MERGE]
               [--w_rounds W_ROUNDS [W_ROUNDS ...]] [--indel INDEL] [-n] [--benchmark] [-f] [--dev] [-v]
               [fastas ...]
 
@@ -64,6 +64,7 @@ optional arguments:
   -w W                  Minimizer window size [1000]
   -t T                  Number of threads [12]
   --fpr FPR             False positive rate for Bloom filter creation [0.025]
+  --hashes HASHES       Number of hash functions for Bloom filter creation [3]
   -b BLOCK_SIZE, --block_size BLOCK_SIZE
                         Minimum synteny block size (bp)
   --merge MERGE         Maximum distance between collinear synteny blocks for merging (bp). 
@@ -108,7 +109,7 @@ conda install -c bioconda -c conda-forge ntsynt
   - [ncls](https://github.com/pyranges/ncls)
   - [python-igraph](https://python.igraph.org/en/stable/)
 - GCC 6+ or Clang 5+ (with OpenMP and C++17 support)
-- [btllib v1.6.2+](https://github.com/BirolLab/btllib)
+- [btllib v1.7.8+](https://github.com/BirolLab/btllib)
 - [meson](https://mesonbuild.com/)
 - [ninja](https://ninja-build.org/)
 - [snakemake](https://snakemake.readthedocs.io/en/stable/)
@@ -173,7 +174,7 @@ optional arguments:
 More information can be found on our [wiki page](https://github.com/BirolLab/ntSynt/wiki/de-novo-statistics-summary)
 
 ### Tips / Visualization <a name=tips></a>
-- To lower the peak memory usage, increase the false positive rate (--fpr) for the constructed Bloom filter
+- To lower the peak memory usage, increase the number of hash functions (--hashes) and/or increase the false positive rate (--fpr) for the constructed Bloom filter
 - Customize parameters such as --merge, --indel, --block_size and --w_rounds for your particular input data and research questions
 - For visualizing the multi-genome output synteny blocks, check out our streamlined visualization tool [ntSynt-viz](https://github.com/BirolLab/ntSynt-viz). Examples of the scripts used to generate the figures in the [ntSynt manuscript](https://link.springer.com/article/10.1186/s12915-025-02455-w) can be found in the sub-directory [visualization_scripts](https://github.com/BirolLab/ntSynt/tree/main/visualization_scripts)
 - If you do not know the approximate sequence divergence between the input assemblies, we recommend using [Mash](https://github.com/marbl/Mash) to estimate the divergences

azure-pipelines.yml

@@ -31,7 +31,7 @@ jobs:
   - script: |
       source activate ntsynt_ci
       mamba install --yes -c bioconda -c conda-forge python=3.12
-      mamba install --yes -c bioconda -c conda-forge intervaltree pybedtools ncls 'python-igraph>=1.0.0' compilers btllib meson ninja snakemake samtools pytest seqtk pylint
+      MAMBA_USE_SHARDED_REPODATA=OFF mamba install --yes -c bioconda -c conda-forge intervaltree pybedtools ncls 'python-igraph>=1.0.0' compilers btllib meson ninja snakemake samtools pytest seqtk pylint
     displayName: Install dependencies
 
   - script: |
@@ -81,7 +81,8 @@ jobs:
 
   - script: |
       mkdir -p ~/miniforge3
-      curl -L https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-MacOSX-x86_64.sh  -o ~/miniforge3/miniforge.sh
+      # May 2026: latest miniforge uses a mamba version which is much slower. Pin to lower version for now.
+      curl -L https://github.com/conda-forge/miniforge/releases/download/26.1.0-0/Miniforge3-MacOSX-x86_64.sh  -o ~/miniforge3/miniforge.sh
       bash ~/miniforge3/miniforge.sh -b -u -p ~/miniforge3
       rm -rf  ~/miniforge3/miniforge.sh
       ~/miniforge3/bin/conda init bash
@@ -95,7 +96,8 @@ jobs:
 
   - script: |
       source activate ntsynt_ci
-      mamba install --yes  -c conda-forge -c bioconda intervaltree pybedtools ncls python-igraph compilers  btllib meson ninja snakemake samtools pytest seqtk pigz
+      mamba --version
+      mamba install --yes -c bioconda -c conda-forge intervaltree pybedtools ncls python-igraph compilers  'btllib>=1.7.8' meson ninja snakemake samtools pytest seqtk pigz
     displayName: Install dependencies
 
   - script: |

bin/ntSynt

@@ -7,6 +7,7 @@ import argparse
 import os
 import shlex
 import subprocess
+import re
 from packaging import version
 
 import snakemake
@@ -30,6 +31,15 @@ def read_fasta_files(filename):
             fastas.append(fasta.strip())
     return fastas
 
+def get_btllib_version():
+    "Return the version of btllib (By checking indexlr) installed"
+    result = subprocess.run(["indexlr", "--version"], text=True, check=True,
+                            stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
+    match = re.search(r"(\d+\.\d+\.\d+)", result.stdout)
+    if not match:
+        raise ValueError(f"Could not parse version from: {result.stdout!r}")
+    return version.Version(match.group(1))
+
 def main():
     "Run ntSynt snakemake file"
 
@@ -59,6 +69,8 @@ def main():
     parser.add_argument("-t", help="Number of threads [12]", type=int, default=12)
     parser.add_argument("--fpr", help="False positive rate for Bloom filter creation [0.025]",
                         default=0.025, type=float)
+    parser.add_argument("--hashes", help="Number of hash functions for Bloom filter creation [3]", default=3,
+                        type=int)
     parser.add_argument("-b", "--block_size", help="Minimum synteny block size (bp)",
                         type=int, required=False)
     parser.add_argument("--merge", help="Maximum distance between collinear synteny blocks for merging (bp). \n"\
@@ -110,6 +122,14 @@ def main():
     if args.fastas and args.fastas_list:
         parser.error("Please supply the input genome fasta files as positional arguments, "\
                     "or specify a single file (one fasta per line) with --fastas_list, NOT both.")
+        
+    if args.hashes > 1:
+        btllib_version = get_btllib_version()
+        if btllib_version < version.Version("1.7.8"):
+            print("WARNING: Using > 1 hash function for Bloom filter creation requires btllib 1.7.8 or higher. " \
+                    f"Current version: {btllib_version}. Please update btllib to use this feature.", flush=True)
+            print("WARNING: Setting --hashes 1 to proceed with ntSynt run.", flush=True)
+            args.hashes = 1
 
     if args.fastas_list:
         fastas_input = read_fasta_files(args.fastas_list)
@@ -134,6 +154,7 @@ def main():
                     f"\t-w {args.w}",
                     f"\t-t {args.t}",
                     f"\t--fpr {args.fpr}",
+                    f"\t--hashes {args.hashes}",
                     ]
     print("\n".join(intro_string), flush=True)
 
@@ -146,7 +167,7 @@ def main():
     command = f"snakemake -s {base_dir}/ntsynt_run_pipeline.smk -p --cores {args.t} " \
                 f"--config references='{fastas_input}' kmer={args.k} window={args.w} threads={args.t} fpr={args.fpr} " \
                 f"prefix={args.prefix} w_rounds='{args.w_rounds}' indel_merge={args.indel} " \
-                f"collinear_merge={args.merge} block_size={args.block_size} "
+                f"collinear_merge={args.merge} block_size={args.block_size} hashes={args.hashes} "
     command += "common=False " if args.no_common else "common=True "
     command += "simplify_graph=False " if args.no_simplify_graph else "simplify_graph=True "
     command += "benchmark=True " if args.benchmark else "benchmark=False "

bin/ntsynt_run_pipeline.smk

@@ -11,6 +11,7 @@ references = config["references"] if "references" in config else "Must specify '
 k = config["kmer"] if "kmer" in config else 24
 w = config["window"] if "window" in config else 1000
 fpr = config["fpr"] if "fpr" in config else 0.025
+hashes = config["hashes"] if "hashes" in config else 3
 max_threads = config["threads"] if "threads" in config else 4
 prefix = config["prefix"] if "prefix" in config else "ntSynt_out"
 common = config["common"] if "common" in config else True
@@ -56,7 +57,7 @@ rule make_common_bf:
     input: refs=references
     output: f"{prefix}.common.bf"
     threads: max_threads
-    params: options=f"-p {prefix}.common --fpr {fpr} -k {k}",
+    params: options=f"-p {prefix}.common --fpr {fpr} --hashes {hashes} -k {k}",
             path_to_script=f"{script_path}/ntsynt_make_common_bf",
             benchmarking=f"{benchmark_path} -o {prefix}.make_common_bf.time" if benchmark else []
     shell: "{params.benchmarking} {params.path_to_script} --genome {input.refs} {params.options} -t {threads}"

src/ntsynt_make_common_bf.cpp

@@ -15,8 +15,6 @@ Creating a common Bloom filter (with k-mers found in all input assemblies) using
 a cascading Bloom filter approach
 */
 
-// Number of hash functions - keep at 1
-const unsigned HASH_FNS = 1;
 // Converting bits to bytes
 const unsigned NUM_BITS_PER_BYTE = 8;
 
@@ -26,7 +24,7 @@ Using the formula found in Broder & Mitzenmacher, 2004, simplified for
 one hash function.
 */
 long long
-approximate_bf_size(std::string genome_file, double fpr, int threads)
+approximate_bf_size(std::string genome_file, double fpr, int hashes, int threads)
 {
   long long genome_size = 0;
   btllib::SeqReader reader(
@@ -35,7 +33,9 @@ approximate_bf_size(std::string genome_file, double fpr, int threads)
     genome_size += record.seq.length();
   }
   std::cout << "Genome size (bp): " << genome_size << std::endl;
-  long long size_bits = ceil(((-1 * genome_size) / log(1 - fpr)));
+  double m_bits = -1.0 * hashes * genome_size / log(1.0 - pow(fpr, 1.0 / hashes));
+
+  long long size_bits = (long long) std::ceil(m_bits);
   return size_bits / NUM_BITS_PER_BYTE;
 }
 
@@ -65,6 +65,11 @@ main(int argc, const char** argv)
   parser.add_argument("--bf")
     .help("Bloom filter size in bytes (optional)")
     .scan<'d', long long>();
+  
+  parser.add_argument("--hashes")
+    .help("Number of hash functions")
+    .default_value(3U)
+    .scan<'u', unsigned>();
 
   parser.add_argument("-t")
     .help("Number of threads")
@@ -85,6 +90,7 @@ main(int argc, const char** argv)
   unsigned num_threads = parser.get<unsigned>("t");
   double fpr = parser.get<double>("fpr");
   unsigned k = parser.get<unsigned>("k");
+  unsigned hashes = parser.get<unsigned>("hashes");
   std::string prefix = parser.get<std::string>("p");
 
   std::cout << "Parameters:" << std::endl;
@@ -96,6 +102,7 @@ main(int argc, const char** argv)
   std::cout << "\t\t-t " << num_threads << std::endl;
   std::cout << "\t\t-k " << k << std::endl;
   std::cout << "\t\t--fpr " << fpr << std::endl;
+  std::cout << "\t\t--hashes " << hashes << std::endl;
   std::cout << "\t\t-p " << prefix << std::endl;
 
   #if _OPENMP
@@ -113,14 +120,14 @@ main(int argc, const char** argv)
     std::cout << "\t\t--bf " << bf_size << std::endl;
   } else {
     std::cout << "Calculating BF size based on input genome size" << std::endl;
-    bf_size = approximate_bf_size(genome_files[0], fpr, num_threads);
+    bf_size = approximate_bf_size(genome_files[0], fpr, hashes, num_threads);
   }
   std::cout << "BF size (bytes): " << bf_size << std::endl;
 
 
   /* Load the initial level 1 BF with k-mers from the first genome  */
   btllib::KmerBloomFilter* bf =
-    new btllib::KmerBloomFilter(bf_size, HASH_FNS, k);
+    new btllib::KmerBloomFilter(bf_size, hashes, k);
   btllib::log_info("Reading " + genome_files[0]);
   btllib::SeqReader reader(
     genome_files[0], btllib::SeqReader::Flag::LONG_MODE, num_threads);
@@ -134,7 +141,7 @@ main(int argc, const char** argv)
   /* For each subsequent genome, only insert into next level BF if the k-mer is
   in the previous level BF. */
   btllib::KmerBloomFilter* new_bf =
-    new btllib::KmerBloomFilter(bf_size, HASH_FNS, k);
+    new btllib::KmerBloomFilter(bf_size, hashes, k);
 
   size_t num_assemblies = genome_files.size();
   for (size_t i = 1; i < num_assemblies; ++i) {
@@ -144,7 +151,7 @@ main(int argc, const char** argv)
       genome, btllib::SeqReader::Flag::LONG_MODE, num_threads);
     #pragma omp parallel
     for (const auto record : reader) {
-      btllib::NtHash nthash(record.seq, HASH_FNS, k);
+      btllib::NtHash nthash(record.seq, hashes, k);
       while (nthash.roll()) {
         if (bf->contains(nthash.hashes())) {
           new_bf->insert(nthash.hashes());
@@ -155,7 +162,7 @@ main(int argc, const char** argv)
     delete bf;
     bf = new_bf;
     if (i < num_assemblies - 1) {
-      new_bf = new btllib::KmerBloomFilter(bf_size, HASH_FNS, k);
+      new_bf = new btllib::KmerBloomFilter(bf_size, hashes, k);
     }
   }