A vector database solution for rational design of CRISPR defense-avoidant phage therapy

A baseline approach to select CRISPR-resistant phages with local BLAST

This section outlines a straightforward approach for selecting phages based on local alignment with CRISPR spacers using BLAST (Altschul et al., 1990). It provides a starting point for researchers aiming to reduce immune responses between phages and bacterial hosts with CRISPR systems. A schematic representation of the pipeline is shown in Fig. 1.

To develop a baseline approach for selecting CRISPR-resistant phages, we employed a combination of bioinformatics tools and techniques. The pipeline involves constructing local BLAST databases from phage genomes and querying them with the CRISPR spacers extracted from a reference genome of the bacterium we are interested in (or CRISPR spacers extracted from a sequenced clinical or metagenomic sample). This enables us to perform targeted searches against the protospacers to identify potential CRISPR-resistant phages with high confidence. This method was chosen for its simplicity and efficiency in comparing DNA sequences, as well as its omnipresence in bioinformatics.

In this section, we describe the detailed methodology for each stage of our approach, including database construction, BLAST search execution, and result analysis. We also discuss the various parameters and considerations that need to be taken into account to ensure accurate and reliable results.

The foundation of our baseline approach is the creation of local BLAST databases containing reference genomes of phages. To achieve this, we first have to obtain genomic sequences from publicly available databases, such as NCBI’s GenBank, PhageDB or INPHARED (Cook et al., 2021). After obtaining the phage sequences, we have to use a string search tool (e.g., CasOffFinder; Bae, Park & Kim, 2014) to find all possible protospacers in the phages. This step is required to ensure that we are comparing spacers with potential protospacers instead of an arbitrary collection of substrings from the phage genomes. A protospacer here is defined as a substring that is located nearby Protospacer-Adjacent Motif—a special set of sequences that CRISPR-associated protein uses to recognize the target DNA sequence (Mojica et al., 2009). The exact length and location of the protospacer in relation to PAM varies highly between different CRISPR systems (Shah et al., 2013), so the exact description of the protospacer-PAM combination is left for the end user, e.g., YCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN for Clostridium difficile where YCN is the PAM (Maikova et al., 2021) and the rest is a mask for the protospacer string. In the latter sections, we provide a case study of designing a phage cocktail targeted to evade a well studied CRISPR-Cas system of C. difficile.

To look for protospacer+PAM combinations, we use CasOffFinder (Bae, Park & Kim, 2014) with an all-permissive string mask. CasOffFinder enables GPU-accelerated search for strings in large genomes. It was initially developed to search for off-target editing sites, but with an all-permissive string mask and allowed number of mismatches equal to the length of protospacer, it can be used to select all possible pairs of protospacer and PAM. After gathering all potential protospacers, we employ makeblastdb tool from BLAST+ suite (Camacho et al., 2009) to create local BLAST databases. We opt for the nucleotide-nucleotide search mode to compare the query spacer sequences against the protospacers from our database. These protospacers represent potential targets for CRISPR defense.

Once the local BLAST databases are constructed, we can proceed with the identification of CRISPR-resistant phage candidates. To achieve this, we first have to extract the spacer sequences from the bacteria of interest. This output provides critical information about the CRISPR-Cas elements within the input sequences, which are essential for understanding potential targets for phage therapy. We can do the extraction using CRISPR-Cas systems annotation tools like CRISPRCasTyper (Russel et al., 2020), CRISPRloci (Alkhnbashi et al., 2021) or CRISPRCasFinder (Couvin et al., 2018) (which we decided to use based on its widespread adoption in the field of CRISPR-Cas research). We choose to run CRISPRCasFinder via the provided Singularity container. While CRISPR-Cas annotation tools provide information beyond CRISPR array content (e.g., location of Cas proteins, type of CRISPR-Cas system and so on), we mainly use only the spacer sequences.

