Merge pull request #38 from jolespin/dev

Refactoring for performance improvements in v2026.7.16

Author: jolespin

CHANGELOG.md

@@ -1,5 +1,9 @@
 #### Daily Change Log:
-* [2026.3.11] - Added support for `-r/--single_reads` in both `profile-taxonomy.py` and `profile-pathway.py`. However, the `profile-pathway.py` needs `oarfish` to use long reads. [issue/#24](https://github.com/jolespin/leviathan/issues/24)
+* [2026.7.16] - Added parallelization support for pathway coverage using `profile_pathway_coverage` from `kegg_pathway_profiler>=2026.7.16` in `leviathan-profile-pathway.py`. The `--n_jobs` flag now applies to both `salmon quant` and pathway coverage computation.
+* [2026.7.16] - Added `--deterministic` flag to `leviathan-profile-pathway.py` which passes `--deterministic` to `salmon quant` for byte-identical results across runs and thread counts.
+* [2026.7.8] - Added `--update_salmon_index` to `leviathan-index.py` and remove `-u` alias on `--update_with_genomes`
+* [2026.7.8] - **BREAKING CHANGE:** `Salmon` is reimplemented with `Rust` from `C++` so any `Leviathan` databases using `Salmon 1.x` will need to be updated.  `Leviathan` now requires `Salmon ≥ 2.x`.
+* [2026.3.11] - Added support for `-r/--single_reads` in both `profile-taxonomy.py` and `profile-pathway.py`. However, the `profile-pathway.py` needs `oarfish` to use long reads but there are no synthetic long-read metatranscriptomics to benchmark. [issue/#24](https://github.com/jolespin/leviathan/issues/24)
 * [2026.3.10] - Added `--veba_major_version` to `compile-manifest-from-veba.py` since `cluster` output directory is changing
 * [2026.3.3] - Added `step_coverage` output files [issue/#22](https://github.com/jolespin/leviathan/issues/22)
 * [2025.12.17] - Added `--table_format parquet|tsv` to `leviathan-merge.py` and merge functions in `leviathan-profile-pathway.py`/`leviathan-profile-taxonomy.py` [Issue #20](https://github.com/jolespin/leviathan/issues)

FAQ.md

@@ -0,0 +1,28 @@
+# Frequently Asked Questions
+## 1. Is it required to run `leviathan-profile-taxonomy` before `leviathan-profile-pathway` or vice versa? 
+No, they are functionally independent.  However, you can always filter the pathway profiling results by taxonomic filters post hoc. 
+
+```python
+# Load taxonomic profiling
+ds_taxonomic = xr.open_dataset("leviathan_output/artifacts/taxonomic_abundances.genome_clusters.nc")
+
+# Load functional profiling
+ds_pathway = xr.open_dataset("leviathan_output/artifacts/pathway.genome_clusters.nc")
+
+# Get (pan)genomes to retain from taxonomic abundances
+minimum_samples_detected = 2
+retained_organisms = (ds_taxonomic["taxonomic_abundances"].to_pandas() > 0).sum(axis=0)[lambda x: x >= minimum_samples_detected].index
+
+# Filter functional profiling with taxonomic abundance gate
+ds_pathway_filtered = ds_pathway.sel(genome_clusters=retained_organisms)
+```
+## 2. How can I estimate the reference coverage (i.e., percent of assigned reads)?
+```python
+# Load taxonomic profiling
+ds_taxonomic = xr.open_dataset("leviathan_output/artifacts/taxonomic_abundances.genome_clusters.nc")
+
+# Sum the sequence abundances (not taxonomic abundances)
+percent_assigned_reads = ds_taxonomic["sequence_abundances"].sum(axis=1).to_pandas()
+percent_unassigned_reads = 100 - percent_assigned_reads
+```
+## 3. Do I need my own genomes or are there reference catalogs?

README.md

@@ -24,69 +24,7 @@ Leviathan: A fast, memory-efficient, and scalable taxonomic and pathway profiler
 Detailed explanation on how to run each module including downloading test data and interpreting output files. 
 
 ## Benchmarking
-### Benchmarking against 10, 100, 1000, and 10000 genomes
-Benchmarking using trimmed `SRR12042303` sample with 4 threads (ml.m5.4xlarge)
-
-| number_of_genomes | number_of_cds_with_features | preprocess | index | profile-taxonomy | profile-pathway |
-|-------------------|-----------------------------|------------|-------|------------------|-----------------|
-| 10                | 1928                        | 0:03       | 0:09  | 0:41             | 2:09            |
-| 100               | 18410                       | 0:31       | 0:26  | 0:41             | 4:29            |
-| 1000              | 191155                      | 5:29       | 3:55  | 0:43             | 12:50           |
-| 10000             | 1684876                     | 46:00      | 39:10 | 0:48             | 18:14           |
-
-### Benchmarking against CAMI-I and CAMI-II using 16 threads
-All benchmarking and analysis was performed using a virtual machine with the following
-specifications: Linux Ubuntu 22.04 64-bit (x86_64), 30 Intel Xeon Platinum 8358 CPU, 222 GB
-memory, and 1 NVIDIA A10 GPU. Benchmarking and analysis was performed using 16 threads
-running 2 jobs simultaneously for *Leviathan* and *HUMAnN*. 
-
-#### Computational Performance
-|                 |          | Leviathan            |                    | HUMAnN               |                    | Fold   Improvement |        |
-|-----------------|----------|----------------------|--------------------|----------------------|--------------------|--------------------|--------|
-|                 |          | Duration   (minutes) | Peak   Memory (GB) | Duration   (minutes) | Peak   Memory (GB) | Duration           | Memory |
-| CAMI_high_toy   | H_S001   | 14.61                | 2.34               | 1083.57              | 32.31              | 74.19              | 13.84  |
-|                 | H_S002   | 14.89                | 2.35               | 949.73               | 32.28              | 63.78              | 13.75  |
-|                 | H_S003   | 14.96                | 2.34               | 875.83               | 32.64              | 58.56              | 13.97  |
-|                 | H_S004   | 15.02                | 2.35               | 852.18               | 32.33              | 56.72              | 13.79  |
-|                 | H_S005   | 15.27                | 2.33               | 826.25               | 32.23              | 54.13              | 13.82  |
-| CAMI_medium_toy | M2_S001  | 5.78                 | 1.62               | 219.25               | 15.95              | 37.96              | 9.83   |
-|                 | M2_S002  | 5.81                 | 1.62               | 174.12               | 16.96              | 29.95              | 10.45  |
-| CAMI_low_toy    | S_S001   | 3.29                 | 1.27               | 76.90                | 10.00              | 23.40              | 7.87   |
-| Marine          | sample_0 | 13.78                | 2.70               | 119.52               | 17.92              | 8.68               | 6.63   |
-|                 | sample_1 | 15.22                | 2.69               | 121.30               | 18.00              | 7.97               | 6.69   |
-|                 | sample_2 | 14.97                | 2.71               | 120.27               | 17.99              | 8.03               | 6.65   |
-|                 | sample_3 | 17.10                | 2.71               | 124.05               | 17.83              | 7.25               | 6.59   |
-|                 | sample_4 | 14.32                | 2.74               | 118.47               | 17.82              | 8.27               | 6.51   |
-|                 | sample_5 | 15.53                | 2.72               | 119.40               | 17.80              | 7.69               | 6.54   |
-|                 | sample_6 | 16.09                | 2.71               | 119.72               | 17.91              | 7.44               | 6.62   |
-|                 | sample_7 | 14.90                | 2.73               | 119.92               | 17.92              | 8.05               | 6.56   |
-|                 | sample_8 | 16.41                | 2.72               | 121.73               | 17.95              | 7.42               | 6.61   |
-|                 | sample_9 | 14.45                | 2.73               | 118.87               | 17.79              | 8.23               | 6.51   |
-
-#### Accuracy Performance
-Ranges from 0.0 - 1.0
-
-|                 | Accuracy | Leviathan |           | HUMAnN |           | Improvement |           |
-|-----------------|----------|-----------|-----------|--------|-----------|-------------|-----------|
-| Dataset         | SampleID | Genome    | Pangenome | Genome | Pangenome | Genome      | Pangenome |
-| CAMI_high_toy   | H_S001   | 0.9492    | 0.9970    | 0.9049 | 0.9610    | 0.0442      | 0.0360    |
-|                 | H_S002   | 0.9551    | 0.9899    | 0.8992 | 0.9591    | 0.0558      | 0.0308    |
-|                 | H_S003   | 0.9556    | 0.9888    | 0.9004 | 0.9598    | 0.0553      | 0.0290    |
-|                 | H_S004   | 0.9496    | 0.9872    | 0.8947 | 0.9588    | 0.0548      | 0.0284    |
-|                 | H_S005   | 0.9420    | 0.9877    | 0.8901 | 0.9573    | 0.0519      | 0.0304    |
-| CAMI_medium_toy | M2_S001  | 0.9692    | 0.9983    | 0.9101 | 0.9620    | 0.0591      | 0.0363    |
-|                 | M2_S002  | 0.9762    | 0.9988    | 0.9177 | 0.9650    | 0.0585      | 0.0338    |
-| CAMI_low_toy    | S_S001   | 1.0000    | 1.0000    | 0.9845 | 0.9845    | 0.0155      | 0.0155    |
-| Marine          | sample_0 | 0.9727    | 0.9933    | 0.8783 | 0.9538    | 0.0944      | 0.0396    |
-|                 | sample_1 | 0.9298    | 0.9922    | 0.8793 | 0.9554    | 0.0505      | 0.0367    |
-|                 | sample_2 | 0.9686    | 0.9817    | 0.8768 | 0.9393    | 0.0918      | 0.0424    |
-|                 | sample_3 | 0.9706    | 0.9842    | 0.8596 | 0.9517    | 0.1110      | 0.0325    |
-|                 | sample_4 | 0.9661    | 0.9880    | 0.8454 | 0.9389    | 0.1207      | 0.0491    |
-|                 | sample_5 | 0.9614    | 0.9856    | 0.8740 | 0.9612    | 0.0874      | 0.0244    |
-|                 | sample_6 | 0.9283    | 0.9869    | 0.8684 | 0.9574    | 0.0599      | 0.0295    |
-|                 | sample_7 | 0.9231    | 0.9942    | 0.8719 | 0.9466    | 0.0512      | 0.0476    |
-|                 | sample_8 | 0.9703    | 0.9889    | 0.8764 | 0.9488    | 0.0940      | 0.0401    |
-|                 | sample_9 | 0.9459    | 0.9859    | 0.8657 | 0.9548    | 0.0802      | 0.0311    |
+Please refer to the publication for benchmarking specs but you should be able to run this easily on a machine with 16GB of RAM. 
 
 ## Modules
 ### `leviathan-preprocess`

WALKTHROUGH.md

@@ -250,16 +250,16 @@ Here are the output files:
 
 ## 7. Reading output files
 ### 7a. Reading Parquet files with Pandas
-
 ```python
-!pip install fastparquet
+!pip install pyarrow
 df = pd.read_parquet("path/to/file.parquet")
 ```
 
 ### 7b. Reading NetCDF files with Xarray
-
+> [!NOTE]
+> To load `.nc` files you must have `h5netcdf` installed before `xarray` is loaded
 ```python
-!pip install xarray
+!pip install xarray h5netcdf
 import xarray as xr
 
 # Taxonomic abundances for genomes
@@ -291,8 +291,34 @@ Data variables:
 
 ```
 
-### 7c. Reformatting Xarray NetCDF files into Pandas DataFrames
+### 7c. Estimating the percent of unassigned reads (i.e., reference coverage)
+Sequence coverage only sums to 100% if every read was assigned so you can subtract from 100% to determine the percent of unassigned reads
+```python
+# Load taxonomic profiling
+ds_taxonomic = xr.open_dataset("leviathan_output/artifacts/taxonomic_abundances.genome_clusters.nc")
+
+# Sum the sequence abundances (not taxonomic abundances)
+percent_assigned_reads = ds_taxonomic["sequence_abundances"].sum(axis=1).to_pandas()
+percent_unassigned_reads = 100 - percent_assigned_reads
+```
+
+### 7d. Filtering functional profiling results with taxonomy gate
+```python
+# Load taxonomic profiling
+ds_taxonomic = xr.open_dataset("leviathan_output/artifacts/taxonomic_abundances.genome_clusters.nc")
+
+# Load functional profiling
+ds_pathway = xr.open_dataset("leviathan_output/artifacts/pathway.genome_clusters.nc")
+
+# Get (pan)genomes to retain from taxonomic abundances
+minimum_samples_detected = 2
+retained_organisms = (ds_taxonomic["taxonomic_abundances"].to_pandas() > 0).sum(axis=0)[lambda x: x >= minimum_samples_detected].index
+
+# Filter functional profiling with taxonomic abundance gate
+ds_pathway_filtered = ds_pathway.sel(genome_clusters=retained_organisms)
+```
 
+### 7e. Reformatting Xarray NetCDF files into Pandas DataFrames
 ```python
 # Load pathway abundances
 ds_pathway = xr.open_dataset("leviathan_output/artifacts/pathway.genome_clusters.nc")
@@ -320,7 +346,7 @@ X_coverage = X_coverage.loc[:,features_passed_qc]
 # Downstream analysis with filtered `X_counts`
 ```
 
-### 7d. Selecting coverage `prevalence` cutoff
+### 7f. Selecting coverage `prevalence` cutoff
 
 ```python
 import compositional as coda

bin/leviathan-index.py

@@ -64,7 +64,8 @@ def main(args=None):
     parser_io.add_argument("-m","--feature_mapping", type=str,  help = "path/to/feature_mapping.tsv [id_gene, feature_set, id_genome, (Optional: id_genome_cluster)] (No header)")
     parser_io.add_argument("-g","--genomes", type=str, help = "path/to/genomes.tsv [id_genome, path/to/genome] (No header)")
     parser_io.add_argument("-d","--index_directory", type=str, required=True, help = "path/to/index_directory/ (Recommended: leviathan_output/index/ if this will only be used for one project or a centralized location if it will be used for multiple projects)")
-    parser_io.add_argument("-u", "--update_with_genomes", action="store_true",  help = "Update databases with genomes for Sylph sketches")
+    parser_io.add_argument("--update_with_genomes", action="store_true",  help = "Update databases with genomes for Sylph sketches")
+    parser_io.add_argument("--update_salmon_index", action="store_true",  help = "Update the Salmon index on an existing database (e.g., Leviathan databases prior to v2026.7.8 were built with Salmon v1.x and need to be updated with Salmon v2.x)")
 
     # Utilities
     parser_utility = parser.add_argument_group('Utility arguments')
@@ -102,7 +103,7 @@ def main(args=None):
     logger.info(f"Command: {sys.argv}")
 
     # Checks
-    if not opts.no_check:
+    if not opts.no_check and not opts.update_salmon_index:
         genomes_with_filepaths = set()
         with open_file_reader(opts.genomes) as f:
             for line in f:
@@ -144,6 +145,16 @@ def main(args=None):
         if not opts.genomes:
             logger.warning("--genomes not provided but can incoporated post hoc by rerunning with --update_with_genomes")
 
+    if opts.update_salmon_index:
+        if opts.update_with_genomes:
+            msg = "--update_salmon_index and --update_with_genomes are mutually exclusive"
+            logger.critical(msg)
+            parser.error(msg)
+        if not opts.fasta:
+            msg = "--fasta is required when --update_salmon_index is specified"
+            logger.critical(msg)
+            parser.error(msg)
+
     # Threads
     if opts.n_jobs == -1:
         from multiprocessing import cpu_count 
@@ -186,6 +197,7 @@ def main(args=None):
     if all([
         os.path.exists(opts.index_directory),
         not opts.update_with_genomes,
+        not opts.update_salmon_index,
         ]):
         msg = f"--index_directory {opts.index_directory} already exists.  If you want to update with genomes, please use --update_with_genomes or remove directory to overwrite"
         logger.critical(msg)
@@ -195,7 +207,22 @@ def main(args=None):
     os.makedirs(os.path.join(opts.index_directory, "logs"), exist_ok=True)
     os.makedirs(os.path.join(opts.index_directory, "tmp"), exist_ok=True)
 
-    if not opts.update_with_genomes:
+    if opts.update_salmon_index:
+        # ==================
+        # Rebuild Salmon Index
+        # ==================
+        logger.info("Rebuilding salmon index")
+        cmd_salmon_indexer = run_salmon_indexer(
+            logger=logger,
+            log_directory=os.path.join(opts.index_directory, "logs"),
+            salmon_executable=opts.salmon_executable,
+            n_jobs=opts.n_jobs,
+            fasta=opts.fasta,
+            index_directory=opts.index_directory,
+            index_options=opts.salmon_index_options,
+        )
+
+    elif not opts.update_with_genomes:
         # Setup config
         logger.info("Setting up config")
 
@@ -295,7 +322,9 @@ def main(args=None):
     # ==================
     # Build Sylph Sketch
     # ==================
-    if config["contains_genome_filepaths"]:
+    if opts.update_salmon_index:
+        logger.info("Skipping Sylph sketching (only updating salmon index)")
+    elif config["contains_genome_filepaths"]:
         logger.info("Writing genome filepaths for Sylph")
 
         # Write filepaths

bin/leviathan-profile-pathway.py

@@ -32,6 +32,8 @@
     check_salmon_index,
 )
 
+from kegg_pathway_profiler.pathways import profile_pathway_coverage
+
 from leviathan.profile_pathway import(
     check_reads_format,
     run_salmon_quant,
@@ -40,7 +42,6 @@
     build_wide_feature_prevalence_matrix,
     build_feature_prevalence_dictionary,
     build_feature_pathway_dictionary,
-    calculate_pathway_coverage,
     aggregate_pathway_abundance_and_append_coverage,
     aggregate_feature_abundance_for_clusters,
 )
@@ -83,6 +84,7 @@ def main(args=None):
     parser_salmon_quant.add_argument("--salmon_include_mappings", action="store_true", help="salmon quant| Include mappings")
     parser_salmon_quant.add_argument("--salmon_gzip", action="store_true", help="salmon quant | Gzip quant.sf")
 
+    parser_salmon_quant.add_argument("--deterministic", action="store_true", help="salmon quant | Deterministic quantification: byte-identical results across runs and thread counts [Default: False]")
     parser_salmon_quant.add_argument("--salmon_quant_options", type=str, default="", help="salmon quant| More options (e.g. --arg=1 ) https://salmon.readthedocs.io/en/latest/ [Default: '']")
 
     # Samtools
@@ -214,7 +216,8 @@ def main(args=None):
         include_mappings=opts.salmon_include_mappings,
         alignment_format=opts.alignment_format,
         salmon_gzip=opts.salmon_gzip,
-        salmon_quant_options=opts.salmon_quant_options, 
+        salmon_quant_options=opts.salmon_quant_options,
+        deterministic=opts.deterministic,
     )
 
     # ===============================
@@ -277,8 +280,13 @@ def main(args=None):
         logger.info(f"[level={level}] Building feature to pathways dictionary")
         feature_to_pathways = build_feature_pathway_dictionary(pathway_to_data)
         logger.info(f"[level={level}] Calculating pathway coverage")
-        coverages, step_coverages = calculate_pathway_coverage(genome_to_features, pathway_to_data)
-        
+        coverages, step_coverages, _ = profile_pathway_coverage(
+            genome_to_kos=genome_to_features,
+            database=pathway_to_data,
+            n_jobs=opts.n_jobs,
+            serialize_output=False,
+        )
+
         # Pathway abundances
         pathway_abundances_filepath = os.path.join(output_directory, "output", f"pathway_abundances.{level}s.{opts.output_format}")
         if opts.output_format != "parquet":
@@ -354,7 +362,12 @@ def main(args=None):
             logger.info(f"[level={level}] Building feature to pathways dictionary")
             feature_to_pathways = build_feature_pathway_dictionary(pathway_to_data)
             logger.info(f"[level={level}] Calculating pathway coverage")
-            coverages, step_coverages = calculate_pathway_coverage(genome_to_features, pathway_to_data)
+            coverages, step_coverages, _ = profile_pathway_coverage(
+                genome_to_kos=genome_to_features,
+                database=pathway_to_data,
+                n_jobs=opts.n_jobs,
+                serialize_output=False,
+            )
 
             # Pathway abundances
             pathway_abundances_filepath = os.path.join(output_directory, "output", f"pathway_abundances.{level}s.{opts.output_format}")

leviathan/__init__.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python
-__version__ = "2026.3.11"
+__version__ = "2026.7.16"
 from . import utils
 from . import index
 from . import profile_taxonomy
-from . import profile_pathway
-
+from . import profile_pathway