Foldseek driver + Structure base Pangenome

.gitignore

@@ -32,4 +32,7 @@ anvio/data/misc/SCG_TAXONOMY/GTDB/SCG_SEARCH_DATABASES/*.dmnd
 anvio/tests/sandbox/update-contigs-and-bams-for-mini-test/output
 anvio/tests/sandbox/test_visualize_split_coverages/TEST_OUTDIR
 anvio/data/misc/KEGG/
+anvio/data/misc/PROSTT5/
 anvio/data/misc/TRNA_TAXONOMY/
+anvio/data/interactive/node_modules
+anvio/data/interactive/package-lock.json

anvio/__init__.py

@@ -3669,7 +3669,27 @@ def TABULATE(table, header, numalign="right", max_width=0):
                      "a comma-separated list. The default stats are 'detection' and "
                      "'mean_coverage_Q2Q3'. To see a list of available stats, use this flag "
                      "and provide an absolutely ridiculous string after it (we suggest 'cattywampus', but you do you)."}
-    )
+    ),
+    'mode': (
+        ['--mode', '-M'],
+        {'default': None,
+         'metavar': 'structure',
+         'type': str,
+         'help': 'Use this flag to set mode to structure or sequence.'}
+    ),
+   'prostt5-weight-dir': (
+        ['--prostt5-weight-dir'],
+        {'default': None,
+         'type': str,
+         'metavar': 'PATH',
+         'help': "The path for the PROSTT5 Weights to be stored. "
+                 "If you leave it as is without specifying anything, anvi'o will set up everything in "
+                 "a pre-defined default directory. The advantage of using "
+                 "the default directory at the time of set up is that every user of anvi'o on a computer "
+                 "system will be using a single data directory, but then you may need to run the setup "
+                 "program with superuser privileges. If you don't have superuser privileges, then you can "
+                 "use this parameter to tell anvi'o the location you wish to use to setup your weights."}
+        ),
 }
 
 # two functions that works with the dictionary above.

anvio/constants.py

@@ -38,6 +38,8 @@
 
 default_interacdome_data_path = os.path.join(os.path.dirname(anvio.__file__), 'data/misc/Interacdome')
 
+default_prostt5_weight_path = os.path.join(os.path.dirname(anvio.__file__), 'data/misc/PROSTT5/weights')
+
 clustering_configs_dir = os.path.join(os.path.dirname(anvio.__file__), 'data/clusterconfigs')
 clustering_configs = {}
 

anvio/data/interactive/package.json

@@ -36,4 +36,4 @@
     "svg-pan-zoom": "3.6.1",
     "toastr": "2.1.2"
   }
-}
+}

anvio/docs/artifacts/prostt5-model-data.md

@@ -0,0 +1,5 @@
+This artifact stores the data downloaded by %(anvi-setup-prostt5)s and is essential for running %(anvi-pan-genome)s. It includes the ProstT5 model weights, which are required by Foldseek to efficiently perform searches for protein structural similarities.
+
+As detailed in the Foldseek documentation, this data consists of the pre-trained ProstT5 model that accelerates protein structure search tasks. The ProstT5 model is crucial for ensuring accurate results when using Foldseek in your anvi'o workflows.
+
+By default, the ProstT5 model weights are stored in anvio/data/misc/PROSTT5/weights, but users can specify a custom path during setup by using the --prostt5-weight-dir parameter in %(anvi-setup-foldseek)s.

anvio/docs/programs/anvi-setup-prostt5.md

@@ -0,0 +1,25 @@
+
+This program, like other anvi-setup commands, prepares your environment by downloading and configuring the ProstT5 model required for %(anvi-pan-genome --mode structure)s. ProstT5 is essential for Foldseek to perform efficient and accurate searches for protein structural similarities. You only need to run this setup once.
+
+By executing this command, the necessary ProstT5 model weights will be downloaded and stored in the %(foldseek-model-data)s artifact, ensuring that Foldseek can function optimally in your anvi'o workflows.
+
+Setting up the ProstT5 model is simple:
+
+{{ codestart }}
+anvi-setup-prostt5
+{{ codestop }}
+
+When running this program, you can provide a path to store your ProstT5 model in. The default path is `anvio/data/misc/PROSTT5/weights`; if you use a custom path, you will have to provide it to %(anvi-pan-genome)s with the same parameter. Here is an example run:
+
+
+{{ codestart }}
+anvi-setup-prostt5 --prostt5-weight-dir path/to/directory 
+{{ codestop }}
+
+If you want to overwrite any data that you have already downloaded (for example if you suspect something went wrong in the download), add the `--reset` flag: 
+
+{{ codestart }}
+anvi-setup-prostt5  --prostt5-weight-dir path/to/directory \ 
+                        --reset
+{{ codestop }}
+

anvio/drivers/foldseek.py

@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# coding: utf-8
+""" Foldseek Driver ,"""
+
+import os
+import argparse
+import pandas as pd
+import tempfile
+
+import anvio
+import anvio.fastalib as f
+import anvio.utils as utils
+import anvio.terminal as terminal
+import anvio.filesnpaths as filesnpaths
+import anvio.constants as constants
+
+from collections import defaultdict
+from anvio.drivers.mcl import MCL
+from anvio.errors import ConfigError
+from anvio.filesnpaths import AppendableFile
+
+
+__copyright__ = "Copyleft 2015-2024, The Anvi'o Project (http://anvio.org/)"
+__credits__ = []
+__license__ = "GPL 3.0"
+__version__ = anvio.__version__
+__maintainer__ = "Metehan Sever"
+__email__ = "[email protected]"
+
+
+run = terminal.Run()
+progress = terminal.Progress()
+pp = terminal.pretty_print
+
+class Foldseek():
+
+    def __init__(self, query_fasta=None, run=run, progress=progress, num_threads=1, weight_dir=None, overwrite_output_destinations=False):
+        self.run = run
+        self.progress = progress
+
+        utils.is_program_exists('foldseek')
+
+        self.query_fasta = query_fasta
+        self.num_threads = num_threads
+        self.weight_dir = weight_dir or constants.default_prostt5_weight_path
+        self.overwrite_output_destinations = overwrite_output_destinations
+        self.tmp_dir = tempfile.gettempdir()
+
+        filesnpaths.is_file_exists(self.weight_dir)
+
+        self.output_file = 'foldseek-search-results'
+
+        if not self.run.log_file_path:
+            self.run.log_file_path = filesnpaths.get_temp_file_path()
+
+        self.names_dict = None
+
+    def create_db(self):
+        self.run.warning(None, header="FOLDSEEK CREATEDB", lc="green")
+        self.progress.new('FOLDSEEK')
+        self.progress.update('creating the search database (using %d thread(s)) ...' % self.num_threads)
+
+        expected_output_dir = os.path.join(self.output_file, "db")
+        expected_output_file = os.path.join(expected_output_dir, "search_db")
+
+        filesnpaths.gen_output_directory(expected_output_dir, delete_if_exists=False)
+
+        cmd_line = ['foldseek',
+                    'createdb',
+                    self.query_fasta,
+                    expected_output_file,
+                    '--prostt5-model', self.weight_dir,
+                    '--threads', self.num_threads
+                    ]
+
+        utils.run_command(cmd_line, self.run.log_file_path)
+
+        self.progress.end()
+        self.run.info('Command line', ' '.join([str(x) for x in cmd_line]), quiet=True)
+        self.run.info('Foldseek search DB', expected_output_file)
+
+    def search(self, query_db, target_db):
+        self.run.warning(None, header="FOLDSEEK EASY SEARCH", lc="green")
+        self.progress.new('FOLDSEEK')
+        self.progress.update('Running search using Foldseek ...')
+
+        query_db = os.path.join(query_db, 'search_db')
+        target_db = os.path.join(target_db, 'search_db')
+
+        result_file_dir = os.path.join(self.output_file, 'result')
+
+        cmd_line = [
+            'foldseek',
+            'easy-search',
+            query_db,
+            target_db,
+            result_file_dir,
+            self.tmp_dir,
+            '--threads', self.num_threads
+        ]
+
+        utils.run_command(cmd_line, self.run.log_file_path)
+
+        self.progress.end()
+
+        self.run.info('Command line', ' '.join([str(x) for x in cmd_line]), quiet=True)
+        self.run.info('Foldseek search Result', result_file_dir)
+
+    def process(self, query_db, target_db):
+
+        self.create_db()
+        self.search(query_db, target_db)
+
+    def get_foldseek_results(self):
+        """ Return result.m8 file """
+        force_makedb, force_search = False, False
+
+        result_dir = os.path.join(self.output_file, 'result')
+
+        return result_dir

anvio/panops.py

@@ -34,6 +34,7 @@
 from anvio.drivers.diamond import Diamond
 from anvio.drivers.mcl import MCL
 from anvio.drivers import Aligners
+from anvio.drivers.foldseek import Foldseek
 
 from anvio.errors import ConfigError, FilesNPathsError
 from anvio.genomestorage import GenomeStorage
@@ -89,12 +90,18 @@ def __init__(self, args=None, run=run, progress=progress):
         self.enforce_hierarchical_clustering = A('enforce_hierarchical_clustering')
         self.enforce_the_analysis_of_excessive_number_of_genomes = anvio.USER_KNOWS_IT_IS_NOT_A_GOOD_IDEA
 
+        self.de_novo_compute_mode = A('mode') or 'sequence'
+        self.prostt5_weight_dir = A('prostt5_weight_dir')
+
         self.additional_params_for_seq_search = A('additional_params_for_seq_search')
         self.additional_params_for_seq_search_processed = False
 
         if not self.project_name:
             raise ConfigError("Please set a project name using --project-name or -n.")
 
+        if self.de_novo_compute_mode == 'structure':
+            self.skip_alignments = True
+
         # when it is time to organize gene_clusters
         self.linkage = A('linkage') or constants.linkage_method_default
         self.distance = A('distance') or constants.distance_metric_default
@@ -245,7 +252,12 @@ def check_params(self):
             filesnpaths.is_file_plain_text(self.description_file_path)
             self.description = open(os.path.abspath(self.description_file_path), 'r').read()
 
-        self.pan_db_path = self.get_output_file_path(self.project_name + '-PAN.db')
+        if self.de_novo_compute_mode == "sequence" or self.user_defined_gene_clusters:
+            self.pan_db_path = self.get_output_file_path(self.project_name + '-PAN.db')
+        elif self.de_novo_compute_mode == "structure":
+            self.pan_db_path = self.get_output_file_path(self.project_name + '-STRUCTURE-PAN.db')
+        else:
+            raise ConfigError("Something is wrong")
 
 
     def process_additional_params(self):
@@ -309,12 +321,33 @@ def run_blast(self, unique_AA_sequences_fasta_path, unique_AA_sequences_names_di
         return blast.get_blast_results()
 
 
+    def run_foldseek(self, unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict):
+        """ Running Foldseek """
+        result_dir = self.get_output_file_path('foldseek-search-results')
+
+        fs = Foldseek(query_fasta=unique_AA_sequences_fasta_path, run=self.run, progress=self.progress,
+                                num_threads=self.num_threads, weight_dir=self.prostt5_weight_dir, overwrite_output_destinations=self.overwrite_output_destinations)
+
+        # It may help downstream analysis in the future but for now have no functionality :/ 
+        fs.names_dict = unique_AA_sequences_names_dict
+        fs.output_file = result_dir
+        fs.log_file_path = self.log_file_path
+
+        db_dir = os.path.join(result_dir, 'db')
+        # FIXME this tmp is necessarry for easy-search instead of giving like that we can set default tempfile.gettempdir() in foldseek driver
+
+        fs.process(db_dir, db_dir)
+
+        return fs.get_foldseek_results()
+
     def run_search(self, unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict):
-        if self.use_ncbi_blast:
-            return self.run_blast(unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict)
+        if not self.de_novo_compute_mode == 'structure':
+            if self.use_ncbi_blast:
+                return self.run_blast(unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict)
+            else:
+                return self.run_diamond(unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict)
         else:
-            return self.run_diamond(unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict)
-
+            return self.run_foldseek(unique_AA_sequences_fasta_path, unique_AA_sequences_names_dict)
 
     def run_mcl(self, mcl_input_file_path):
         mcl = MCL(mcl_input_file_path, run=self.run, progress=self.progress, num_threads=self.num_threads)