We then need to perform BLAST searches against the local phage databases using the blastn tool from the BLAST+ suite. The choice of parameters is critical to achieving accurate and reliable results. Our BLAST search is unusual—while in a typical search, we would like to look for the most similar sequences to the query (Altschul et al., 1994), here we are looking for the sequences that are not similar since we are interested in the phages that do not contain targetable protospacers. We want to get as many sequences as possiblie to reduce the probability of selecting the wrong phages. Therefore, our parameter choice is as follows—evalue is 10⁵ (unlike typical BLAST search that would use evalues like ${10}^{-10}$), max_target_seq is 10⁶ (while BLAST default is 500). We also use a word size of 11 to balance sensitivity and efficiency.

From that, the selection of phage cocktail components is a simple rearrangement of the BLAST search result according to their similarity to the query. Let $PS(n,i,q)$ be the i-th protospacer (of I protospacers possible) from the n-th phage (of N total phages) that corresponds to the query spacer $q$ (of Q total queries). Such a protospacer-query pair is represented by the alignment $A(PS(n,i,q))$ which can be characterized by the length of the largest aligned portion, $L(A(PS(n,i,q)))$. In order to find the most dissimilar phage, we have to solve the following optimization problem:

(1) $$N_{phage}=\underset{n}{argmin}(\mathrm{\forall }n\in N,i\in I,q\in Q:maxL(A(PS(n,i,q))).$$

We are interested in a phage (or phages) that has the smallest maximal alignment length among its protospacer-query pairs. By sorting and selecting the top $m$ phages with the smallest maximum alignment lengths, we can identify a promising composition for a CRISPR-resistant phage cocktail of size $m$.

Vector database with hierarchical navigable small world structure

To address the shortcomings of the baseline BLAST approach, we introduce VDPhage—a novel vector-based database solution designed with a HNSW structure. The core concept behind this approach is to transform protospacer sequences into high-dimensional vectors while preserving their sequential similarity and then navigate this vector space efficiently following the small world structure. A schematic representation of the VDPhage pipeline is shown in Fig. 2. VDPhage pipeline follows closely the baseline described in previous section, only substitutes the local BLAST database for vector and ordering by maximal alignment length for ordering by minimal cosine distance.

The incorporation of a layered navigable structure further enhances the database’s functionality by enabling efficient navigation through interconnected nodes. By leveraging these structural features, the vector database offers significant advantages over traditional database systems like BLAST.

The HNSW graph (Malkov & Yashunin, 2018) is an extension of k-nearest neighbors (k-NN) proximity graphs. A k-NN proximity graph is an undirected graph weighted by distance between its nodes (e.g., alignment score or other forms of distance between sequences) where the only links that exist are the connections between the closest nodes. Search in an exact k-nearest neighbor proximity graph involves greedy traversal over the links starting from a (generally) random node. Such a search scales linearly with the size of the graph, therefore for performance’s sake, various approximations to the exact k-NN algorithm had been developed.

HNSW introduces two main improvements to the exact k-NN search—navigability and hierarchy. Navigability follows from the Kleinberg’s small world model (Kleinberg, 2000), which itself is an extension of Watts-Strogatz graph (Watts & Strogatz, 1998). In Kleinberg’s model, there are two types of connections between nodes—short connections that link closest neighbors and long connections that link nodes at random with probability inversely proportional to their weight (distance between sequences). Addition of long connections allows for $O((logn{)}^2)$ complexity of search (Kleinberg, 2000).

Hierarchy in HNSW involves separation of the navigable proximity graph into layers according to their distance. The greedy search commences starting from the top layer to bottom layer, and on each iteration only a small portion of connected nodes are thus traversed. This allows for reduction of search complexity from $O((logn{)}^2)$ to $O(logn)$ (Malkov & Yashunin, 2018).

For our HNSW engine, we use Usearch library (Vardanian, 2023). Usearch is an advanced, open-source vector database solution designed for efficient searching and clustering within diverse applications involving vectors and strings. Utilizing high-performance algorithms such as HNSW, Usearch offers customizable metrics support and low-latency performance through native language bindings. Its compact codebase and minimal dependencies make it a lightweight alternative to other tools like FAISS (Douze et al., 2024). With built-in capabilities for clustering, dimensionality reduction, and multi-language support, Usearch excels in scalability and deployment ease. To work with Usearch, we use the vicinity (van Dongen et al., 2024) library that provides a simple and uniform interface for different vector database engines.

For the construction of the vector database, the sequences undergo feature extraction, which involves transforming the raw genetic information into high-dimensional vectors. This transformation could be in principle achieved using advanced machine learning techniques, such as word embeddings or sequence-to-vector models like dna2vec (Ng, 2017), kmer2vec (Ren, Yin & Yau, 2022) or kmer-node2vec (Yu et al., 2022); however, we strive for simplicity and use the following simple transformation (Algorithm 1).

This is a classical method (consult Moeckel et al., 2024) for a thorough modern review of k-mer methods, we use the Rust implementation of Weber (2024)) of alignment-free sequence comparison and it fits perfectly within the framework of vector databases. We compute k-mer frequencies and use them as our embeddings for the vector database. This allows us to introduce cosine distance between two embeddings as our distance between spacer and protospacer.

(2) $$cos(S,P)=1-\frac{Algorithm1(S)\cdot Algorithm1(P)}{||Algorithm1(P)||||Algorithm1(P)||},$$where S and P are spacer and protospacer strings respectively, $\cdot$ represents dot product and ||X|| represents the magnitude of the vector.

The next stage involves the creation of the hierarchical structure of the database. This hierarchy organizes the vectors into clusters based on their similarity, facilitating efficient data retrieval during query processing. The navigable small-world property allows for faster traversal through the network, enabling users to quickly access related sequences and associated labels. This results in a robust and efficient database system that ensures that each genetic sequence is accurately represented and organized in a manner that maximizes retrieval speed and usability. Construction of the HNSW is $O(nlogn)$ because it is constructed by a consecutive search of $n$ elements of the data (Malkov & Yashunin, 2018).

When a user submits a query to the vector database, the system utilizes its hierarchical structure to quickly narrow down potential matches. The small-world organization allows the database to efficiently partition the high-dimensional vector space into clusters, enabling the search algorithm to focus on relevant regions of the space rather than scanning the entire dataset. This drastically reduces the complexity of the search process and accelerates query execution times. With the small-world network, users can navigate directly to related sequences through interconnected nodes, significantly reducing query response times and improving overall usability, unlike the traditional relational database that requires navigating of a series of linear lists that result from specific queries.

A design for comparison study between the baseline and VDPhage

We test our approach for CRISPR-resistant phage cocktail design on a simulated scenario of C. difficile infection. We use reference genome of C. difficile strain S-0253 (NCBI, 2025a) as a source of CRISPR spacers and a selection of Clostridioides phages from NCBI virus database (NCBI, 2025b) as a source of protospacers (70 phages total). Protospacer-Adjacent Motif for CRISPR system of C. difficile is “YCN” (Maikova et al., 2021). In the worst-case scenario, one might expect that the majority of spacers from the native CRISPR-Cas system of the bacteria would exhibit high similarity to protospacers in the vector or BLAST database, leading to potential CRISPR interference and reduced therapeutic efficacy. To assess our approach against this scenario, we have to evaluate the distribution of sequence similarities for best protospacer-spacer alignments. To compute such alignments, we use local BLAST database built on the protospacer sequences found by CasOffFinder with rather standard parameters (e-value of ${10}^{-6}$ and word size of 11). These alignments would also give us the worst possible selection of phages—a set of phages, all of which do indeed contain protospacers highly similar to ones in the CRISPR array.

This comparison should reveal that our proposed vector database solution effectively reduces the likelihood of selecting phages that have protospacers with high similarity to known spacers, thereby minimizing CRISPR interference. An evidence for that would be a significant difference in distributions of sequence similarities from the best alignments and from the protospacer-spacer pairs that motivate the phage’s ranking in VDPhage or baseline outputs.

To measure sequence similarity, we use Levenshtein distance—a minimal number of insertions, deletions and substitutions required to turn one string into another (Levenshtein, 1966). To test whether two samples (e.g., best alignments and VDPhage) came from the same distribution, we use Kolmogorov–Smirnov statistics (Smirnov, 1944). We select p-value of $\frac{0.05}{N}$ as our threshold for rejecting the null hypothesis (two samples came out of the same distribution, N is the number of comparisons, as the Bonferroni correction for multiple comparisons suggests). In addition to that, we visualize the diverging distributions with violin plots.

The top $m$ phages that constitute the proposed phage cocktail of size $m$ should not overlap with the worst possible phages. We can check it by computing set intersections of phage lists and visualize it with a Venn diagram.

Performance evaluation

This section describes the evaluation of two proposed algorithms to select CRISPR-resistant phage cocktails. The comparison between two algorithms is based on a simulated scenario of C. difficile infection described in ‘A design for comparison study between the baseline and VDPhage’ of the ‘Materials and Methods’ section. We compare the quality of resulting phage cocktails measured in deviation of sequence similarity distributions from the sequence similarity distribution of best available protospacer-spacer alignments.

In the case of phage therapy design, our objective is to identify a set of phage genomes that could evade the CRISPR-Cas system of a specific bacterial species, thereby offering an effective treatment strategy. To do so, we have to conduct a sequence similarity search using a query representing the CRISPR array of the bacterial species. The HNSW structure of the vector database enables rapid identification and retrieval of phage genomes that share high similarity with the query. The results of Levenshtein distance distribution comparison is shown in Fig. 3.

All of the selected phages distributions differ highly from the distribution for best available alignments. Cocktails suggested by baseline and VDPhage overlap significantly (six out of 10 for top 10 cocktail), but neither of these methods select any of the worst possible phages. Distributions for BLAST have more variance in Levenshtein distance than the ones for VDPhage which suggests that VDPhage is more robust and consistent at locating the global minimum of sequence similarity. Considering the statistical significance of differences between the distributions (best alignments, VDPhage all and BLAST all), Kolmogorov-Smirnov test (corrected for three comparisons by Bonferroni correction) gives us the following p-values: $0.0427$ for BLAST vs VDPhage, $0$ for VDPhage vs best alignments and $0$ for BLAST vs best alignments. Since our threshold for statictical significance is $\frac{0.05}{3}\approx 0.0167$, we consider significant only differences between best alignments and the rest. In the case of BLAST and VDPhage, we do not have enough evidence to reject the null hypothesis of the distributions being the same.

Runtime evaluation

The runtime evaluation of our proposed vector database solution for CRISPR defense avoidant phage therapy was conducted to assess its efficiency, scalability, and accuracy compared to baseline methods (local BLAST). We have performed 100 of each operations (VDPhage build, VDPhage query, baseline build, baseline query) and computed mean computation time and its standard deviation. Building the baseline database runs much faster than vector database – $3.87\pm 0.17$ s vs $22.33\pm 0.54$ s respectively. However, the querying is much faster for VDPhage – $11.53\pm 0.38$ for $\text{770,000}$ extracted records vs $13.67\pm 0.25$ for $\text{35,689}$ extracted records, which means on average $\text{66,782.31}$ records per second vs $\text{2,610.75}$ records per second. VDPhage extracts a single record almost 26 times faster than local BLAST, which is expected because of logarithmic complexity of the search in HNSW. All experiments were performed on an Alienware 17 R4 laptop with 24Gb of RAM and 2.9 GHz Intel Core i7-7820HK Quad-Core CPU.

Overview of VDPhage CLI application

The VDPhage software is available as a command line interface tool. It supports both vector database and local BLAST methods, outputting a sorted list of phages and the accompanying information on its decisions. The script relies on several external libraries, including “os”, “sys”, “json”, “vicinity”, and “argparse”, which facilitate tasks such as parsing command-line arguments, reading configuration files, handling file input/output operations, executing bioinformatics analyses (e.g., BLAST queries), and managing phage data. VDPhage, as well as the local BLAST baseline, provides two modes of operation:

• vdphage build—to construct the vector database (or BLAST database) from a FASTA file of protospacers. The build command constructs a searchable database containing information about protospacers, spacer sequences, and their respective alignments. Additionally, the user has to specify the k-mer length. In the code, build mode is enabled by build_baseline and build_vdphage functions that read the configuration file, create necessary output directories, copy the protospacers file to the output directory, build a BLAST or VDPhage database and save metadata about the build;

• vdphage query—to query the database for spacers from bacterial genome (in FASTA file). The query command searches the database for potential protospacer-spacer interactions based on input spacer sequences. Additionally, the user has to specify the number of protospacer-spacer pairs to return (default is 10,000) The query mode is represented in the code with query_baseline and query_vdphage functions that read the configuration from a JSON file, create a temporary directory if it doesn’t exist, perform a BLASTN search or a vector database query with the new spacers against the respective database, parse the results, and save a summary as a TSV file.

The tool can be configured using a JSON configuration file config.json. The configuration file allows users to customize paths to temporary directories and external tools (e.g., k-mer counter), parameters like makeblastdb and blastn options, and other settings specific to the chosen method. The configuration file has the following structure:

1. TEMPDIR—a path to temporary directory (/tmp/vdphage/ by default);

2. PROTOSPACER KMER COUNTS—a name of file with k-mer counts (protospacers.kmer_counts.npy by default);

3. PROTOSPACER KMER COUNT IDS—a name of file with the corresponding identification information (protospacers.kmer_counts.ids.txt by default);

4. COUNT KMERS—a command that runs kmer-counter (the default one is ../kmer-counter/target/release/kmer-counter –collapse 0 –klength KKKK –file INPUTFILE –ids OUTIDS –out OUTFILE). COUNT KMERS takes INPUTFILE with the protospacers or spacers as input and outputs PROTOSPACER KMER COUNTS and PROTOSPACER KMER COUNT IDS;

5. MAKEBLASTDB—a command to run makeblastdb (the default one is makeblastdb -in INPUTFILE -title “Phages local BLAST” -dbtype nucl). MAKEBLASTDB takes the INPUTFILE with protospacers as the input and creates a local BLAST database within the same directory as the input file;

6. BLAST OUT—a filename of BLAST query result (BLAST_QUERY_RESULT.xml by default);

7. BLASTN—a command to run BLASTN (the default one is blastn -query INPUTFILE -db DATABASE -evalue 10000 -num_threads 1 -out OUTFILE -outfmt 5 -max_target_seqs 1000000 -word_size 11). BLASTN takes INPUTFILE with the spacers and the database name DATABASE as the input and outputs the query result into BLAST OUT.

The example commands are as follows:

The functionality of vicinity library for NHSW computations is encapsulated in a VicinityWrapper class that exposes __init__ and query methods. The initialization method __init__ creates the VicinityWrapper object with protospacers, spacers, and a database (V), loads the spacers from the provided list spacers, loads the protospacers from the provided list protospacers, which includes both protospacers and associated phages, stores the loaded database in self.db attribute, stores the loaded spacers in self.spacers attribute after loading them using a custom function load_spacers(), extracts unique phages from the protospacers with np.unique() and stores them in self.phages attribute after converting them to a numpy array.

The query method queries the database with a given vector qv and returns a summary dataframe. It performs a query on the loaded database using self.db.query() method with the provided query vector qv and keyword argument k. Then it concatenate dataframes for each spacer in the spacers dictionary. Each dataframe contains information about spacers, protospacers, distances, queries, and hits. After that a new column phage is added to the dataframe by splitting the protospacer column and extracting the phage ID. Then any phages that are not present in the queried results are identified and added to the dataframe with default values. The modified dataframe is then concatenated with the additional rows for absent phages. After that, a summary dataframe with minimum distances, spacers, protospacers, queries, hits, and phage IDs grouped by phage is created and sorted in descending order based on the distance column to prioritize longer minimal distances. Then the query method calculates the Levenshtein distance between query strings and their corresponding hit strings and adds it as a new column levenshtein to the summary dataframe. The resulting dataframe is then returned from the method and saved to the output file.

The example phage cocktail of size 5 (top five results) selected by VDPhage is shown in the Table 1. Here “phage” column denotes the accession number of a phage, “spacer” column denotes the spacer with the minimal cosine distance of a spacer-protospacer pair matched to that phage. Each spacer is represented as the accession number of the bacterial genome followed by chromosome, start and end position. The cosine distance itself is shown in “distance” column and “levenshtein” column shows the Levenshtein distance between spacer and protospacer in the pair with minimal cosine distance.