Methods and systems for discovery of embedded target genes in biosynthetic gene clusters

Abstract

The present disclosure relates to computer-based methods and systems for using comparative genomics techniques and machine learning models to identify genes associated with biosynthetic gene clusters (BGCs), including embedded target genes (ETaGs) that are homologs of potential therapeutic targets.

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of priority to U.S. Provisional Patent Application No. 63 / 263,638, filed November 5, 2021, and to U.S. Provisional Patent Application No. 63 / 278,065, filed November 10, 2021, the contents of each of which are incorporated by reference herein in their entireties.

[0002] The present disclosure relates generally to methods and systems for identifying genes associated with a gene cluster (e.g., a biosynthetic gene cluster), and applications thereof, including methods for determining boundaries of a gene cluster (e.g., boundaries of a biosynthetic gene cluster), methods for identifying therapeutic targets, and methods for drug discovery. [Background technology]

[0003] Microorganisms produce a wide variety of small molecule compounds known as secondary metabolites or natural products with diverse chemical structures and functions. Some secondary metabolites enable microorganisms to tolerate hostile environments, while others serve as weapons of inter- and intra-species competition. See, e.g., Piel, J. Nat. Prod. Rep., 26:338-362, 2009. Many human pharmaceuticals (including, e.g., antibacterial agents, antitumor agents, and insecticides) are derived from secondary metabolites. See, e.g., Newman DJ and Cragg GM, J. Nat. Prod., 79:629-661, 2016.

[0004] Microorganisms synthesize secondary metabolites using enzymatic proteins encoded by clusters of co-localized genes called biosynthetic gene clusters (BGCs). Evidence is emerging that some microbial biosynthetic gene clusters contain genes that do not appear to be involved in the synthesis of the associated biosynthetic products produced by the enzymes encoded by the clusters. In some cases, such non-biosynthetic genes have been described as "self-protecting" because they code for proteins that can apparently render the host organism resistant to the associated biosynthetic product. For example, in some cases, resistance mutants of non-biosynthetic genes that code for transporters of biosynthetic products, detoxifying enzymes that act on biosynthetic products, or proteins whose activity is targeted by biosynthetic products have been reported. See, for example, Cimermancic et al., Cell 158:412, 2014; Keller, Nat. Chem. Biol. 11:671, 2015. Researchers have proposed that the identification of such genes and the determination of their functions may be useful in determining the role of the biosynthetic products synthesized by the enzymes of the clusters. See, e.g., Yeh et al., ACS Chem. Biol. 11:2275, 2016; Tang et al., ACS Chem. Biol. 10:2841, 2015; Regueira et al., Appl. Environ. Microbiol. 77:3035, 2011; Kennedy et al., Science 284:1368, 1999; Lowther et al., Proc. Natl. Acad. Sci. USA 95:12153, 1998; Abe et al., Mol. Genet. Genomics 268:130, 2002. US Patent Application Publication No. 2020 / 0211673 provides the insight that certain non-biosynthetic genes present in biosynthetic gene clusters or in close proximity to biosynthetic genes in clusters (particularly eukaryotic, e.g., fungal, biosynthetic gene clusters, as opposed to bacterial, biosynthetic gene clusters) may represent homologs of human genes that are targets for therapeutic purposes. Such non-biosynthetic genes are referred to as "embedded target genes" or "ETaGs."

[0005] Traditionally, secondary metabolites have been identified from microbial cultures and screened for therapeutic activity against human targets of interest. However, most microorganisms are not culturable, and even BGCs in culturable microorganisms may remain transcriptionally silent under laboratory conditions. Recent developments in nucleic acid and protein sequencing technologies as well as bioinformatics pipelines have made it possible to rapidly identify large numbers of BGCs from environmental microorganisms without the need to culture the microorganism and test the biological activity of the BGC. See, for example, Palazzotto E. and Weber T, Curr. Opin. Microbiol., 45:109-116, 2018. However, it remains a challenge to precisely define the genomic boundaries of BGCs using purely computational methods. There is also no computational pipeline available to identify genes embedded in BGCs that confer self-protection against secondary metabolites produced by the BGC. Summary of the Invention

[0006] Exemplary methods and systems are disclosed herein for identifying genes associated with biosynthetic gene clusters (BGCs) in a target genome. The disclosed methods and systems provide genome database search and analysis tools that can be used to identify putative embedded target gene sequences (pETaGs) in one or more target genomes that are homologs of known ETaG sequences in a first genome. Methods and systems are also disclosed for evaluating the likelihood that a given pETaG is an actual ETaG (e.g., a resistance gene for secondary metabolites produced by BGCs). The disclosed methods and systems can also be used to determine the boundaries of BGCs in one or more target genomes, or to aid in the identification of small molecule modulators of target genes of interest.

[0007] Disclosed herein is a computer-implemented method for identifying embedded target genes (ETaGs), the method comprising: specifying one or more query sequences or proxies thereof; selecting one or more target genomes; performing a search of the one or more target genomes using the one or more query sequences or proxies thereof to identify putative embedded target gene (pETaG) sequences that are homologs of the one or more query sequences based on a comparison of one or more homology-based metrics for candidate pETaGs to one or more predetermined homology-based metric thresholds; and determining whether a given pETaG is an ETaG based on a comparative genomics analysis of multiple genomes.

[0008] In some embodiments, the comparative genomics analysis comprises generating a comparative genomics heat map based on the plurality of genomes. In some embodiments, the plurality of genomes comprises a plurality of positive genomes and a plurality of negative genomes. In some embodiments, the comparative genomics analysis comprises determining phylogenetic characteristics, co-occurrence characteristics, co-evolution characteristics, or any combination thereof for a given pETaG based on the plurality of genomes.

[0009] In some embodiments, the comparative genomics analysis involves analysis of an input dataset that includes phylogenetic features, co-occurrence features, co-evolution features, comparative genomics heat maps, data derived from comparative genomics heat maps, or any combination thereof for pETaG using a machine learning model or an empirical algorithm to predict the probability that pETaG is an ETaG.

[0010] In some embodiments, the computer-implemented method further includes determining that the identified pETaG is associated with the resistance mechanism based on determining a copy number of the identified pETaG. In some embodiments, the computer-implemented method further includes determining that the pETaG is associated with the resistance mechanism based on determining a copy number difference between a positive genome containing pETaG and a negative genome not containing pETaG.

[0011] In some embodiments, the one or more query sequences, or proxies thereof, comprise one or more protein sequences, one or more nucleic acid sequences, one or more Universal Protein Resource (Uniprot) identification numbers, one or more profile hidden Markov models (pHMMs), a specified set of protein sequence domains, or any combination thereof. In some embodiments, the one or more query sequences, or proxies thereof, are selected from a bacterial genome, an archaeal genome, a fungal genome, a plant genome, an animal genome, a human genome, or any combination thereof.

[0012] In some embodiments, one or more target genomes are selected from bacterial genomes, fungal genomes, plant genomes, or any combination thereof. In some embodiments, two or more target genomes are selected based on pairwise similarity scores, pairwise phylogenetic distances, or any combination thereof.

[0013] In some embodiments, the computer-implemented method further comprises filtering the two or more selected target genomes to retain only those target genomes whose (i) pairwise similarity scores are greater than a specified pairwise similarity threshold, or (ii) pairwise phylogenetic distances are less than a specified phylogenetic distance threshold.

[0014] In some embodiments, the computer-implemented method further comprises clustering the retained target genomes into a set using a clustering algorithm, and performing a search using one or more of the sets of clustered target genomes. In some embodiments, the clustering algorithm comprises a Markov cluster algorithm.

[0015] In some embodiments, the search is performed using BLAST, DIAMOND, HMMER, Exonerate, or ggsearch. In some embodiments, the search is limited to one or more specific regions of the one or more target genomes. In some embodiments, the one or more specific regions include one or more biosynthetic gene clusters (BGCs). In some embodiments, the one or more BGCs in the one or more target genomes are predicted using a BGC search algorithm. In some embodiments, the BGC search algorithm includes antiSMASH, SMURF, TOUCAN, or deepBGC.

[0016] In some embodiments, one or more BGCs are predicted for one or more target genomes by extracting sequence regions of a specified length proximal to gene sequences that match known biosynthetic core synthases determined using a sequence search tool. In some embodiments, the sequence search tool includes BLAST, DIAMOND, HMMER, Exonerate, or ggsearch. In some embodiments, one or more BGCs are predicted for one or more query genomes using a hidden Markov model (HMM) of known core synthases. In some embodiments, one or more BGCs are predicted for one or more target genomes based on the co-localization of protein sequence domains associated with known core synthases.

[0017] In some embodiments, the one or more homologous sequence-based metrics include a sequence identity percentage, a sequence coverage percentage, an E-value, a bit score, an HMM score, or any combination thereof. In some embodiments, the one or more predetermined homologous sequence-based metric thresholds include a sequence identity percentage threshold, a sequence coverage percentage threshold, an E-value threshold, a bit score threshold, an HMM score threshold, or any combination thereof. In some embodiments, the one or more predetermined homologous sequence-based metric thresholds include a sequence identity percentage threshold having a value of at least 20%, at least 30%, at least 40%, at least 50%, at least 60%, at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least 98%. In some embodiments, the one or more predetermined homologous sequence-based metric thresholds include a sequence coverage percentage threshold having a value of at least 20%, at least 30%, at least 40%, at least 50%, at least 60%, at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least 98%. ...10, at least 9, at least 8, at least 7, at least 6, at least 5, at least 4, at least 3, at least 2, at least 1, at least 0.01, at least 0.001, at least 1e -10 Less than 1e -20 Less than 1e -30 Less than 1e -40 Less than 1e -50 Less than 1e -60 Less than 1e -70 Less than 1e -80 Less than 1e -90 Less than or equal to 1e -100In some embodiments, the one or more predetermined homologous sequence-based metric thresholds include a bit score threshold having a value of at least 40, at least 50, at least 60, at least 70, at least 80, at least 90, at least 100, at least 250, at least 500, at least 1000, or at least 5000. In some embodiments, the one or more predetermined homologous sequence-based metric thresholds include an HMM score threshold having a value of at least 10, at least 25, at least 50, at least 100, at least 250, at least 500, at least 1000, or at least 5000.

[0018] In some embodiments, performing the search includes converting one or more query sequences, including protein sequences, into nucleic acid sequences, and performing a search of one or more target genomes using the one or more query sequences converted into nucleic acid sequences to identify homologous nucleic acid sequences based on a comparison of one or more homology-based metrics of the candidate pETaGs to one or more predetermined homology-based metric thresholds, and comparing the genomic coordinates of the homologous nucleic acid sequences to the genomic coordinates corresponding to the predicted protein sequences in the one or more target genomes. In some embodiments, if a homologous nucleic acid sequence overlaps with a nucleic acid sequence corresponding to a single predicted protein sequence and the overlap is greater than a specified nucleic acid sequence overlap threshold, the predicted protein sequence is reported as pETaG. In some embodiments, if a homologous nucleic acid sequence overlaps with a nucleic acid sequence corresponding to multiple predicted protein sequences and the respective overlaps are greater than a specified nucleic acid sequence overlap threshold, only one of the predicted protein sequences is reported as pETaG.

[0019] In some embodiments, the predicted protein sequence reported as pETaG is the predicted protein sequence in which the homologous nucleic acid sequence and the nucleic acid sequence corresponding to the predicted protein sequence show the highest sequence identity percentage, sequence coverage percentage, E value, or bit score value. In some embodiments, the predicted protein sequence reported as pETaG is the predicted protein sequence in which the homologous nucleic acid sequence and the nucleic acid sequence corresponding to the predicted protein sequence show the longest overlapping sequence. In some embodiments, if the homologous nucleic acid sequence overlaps with one or more nucleic acid sequences corresponding to one or more predicted protein sequences, but the respective overlaps are less than a specified nucleic acid sequence overlap threshold, the longest predicted protein sequence is reported as pETaG. In some embodiments, if the homologous nucleic acid sequence does not overlap with the nucleic acid sequence corresponding to the predicted protein sequence, the genomic coordinates of the homologous nucleic acid sequence are reported as pETaG. In some embodiments, the specified nucleic acid sequence overlap threshold has a value of at least 20%, 30%, 40%, 50%, 60%, 70%, 75%, 80%, 85%, 90%, 95%, or 98%.

[0020] In some embodiments, the comparative genomics heat map generated for a given pETaG comprises a plurality of cells arranged in a grid according to a first axis and a second axis, the first axis corresponds to a plurality of different target genomes, the plurality of different target genomes each comprising a plurality of positive genomes having an ortholog of an anchor gene sequence of one known BGC of the target genome and a plurality of negative genomes having no ortholog of the anchor gene sequence, the second axis corresponds to a plurality of query gene sequences or their orthologs co-localized with the anchor gene sequence of the known BGC, the putative embedded target gene (pETaG) is one of the plurality of co-localized query gene sequences, and the numerical value of each cell is based on (i) the presence or absence of the ortholog of each co-localized query gene sequence in the respective target genome, (ii) the sequence similarity of the ortholog to each co-localized query gene sequence, and (iii) whether the ortholog of each query gene sequence co-localizes with the ortholog of the anchor gene sequence in the respective genome.

[0021] In some embodiments, the computer-implemented method further includes analyzing the comparative genomics heatmap or underlying data thereof using a trained machine learning model, where the machine learning model is a classification model configured to output a probability for each of a plurality of predefined likelihood categories (e.g., categories of likelihood that a putative embedded gene is embedded in a gene cluster (e.g., BGC)). In some embodiments, the classification model is a long short-term memory (LSTM) model. In some embodiments, the classification model is a convolutional neural network (CNN) model. In some embodiments, the classification model is a vision transformer model, a generative adversarial network model, a variational autoencoder model, or a latent diffusion model. In some embodiments, for example, the plurality of predefined likelihood categories include (1) high likelihood, (2) somewhat high likelihood, (3) somewhat low likelihood, and (4) low likelihood.

[0022] 1. A computer-implemented method for determining the likelihood that a putative embedded target gene (pETaG) is a resistance gene for a secondary metabolite produced by a biosynthetic gene cluster (BGC) in a query genome, comprising: a) determining a likelihood that pETaG is associated with a BGC based on the presence or absence of orthologs of each of a plurality of query genes that are co-localized with the BGC in a plurality of different genomes, the plurality of genomes including a plurality of positive genomes that contain an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not contain an ortholog of the anchor gene of the BGC, the anchor gene being known to be associated with the BGC; ii) one or more phylogenetic characteristics of a last common ancestor (LCA) of homologs of pETaG in a phylogenetic tree of the plurality of genomes; and iii) determining a likelihood that pETaG is associated with a BGC based on the presence or absence of orthologs of each of a plurality of query genes that are co-localized with the BGC in a plurality of different genomes, the plurality of genomes including a plurality of positive genomes that contain an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not contain an ortholog of the anchor gene of the BGC, the anchor gene being known to be associated with the BGC. Also disclosed herein is a computer-implemented method comprising: determining one or more parameters selected from one or more scores indicative of co-occurrence of orthologs of pETaG and orthologs of anchor genes among a plurality of positive genomes; iv) one or more scores indicative of co-evolution of sequence diversity among orthologs of pETaG with respect to sequence diversity among orthologs of anchor genes in positive genomes including both orthologs of pETaG and orthologs of anchor genes; and v) one or more scores indicative of copy numbers of homologs of pETaG in a plurality of positive genomes and copy numbers of homologs of pETaG in a plurality of negative genomes; and b) determining a likelihood that pETaG is a resistance gene for secondary metabolites produced by BGC based on the one or more parameters. In some embodiments, the computer implementation further comprises predicting the probability that pETaG is an actual ETaG using one or more of the above parameters as inputs for a trained machine learning model. In some embodiments, the trained machine learning model is a trained neural network (e.g., an artificial neural network (ANN), a multi-layer perceptron (MLP), a deep neural network (DNN), a convolutional neural network (CNN), etc.). In some embodiments, the data can be utilized to train other types of machine learning models (e.g., Bayesian inference, decision tree-based methods such as XGBoost or random forests, etc.).In some embodiments, the data may be utilized to train a logistic regression model or other type of supervised model.

[0023] In some embodiments, pETaG co-localizes with the BGC in the query genome. In some embodiments, pETaG is not involved in the production of secondary metabolites by the BGC. In some embodiments, the anchor gene is the core synthase gene of the BGC.

[0024] In some embodiments, the computer-implemented method includes determining, for each of a plurality of pETaGs, a likelihood that the pETaG is a resistance gene for a secondary metabolite produced by a BGC in a target genome. In some embodiments, the computer implementation includes: a) identifying putative BGCs in the plurality of genomes having pairwise sequence similarity above a threshold; b) identifying a non-biosynthetic gene that co-localizes with an ortholog of an anchor gene in the putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in the organism of interest, and the non-biosynthetic gene is not involved in the production of a secondary metabolite by the BGC; c) for each of the plurality of query genes, identifying as pETaG the non-biosynthetic gene that encodes a protein having the highest sequence similarity to the protein of the respective target gene, and identifying the genome encoding the non-biosynthetic gene as the target genome; and d) determining, for each of the plurality of query genes, a likelihood that the respective pETaG is a resistance gene for a secondary metabolite produced by the respective BGC in the respective target genome.In some embodiments, the computer-implemented method includes: a) clustering the genomes in the database into a plurality of clusters, each cluster including genomes with pairwise sequence similarity above a threshold; and b) for each of the plurality of clusters, i) identifying a non-biosynthetic gene that co-localizes with an ortholog of an anchor gene in the putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in the organism of interest, and wherein the non-biosynthetic gene is not involved in the production of a secondary metabolite by the BGC; and ii) for each of the plurality of query genes, identifying an ortholog of an anchor gene that co-localizes with an ortholog of an anchor gene in the putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in the organism of interest, and wherein the non-biosynthetic gene is not involved in the production of a secondary metabolite by the BGC; The method includes identifying non-biosynthetic genes encoding proteins having a resistance to a secondary metabolite as candidate pETaGs; c) clustering the candidate pETaGs into a plurality of clusters based on sequence similarity among the pETaGs, and identifying candidate pETaGs encoding proteins having the highest sequence similarity to the protein of the respective query gene in each cluster as pETaGs and each genome encoding the pETaGs as a target genome; and d) determining, for each of the plurality of query genes, a likelihood that each pETaG is a resistance gene for a secondary metabolite produced by each BGC in each target genome.

[0025] In some embodiments, the threshold is at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least at least 98% pairwise sequence similarity. In some embodiments, each of the identified non-biosynthetic genes encodes a protein having at least about 30% sequence identity to the protein encoded by the respective query gene.

[0026] In some embodiments, the plurality of query genes is all protein-coding genes in an organism of interest. In some embodiments, the organism of interest is a mammal. In some embodiments, the organism of interest is a human. In some embodiments, the organism of interest is a reptile, bird, amphibian, animal, plant, fungus, or bacterium.

[0027] In some embodiments, the plurality of genomes are fungal genomes. In some embodiments, the plurality of genomes are bacterial genomes. In some embodiments, the plurality of genomes are plant genomes.

[0028] In some embodiments, each of the plurality of clusters comprises about 10 to about 100 genomes.

[0029] Disclosed herein is a computer-implemented method for identifying a druggable target in an organism of interest, comprising performing any one of the methods described herein and identifying the query gene as a druggable target based on the likelihood that the respective pETaG of the query gene is a resistance gene to a secondary metabolite produced by a BGC in the target genome. In some embodiments, the computer-implemented method further comprises identifying a secondary metabolite or an analog thereof as a small molecule modulator of the query gene or a protein encoded by the query gene. In some embodiments, the computer-implemented method of claim 62 further comprises contacting the secondary metabolite or an analog thereof with the protein encoded by the query gene and detecting the activity of the protein encoded by the query gene. In some embodiments, the number of positive genomes is equal to the number of negative genomes. In some embodiments, the computer-implemented method comprises selecting a plurality of positive genomes and a plurality of negative genomes from a database of genomes. In some embodiments, the computer-implemented method comprises clustering the database of genomes into a plurality of clusters based on sequence similarity and selecting one positive genome per cluster to provide a plurality of positive genomes. In some embodiments, the computer-implemented method comprises selecting the negative genome with the highest sequence similarity to each positive genome in the cluster.In some embodiments, the average pairwise sequence identity percentage of the orthologs of one or more single-copy genes in the positive genome is about 95% or less, and / or the average pairwise sequence identity percentage of the orthologs of one or more single-copy genes in the negative genome is about 95% or less.In some embodiments, the number of positive genomes is at least 5.

[0030] In some embodiments, the one or more parameters include a likelihood that pETaG is associated with the BGC based on the presence or absence of each of a plurality of query genes in the BGC in a plurality of different genomes. In some embodiments, determining the likelihood that pETaG is associated with the BGC includes: a) receiving a grid representation including a plurality of cells arranged according to a first axis and a second axis, the first axis corresponding to a plurality of genomes and the second axis corresponding to a plurality of query genes in the BGC in the query genome, each cell being based on i) the presence or absence of an ortholog of each query gene in the respective genome, ii) sequence similarity of the ortholog to the respective query gene, and iii) whether the ortholog of each query gene co-localizes with an ortholog of an anchor gene in the respective genome; and b) inputting the grid representation into a machine learning model, the machine learning model being trained to determine the likelihood that pETaG is associated with the BGC based on the values ​​of the plurality of cells in the grid representation, thereby providing the likelihood that pETaG is associated with the BGC.

[0031] In some embodiments, the computer-implemented method includes: a) identifying a putative BGC comprising pETaG from a library of putative BGCs from a plurality of genomes, and identifying the longest biosynthetic gene in the putative BGC as a core synthase gene; b) obtaining a plurality of positive genomes comprising an ortholog of the core synthase gene and a plurality of negative genomes not comprising an ortholog of the core synthase gene, wherein the plurality of positive genomes have pairwise sequence similarity below a threshold, and the plurality of negative genomes are selected based on sequence similarity to the plurality of positive genomes; and c) obtaining a grid comprising a plurality of cells arranged according to a first axis and a second axis. and creating a representation, wherein a first axis corresponds to all protein-coding genes that co-localize with the core synthase gene in the putative BGC in the query genome, and a second axis corresponds to a plurality of positive genomes and a plurality of negative genomes, and each cell is calculated based on i) the presence or absence of an ortholog of each protein-coding gene in the respective genome, ii) the sequence similarity of the ortholog to the respective protein-coding gene, and iii) whether the ortholog of each protein-coding gene co-localizes with an ortholog of the core synthase gene in the respective genome.

[0032] In some embodiments, the machine learning model is a classification model configured to output a probability for each of a plurality of predefined likelihood categories. In some embodiments, the classification model is a long short-term memory (LSTM) model. In some embodiments, the classification model is a convolutional neural network (CNN). In some embodiments, the classification model is a vision transformer model, a generative adversarial network model, a variational autoencoder model, or a latent diffusion model. In some embodiments, the plurality of predefined likelihood categories include (1) high likelihood, (2) somewhat high likelihood, (3) somewhat low likelihood, and (4) low likelihood.

[0033] In some embodiments, the one or more parameters include one or more phylogenetic features of the last common ancestor (LCA) of the homologs of pETaG in the phylogenetic trees of the multiple positive genomes and the negative genomes. In some embodiments, the one or more phylogenetic features are selected from the group consisting of the average copy number difference (CND) between the genes in the multiple positive genomes and the genes in the multiple negative genomes and the value determined from the multiple positive genomes, the ratio of the average to the LCA, the ratio of the standard deviation to the LCA, the ratio of the average of the adjacent distances, the standard deviation of the ratio of the adjacent distances, and the sum of the clade ratios. In some embodiments, the one or more parameters include one or more scores indicative of the co-occurrence of the orthologues of pETaG and the orthologues of the anchor gene in the multiple positive genomes. In some embodiments, the one or more scores indicative of the co-occurrence are selected from the group consisting of the co-occurrence pETaG distance, the co-occurrence pETaG rank, the co-occurrence core distance, and the co-occurrence core rank. In some embodiments, the one or more parameters include one or more scores indicative of co-evolution of sequence diversity between orthologs of pETaG with respect to sequence diversity between orthologs of the anchor gene in positive genomes that include both orthologs of pETaG and orthologs of the anchor gene. In some embodiments, the one or more scores indicative of co-evolution are selected from the group consisting of co-evolution correlation, co-evolution rank, and co-evolution slope. In some embodiments, the one or more parameters further include one or more features of the plurality of positive genomes and the plurality of negative genomes. In some embodiments, the one or more features are selected from the group consisting of the number of positive genomes, the average pairwise genome identity (PGI) between the positive genomes, the standard deviation of PGI between the positive genomes, the number of negative genomes, the average PGI between the negative genomes, and the standard deviation of PGI between the negative genomes.

[0034] In some embodiments, determining the likelihood based on the one or more parameters includes inputting the one or more features into a machine learning model, and the machine learning model is trained to determine the likelihood that pETaG is a resistance gene. In some embodiments, the machine learning model is a deep learning model. In some embodiments, the machine learning model is a decision tree model. In some embodiments, the machine learning model is a Bayesian inference model. In some embodiments, the one or more features can be utilized to train a logistic regression model or other type of supervised model. In some embodiments, whether a gene co-localizes with an anchor gene of a BGC is determined using antiSMASH. In some embodiments, whether a gene co-localizes with an anchor gene of a BGC is determined based on whether the gene is located within a proximal distance from the anchor gene. In some embodiments, the proximal zone is about 50 kb or less. In some embodiments, the proximal zone is about 20 kb.

[0035] Disclosed herein is a system comprising one or more processors and a memory, the memory being communicatively coupled to the one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to: i) receive as input one or more query sequences, or proxies thereof; ii) receive as input a selection of one or more target genomes; iii) perform a search of one or more target genomes using the one or more query sequences or proxies thereof to identify putative embedded target gene (pETaG) sequences that are homologs of the one or more query sequences based on a comparison of one or more homology-based metrics for candidate pETaGs to one or more predetermined homology-based metric thresholds; and iv) determine whether a given pETaG is an actual ETaG based on a comparative genomics analysis of multiple genomes related to the one or more target genomes.

[0036] 1. A system comprising: one or more processors; and a memory, the memory being communicatively coupled to the one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to execute a method for determining a likelihood that a putative embedded target gene (pETaG) is a resistance gene for a secondary metabolite produced by a biosynthetic gene cluster (BGC) in a query genome, the method comprising: a) determining a likelihood that pETaG is associated with a BGC based on the presence or absence of an ortholog of each of a plurality of query genes that are co-localized with the BGC in a plurality of different genomes, the plurality of genomes including a plurality of positive genomes that include an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not include an ortholog of the anchor gene of the BGC, the anchor gene being known to be associated with the BGC; ii) determining a likelihood that a putative embedded target gene (pETaG) is associated with a BGC based on the presence or absence of an ortholog of each of a plurality of query genes that are co-localized with the BGC in a plurality of different genomes, the plurality of genomes including a plurality of positive genomes that include an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not include an ortholog of the anchor gene of the BGC, the anchor gene being known to be associated with the BGC; Also disclosed herein is a system that includes determining one or more parameters selected from one or more phylogenetic features of the last common ancestor (LCA) of pETaG homologs in a phylogenetic tree of genomes, iii) one or more scores indicative of co-occurrence of pETaG orthologs and anchor gene orthologs among a plurality of positive genomes, iv) one or more scores indicative of co-evolution of sequence diversity between pETaG orthologs with respect to sequence diversity between anchor gene orthologs in a positive genome including both pETaG orthologs and anchor gene orthologs, and v) one or more scores indicative of copy numbers of pETaG homologs in a plurality of positive genomes and copy numbers of pETaG homologs in a plurality of negative genomes, and b) determining a likelihood that pETaG is a resistance gene for secondary metabolites produced by BGC based on the one or more parameters. In some embodiments, the method performed by the system further includes predicting the probability that pETaG is an actual ETaG using one or more of the above parameters as inputs to a trained machine learning model. In some embodiments, the trained machine learning model is a trained neural network (e.g., an artificial neural network (ANN), a multi-layer perceptron (MLP), a deep neural network (DNN), a convolutional neural network (CNN), etc.).In some embodiments, the data may be utilized to train other types of machine learning models (e.g., Bayesian inference, decision tree-based methods such as XGBoost or random forests, etc.) In some embodiments, the data may be utilized to train logistic regression models or other types of supervised models.

[0037] Disclosed herein is a system comprising one or more processors and a memory, the memory being communicatively coupled to the one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to perform any of the methods described herein.

[0038] Disclosed herein is a non-transitory computer-readable storage medium that stores one or more programs and includes instructions that, when executed by one or more processors of an electronic device, cause the electronic device to perform any of the methods described herein.

[0039] It should be understood that all combinations of the foregoing concepts, and additional concepts described in more detail below, are contemplated as part of the inventive subject matter disclosed herein (unless such concepts are mutually inconsistent). In particular, all combinations of claimed subject matter appearing at the end of this disclosure are contemplated as part of the inventive subject matter disclosed herein.

[0040] Incorporation by Reference All publications, patents, and patent applications mentioned in this specification are incorporated herein by reference in their entirety to the same extent as if each individual publication, patent, or patent application was specifically and individually indicated to be incorporated by reference in its entirety. In the event of a conflict between a term in this specification and a term in an incorporated reference, the term in this specification shall control.

[0041] Various aspects of the disclosed methods, devices, and systems are set forth with particularity in the appended claims. A better understanding of the features and advantages of the disclosed methods, devices, and systems will be obtained by reference to the following detailed description of exemplary embodiments and the accompanying drawings. [Brief description of the drawings]

[0042] [Figure 1] FIG. 1 shows exemplary putative biosynthetic gene clusters (BGCs) predicted by antiSMASH.

[0043] [Diagram 2] FIG. 1 provides a non-limiting example of a process flow chart for identifying and evaluating putative embedded target genes (pETaG).

[0044] [Diagram 3] FIG. 1 is an exemplary diagram of a positive genome (i.e., a genome that contains the core synthase gene sequence) and a negative genome (i.e., a genome that does not contain the core synthase gene sequence).

[0045] [Figure 4] FIG. 1 provides a non-limiting example of a process flow chart for identifying and evaluating putative embedded target genes (pETaG).

[0046] [Diagram 5] FIG. 13 is an exemplary diagram of different outcomes for evaluation of target pETaG search results.

[0047] [Figure 6] FIG. 1 illustrates an exemplary method for generating a grid representation (e.g., a heat map) of comparative genomics data.

[0048] [Figure 7]FIG. 1 shows an exemplary method for determining the likelihood that a putative embedded gene that co-localizes with a core synthase gene of a BGC in a query genome is associated with a BGC, according to some examples.

[0049] [Figure 8A-1] FIG. 1 shows an exemplary comparative genomics heat map. [Figure 8A-2] FIG. 1 shows an exemplary comparative genomics heat map. [Figure 8A-3] FIG. 1 shows an exemplary comparative genomics heat map.

[0050] [Figure 8B-1] FIG. 1 illustrates an exemplary long-short-term memory (LSTM) model used to classify input comparative genomics heatmaps into one of four likelihood categories (i.e., hierarchies). [Figure 8B-2] FIG. 1 illustrates an exemplary long-short-term memory (LSTM) model used to classify input comparative genomics heatmaps into one of four likelihood categories (i.e., hierarchies).

[0051] [Figure 8C-1] FIG. 8B-1 is a diagram showing an example of an output of the LSTM model shown in FIG. 8B-2. [Figure 8C-2] FIG. 8B-1 is a diagram showing an example of an output of the LSTM model shown in FIG. 8B-2.

[0052] [Figure 9A] 13 is a table comparing manual and machine learning-based classification of comparative genomics heat maps for "Tier A+", "Tier 1", "Tier 2", and "Tier 3" categories, respectively.

[0053] [Figure 9B]FIG. 1 shows a table comparing manual and machine learning-based classification of comparative genomics heat maps for "Tier A+", "Tier 1", "Tier 2", and "Tier 3", including positive predictive value, negative predictive value, sensitivity value, and specificity value.

[0054] [Figure 10A] FIG. 13 shows an exemplary comparative genomics heat map of the BGC of lovastatin, identifying the true boundaries of the BGC of lovastatin compared to the BGC predicted by antiSMASH. [Figure 10B] FIG. 13 shows an exemplary comparative genomics heat map of the BGC of lovastatin, identifying the true boundaries of the BGC of lovastatin compared to the BGC predicted by antiSMASH.

[0055] [Figure 11A-1] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier A+." [Figure 11A-2] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier A+."

[0056] [Figure 11B-1] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier 1." [Figure 11B-2] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier 1."

[0057] [Figure 11C-1] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier 2." [Figure 11C-2] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier 2."

[0058] [Figure 11D-1]FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier 3." [Figure 11D-2] FIG. 1 shows an exemplary comparative genomics heat map that was manually reviewed and classified as "Tier 3."

[0059] [Figure 12A] FIG. 1 illustrates a data table of a set of features (e.g., up to 27 or more features) that may be organized into a data table and utilized to train a machine learning model.

[0060] [Figure 12B] FIG. 1 illustrates the initial training stage of a neural network trained to output a probability value (i.e., an "embedded gene probability value" (e.g., "ETaG probability value")) that a putative embedded gene (e.g., pETaG) is associated with a BGC.

[0061] [Figure 12C] FIG. 13 illustrates a further training stage of a neural network trained to output probability values ​​for embedded genes (e.g., “ETaG probability values”).

[0062] [Figure 12D] FIG. 1 illustrates the inference stage of a neural network trained to output ETaG probability values.

[0063] [Figure 12E] FIG. 1 illustrates an exemplary data table containing an unknown input data set of features and the corresponding output of ETaG probability values ​​for each of the features.

[0064] [Figure 13] FIG. 1 is an exemplary diagram of a workflow for calculating phylogenetic features using a custom software algorithm (PhyloCCC).

[0065] [Figure 14]FIG. 1 provides non-limiting examples of phylogenetic signatures calculated from the lovastatin ETaG and a selected set of positive and negative genomes.

[0066] [Figure 15] FIG. 1 provides a non-limiting example of a workflow for assessing coevolution when using percent identity to compare coevolution between pairs of COGs.

[0067] [Figure 16] FIG. 1 provides a non-limiting example of the number of units used per hidden layer of a deep learning model trained to assess the likelihood that a putative ETaG is an ETaG.

[0068] [Figure 17] FIG. 1 provides a non-limiting example of performance data (test loss) for a deep learning model trained to assess the likelihood that a putative ETaG is an ETaG.

[0069] [Figure 18] FIG. 1 provides a non-limiting example of performance data (trial specificity) of a deep learning model trained to assess the likelihood that a putative ETaG is an ETaG.

[0070] [Figure 19] FIG. 1 provides a non-limiting example of performance data (test sensitivity) of a deep learning model trained to assess the likelihood that a putative ETaG is an ETaG.

[0071] [Figure 20] FIG. 1 provides a non-limiting example of performance data (test accuracy) for a deep learning model trained to assess the likelihood that a putative ETaG is an ETaG. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0072] The present disclosure provides a comparative genomics method and system for identifying orthologs of genes associated with biosynthetic gene clusters (BGCs) in a first genome present in one or more target genomes. In particular, the disclosed method and system can be used to identify putative embedded target genes (pETaGs) associated with, for example, resistance mechanisms in one or more target genomes based on a specified set of sequence homology search criteria. The pETaGs thus identified can then be filtered and evaluated for the likelihood that they are actual ETaGs associated with, for example, resistance mechanisms in the host organism from which the one or more target genomes are derived.

[0073] In some examples, a grid representation (e.g., a heat map) of the genomic data associated with a given pETaG can be used to filter and evaluate pETaGs based on comparative genomics analysis of each pETaG. The grid representation allows for visual and / or machine learning-based evaluation of co-occurrence and co-localization between genes found in close proximity to biosynthetic genes, biosynthetic gene clusters (BGCs), or other genes of interest, as described in more detail below.

[0074] In some examples, pETaGs may be filtered and evaluated based on evolutionary metrics based on various comparative genomics analyses. Examples of such evolutionary metrics include, but are not limited to, phylogenetic features, co-occurrence features, co-evolution features, and genomic dataset features, which may be determined for multiple genomes, including both positive genomes (i.e., genomes that contain core synthase gene sequences) and negative genomes (i.e., genomes that do not contain core synthase gene sequences), as described in more detail below.

[0075] In some examples, as described in more detail below, data derived from comparative genomics-based heatmaps and / or evolutionary metrics for pETAGs can be further processed, for example, using empirical algorithms and / or machine learning-based models, to determine a likelihood score or probability that a given pETaG is a real ETaG. Output from such analysis can then be used to compile a target lookup table (e.g., an array of data that collects sequence features and values ​​or ranges of evolutionary metrics for each pETaG, as well as the probability that it is a real ETaG).

[0076] definition Unless otherwise defined, all technical terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this disclosure belongs.

[0077] As used in this specification and the appended claims, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. Any reference herein to "or" is intended to include "and / or" unless specifically stated otherwise and includes any and all possible combinations of one or more of the associated listed items.

[0078] As used herein, the terms "includes," "including," "comprises," and / or "comprising" specify the presence of stated features, integers, steps, operations, elements, components, and / or units, but do not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components, units, and / or groups thereof.

[0079] As used herein, the term "about" a number refers to ±10% of that number. When used in the context of a range, the term "about" refers to the range minus 10% of its lowest value and plus 10% of its highest value.

[0080] As used herein, "secondary metabolite" refers to an organic small molecule compound produced by archaea, bacteria, fungi, or plants that is not directly involved in the normal growth, development, or reproduction of the host organism, but is required for the host organism's interaction with its environment. Secondary metabolites are also known as natural products or genetically encoded small molecules. The term "secondary metabolite" is used interchangeably herein with "biosynthetic product" when referring to the product of a biosynthetic gene cluster.

[0081] The term "biosynthetic gene cluster" or "BGC" is used interchangeably herein and refers to a locally clustered group of one or more genes that together code for a biosynthetic pathway for the production of a secondary metabolite. A BGC contains genes that encode signature biosynthetic proteins characteristic of each type of BGC. The longest biosynthetic gene in a BGC is referred to herein as the "core synthase gene" of the BGC. In addition to genes involved in the biosynthesis of secondary metabolites, a BGC may also contain non-biosynthetic genes, i.e., genes that encode products that are not involved in the biosynthesis of secondary metabolites, interspersed among the biosynthetic genes. Non-biosynthetic genes are referred to herein as "associated" or "embedded" in a BGC if their products are functionally related to the secondary metabolites of the BGC. The term "anchor gene" as used herein refers to biosynthetic or non-biosynthetic genes that are known to be co-localized and functionally related (i.e., associated) with a BGC.

[0082] The term "co-localized" refers to the presence of two or more genes closely spaced in a genome, not more than about 200 kb apart, not more than about 100 kb apart, not more than about 50 kb apart, not more than about 40 kb apart, not more than about 30 kb apart, not more than about 20 kb apart, not more than about 10 kb apart, not more than about 5 kb apart, or less apart.

[0083] The term "homolog" refers to a gene that is part of a group of genes related by descent from a common ancestor (i.e., the gene sequences (i.e., nucleic acid sequences) of the group of genes and / or the sequences of their protein products are inherited from a common origin). Homologs can arise through speciation events (giving rise to "orthologs"), or through gene duplication events, or through horizontal gene transfer events. Homologs can be identified by phylogenetic methods, through the identification of common functional domains in aligned nucleic acid or protein sequences, or through sequence comparison.

[0084] The term "ortholog" refers to a gene that is part of a group of genes predicted to have evolved from a common ancestral gene by speciation.

[0085] The terms "bidirectional best hit" and "BBH" are used interchangeably herein and refer to the relationship between a pair of genes in two genomes (i.e., a first gene in a first genome and a second gene in a second genome), where the first gene or its protein product is identified as having the most similar sequence in the first genome compared to the second gene or its protein product in the second genome, and the second gene or its protein product is identified as having the most similar sequence in the second genome compared to the first gene or its protein product in the first genome. The first gene is the bidirectional best hit (BBH) of the second gene, and the second gene is the bidirectional best hit (BBH) of the first gene. BBH is a commonly used method to infer orthology.

[0086] As used herein, "sequence similarity" between two genes means similarity in either the nucleic acid (eg, mRNA) sequences encoded by the genes or the amino acid sequences of the gene products.

[0087] "Percent sequence identity" or "percent sequence homology" with respect to nucleic acid sequences (or protein sequences) described herein is defined as the percentage of nucleotide residues (or amino acid residues) in a candidate sequence that are identical or homologous to the nucleotide residues (or amino acid residues) in an oligonucleotide (or polypeptide) to which the candidate sequence is compared, after aligning the sequences and considering any conservative substitutions as part of the sequence identity. The homology between different amino acid residues in a polypeptide is determined based on a substitution matrix such as BLOSUM (BLOcks Substitution MAtrixAlignment). Methods for aligning sequences and determining percent sequence identity or percent sequence homology of nucleic acid or protein sequences are well known to those skilled in the art. Examples of publicly available computer software that can be used include, but are not limited to, BLAST (Basic Local Alignment Search Tool; software for comparing amino acid sequences of proteins or nucleotide sequences of DNA and / or RNA molecules), BLAST-2, ALIGN or Megalign (DNASTAR) software. Any of a variety of appropriate parameters for measuring sequence alignment and determining percent sequence identity or homology can be determined by those skilled in the art, including the use of algorithms necessary to achieve maximal alignment over the full length of the sequences being compared.

[0088] Certain aspects of the present disclosure include process steps and instructions described herein in the form of an algorithm. It should be noted that the process steps and instructions of the present disclosure may be embodied in software, firmware, and / or hardware, and when embodied in software, may be downloaded to reside and operate on different platforms used by various operating systems. Unless otherwise stated in the following disclosure, descriptions utilizing terms such as "processing," "calculating," "computing," "determining," "displaying," "generating," and the like, will be understood to refer to the operations and processes of a computer system or similar electronic computing device that manipulates and transforms data represented as physical (electronic) quantities in the memory or registers of the computer system or other such information storage, transmission, or display devices.

[0089] Methods for discovering target genes embedded in biosynthetic gene clusters The systems and methods described herein relate to the identification of genes associated with biosynthetic gene clusters (BGCs). Figure 1 shows exemplary putative biosynthetic gene clusters (BGCs) predicted by antiSMASH, a genomics database search tool that enables rapid genome-wide identification, annotation and analysis of secondary metabolite biosynthetic gene clusters in, for example, bacterial, plant and fungal genomes.

[0090] As shown in Figure 1, a putative BGC may contain a set of genes encoding biosynthetic enzymes (the longest of which are referred to herein as "core biosynthetic proteins" or "core synthases") in a biosynthetic pathway for the production of a secondary metabolite, as well as a set of interspersed non-biosynthetic genes. Exemplary BGCs include, but are not limited to, biosynthetic gene clusters for synthesizing non-ribosomal peptide synthetases (NRPSs), polyketide synthases (PKSs), terpenes and bacteriocins. See, e.g., Keller N, "Fungal secondary metabolism: regulation, function and drug discovery," Nature Reviews Microbiology 17:3 (2019):167-180; and Fischbach M. and Voigt CA, PROKARYOTIC GENE CLUSTERS: A RICH TOOLBOX FOR SYNTHETIC BIOLOGY.: Institute of Medicine (US) Forum on Microbial Threats. The Science and Applications of Synthetic and Systems Biology: Workshop Summary. Washington (DC): National Academies Press (US); 2011, A21.

[0091] In some cases, the non-biosynthetic gene co-localized with the biosynthetic gene in the BGC can be a homolog of a human protein that includes a therapeutic target of interest. In some examples, the non-biosynthetic gene co-localized with the BGC can include a homolog of a mammalian gene, a reptilian gene, an avian gene, an amphibian gene, or a gene of interest from any other organism that includes a gene encoding a veterinary therapeutic target of interest. In some embodiments, the non-biosynthetic gene co-localized with the BGC can be a homolog of a fungal protein that includes a fungicidal target of interest. In some embodiments, the non-biosynthetic gene co-localized with the BGC can be a homolog of a plant protein that includes a herbicide target of interest. In some embodiments, the non-biosynthetic gene co-localized with the BGC can be a homolog of a bacterial protein that includes a microbial target of interest. In some examples, the non-biosynthetic gene can be functionally related to a secondary metabolite produced by the BGC. In some examples, the non-biosynthetic gene functionally related to a secondary metabolite produced by the BGC can be functionally related to a resistance mechanism that protects the host organism from the effects of the secondary metabolite. In some instances, the non-biosynthetic gene may encode a functionally unrelated protein product. The non-biosynthetic gene in the BGC that is a homolog of the human protein of interest is a putative embedded target gene (pETaG). The methods described herein allow for the identification of pETaG in one or more target genomes based on a known (or query) ETaG sequence in a first (or reference) genome (e.g., a human genome, a mammalian genome, a reptile genome, an avian genome, an amphibian genome, a plant genome, a bacterial genome, a fungal genome, or any other genome of interest), and utilize comparative genomics to determine the likelihood that a given pETaG is in fact a true ETaG associated with the BGC.

[0092] FIG. 2 provides a non-limiting example of a process 200 in which a query may be formulated and a search of one or more query genomes (or target genomes) may be performed in one or more genome databases to identify pETaGs, and then a given pETaG may be evaluated for likelihood of being a true pETaG. The process 200 may be performed, for example, as a computer-implemented method using software executing on one or more processors of one or more electronic devices, computers, or computer platforms. In some examples, the process 100 is performed using a client-server system, with blocks of the process 100 being divided in any manner between a server and a client device. In other examples, blocks of the process 100 are divided between a server and multiple client devices. Thus, although portions of the process 100 are described herein as being performed by a particular device of a client-server system, it will be understood that the process 100 is not so limited. In other examples, the process 100 is performed using only one client device or only multiple client devices. In process 100, some blocks are optionally combined, the order of some blocks is optionally changed, and some blocks are optionally omitted. In some examples, additional steps may be performed in combination with process 100. Thus, the operations shown (and described in more detail below) are exemplary in nature and, therefore, should not be considered limiting.

[0093] In step 202 of FIG. 2, a search query is formulated that includes one or more known query sequences (e.g., one or more gene sequences of interest). The search query may include one or more query sequences (or target sequences, e.g., one or more known ETaG sequences) or proxies thereof. For example, the search query may include one or more protein (amino acid) sequences, one or more nucleic acid (nucleotide) sequences, one or more Universal Protein Resource (Uniprot) identification numbers, one or more Profile Hidden Markov Models (pHMM; i.e., a probabilistic model that captures biological diversity of protein domains based on position-dependent scoring of multiple sequence alignments for corresponding proteins or nucleic acid sequences), a specified set of protein or nucleic acid sequence domains, or any combination thereof. The search query may include one or more query sequences (e.g., one or more gene sequences of interest) or proxies thereof selected from, for example, archaeal genomes, bacterial genomes, fungal genomes, plant genomes, animal genomes, human genomes, or any combination thereof. The search query may include one or more query sequences (e.g., one or more gene sequences of interest) selected according to the involvement of their corresponding protein products in, for example, a cellular biosynthetic pathway, a cell signaling pathway, a disease state (e.g., cancer, immunological and infectious diseases, or a rare disease), or any combination thereof.

[0094] In some examples, for example, in a target agnostic search, all protein sequences (as an example) from one or more genomes (e.g., from one or more organisms) can be selected as a query sequence for the query. In some examples, for example, in a targeted search, one or more specific protein sequences (as an example) from one or more genomes can be selected as a query sequence for the query. In some examples, the search query can include 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1000, 10,000, 20,000, 30,000, 40,000, 50,000, 60,000, 70,000, 80,000, 90,000, 100,000, or more than 100,000 query sequences (or any number of query sequences within this range). The formulation of the search query is described in more detail below.

[0095] In step 204 of FIG. 2, one or more target genomes are selected to be searched for homologous sequences. In some examples, the target genome can be selected from organisms of any kingdom of life (e.g., archaea, bacteria, fungi, plants, etc.) in which secondary metabolite molecules are known to be produced. In some examples, the target genome can be selected from organisms of any kingdom of life, e.g., archaea, bacteria, fungi, plants, animals, humans, or any combination thereof. In some examples, an individual genome can be used as the target genome. In some examples, multiple genomes obtained from the same or different kingdoms of life can be used as the target genome. In some examples, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 150, 200, 250, 500, 1,000, 5,000, 10,000, 20,000, 30,000, 40,000, 50,000, 60,000, 70,000, 80,000, 90,000, 100,000, or more than 100,000 genomes (or any number of genomes within this range) can be selected for search. In some examples, genomes can be grouped according to, for example, pairwise sequence identity percentage or phylogenetic distance before performing the search. The selection and grouping of target genomes to be searched are discussed in more detail below.

[0096] In step 206 of FIG. 2, a search is performed using one or more genomic and / or proteomic databases to identify sequences homologous to one or more query sequences (e.g., putative embedded target genes (pETaG)) in one or more target genomes. Identifying homologous sequences may include aligning the query sequence to sequences in one or more target genomes and determining one or more homology-based metrics. In some examples, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 500, 1000, 5,000, 10,000, 20,000, 30,000, 40,000, 50,000, 60,000, 70,000, 80,000, 90,000, 100,000, or more than 100,000 genomes (or any number of genomes within this range) may be selected as target genomes. Examples of suitable homologous sequence-based metrics include, but are not limited to, percent sequence identity, percent sequence coverage, E-value (a parameter indicating the expected number of hits of equivalent sequence similarity scores that can be found by chance), bit score (a measure of sequence similarity that is independent of query sequence length and database size and is normalized based on raw pairwise alignment scores), HMM score (a measure from a hidden Markov model that captures the biological diversity of protein domains based on position-dependent scoring of multiple sequence alignments of corresponding protein or nucleic acid sequences), or any combination thereof. Suitable search methods are described in more detail below.

[0097] In step 208 of FIG. 2, the pETaGs identified by the search are evaluated to determine the likelihood that a given pETaG is a true ETaG. The evaluation of the pETaG may be based, for example, on the use of comparative genomics heat maps and / or calculation of genome-derived metrics to determine whether a candidate pETaG is related to a candidate core synthase (or candidate BGC). In some examples, the comparative genomics analysis used to evaluate the pETaG may be based on identifying a set of positive and negative genomes. Methods for evaluating pETaGs to identify true ETaGs are described in more detail below.

[0098] FIG. 3 provides a schematic diagram of positive genomes (i.e., genomes of core synthase gene sequences) and negative genomes (i.e., genomes not containing core synthase gene sequences). In this example, genomes 1, 2, 3, 4, ......, N were aligned and searched to identify pETaG as well as other genes in a putative BGC. Genomes that contain embedded copies of pETaG sequences as well as other genes in a BGC are considered to be "positive" genomes. In addition to copies of pETaG sequences embedded in the BGC region, some positive genomes may contain additional copies of pETaG sequences (e.g., housekeeping copies). In some cases, the presence of one or more additional copies of pETaG sequences (e.g., housekeeping copies) may indicate pETaG involved in the resistance mechanism. Genomes lacking a BGC are considered to be "negative" genomes. As shown, some negative genomes may contain non-embedded pETaG sequences (e.g., only housekeeping copies).

[0099] FIG. 4 provides another non-limiting example of a flow chart of a process 400 for identifying and evaluating putative embedded target genes (pETaG).

[0100] Formulation of a search query: In step 402 of FIG. 4, a search query for identifying putative ETaGs is formulated. In some examples, the search may include, for example, an untargeted search, where all of the proteins in one or more organisms are selected as the query sequence. In some examples, the search may include, for example, a targeted search, where one or more specific proteins from one or more organisms are selected as the query sequence. In some examples, the search query may include 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1000, or more than 1000 known query sequences. In some examples, the query target (either an untargeted search or a targeted search) may be selected from 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or more than 10 organisms. The one or more organisms selected as the source of the query target may be obtained from any kingdom of life, for example, animals, plants, fungi, bacteria, archaea, etc.

[0101] Query targets for either target-agnostic or targeted searches can be specified as protein sequences (amino acids), gene sequences (nucleotides), Uniprot IDs, profile hidden Markov models (pHMMs), sets of specified protein or nucleic acid sequence domains, or any combination thereof. In some examples, queries can be grouped by the involvement of proteins of interest in a particular pathway, sequence similarity between proteins of interest, involvement in a particular disease, etc.

[0102] Selection of target genome: In step 404 of FIG. 4, a target genome to be searched for homologs of the protein of interest is selected from a genome database. The target genome can be selected for an organism from any kingdom of life (e.g., animals, plants, fungi, bacteria, archaea, etc.) where secondary metabolite molecules are known to be produced. In some examples, an individual genome can be used as the target genome. In some examples, all genomes in a given database can be used as the target genome. In some examples, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1000, 5,000, 10,000, 20,000, 30,000, 40,000, 50,000, 60,000, 70,000, 80,000, 90,000, 100,000, or more than 100,000 genomes (or any number of genomes within this range) can be selected as the target genome. Examples of suitable genomic databases include, but are not limited to, UniProt Knowledge Base (protein sequence and function database), Swiss-Prot (curated protein sequence database), GenBank (annotated collection of all publicly available DNA sequences), MycoCosm (Joint Genome Institute), PhycoCosm (Joint Genome Institute), Phytozome (Joint Genome Institute) and RefSeq (annotated set of sequences including genomic DNA, transcripts, and proteins). In some examples, the genomic database may include a publicly available database. In some examples, the genomic database may include a proprietary or private database.

[0103] In some examples, before performing the search candidate query (or target), the genomes selected for target-free or targeted search can be grouped into a set of two or more target genomes, for example, to maximize genome diversity, focus the search on a given genus, or focus the search on closely related genomes.Non-limiting examples of criteria that can be used to group two or more genomes include pairwise sequence identity, pairwise sequence similarity, and phylogenetic distance.In some examples, the group of query genomes can include 2, 4, 6, 8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or more than 100 genomes.

[0104] Pairwise sequence alignment methods are used to find the best pairwise (local or global) alignment of two query sequences. In some examples, pairwise sequence identity or sequence similarity can be calculated, for example, by comparing single copy protein (amino acid) sequences or gene (nucleotide) sequences shared between two genomes. Single copy proteins or genes can be annotated using tools such as BUSCO (Manni et al. (2021), BUSCO Update: Novel and Streamlined Workflows along with Broader and Deeper Phylogenetic Coverage for Scoring of Eukaryotic, Prokaryotic, and Viral Genomes, Mol. Biol. Evol. 38(10):4647-4654), which estimates the completeness and redundancy of processed genomic data based on universal single copy orthologs. In some examples, a predefined subset of proteins or genes can also be used. The selected proteins or genes are individually aligned, trimmed and concatenated to calculate pairwise sequence identity, i.e., the number of identical amino acid residues or nucleotides in the alignment. In some examples, pairwise sequence identity can be calculated by performing an alignment of the entire genome and calculating the percent sequence identity between the pair of genomes.

[0105] Alternatively, sequence similarity scores can be calculated based on substitution matrices such as BLOSUM (used to score local alignments between evolutionarily diverse protein sequences) and PAM (used to score global alignments between closely related protein sequences). See, e.g., Trivedi et al. (2020), "Substitution Scoring Matrices for Proteins-An Overview," Protein Science 29:2150-2163.

[0106] Genome grouping: In some examples, two or more target genomes can be grouped according to phylogenetic distance. Specific sequences such as single copy proteins or genes, or internal transcribed spacer (ITS) regions (i.e., DNA spacers located between small and large subunit ribosomal RNA (rRNA) genes in chromosomes (or corresponding transcribed regions in polycistronic rRNA precursor transcripts) can be used to create a phylogenetic tree from a set of genomes. The resulting phylogenetic tree is analyzed to identify the phylogenetic distance (or degree of genome divergence as indicated by branch lengths) between all pairs of genomes.

[0107] In some examples, when a genome grouping value such as a pairwise sequence identity percentage score, a pairwise sequence similarity score, or a phylogenetic distance score is calculated, the target genomes can be grouped by first applying a specified threshold value to the given grouping value to remove the candidate target genomes (or pairs of candidate target genomes) that do not meet the threshold, and then using a Markov Cluster (MCL) algorithm or another clustering algorithm to cluster the remaining target genomes into a set. In some examples, the pairs of candidate target genomes can be retained and grouped if their pairwise identity percentage or pairwise sequence similarity is greater than 50%, greater than 60%, greater than 70%, greater than 75%, greater than 80%, greater than 85%, greater than 90% or greater than 95%. Examples of other clustering algorithms that can be used to group the target genomes include, but are not limited to, k-means clustering methods, hierarchical clustering methods, mixed model methods, or any combination thereof.

[0108] In some examples, taxonomic relationships from known or curated databases can be used to group genomes based on, for example, phylum, class, order, family, genus, or species.

[0109] Target-agnostic search: In step 406 of Figure 4, a target-agnostic search is performed to determine the absence or presence of homologs of multiple query targets in one or more selected target genomes. The results of the target-agnostic search are then optionally de-replicated in step 410 by selecting one representative genome when multiple target genomes of the same species are identified as positive genomes.

[0110] Targeted Search: In step 408 of Figure 4, a targeted search is performed to determine the absence or presence of homologs of one or more query targets in one or more selected target genomes. The results of the untargeted search are then optionally de-replicated in step 412 by selecting one representative genome if multiple target genomes of the same species are identified as positive genomes.

[0111] For either target-free or target-targeted search, the search can be carried out using any of a variety of genomics search tools.Examples include but are not limited to BLAST, Diamond, HMMER (software for searching sequence databases for sequence alignment and homology using a probabilistic model called profile hidden Markov model (pHMM)), or ggsearch.

[0112] For either target-agnostic or targeted searches, the search may be performed across the entire genome of one or more selected target genomes, or may be limited to specific locations or regions of one or more selected target genomes (e.g., promoter regions, coding regions, introns, exons, termination sequences, etc.). In some examples, the search may be limited to regions of one or more selected target genomes that are known to be or predicted to be biosynthetic gene clusters (BGCs).In some instances, the BGC is integrated with antiSMASH (software for identifying, annotating, and comparing gene clusters encoding the biosynthesis of secondary metabolites in bacterial and fungal genomes; see, e.g., Medema et al. (2011), "antiSMASH: Rapid Identification, Annotation and Analysis of Secondary Metabolite Biosynthesis Gene Clusters in Bacterial and Fungal Genome Sequences," Nucleic Acids Research vol. 39, web server issue W339-W346), SMURF (web-based software for predicting clustered secondary metabolite genes based on their genomic context and domain content; see, e.g., Khaldi et al. (2010), "SMURF: Genomic Mapping of Fungal Secondary Metabolite Clusters," Fungal Genetics vol. 39, web server issue W339-W346), and / or SMURF (web-based software for predicting clustered secondary metabolite genes based on their genomic context and domain content; see, e.g., Khaldi et al. (2010), "SMURF: Genomic Mapping of Fungal Secondary Metabolite Clusters," Fungal Genetics vol. 39, web server issue W339-W346). Biol. 47(9):736-741), TOUCAN (a supervised learning framework for identifying fungal biosynthetic gene clusters based on protein or nucleotide sequences) e.g., Almeida et al. (2020), "TOUCAN: a framework for fungal biosynthetic gene cluster discovery", NAR Genomics and Bioinformatics 2(4):lqaa098, Deep-BGC (deep learning-based software for biosynthetic gene cluster prediction; see e.g., Hannigan et al. (2019), "A Deep Learning Genome-Mining Strategy for Biosynthetic Gene Cluster Prediction", Nucleic Acids Research 47(18):e110), or custom search algorithms.

[0113] BGC Prediction: BGCs can be predicted by extracting regions (of a specified length) from the target genome adjacent to gene sequences that are likely to be homologs of known biosynthetic core synthase genes based on sequence searches (using tools such as BLAST, Diamond, or ggsearch), HMMs of known core synthases (using tools such as HMMER), or co-localization of protein domains associated with core synthases. If a match with the target gene sequence is found in the predicted BGC or in a specified sequence region of the target genome, a candidate putative embedded target gene (pETaG) is identified. If the search results exceed a specified threshold for sequence similarity according to any of a variety of metrics known to those skilled in the art, the search results are considered a match with the query target. Examples include percent sequence identity (e.g., at least 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, 96%, 97%, 98%, 99%, or more), percent sequence coverage (e.g., at least 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, 96%, 97%, 98%, 99%, or more), E-value (e.g., 10, 1, 0.1, 0.001, 0.0001, 1e -10 , 1e -20 , 1e -100 , or less), bit score (e.g., 5, 10, 25, 50, 100, 250, 500, 1000, 5000 or more), or HMM score (e.g., 5, 10, 25, 50, 100, 250, 500, 1000, 5000 or more).

[0114] Highly sensitive targeted search techniques: In some examples of targeted searches, additional search methods can be utilized to increase the sensitivity of the search and capture the hard-to-find signal of the presence of pETaG. Searches of query targets, including protein or nucleotide sequences, or pHMMs of query targets, are often performed using sequence alignment tools such as BLAST, Exonerate, and HMMER. When the query target is a protein sequence, a protein-DNA search tool such as TBLASTN or Exonerate can be used to convert the protein sequence target to a nucleotide sequence, and then a search of the nucleotide sequence can be performed against one or more query genomes (or target genomes). Alternatively, a probability model such as PFAM or TIGRFAM can also be used with an HMM search tool such as HMMER. The resulting search hits can be filtered by genome location and / or by comparison of sequence metrics such as E-value, sequence identity percentage, query coverage (e.g., the degree of coverage of the query sequence by the identified homologous sequences), target coverage (e.g., the degree of coverage of the identified homologous sequences by the query sequence of interest), or bit score with a corresponding cutoff threshold. The cutoff threshold for percent sequence identity may be at least 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, 96%, 97%, 98%, 99%, or more. The cutoff threshold for target coverage and / or query coverage may be at least 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, 96%, 97%, 98%, 99%, or more. The cutoff threshold for E-value may be 10, 1, 0.1, 0.001, 0.0001, 1e -10 , 1e -20 , 1e -100 , or lower. The bit score cutoff threshold may be 5, 10, 25, 50, 100, 250, 500, 1000, 5000 or more.

[0115] DNA alignment regions from the search results can be evaluated by comparing the genomic coordinates of the alignment region with a given protein sequence prediction in the query genome. If the DNA alignment region overlaps with a single predicted protein sequence and the corresponding sequence similarity metrics (e.g., E-value, sequence identity percentage, query coverage, target coverage, or bit score) of the DNA alignment region and the predicted protein sequence exceed a specified threshold of sequence similarity metrics (as above), the corresponding nucleic acid sequence is reported as pETaG (region case 1).

[0116] If a DNA alignment region overlaps with multiple predicted protein sequences and the corresponding sequence similarity metric (e.g., E-value, percent sequence identity, query coverage, target coverage, or bit score) for the DNA alignment region and each of the multiple predicted protein sequences exceeds a specified threshold, only one of the predicted protein sequences (or its corresponding nucleic acid sequence) is reported as pETaG (Region Case 2). In some examples, the decision of which overlapping protein sequence to report as pETaG may be based on which has the highest corresponding sequence metric (i.e., the protein sequence (or its corresponding nucleic acid sequence) with the highest E-value, percent sequence identity, query coverage, target coverage, or bit score is reported as pETaG). In some examples, the decision of which overlapping protein sequence to report as pETaG may be based on protein sequence length (i.e., the longest overlapping protein sequence (or its corresponding nucleic acid sequence) is reported as pETaG).

[0117] If a DNA alignment region overlaps with a single predicted protein sequence or multiple predicted protein sequences, but none of the corresponding sequence similarity metrics exceed a specified threshold, the longest protein sequence or the protein sequence with the highest corresponding sequence similarity metric value is reported as pETaG (region case 2).

[0118] If the DNA alignment region does not overlap with the predicted protein sequence, the coordinates of the DNA alignment region are reported as pETaG (region case 3).

[0119] Reviewing the results of a targeted search and reporting pETaG: Figure 5 shows the different cases described for the results of a targeted search. The protein case shown at the top of the figure illustrates a scenario where a targeted search method (e.g., protein search, HMM search, etc.) identifies a complete protein in a target genome of interest as a direct hit to the search query (e.g., human gene) and the protein ID is reported to correspond to pETaG. Region case 1 illustrates a scenario where a DNA sequence search identifies a DNA alignment region that overlaps with a portion of a predicted protein sequence with high sequence similarity (e.g., >70% sequence identity in this example) and the protein ID is reported to correspond to pETaG. Region Case 2 illustrates a scenario where the DNA sequence search identifies a DNA alignment region that overlaps with a predicted gene sequence, a stretch of a predicted gene sequence, or portions of multiple predicted gene sequences, none of the overlapping sequences exhibit a sequence similarity metric above a specified threshold (e.g., >70% sequence identity), and the protein ID of the protein sequence exhibiting the greatest overlap and the coordinates of the DNA alignment region are reported as pETaG. Region Case 3 illustrates a scenario where the DNA sequence search identifies a DNA alignment region that does not overlap with a predicted protein sequence in the target genome of interest, and the coordinates of the DNA alignment region are reported as pETaG.

[0120] Generation of Comparative Genomics Heat Maps: Returning to FIG. 4, in step 414, a comparative genomics heat map is generated for each pETaG identified by the untargeted or targeted search, and then used to evaluate the pETaG. The heat map, or the underlying data from the heat map, can be used, for example, to evaluate the degree of "embeddedness" (i.e., the degree of association with the BGC) of a given pETaG. In some examples, a vector representation of the data contained in the heat map can be provided to one or more trained neural networks (e.g., trained long-short-term memory (LSTM) models) to perform embeddedness classification of the corresponding pETaG. Methods for generating comparative genomics heat maps and using them to evaluate pETaGs are described in more detail below.

[0121] Calculation of grouped genome-derived metrics: In step 416 of FIG. 4, when two or more target genomes are optionally grouped into a set, a grouped genome-derived metric is calculated. The grouped genome-derived metric can then be used to evaluate the pETaGs. In some examples, the target genomes can be grouped before performing a query search, as described above. In some examples, the target genomes can be grouped after performing a query search of all target genomes for a query target of interest. For each set of target genomes and each pETaG in the set of genomes, a number of search features can be calculated. Examples of grouped genome-derived metrics (or search features) include, but are not limited to, (i) the number of "positive" genomes in the genome set (i.e., the number of genomes that contain a core synthase and for which a BGC is predicted, e.g., by using antiSMASH, SMURF, TOUCAN, deepBGC, or a custom BGC identification method), (ii) the number of "negative" genomes in the genome set (i.e., the number of genomes that do not contain a core synthase), and (iii) the copy number difference (CND) of a given gene between the positive and negative genomes (pETaG is assessed, for example, as an indicator of a gene's involvement in the resistance mechanism is duplication, e.g., an extra copy of a housekeeping gene). In some examples, a positive genome may be defined as a genome that contains a candidate BGC independent of the core synthase gene. In some examples, a negative genome may be defined as a genome that does not contain a candidate BGC independent of the core synthase gene.

[0122] Clusters of Orthologous Groups (COGs): In some examples, clusters of orthologous groups (COGs) of proteins or genes can be used to determine copy number differences ("COG CNDs"). COGs can be rapidly generated and provide a somewhat orthogonal / alternative method for calculating copy number differences than tree-based approaches. COGs can be identified by performing an all-versus-all protein (amino acid) sequence search using all of the target genomes in a given genome set using tools such as BLASTp, ggsearch or Diamond-BLASTp. Members of the same COG are presumed to have orthologous functions. COGs can also be defined as gene or protein families (or orthogroups, i.e., a set of genes or proteins derived from a single gene in the last common ancestor of the set of target species / genomes under consideration). COGs can also be identified using tools such as OrthoMCL, OrthoFinder, PanX, or other orthogroup / pan-genome generating tools, or using protein clustering tools such as USEARCH, CD-HIT, and MMseqs. Pairwise associations are established between each query target and the corresponding target protein sequences identified in the search, for example, based on either percent sequence identity or E-value, and then clustered using MCL, SiLiX, or other clustering algorithms. For a COG of interest, the copy number difference is calculated by subtracting the average number of target gene homologs in the negative genome from the average number of target gene homologs in the positive genome present in the COG.

[0123] Agnostic tree CND: In some examples, copy number differences can be calculated from a phylogenetic tree created for a set of grouped target genomes (i.e., "agnostic tree CND"). COGs are identified as described above, or in some examples, the same COGs used to determine the COG CND can also be used to determine the agnostic tree CND. For COGs containing pETaG, a multiple sequence alignment is created using any multiple sequence alignment software tool (e.g., MAFFT, MUSCLE, ClustalW, etc.), and then trimmed using any sequence trimming software tool (e.g., trimAI, GBlocks, ClipKIT, etc.). The resulting trimmed sequence alignment is then used to create a phylogenetic tree using any phylogenetic tree reconstruction software (e.g., FastTree, IQ-TREE, RAxML, MEGA, MrBayes, BEAST, PAUP, etc.). Additionally, the phylogenetic tree can also be constructed using various algorithms, such as maximum likelihood algorithms, maximum parsimony algorithms, neighbor-joining algorithms, distance matrix algorithms, Bayesian estimation algorithms, or any combination thereof. The resulting phylogenetic tree is then computationally analyzed to identify the last common ancestor clade of pETaG and the housekeeping version of pETaG. The last common ancestor can be determined by first identifying a clade (i.e., the pETaG clade) that contains all or a defined subset of pETaG from the positive genome. The tree can then be traversed from the pETaG clade back toward the root, checking whether each new clade contains genes from all or a defined subset of the negative genome. Once the correct clade is identified (i.e., a clade that contains all or a defined subset of pETaG and genes from all or a defined subset of the negative genome), it is identified as the last common ancestor clade.The independent tree CND is then calculated for a given pETaG by subtracting the average number of gene homologs in the negative genome from the average number of gene homologs in the positive genome that are present in the last common ancestor clade.

[0124] Evolutionary Metrics: In some examples, the candidate pETaGs may also be evaluated based on one or more evolutionary metrics (e.g., phylogenetic features, co-occurrence features, or co-evolution features) determined in step 416 of FIG. 4 using a custom software algorithm that utilizes a set of positive and negative genomes to evaluate whether the candidate pETaGs are associated with a candidate core synthase or BGC. A positive genome is defined as a genome that contains a candidate core synthase in a candidate BGC, and at least one of the positive genomes contains a candidate pETaG. In some embodiments, a positive genome may be defined as a genome that contains a candidate BGC independent of the core synthase gene. A negative genome is defined as a genome that does not contain a candidate core synthase in a candidate BGC. In some embodiments, a negative genome may be defined as a genome that does not contain a candidate BGC independent of the core synthase gene. The positive and negative genomes used to calculate the evolutionary metrics may be identified using any of a variety of methods known to those of skill in the art. For example, the comparative genomics heatmaps described elsewhere herein or the data underlying them may be fed into the custom software algorithm. Alternatively, any method may be used that utilizes protein or DNA sequence searches to determine whether a query sequence of interest co-localizes with a core synthase and / or is located in a predicted BGC.

[0125] To evaluate the evolutionary metrics of candidate pETaGs, a minimum number of positive and / or negative genomes are required to perform a particular analysis. For example, phylogenetic features require at least one positive genome and at least one negative genome for analysis. Co-occurrence features require at least one positive genome and at least one negative genome. Co-evolution features require at least three positive genomes, but do not require identifying negative genomes. The actual number of positive and negative genomes used can be any number greater than the minimum number required for a given analysis. In some examples, the number of positive and / or negative genomes used to evaluate evolutionary metrics can be at least 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, or more than 60 genomes. In some examples, the number of combined positive and / or negative genomes used does not exceed 20 genomes to maintain computational efficiency.

[0126] In some examples, positive genomes may be optionally de-replicated by species (i.e., if multiple genomes of the same species are identified as positive genomes, they are filtered so that only one representative genome is kept). Genome de-replicated in several ways, for example, based on taxonomy if the species name is known, based on pairwise sequence identity or similarity as described above for the search methodology (species level distinction can be made at or below a specified threshold for pairwise sequence identity (e.g., 99.9%, 99.8%, 99.7%, 99.6%, 99.5%, 99.4%, 99.3%, 99.2%, 99.1%, 99%, 98%, 97%, 96%, 95%, 94%, 93%, 92%, 91%, 90%, 85%, 80% or less) or a specified threshold for similarity (threshold is , which varies depending on sequence length, number of sequences used, and similarity matrix used), or based on phylogenetic distance as described above for the search methodology (species level distinction can be determined based on a specified threshold for pairwise phylogenetic distance (e.g., 0.0001, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.02, 0.03, 0.04, 0.05, 0.1, or higher)). In some embodiments, positive genomes may optionally be de-replicated by other taxonomic rankings, such as, but not limited to, genus, family, order, class, or phylum.

[0127] In some examples, the genome selected for use as the negative genome may be the genome that is closest to the positive genome but is classified as negative. The closeness of the negative genome to the positive genome can be determined in any of several ways, for example, based on taxonomy if the species name is known, based on pairwise sequence identity or similarity as described above for the search methodology (e.g., the negative genome with the highest identity or similarity percentage to the positive genome is selected), or based on phylogenetic distance as described above for the search methodology (e.g., the negative genome with the lowest phylogenetic distance to the positive genome is selected).

[0128] Phylogenetic features: Phylogenetic features can be calculated from the phylogenetic tree of the last common ancestor clade of the candidate pETaG and its housekeeping copy. The process of creating a phylogenetic tree and identifying the last common ancestor is described below. In the first step of the process, pETaG is used as a query sequence to identify all homologs in the set of selected positive and negative genomes using any sequence search and alignment tool (e.g., BLASTp, ggsearch, or Diamond-BLASTp, etc.). The homolog sequences are then aligned using any alignment software tool (e.g., MAFFT, MUSCLE, ClustalW, etc.) and trimmed using any sequence trimming software (e.g., trimAI, GBlocks, ClipKIT, etc.). The resulting trimmed sequence alignment is then used to create a phylogenetic tree using any phylogenetic reconstruction software tool (e.g., FastTree, IQ-TREE, RAxML, MEGA, MrBayes, BEAST, PAUP, etc.). Alternatively, the phylogenetic tree can be constructed using various algorithms, such as maximum likelihood algorithms, maximum parsimony algorithms, neighbor-joining algorithms, distance matrix algorithms, Bayesian inference algorithms, or any combination thereof. The resulting phylogenetic tree is then computationally analyzed to identify the last common ancestor clade of pETaG and the housekeeping version of pETaG. The last common ancestor is determined by first identifying a clade (i.e., a pETaG clade) that contains all or a subset of pETaG from the positive genome. The tree can then be traversed by working back from the pETaG clade toward the root and checking whether each new clade contains genes from all or a subset of the negative genome. Once the correct clade is found (i.e., a clade that contains all or a subset of pETaG and genes from all or a defined subset of the negative genome), it is identified as the last common ancestor clade (LCA).

[0129] Last Common Ancestor (LCA) clades can be used to calculate various metrics from phylogenetic trees, including but not limited to:

[0130] Phylogenetic CND - calculated for a given gene (e.g., pETaG) by subtracting the average number of gene copies in the negative genomes from the average number of gene copies in the positive genomes present within the last common ancestor clade.

[0131] Average ratio to LCA - calculated as the average branch length of all pETaG genes to an LCA node divided by the average branch length of all house copy genes to an LCA node.

[0132] Ratio of standard deviation (stdev) for LCA - calculated as the standard deviation of branch lengths of all pETaG genes for an LCA node divided by the standard deviation of branch lengths of all house copy genes for an LCA node.

[0133] Average neighbor distance ratio - calculated as the average branch length of each pETaG gene to all other pETaGs divided by the average branch length of each house copy gene to all other house copy genes.

[0134] Ratio of standard deviation (stdev) for adjacent distance - calculated as the standard deviation of the branch lengths of each pETaG gene relative to all other pETaGs divided by the standard deviation of the branch lengths of each house copy gene relative to all other house copy genes.

[0135] Clade sum ratio - calculated as the sum of all branch lengths in pETaG clades (i.e., clades in the phylogenetic tree that contain all or a subset of pETaGs) divided by the sum of all branch lengths in house copy clades (i.e., clades in the phylogenetic tree that contain all or a subset of the house copy genes).

[0136] Alternative phylogenetic metrics: In some cases, alternative metrics can be calculated from the phylogenetic tree and used as additional evidence that pETaG is likely to be an actual ETaG. Examples include, but are not limited to, the Robinson-Foulds (RF) distance metric used to calculate the distance between two phylogenetic trees (in this case, the distance between the pETaG clade and the House copy clade). Similar RF distance metrics or adjustment-based distance metrics can be calculated using tools such as treeKO to compare the topology and distance between pETaG and the House copy clade.

[0137] Co-occurrence features: Co-occurrence features can be calculated from clusters of orthologous groups (COGs) of protein or gene sequences. COGs can be identified by performing all protein (amino acid) to all nucleotide to all sequence searches using software tools such as BLAST, ggsearch, or Diamond, using all positive and negative genomes. The mutual best hits (i.e., the best protein match between genome A and genome B must also be the best protein match between genome B and genome A) are used to establish associations between protein sequences (or their nucleotide sequence counterparts), and then clustered using, for example, MCL, SiLiX, or other clustering algorithms. Alternatively, in some instances, unidirectional search results can be used (instead of mutual searches) to create associations prior to clustering. COGs can also be identified using tools such as OrthoMCL, OrthoFinder, PanX, or other orthogroup / pan-genome generation tools, or using protein clustering tools such as USEARCH, CD-HIT, and MMseqs.

[0138] Various co-occurrence metrics can be calculated from the identified COGs. These include, but are not limited to:

[0139] Normalized Co-occurrence Distance Metric: The normalized co-occurrence distance metric is calculated using the following formula:

number

[0140] Co-occurrence pETaG distance - calculated as the distance score of COGs containing pETaG.

[0141] Co-occurrence pETaG rank - the ascending rank of the distance score of the COG containing the pETaG relative to the distance scores of the other COGs. In the event of a tie in distance scores, the rank assigned to all COGs in the tie is the lowest rank in the group.

[0142] Co-occurrence core distance - distance score of COGs containing the core synthase.

[0143] Co-occurrence CoreRank - Rank of the COGs containing the core synthase in ascending order of distance score relative to the distance scores of other COGs. In the event of a tie in distance scores, the rank assigned to all COGs in the tie is the lowest rank in the group.

[0144] Co-evolutionary features: Co-evolutionary features can also be calculated from clusters of orthologous groups (COGs) of protein or gene sequences. COGs can be identified in a similar manner as above, for example using tools such as BLASTp, ggsearch or Diamond-BLASTp to perform a sequence search of all proteins (amino acids) using all of the target genomes in a given genome set, followed by clustering the pairwise sequence similarity results. Only COGs that contain at least three genes that are single copy genes are passed to the co-evolutionary analysis. Single copy means that the COG consists of only one gene / protein per genome of the genomes present in the COG.

[0145] Coevolutionary analysis involves the determination of COG-COG comparison values, which can be done using several different techniques, such as multiple sequence alignment (MSA) or gene-tree comparison.

[0146] A multiple sequence alignment (MSA) can be created for each pair of COGs using any of a variety of alignment software tools, including but not limited to MAFFT, MUSCLE, and ClustalW. The MSA can then be trimmed based on specified parameters (e.g., removing all gaps, removing gaps where the number of consecutive gaps is greater than a specified threshold, keeping all gaps, etc.), and the percent sequence identity (the number of identical residues in the alignment) is calculated. Alternatively, sequence similarity scores can be calculated based on alternative matrices such as BLOSUM and PAM.

[0147] Gene phylogenetic tree comparisons can be calculated based on phylogenetic trees generated from amino acid or nucleotide sequences within each pair of COGs by sequence alignment, trimming, and phylogenetic reconstruction. The phylogenetic trees must be constrained to the topology of the species tree of the target genome, and genes present in both COGs are compared. Pairwise branch lengths between the two resulting COG phylogenetic trees are then calculated. This analysis can be performed using custom scripts or using software tools such as the Co-Variance algorithm in PhyKIT (https: / / github.com / JLSteenwyk / PhyKIT).

[0148] Correlations of all pairwise COG combinations are then calculated using either percent sequence identity, sequence similarity, or branch length comparisons using Pearson R, or any other correlation metric known to those skilled in the art. Correlations can only be calculated between pairs of COGs that share at least three genomes.

[0149] Examples of co-evolutionary metrics that may be calculated include, but are not limited to, the following:

[0150] Coevolutionary Correlations - COG X Pairwise sequence identity percentage and COGs Y Correlation with percent pairwise sequence identity.

[0151] Coevolutionary rank - the ascending rank of the correlation between COGx and COGy relative to all other pairwise COG correlations. In the case of a tie, the rank assigned to all pairwise COG correlations in the tie is the lowest rank in the group.

[0152] Coevolutionary Gradient - COG X Pairwise identity percentage and COGs Y Orthogonal regression with pairwise identity percentages.

[0153] In some instances, COGx is a COG that contains a core synthase and COGy is a COG that contains a pETaG.

[0154] Genomic Dataset Features: For a group of target genomes clustered as described above, various genomic dataset features can be calculated. These include, but are not limited to:

[0155] Number of positive genomes - The final number of positive genomes used in any and all metrics above. In some examples, 10, 15, 20, 25, 30, or more positive genomes may be used as input.

[0156] PGI positive average - the average pairwise genomic identity (PGI) between the final positive genome counts. Pairwise genomic identity can be calculated, for example, by comparing single copy protein (amino acid) or gene (nucleotide) sequences shared between two genomes.

[0157] PGI positive standard deviation (stdev) - the standard deviation of pairwise genome identities between the final positive genome counts.

[0158] Number of negative genomes - The final number of negative genomes used in each and every metric above. In some examples, 10, 15, 20, 25, 30, or more negative genomes may be used as input.

[0159] PGI Negative Mean - The average of pairwise genome identities between the final negative genome counts.

[0160] PGI negative standard deviation (stdev) - the standard deviation of pairwise genome identities between the final negative genome counts.

[0161] Returning to FIG. 4, in step 418, the comparative genomics heatmaps generated for each pETaG identified in the search, or its underlying data, as described above, can be used to evaluate the degree of "embeddedness" of a given pETaG. In some examples, the vector representation of the data contained in the heatmaps may be provided to one or more trained neural networks (e.g., trained long-short-term memory (LSTM) models, convolutional neural networks (CNNs)) to perform embeddedness classification of the corresponding pETaG. For example, the comparative genomics heatmaps, or the underlying data therefrom, can be analyzed using a trained LSTM model to predict the probability that a pETaG is associated with a BGC. Methods for generating comparative genomics heatmaps and using them to evaluate pETaGs are described in more detail below.

[0162] In step 420 of Figure 4, a "feature table" is compiled that summarizes the grouped genome-derived metrics and / or embeddedness classification data for the pETAGs identified in the search. In some examples, the feature table, or the data contained therein, can be used as input for a machine learning-based analysis (e.g., a deep learning-based analysis) to evaluate the pETaGs.

[0163] Machine learning-based (e.g., deep learning-based) evaluation of pETaG: In step 422 of FIG. 4, a machine learning-based analysis (e.g., deep learning-based analysis) can be performed to evaluate the likelihood (or probability) that a given pETaG is an actual ETaG. The input for the machine learning-based analysis is the feature table compiled in step 420 of FIG. 4 or the data contained therein. In some examples, predicted BGC regions of the target genome identified using a BGC search tool as shown in step 424 can also be used as input for the machine learning-based evaluation. As shown, in some examples, the predicted BGC regions of the target genome can also be used as part of the selection of the target genome to be used for the search (e.g., in step 404 of FIG. 4). As mentioned above, BGCs can be predicted by extracting regions (of a specified length) from the target genome adjacent to gene sequences that are likely to be homologs of known biosynthetic core synthase genes based on sequence searches (using tools such as BLAST, Diamond, or ggsearch), HMMs of known core synthases (using tools such as HMMER), or co-localization of protein domains associated with core synthases.

[0164] Any of a variety of machine learning algorithms can be used in implementing the disclosed pETaG evaluation methods. For example, the machine learning algorithms employed may be supervised learning algorithms (e.g., algorithms that rely on the use of a set of labeled training data to infer a relationship between a set of one or more pETaG features and a prediction of the probability that a given pETaG is an actual ETaG), unsupervised learning algorithms (e.g., algorithms used to draw inferences from a training data set consisting of pETaG features that are not paired with labeled ETaG classification or probability data), semi-supervised learning algorithms (e.g., algorithms that utilize both labeled and unlabeled pETaG feature data for training (typically using a relatively small amount of labeled data along with a large amount of unlabeled data), or the like. The learning algorithms may include algorithms that model the learning task as a series of individual decisions, such as gradient boosted trees where new trees are created to model errors or residuals to complement and strengthen existing trees that may be used to split or classify pETaGs based on feature data, deep learning algorithms, such as algorithms inspired by the structure and function of the human brain, such as artificial neural networks (ANNs), specifically large neural networks containing many hidden layers of connected "nodes" that can be used to map pETaG feature data to a probabilistic prediction or classification decision, or any combination thereof.

[0165] A deep learning model (i.e., a trained deep learning algorithm) can include any total number of layers, and any number of hidden layers, where the hidden layers act as trainable feature extractors that allow mapping a set of input data to a preferred output value or set of output values. Each layer of the neural network includes several nodes (or units). A node receives inputs coming directly from the input data (e.g., pETaG feature data derived using the method described above) or from the output of a node of a previous layer, and performs a certain operation, e.g., an addition operation. In some cases, the connection from the input to the node is associated with a weight (or weight coefficient). In some cases, a node may be connected to, for example, an input X i and their associated weights W i All pairwise products with may be summed. In some cases, the weighted sum is offset by a bias b. In some cases, the output of the node may be gated using an activation function f, which may be, for example, a threshold, or a linear or nonlinear function. The activation function may be, for example, a rectified linear unit (ReLU) activation function, or other functions such as a saturating hyperbolic tangent, identity, binary step, logistic, arcTan, soft sine, parametric rectified linear unit, exponential linear unit, softPlus, bent identity, softExponential, Sinusoid, sine, Gaussian, or sigmoid function, or any combination thereof.

[0166] The weighting coefficients, bias values, and thresholds, or other computational parameters of a neural network can be "taught" or "learned" in a training phase that uses one or more sets of training data. For example, the parameters can be trained using input data from a training data set and gradient descent or backpropagation techniques such that the output values ​​that the deep learning model computes (e.g., predictions of the probability that a given pETaG is an actual ETaG) match the examples contained in the training data set.

[0167] For example, in some examples, a deep learning model trained to predict the probability that a given pETaG is an actual ETaG can include a fully connected neural network with a specified number of hidden layers (e.g., 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 100, 1000, or more than 1000 hidden layers) and a specified number of units per hidden layer (e.g., 1, 2, 4, 8, 16, 32, 64, 126, 256, or more than 256 units per hidden layer). The output layer can comprise a single unit configured to predict the similarity of a given pETaG to known ETaGs used to train the model (e.g., train the model to predict the probability that a given pETaG is an actual ETaG). In the case of supervised learning, for example, the training dataset can include feature data of a set of known (or positive) ETaGs to create a set of true positives. The training dataset may also include feature data for a set of negative ETaGs to create a set of true negatives (i.e., negative ETaGs are gene or protein sequences that are not ETaGs). Known examples of positive and negative ETaGs may be identified, for example, from the scientific literature and / or through in-house research. Examples of pETaG feature data that may be used to train a model (or to perform an empirical evaluation of pETaGs, as described below) are summarized in Table 1. [Table 1]

[0168] In some examples, training data including true positives and true negatives can be split into training test datasets (e.g., using a 90 / 10, 80 / 20, 70 / 30, etc. split) and used to first train and then test the deep learning model. In some examples, if the initial dataset is imbalanced, positive or negative weights can be used to balance the true positive and true negative training datasets.

[0169] In one non-limiting example, a deep learning model for predicting the probability that a given pETaG is an actual ETaG was trained using a stochastic gradient descent algorithm and data from an initial dataset with an 80 / 20 train-test split, using the following training parameters: gradient clip max norm = 1.0 (i.e. gradients are "clipped" if their normalized value exceeds the maximum of 1.0) Mini-batch size = 512 (i.e., 512 training data pairs were used per training iteration) Learning rate = 1e-5 (i.e., weights were updated by 1e-5 increments at each training iteration) Max epochs = 3,000 (i.e., the number of complete passes through the training dataset during training) Number of hidden layers: We used models with different numbers of hidden layers and later evaluated their performance. Number of units per hidden layer: We used models with different numbers of units per hidden layer and later evaluated their performance.

[0170] Once a deep learning model has been trained, it can be tested, for example, using a test dataset from the initial training-test dataset split. Examples of test results from models with two hidden layers and various numbers of units per hidden layer are shown below.

[0171] Use of Empirical Scoring to Evaluate pETaGs: In some examples, empirical scoring methods can be used instead of or in addition to deep learning-based methods to evaluate pETaGs. Similar to deep learning-based techniques, the analysis is performed using a pETaG dataset that includes examples of data for positive (known) and negative (known to not exist) ETaGs. A set of positive ETaGs can be identified to create a set of true positives. A set of negative ETaGs (i.e., gene or protein sequences that are known not to be ETaGs) can be identified to create a set of true negatives. The positive and negative ETaGs can be identified, for example, from the scientific literature and / or through in-house research. Candidate pETaGs selected from the positive and negative ETaG datasets can be analyzed using the search methods described above (e.g., target-agnostic search methods) and evaluated using a specified combination of the metrics described above.

[0172] The results from the analysis of positive and negative ETaGs can be used to determine a set of rules, weights, scores, and / or thresholds based on the values ​​of the different metrics summarized in Table 1. The resulting set of rules, weights, scores, and / or thresholds can then be used to develop an algorithm for filtering pETaG results and classifying candidate pETaGs as positive or negative ETaGs.

[0173] The set of rules, weights, scores, and / or thresholds assigned to the different metrics can be determined based on an analysis to identify whether the corresponding values ​​of the candidate pETaGs are more similar to the values ​​of the positive ETaGs or the values ​​of the negative ETaGs. Decision points and metric value thresholds can be based on commonly used statistical measurements and analyses, including, but not limited to, for example, the mean, median, standard deviation, standard error, quartiles, confidence intervals, bootstrap (or any variation thereof), jackknife (or any variation thereof), or any combination thereof.

[0174] Algorithm testing may include testing all combinations of metric values ​​of positive or negative ETaG datasets and corresponding rules, weights, scores, and / or thresholds for their performance in predicting the likelihood that a pETaG is a real ETaG. In some examples, the best algorithm may be selected based on a given target statistic to be maximized. For example, if one seeks to maximize accuracy, a first empirical model may be selected, if one seeks to maximize sensitivity, a second empirical model may be selected, and if one seeks to maximize, for example, an F1 score (a statistical measure of the accuracy of the test), a third empirical model may be selected. The selected model may then be used to, for example, analyze results from a target evaluation for a given candidate pETaG and provide an empirical score that may be used to assess the likelihood or probability (and / or confidence level) that a given pETaG is a true ETaG.

[0175] Returning to FIG. 4, in step 426, a target lookup table is compiled based on deep learning-based and / or empirical pETaG evaluation methods that includes an array of data that maps input values ​​or ranges of an evolutionary metric of a pETaG to output values ​​of the probability that the pETaG is an actual ETaG.

[0176] Grid representation (heat map) analysis method In some examples, the methods and systems described herein relate to identifying genes associated with gene clusters, e.g., biosynthetic gene clusters (BGCs), using machine learning algorithms to evaluate grid representations (e.g., heat maps, or data matrices) of orthologs of genes that co-localize with anchor genes (e.g., core synthase genes) of gene clusters, e.g., BGCs, across diverse genomes.

[0177] The most commonly used bioinformatics tool to identify biosynthetic gene clusters is antiSMASH, which annotates over 40 types of BGCs based on the presence of certain important protein domains. Currently, there are no good methods available to accurately predict the genomic boundaries of gene clusters, e.g., BGCs. antiSMASH defines BGC regions using the following algorithm: In the first step, all gene products of the analyzed sequence are searched against a database of hidden Markov model (HMM) profiles for highly conserved enzymes (e.g., core enzymes) that are indicative of a particular BGC type. In the second step, predefined cluster rules are employed to define individual "clusters" encoded in the analyzed sequence region. Each identified cluster contains a core gene product or a core synthase gene that triggers the cluster rule. antiSMASH defines BGC regions by extending a predefined length, e.g., 20 kb, upstream and downstream of the core synthase gene. The predefined lengths for the different cluster types are empirically determined and generally tend to over-include neighboring genes as part of the BGC. See, e.g., Blink K. et al. (2017), Nucleic Acids Res., vol. 45, pp. W36-W41, and Weber T. et al., antiSMASH5, antiSMASH Database Manual (2019). Thus, genes identified as part of a BGC using antiSMASH or based on proximity to the core synthase genes of the BGC may not be functionally related to the secondary metabolites produced by the BGC.

[0178] To solve this problem, the methods and systems described herein leverage comparative genomics and machine learning algorithms to evaluate heatmaps representing the distribution of orthologs of genes, e.g., bidirectional best hits (BBHs), that co-localize with anchor genes (e.g., core synthase genes) known to be associated with BGCs across multiple diverse genomes to determine the likelihood that an ortholog of a query gene (e.g., a gene in a reference genome) that co-localizes with an anchor gene (e.g., core synthase gene) of a gene cluster (e.g., BGC) in a query (or target) genome is associated with a BGC in the query genome. The machine learning model can be trained using manually curated heatmaps, including manually curated heatmaps representing co-localized non-biosynthetic genes known or experimentally validated to be associated with gene clusters (e.g., BGCs), co-localized non-biosynthetic genes known or experimentally validated to have no functional association with gene clusters (e.g., BGCs), and genes including borderline cases. The trained machine learning algorithm uses sequence information from multiple genomes to enable rapid evaluation of a large number of putative embedded genes that colocalize with anchor genes (e.g., core synthase genes) in gene clusters such as BGCs, greatly improving the accuracy of delineating gene cluster boundaries. Furthermore, the method enables prioritization of putative embedded genes in gene clusters such as BGCs for downstream evaluation, including evaluation by a time-consuming and costly experimental validation process.

[0179] The methods and systems described herein can be used to define gene cluster boundaries, i.e., to identify functionally related genes that are co-localized on a chromosome. Gene clusters or their protein products can be involved in various cellular functions, such as biosynthesis (e.g., secondary and primary metabolism), immunity, cell structure, scavenging, energy, and sensing. In particular, the methods described herein can be used to define the boundaries of gene clusters (e.g., BGCs) in different genomes by identifying genes associated with the gene clusters.

[0180] Furthermore, the methods and systems can be used to identify resistance genes embedded in BGCs that are necessary for the production of secondary metabolites (e.g., genes that confer resistance to the host organism against the action of secondary metabolites produced by the BGC). Identification of resistance genes embedded in BGCs allows for deorphanization of protein targets (ETaG products) of BGC-encoded small molecules that may have homologs in mammalian genomes. Mammalian homologs of ETaG may serve as candidate therapeutic targets, and secondary metabolites may provide small molecule scaffolds for developing modulators against such mammalian homologs.

[0181] FIG. 6 illustrates an exemplary method 600 for generating a grid representation (e.g., a heat map) of genomic data that can be input into a machine learning algorithm (such as an artificial neural network (ANN), a convolutional neural network (CNN), a multi-layer perceptron (MLP), a deep neural network (DNN), a LSTM, a vision transformer model, a generative adversarial network (GAN) model, a variational autoencoder model, a latent diffusion model, etc.) to determine the likelihood that a putative embedded gene (e.g., pETaG) is associated with a gene cluster (e.g., BGC) in the genome.

[0182] FIG. 7 illustrates an exemplary method 700 for determining the likelihood that a putative embedded gene is associated with a gene cluster (e.g., BGC). Process 600 and process 700 are performed, for example, using one or more electronic devices implementing a software platform. In some examples, process 600 and / or process 700 are performed using a client-server system, with blocks of process 600 and / or process 700 being divided in any manner between a server and one or more client devices. In some examples, process 600 and / or process 700 are performed using only one client device or only multiple client devices. In process 600 and / or 700, some blocks are optionally combined, the order of some blocks is optionally changed, and some blocks are optionally omitted. In some examples, additional steps can be performed in combination with process 600 and / or process 700. Thus, the operations shown (and described in more detail below) are exemplary in nature and, therefore, should not be considered limiting.

[0183] Grid Representation In block 702 of FIG. 7, an exemplary system (e.g., comprising one or more electronic devices) receives a grid representation (e.g., a heat map representation) of genomic data including a plurality of cells arranged according to a first axis and a second axis, where the first axis corresponds to a plurality of different genomes (e.g., non-mammalian genomes), the second axis corresponds to a plurality of query genes that co-localize with anchor genes (e.g., core synthase genes) of a gene cluster (e.g., BGC) in a query (or reference) genome, and the putative embedded gene is one of the plurality of query genes. Each cell in the grid representation has a value based on: (i) whether an ortholog of the respective query gene (i.e., the query gene corresponding to the cell) is present or absent in the respective genome (i.e., the genome corresponding to the cell), (ii) the sequence similarity of the ortholog to the respective query gene, and (iii) whether the ortholog of the respective query gene co-localizes with an ortholog of the anchor gene (e.g., core synthase gene) in the respective genome.

[0184] The grid representations described herein can take any of a variety of forms known to those of skill in the art. For example, the grid representation can be a data matrix, such as a two-dimensional data matrix (e.g., a table or array) arranged according to a first axis and a second axis, with values ​​as described herein for each of the cells in the data matrix. In some examples, the grid representation includes one or more matrices (e.g., a table) of data. For example, each set of values ​​of the cells in the grid representation (i.e., (i) whether an ortholog of the respective query gene (i.e., the query gene corresponding to the cell) is present or absent in the respective genome (i.e., the genome corresponding to the cell), (ii) the sequence similarity of the ortholog to the respective query gene, and (iii) whether the ortholog of the respective query gene co-localizes with an ortholog of an anchor gene (e.g., a core synthase gene) in the respective genome), or a combination thereof, can be stored in a separate table and used as input for the machine learning-based methods described herein. In some examples, the grid representation can be a physical representation of the underlying data matrix, such as a heat map, that facilitates visualization of the data.

[0185] For each query gene in a query (or reference) genome, orthologs in any given genome (e.g., target genome) can be identified based on the coding sequence of the query gene or the protein sequence encoded by the query gene, or based on phylogenetic relationships using methods known in the art. For example, an ortholog of a query gene in a given genome can be a gene in the given genome that has the highest sequence similarity to the query gene or encodes a protein whose sequence similarity exceeds a predetermined threshold. Sequence similarity can be quantified by any of a variety of parameters known to those skilled in the art, including percent sequence identity, percent sequence homology, bit score, and e-value. The predetermined threshold can be, for example, at least about any one of 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, 96%, 97%, 98%, 99% or higher percent sequence identity or percent sequence homology.

[0186] In some examples, an ortholog of a query gene in a given genome may be a bidirectional best hit (BBH) of the query gene in the given genome. Methods for identifying BBHs are described, for example, in Moreno-Hagelsieb G, Latimer, K., Bioinformatics. 2008 Feb. 1; vol. 24(3): pp. 319-24. For example, to identify a BBH of a query gene in a given genome, the given genome is first searched for genes that code for proteins with the highest sequence similarity to the protein encoded by the query gene ("putative BBHs"). This search is followed by a reciprocal search, in which the query genome is searched for genes that code for proteins with the highest sequence similarity to the putative BBHs identified in the query genome. If the gene identified in the reciprocal search is the original query gene, then the putative BBH is a true BBH. Alternatively, orthologs of a query gene in a given genome can be identified using the mutual minimum distance method as described, for example, in Wall DP, Deluca T, Methods Mol Biol. 2007, 396:95-110.

[0187] A query gene, including a putative embedded gene in a query (or reference) genome, co-localizes with an anchor gene (e.g., a core synthase gene) of a gene cluster, e.g., a BGC. Two genes may be considered to be co-localized if one gene is within a specified distance or proximity zone of the other. In some examples, the distance between two genes may be considered, for example, as the shortest distance between the genomic coordinates of the two genes. For example, if gene A is present on the + strand and contains a sequence ranging from position 1 to 100, and gene B is present on the - strand and contains a sequence ranging from position 300 to 200 (i.e., position 300 is the start of gene sequence B due to its position on the - strand), then the distance between the two genes is 200-100=100 bp. In some examples, the distance between two genes may be considered as the longest distance between the genomic coordinates of the two genes. In some examples, the distance between two genes may be considered as the distance between the genomic coordinates of the midpoints of the two genes. A putative embedded gene in a target genome co-localizes with an anchor gene (e.g., a core synthase gene) of a BGC in a target genome if the putative embedded gene is within a designated proximal zone relative to an anchor gene (e.g., a core synthase gene) in a BGC of a target genome. In some examples, the proximal zone is about 1-100 kb, e.g., about 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, or 100 kb or less upstream or downstream of an anchor gene (e.g., a core synthase gene) in a BGC. In some examples, the proximal zone is about 1-10 kb, e.g., 1, 2, 3, 4, 5, 6, 7, 8, 9, or 10 kb or less upstream or downstream of an anchor gene (e.g., a core synthase gene) in a gene cluster (e.g., a BGC). In some examples, the proximal zone is 5 kb or less upstream or downstream of a gene. In some examples, the proximal zone is no more than 10 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 15 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 20 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 25 kb upstream or downstream of the gene.In some examples, the proximal zone is no more than 30 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 35 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 40 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 45 kb upstream or downstream of the gene. In some examples, the proximal zone is no more than 50 kb upstream or downstream of the gene.

[0188] A putative BGC may, for example, contain both biosynthetic and non-biosynthetic genes and may further contain pseudogenes (e.g., non-functional segments of DNA that resemble functional genes). Curated libraries of profile hidden Markov models (pHMMs) of biosynthetic domains of BGCs (i.e., probabilistic models that capture the biological diversity of biosynthetic domains based on position-dependent scoring of multiple sequence alignments for corresponding protein or nucleic acid sequences) are known in the art and can be used to identify biosynthetic genes clustered in a genome. The anchor gene may be any one of the signature genes associated with the BGC. For example, the anchor gene may be a core synthase gene in the BGC that is the largest biosynthetic gene in the BGC, or the anchor gene may be the closest core synthase gene to the query gene. Alternatively, the anchor gene may be a non-biosynthetic gene known to be associated with the BGC, for example, a gene encoding a transporter of a secondary metabolite produced by the BGC. In some embodiments, multiple genes in a putative BGC can be identified by extending a window of a predetermined length upstream and downstream of an anchor gene (e.g., a core synthase gene). In some embodiments, multiple genes in a putative BGC of a given genome can be identified using bioinformatics methods such as antiSMASH.

[0189] For example, the grid representation is genome Q1 to Q n The first axis (e.g., Y axis) corresponds to the gene G1 to G mand a second axis (e.g., X-axis) corresponding to the genome Q. i (1≦i≦n) and a cell corresponding to gene Gj(1≦i≦n) has a value and a color selected from a first color, a second color and a third color according to the following: (i)Q i G j If the color does not have a BBH, then the color is the first color and the value is 0, or (ii)Q i G j If you have BBH in the Q i G in j G against BBH j Based on the percentage of sequence identity of (ii-1)G j BBH in the middle is Q i If it colocalizes with the BBH of the core synthase gene in the (ii-2)G j BBH in the middle is Q i If the core synthase genes in the BBH do not colocalize with the BBH, the color is a third color.

[0190] The grid representation may be hierarchically clustered to aid in visualization and manual annotation. For example, the grid representation may be clustered based on pairwise sequence identity or homology between genomes, phylogeny of genomes, or the presence or absence of orthologs corresponding to all query genes in the grid representation. In some examples, the first axis of the grid representation (e.g., heat map) is organized according to the phylogenetic tree of the multiple genomes represented in the grid representation.

[0191] In some examples, the putative embedded gene is a putative embedded target gene (pETaG) that encodes a homolog of a mammalian protein of interest, e.g., a homolog of a human protein of interest. In some examples, pETaG is homologous to an expressed mammalian nucleic acid sequence. In some examples, the mammalian nucleic acid sequence is an expressed mammalian nucleic acid sequence. In some examples, the mammalian nucleic acid sequence is a mammalian gene. In some examples, the mammalian nucleic acid sequence is an expressed mammalian gene. In some examples, the mammalian nucleic acid is a human nucleic acid sequence. In some examples, the human nucleic acid sequence is an expressed human nucleic acid sequence. In some examples, the human nucleic acid sequence is a human gene. In some examples, the human nucleic acid sequence is an expressed human gene.

[0192] An example of a genome heat map is shown in Figure 8A-1 to Figure 8A-3. This example shows a heat map of pETaG in a BGC identified by antiSMASH in a genome marked with an asterisk (*). Each column along the X-axis represents a protein ("protein X") encoded by a query gene in a BGC identified in a query (or reference) genome marked with an asterisk (*). In some examples, the BGC is identified by antiSMASH. In some examples, the gene in the BGC is within a 20 kb proximal zone from the core synthase gene of the BGC (i.e., within ±20 kb of the core synthase gene). The columns corresponding to pETaG and the core synthase gene are indicated by arrows. Each row along the Y-axis represents a unique genome ("genome Y") selected from a genome database. Half of the genomes contain the BBH of the core synthase gene and are referred to as positive genomes for the purpose of identifying genes in putative BGCs. Half of the genomes do not contain the BBH of the core synthase gene and are designated negative genomes for the purpose of identifying genes in putative BGCs. Each cell is colored or shaded according to the presence or absence of the BBH of the respective query gene and the percentage sequence identity (number in the cell) of the BBH to the respective query gene. For example, if the BBH of protein X is absent in genome Y, the cell (X, Y) is blank, if the BBH of protein X is present in genome Y and the BBH is in the same antiSMASH BGC cluster as the BBH of the core synthase gene in genome Y, the cell (X, Y) is e.g. blue or positive, or if the BBH of protein X is present in genome Y and the BBH is not in the same antiSMASH BGC cluster as the BBH of the core synthase gene in genome Y, the cell (X, Y) is e.g. red or negative. The intensity of the red or blue (or grayscale shading) of the cell (X, Y) is based on the percentage sequence identity of the BBH of protein X in genome Y to protein X. The heatmap is hierarchically clustered based on pairwise sequence identity between the genomes.

[0193] Each genome shown in the grid representation may correspond to an assembled genome, or multiple genome fragments obtained from genome sequencing. In some examples, the genome is annotated using a bioinformatics tool such as antiSMASH prior to analysis using any one of the methods described herein. For example, a database of genomes can be constructed such that all putative biosynthetic gene clusters are identified and annotated. For example, genome fragments containing putative BGCs instead of the entire genome can be queried in one or more steps of the methods described herein. For example, the co-localization of a gene (e.g., an ortholog of a query gene) with an anchor gene (e.g., a core synthase gene) can be determined based on the putative BGC annotation in the genome.

[0194] The methods described herein are suitable for any genome that contains a gene cluster, e.g., a BGC. Bacterial, plant, and fungal genomes are known to encode biosynthetic gene clusters. In some embodiments, the query genome and the multiple queried genomes used to generate the grid representation belong to the same kingdom. In some embodiments, the query (or reference) genome and the multiple queried (or target) genomes used to generate the grid representation belong to different kingdoms. Suitable genomes include, but are not limited to, genomes from Archaea, Protozoa, Chromista (e.g., brown algae, diatoms, cryptophytes, etc.), Plantae (e.g., green algae and plants), Fungi, and Animalia. In some embodiments, the query genome and the multiple queried genomes are fungal genomes, such as genomes of different fungal strains. In some embodiments, the query genome and the multiple queried genomes are bacterial genomes, such as genomes of different bacterial strains. In some embodiments, the query genome and the multiple queried genomes are plant genomes, such as genomes of different plant strains. Without wishing to be bound by any theory or hypothesis, fungal genomes are eukaryotic genomes that are phylogenetically related to mammalian genomes rather than to bacterial or plant genomes. Therefore, fungal genomes may be preferred for identifying ETaGs that correspond to human target genes for secondary metabolites produced by ETaG-carrying BGCs.

[0195] At least two genomes are required to construct the grid representation. In some examples, the first axis of the grid representation corresponds to any one of at least about 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 150, 200, 250 or more genomes. In some examples, the first axis of the grid representation corresponds to at least 20 genomes. In some examples, the first axis of the grid representation corresponds to about 50 genomes. Although a large number of genomes may provide more comparative genomics information, it also requires a large amount of computational power and time. Therefore, it may be desirable to sample a representative set of genomes that are diverse with respect to their sequence similarity and / or phylogenetic relationship to generate the grid representation to balance the performance of the method (e.g., accuracy of prediction) and computational resources.

[0196] The genomes represented in the grid representation may include a "positive genome" and a "negative genome." A positive genome is a genome that has an ortholog of an anchor gene, such as a core synthase gene, in the query (or reference) genome. A negative genome is a genome that does not have an ortholog of an anchor gene, such as a core synthase gene, in the query (or reference) genome. In some examples, the positive genome and the negative genome are selected from a database of genomes. In some examples, the multiple genomes used to construct the grid representation include multiple positive genomes that each have an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene) and multiple negative genomes that do not have an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene). In some examples, the positive genome is selected from a genome database by identifying a genome that has an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene), and the negative genome is selected from a genome database by identifying a genome that does not have an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene). The negative genomes may be phylogenetically adjacent to the selected positive genome. In some embodiments, the number of positive genomes and the number of negative genomes are equal to each other.

[0197] The positive genome and the negative genome can be selected from a database with a large number of genomes. For example, the database can contain at least 2, 10, 100, 500, 1000, 5000, 10000, 15000, 20000, 25000, 30000, 35000, 40000, 45000, 50000, 100000, 200000, 500000, 1000000, or more than 1000000 genomes. The selection of the positive genome and the negative genome from a large genome database may require clustering of the genome to allow sampling of a variety of genomes including the positive genome and the negative genome from the database. For example, in some examples, the database genomes can be clustered according to the sequence similarity of one or more single copy genes in the genome, or the sequence similarity of the orthologs of the anchor gene (e.g., core synthase gene), or the putative embedded gene (e.g., pETaG) in the genome. Clustering may be performed using unsupervised clustering methods. Unsupervised clustering methods may include, for example, the use of Markov Cluster Algorithm (MCL), Restricted Neighborhood Search Cluster (RNSC) algorithm, affinity propagation clustering algorithm, spectral clustering algorithm, k-means clustering algorithm, or any other method known in the art. Alternatively, clustering may include the use of supervised clustering methods known in the art, such as supervised k-means clustering or semi-supervised spectral clustering. The threshold for clustering may be determined by a predetermined goal for the number of clusters. For example, the threshold for clustering may be a predetermined sequence similarity level between genome groups, for example, requiring that the sequence similarity between different genome groups is equal to or less than about any one of 99.5%, 99%, 98%, 95%, 90%, 85%, 80%, 75%, 70%, 65%, 60%, 50%, 40%, 30%, or less.In some embodiments, it may be desirable to select a positive genome in which the pairwise sequence similarity (e.g., sequence identity) percentage of the orthologs of one or more single-copy genes in the positive genome exceeds any one of about 99.5%, 99%, 98%, 97%, 96%, 95%, 94%, 93%, 92%, 91%, 90%, 85%, 80%, 70%, 60%, 50%, 40%, 30%, or less. In some embodiments, it may be desirable to select a negative genome in which the pairwise sequence similarity (e.g., sequence identity) percentage of the orthologs of one or more single-copy genes in the negative genome exceeds any one of about 99.5%, 99%, 98%, 97%, 96%, 95%, 94%, 93%, 92%, 91%, 90%, 85%, 80%, 70%, 60%, 50%, 40%, 30%, or less. Representative genomes from each cluster may be further selected for use in the analysis steps described herein. In some embodiments, negative genomes are selected from the database by identifying genomes that have the highest sequence similarity to the positive genome but lack an ortholog of an anchor gene (e.g., a core synthase gene).

[0198] For example, the grid representation can be constructed based on a number n of positive genomes selected from a database of m genomes. As a first step, from the m genomes in the database, a number m1 of positive genomes that have an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene) and a number (m-m1) of negative genomes that do not have an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene) are identified. In a situation where m1 is greater than n, the m1 positive genomes are clustered into n clusters using an MCL based on the average sequence similarity of one or more single copy genes between the positive genomes (e.g., identified using the BUSCO tool, see busco.ezlab.org). Then, one positive genome is selected from each of the n clusters to provide n positive genomes for constructing the grid representation. Each of the n negative genomes is selected by identifying the genome that is most similar to the selected positive genome among the (m-m1) negative genomes (e.g., the genome with the highest sequence similarity or the genome with the shortest phylogenetic distance). This genome selection method results in a grid representation constructed from n positive genomes and n negative genomes.

[0199] In the situation where m1 is smaller than n, more negative genomes than positive genomes may be selected to construct a grid representation with a total of 2n genomes. In this case, the (m-m1) negative genomes may be clustered into (2n-m1) clusters, and one negative genome is selected from each cluster. Alternatively, for each of the m1 positive genomes, two or more negative genomes that are closely related to the positive genome are selected based on their sequence similarity or phylogenetic distance to the positive genome, resulting in a total of 2n-m1 negative genomes being selected.

[0200] The grid representation can be generated, for example, using the method illustrated in FIG.

[0201] As a first optional step, resource files can be prepared for downstream computational processes. Exemplary resource files include pairwise genome comparison files, files (e.g., FASTA files) containing related proteins or genes from the target genome (i.e., the genome selected for analysis), and, optionally, resource files containing clusters of orthologous groups (COGs) of proteins or genes.

[0202] A pairwise genome comparison file is a file created to show the homology relationships between all pairs of genomes in a database. Genome similarity can be determined based on either pairwise genome sequence similarity or pairwise phylogenetic distance between genomes.

[0203] In some examples, pairwise identity or similarity between genomes can be determined by comparing the entire genome sequence or by comparing the sequences of a subset of proteins or genes. For example, to determine the entire genome sequence identity, the entire genome can be aligned and the pairwise identity between the alignments is calculated. Alternatively, pairwise genome identity can be calculated by comparing single copy proteins (i.e., amino acid sequences) or genes (i.e., nucleotide sequences) shared between pairs of genomes. In some preferred embodiments, single copy proteins or genes are used as duplicated or fragmented proteins, which may provide erroneous estimates of genome homology. Single copy proteins or genes in genomes can be annotated using BUSCO (doi.org / 10.1093 / molbev / msab199), or specific known single copy proteins or genes can be used to determine genome sequence similarity. In some examples, single copy proteins can be identified using known bioinformatics tools such as OrthoMCL, OrthoFinder, or PanX. A subset or all single copy proteins or genes shared between genomes are aligned separately, trimmed, and concatenated to form a super-alignment. Pairwise identity is the number of identical residues in the super-alignment. Alternatively, similarity scores can be calculated based on alternative matrices such as BLOSUM and PAM when protein sequences are used to determine sequence similarity.

[0204] In other examples, genome similarity is determined based on the phylogenetic distance between genomes. To determine phylogenetic distance, a set of single copy proteins (i.e., amino acid sequences) or genes (i.e., nucleotide sequences) can be individually aligned using any alignment software such as MAFFT, MUSCLE, or ClustalW, trimmed using any sequence trimming software such as trimAI, GBlocks, or ClipKIT, and concatenated to create a superalignment. The superalignment can be used by any phylogenetic tree building software such as FastTree, IQ-TREE, RAxML, MEGA, MrBayes, BEAST, or PAUP to provide a phylogenetic tree of the genomes. The tree can be constructed using different algorithms such as maximum likelihood algorithm, maximum parsimony algorithm, neighbor-joining algorithm, distance matrix algorithm, or Bayesian inference algorithm. Alternatively, instead of the superalignment approach, a gene-coalescent phylogenetic model approach can be used to reconstruct the phylogeny.

[0205] In some examples, a FASTA file containing all protein or gene sequences from the target genome is created as an input resource file. Alternatively, a smaller subset of proteins, such as proteins in putative gene clusters (e.g., putative BGCs) based on bioinformatics predictions, can be provided as input instead of all proteins from all genomes. Putative BGCs can be predicted using publicly available BGC prediction tools, such as, for example, antiSMASH, SMURF, TOUCAN, deepBGC, or using custom BGC prediction tools. The FASTA file can contain either protein or nucleic acid sequences, since sequence similarity (e.g., homology) can be determined using either protein or nucleic acid sequences.

[0206] Gene clusters (e.g., BGCs) protein or gene clusters of orthologous groups (COGs) may be provided as optional resource files. Members of the same COG are predicted to have orthologous functions. COGs can be created using protein clustering tools such as USEARCH, CD-HIT and MMseqs. Alternatively, COGs can be generated using clustering algorithms such as MCL or SiLiX after sequence (e.g., amino acid or nucleotide) alignment searches using BLAST, ggsearch, or Diamond. COGs can also be identified using custom-developed scripts or using known bioinformatics tools such as OrthoMCL, OrthoFinder, or PanX. COGs can also be defined as gene families or protein families. A COG file can contain, for example, protein or gene IDs from the same COG in each row of a table.

[0207] In block 602 of FIG. 6, the method for generating a grid representation includes identifying a putative gene cluster (e.g., a putative biosynthetic gene cluster (BGC)) containing a putative embedded gene in a query (or reference) genome from a plurality of genomes, the putative gene cluster containing an anchor gene (e.g., a core synthase gene) known to be associated with the gene cluster (e.g., BGC), the anchor gene co-localizing with the putative embedded gene. In some examples, the method includes identifying the longest biosynthetic or structural gene in the putative gene cluster as the anchor gene. A resource file can be used to perform the steps of block 602. This step establishes a first axis of the grid representation, e.g., the X-axis of a heat map, corresponding to the plurality of query genes that co-localize with the putative gene cluster (e.g., BGC) in the query genome and contain the putative embedded gene. The co-localization can be determined based on the putative BGC annotation or based on the distance between the two genes, e.g., the two genes within a specified proximity zone of about 50 kb or less, or about 20 kb or less.

[0208] For example, the X-axis of the heatmap can be established based on a single protein or gene ID of interest (e.g., the gene ID corresponding to pETaG) or based on multiple protein or gene IDs of interest as input. If multiple protein or gene IDs are input, the correlation and co-localization of these genes (along with adjacent genes) to each other across multiple genomes is determined. For example, the multiple protein or gene IDs of interest may correspond to ETaG and core synthase genes in the query (or reference) genome. If a single protein or gene ID is used as input, it can be determined whether the adjacent genes surrounding it co-localize across multiple genomes. The X-axis corresponds to the search region that contains the adjacent genes of the protein or gene of interest (e.g., pETaG). The search region can be defined by predicted BGC or based on coordinate positions on the genome. Putative BGCs can be predicted using tools such as antiSMASH, SMURF (dx.doi.org / 10.1016 / j.fgb.2010.06.003), TOUCAN (doi.org / 10.1093 / nargab / lqaa098), deepBGC (doi.org / 10.1093 / nar / gkz654), or other custom search algorithms.

[0209] If a predefined region is used as input, it is assumed that the specified input protein or corresponding gene ID is present in that region. Alternatively, a custom neighbor distance (i.e., proximal zone distance), specified, for example, in units of base pairs (bp), can be used to identify a genomic window region that contains a certain number of bp upstream and downstream adjacent to either side of the specified protein or gene ID. All proteins (or genes) located in the defined region are assigned as the X-axis protein (or gene). The proteins (or genes) in the region are labeled. If the input protein ID is a single protein, the input label is used as the label for the protein ID. For example, a core synthase gene can be input and genes that co-localize with core synthase can be determined. If the input protein ID is a single protein and it is desired to determine its correlation with another gene (for which no ID was input), another specified input can be used to identify genes in the region that will be labeled as such. For example, if a gene of interest corresponding to ETaG is used as input and it is desired to determine whether it correlates with core synthase, a search for core synthase can be used to identify and the identified protein can be labeled as such. Proteins can be searched based on gene annotation. When there are multiple target proteins in a region matching the requested search criteria, several options are available. For example, a heatmap for each target protein can be created. Target proteins can be selected based on the length or proximity of the target protein to the input protein ID. When multiple protein or gene IDs are used as input, the protein IDs are labeled based on their input labels. For example, one input protein can be labeled as an ETaG and the other as a core synthase.

[0210] In block 604 of Figure 6, the method includes obtaining a plurality of positive genomes that include an ortholog of the core synthase gene and a plurality of negative genomes that do not include an ortholog of the core synthase gene, where the plurality of negative genomes are selected based on sequence similarity or phylogenetic distance to the plurality of positive genomes. The resource file can be used to perform the step of block 404. This step establishes a second axis of the grid representation, e.g., the Y axis of a heat map, that corresponds to the target genome that includes the plurality of positive genomes and the plurality of negative genomes.

[0211] For example, to establish the Y-axis of the heatmap, positive and negative genome IDs can be obtained as follows: To establish positive genome IDs, protein homologs of the core synthase IDs are searched for by running a protein sequence alignment tool such as ggsearch, BLASTp or Diamond-blastp against the set of genomes. Alternatively, the gene (i.e., nucleotide sequence) can be used to find homologs of the gene of the core synthase in the set of genomes using tools such as ggsearch or BLASTn. Protein IDs that fall within a specified range of minimum and maximum sequence identity are identified. A specified number of protein homologs can be selected from those with the highest sequence identity to the query core synthase. Core synthase homologs can be de-replicated by selecting representatives from protein clusters created by using a protein clustering tool with a specified cutoff. Alternatively, de-replicating can be done using BUSCO pairwise cutoff, phylogenetic distance cutoff, or taxonomic classification. If a phylogenetic tree is used as the genome homology resource file, it can be traversed to select a diverse set of positive genomes that meet the protein homolog presence criteria. Alternatively, positive genomes can be selected based on pairwise BUSCO identity cutoffs, phylogenetic distance or clade, or taxonomic classification (e.g., one isolate per species from genomes with protein homologs or one species per genus or family). These methods can be utilized to help ensure the selection of a diverse set of genomes, which provides greater accuracy in identifying co-localized genes and avoids confounding results with multiple genomes of the same species, which can bias results toward co-localization and increase false positive rates. Obtain each genome ID from the selected core synthase homologs and assign them as positive genomes.

[0212] To obtain the negative genome IDs, the pairwise genome homology file can be used as input to select the genome ID with the highest sequence identity or closest phylogenetic distance to each positive genome. If a genome contains a core synthase homolog within the specified sequence identity range, it is removed from the list of candidates and the search skips to the next candidate genome. If a phylogenetic tree is used as the homology resource file, the tree can be traversed to find the closest negative genome to each positive genome. The positive and negative genome IDs are combined and assigned as the Y-axis genome ID.

[0213] Optionally, a file containing all protein (or nucleotide) sequences and gene annotation files (GFF, GTF, GenBank, or similar) related to the Y-axis genome is obtained. When gene cluster (e.g., BGC) prediction is used to define the search region, gene cluster (e.g., BGC) information of the Y-axis genome, including gene cluster (e.g., BGC) ID, cluster number, and protein (or gene) ID located in the gene cluster (e.g., BGC), is stored as a resource file.

[0214] In some examples, a linkage matrix (also known as a cladogram) is constructed to visually show the distance between all genomes on the Y-axis. The linkage matrix can be created using a pairwise homology matrix of the Y-axis genomes from a genome homology resource file (e.g., pairwise identity, similarity, or phylogenetic distance). A hierarchical clustering method can be used with the pairwise homology matrix to create the linkage matrix. Alternatively, the BBH presence / absence or forward alignment results of the X-axis protein homologs can be used with a hierarchical clustering method to create the linkage matrix. Alternatively, a phylogenetic tree can be used. The phylogenetic tree can be created from either a set of proteins (i.e., amino acid sequences) or transcripts (i.e., nucleotide sequences).

[0215] In block 606 of FIG. 6, the method includes creating a grid representation including a plurality of cells arranged according to a first axis and a second axis, where the first axis corresponds to all protein-coding genes (e.g., the query gene) that co-localize with an anchor gene in a putative gene cluster (e.g., a putative BGC) in the query (or reference) genome, and the second axis corresponds to a plurality of positive genomes and a plurality of negative genomes, and each cell is based on (1) the presence or absence of an ortholog of the respective protein-coding gene in the respective genome, (2) the sequence similarity of the ortholog to the respective protein-coding gene, and (3) whether the ortholog of the respective protein-coding gene co-localizes with an ortholog of the anchor gene in the respective genome.

[0216] For example, the following steps can be used to create a heatmap matrix. First, obtain bidirectional best hit (BBH) results for the X-axis proteins (or genes) in all Y-axis genomes to provide a BBH table. BBHs are identified as pairs of proteins (or genes) from two different genomes that are more similar to each other than either one is to any gene in the other genome. BBHs can be useful for identifying true orthologs of a given protein (or gene), which is particularly useful for identifying current orthologs of genes (such as ETaGs) that have had duplication events. Identifying BBHs involves a forward alignment step and a reverse alignment step. In the forward alignment step, a sequence alignment tool is run using each X-axis protein as a query against the protein FASTA file of each Y-axis genome with a specified cutoff. Sequence alignment tools such as ggsearch, BLASTp or Diamond-blastp can be used. Alternatively, genes or transcripts (i.e., nucleotide sequences) can be used in place of protein sequences for the forward alignment step using tools such as ggsearch or BLASTn. The protein (or gene) IDs of the best matches from each alignment are stored in a table with the X-axis protein (or gene) as column and the Y-axis genome as index. The sequence identities of the best matches from each alignment are stored in a table with the X-axis protein as column and the Y-axis genome as index. In the reverse alignment step, reverse alignment is performed using a protein sequence alignment tool such as ggsearch, BLASTp or Diamond-blastp. Each protein stored as the best hit from the forward alignment step is used as a query protein (or query gene). The protein alignment is performed against the protein FASTA file of the Y-axis genome. Alternatively, genes or transcripts (i.e., nucleotide sequences) can be used in place of protein sequences for reverse alignment using tools such as ggsearch or BLASTn.If the best hit from the reverse alignment is the same as the query protein used in the forward alignment, the X-axis protein (or gene) and its forward alignment hit are BBH and the BBH value in the table is defined as true. If the best hit from the reverse alignment is different from the query protein used in the forward alignment, the X-axis protein (or gene) and its forward alignment hit are not BBH and the BBH value in the table is defined as false. This binary data is stored in the BBH table with the X-axis protein as the column and the Y-axis genome as the index. Alternatively, only the forward alignment results are used instead of the complete BBH results.

[0217] Also, a colocalization table is created. For example, when the search region is defined using BGC prediction, the cluster number is obtained and the cluster number of each forward alignment hit from the table containing the BGC information of each genome is obtained. The values ​​are stored in the table with the X-axis protein as the column and the Y-axis genome as the index. If the cluster number of each forward alignment hit is the same as the core synthase homolog of the genome, the forward alignment hit colocalizes with the core synthase homolog and the colocalization value is defined as true. Otherwise, the forward alignment hit does not colocalize with the core synthase homolog and the colocalization value is defined as false. This binary information of colocalization is stored in the colocalization table.

[0218] Alternatively, if a custom neighbor distance (i.e., proximity zone distance) is used to define the search region, the genomic location is obtained and a table containing the genomic location of each forward alignment hit is created. The scaffold ID and the coordinates of the start and end positions of each protein are stored. A colocalization table is created containing binary information of the colocalization between each forward alignment hit and the core synthase homolog of the corresponding genome. The colocalization value is defined as true if the query protein and the core synthase homolog are located within the specified distance (i.e., proximity zone) in the same scaffold.

[0219] A final table that stores the values ​​of the cells of the heatmap is created based on the BBH and colocalization tables as described above. The following conversions are applied: convert true to 1 and false to 0 from the binary data of the table containing BBH information; and convert true to 1 and false to -1 from the binary data of the table containing colocalization information. To calculate the value of each cell of the heatmap, the multiple sequence identity, BBH information, and colocalization information value of all combinations of X-axis proteins and Y-axis genomes are calculated. For example, a cell may have a value of sequence identity (96.28)*BBH (1 or 0)*colocalization (1 or -1). If there is no BBH, the cell value will be 0. If the cell corresponds to a gene that is not colocalized with the core synthase gene, the cell value will be negative sequence identity. For heatmaps based on forward alignment hits, the value of each cell is calculated as the product of sequence identity and colocalization information value of all combinations of X-axis proteins and Y-axis genomes. The heatmap is plotted using the final table and the linkage matrix. A divergent color map can be used to visualize values ​​ranging from -100 to 100. The Y-axis can be sorted based on a hierarchical clustering of the connectivity matrix.

[0220] The generated grid representation (e.g., a heatmap, or a data matrix) and / or a subset of the generated grid representation can be input into a machine learning model, such as an LSTM or CNN, to provide the likelihood that a putative embedded gene (e.g., pETaG) is associated with a BGC.

[0221] A machine learning model for performing embeddedness classification In block 704 of Figure 7, the system inputs the grid representation into a machine learning model, which is trained to determine the likelihood that a putative embedded gene is associated with a gene cluster (e.g., BGC) based on the values ​​of a number of cells in the grid representation. Non-limiting examples of machine learning models that may be used include, but are not limited to, artificial neural networks (ANNs), convolutional neural networks (CNNs) or other recurrent neural networks (RNNs), multi-layer perceptrons (MLPs), deep neural networks (DNNs), LSTMs, vision transformer models, generative adversarial networks (GANs), variational autoencoders, latent diffusion models, and the like.

[0222] The likelihood may be the probability that the putative embedded gene is associated with a gene cluster (e.g., BGC), or the probability that the putative embedded gene falls into one of a number of predefined likelihood categories. In some examples, the likelihood may fall into one of the following four categories: (1) the putative embedded gene is highly likely to be associated with a gene cluster (e.g., BGC) (also referred to herein as "tier A+"); (2) the putative embedded gene is moderately likely to be associated with a gene cluster (e.g., BGC) (also referred to herein as "tier 1"); (3) the putative embedded gene is moderately likely to be associated with a gene cluster (e.g., BGC) (also referred to herein as "tier 2"); and (4) the putative embedded gene is low likely to be associated with a gene cluster (e.g., BGC) (also referred to herein as "tier 3"). A tier A+ heatmap has clearly defined gene cluster (e.g., BGC) boundaries, and the putative embedded gene (e.g., putative resistance gene or pETaG) is within those boundaries. See, for example, Figures 11A-1 and 11A-2. Tier 1 heat maps have poorly defined gene cluster (e.g., BGC) boundaries, but putative embedded genes (e.g., putative resistance genes or pETaG) are correlated with anchor genes (e.g., core synthase genes). See, for example, Figures 11B-1 and 11B-2. Tier 2 heat maps tend to provide insufficient information to identify gene cluster (e.g., BGC) boundaries, or allow for the acceptance or rejection of correlations between putative embedded genes (e.g., putative resistance genes or pETaG) and anchor genes (e.g., core synthase genes). See, for example, Figures 11C-1 and 11C-2.The tier 3 heatmap shows that the putative embedded gene is a false positive because (1) the boundaries of the gene cluster (e.g., BGC) are well defined and the putative embedded gene (e.g., putative resistance gene or pETaG) is not within the boundaries, or (2) there is no correlation or co-localization between the putative embedded gene (e.g., putative resistance gene or pETaG) and the anchor gene (e.g., core synthase gene). See, e.g., Figure 11D-1 and Figure 11D-2. The method can output a probability associated with each of these likelihood categories. For example, the method can provide a probability associated with each of four categories, where the sum of the four probabilities is equal to 100%.

[0223] In block 706 of FIG. 7, the system obtains the likelihood that the putative embedded gene is associated with a gene cluster (e.g., a BGC) from the machine learning model.

[0224] In block 708 of FIG. 7, the system displays on the display a grid representation (eg, a heat map) and the likelihood that the putative embedded genes are associated with a gene cluster (eg, a BGC).

[0225] 8A-1-8A-3 show heatmaps depicting the degree of "embeddedness" (i.e., association with BGCs) of pETaG in a calculated heatmap, according to an example of the present disclosure. In the heatmaps, putative BGCs are identified by antiSMASH in genomes marked with an asterisk (*). Each column along the X-axis represents a protein ("protein X") encoded by a query gene in a BGC identified in a query genome marked with an asterisk (*). Genes in a BGC can be identified by antiSMASH or based on the proximity (e.g., within 20 kb) of the BGC to a core synthase gene. Columns corresponding to pETaG and the core synthase gene are indicated by arrows. Each row along the Y-axis represents a unique genome ("genome Y") selected from a genome database. Half of the genomes contain a BBH of the anchor gene and are referred to herein as "positive genomes." Half of the genomes do not contain a BBH of the anchor gene and are referred to herein as "negative genomes." Each cell is colored or shaded according to the presence or absence of the BBH of the respective query gene and the percentage sequence identity (number in the cell) of the BBH to the respective query gene. For example, if the BBH of protein X is absent in genome Y, the cell (X, Y) is blank, if the BBH of protein X is present in genome Y and the BBH is in the same antiSMASH BGC cluster as the BBH of the core synthase gene in genome Y, the cell (X, Y) is, for example, blue or positive, or if the BBH of protein X is present in genome Y and the BBH is not in the same antiSMASH BGC cluster as the BBH of the core synthase gene in genome Y, or if the BBH of the core synthase gene is absent in genome Y, the cell (X, Y) is, for example, red or negative. The intensity of the red or blue (or grayscale shading) of the cell (X, Y) is based on the percentage sequence identity of the BBH of protein X in genome Y to protein X. The heatmap is hierarchically clustered based on pairwise sequence identity between genomes.

[0226] 8B-1 and 8B-2 and 8C-1 and 8C-2 show an exemplary long short-term memory (LSTM) model used to classify an input heatmap into one of four hierarchies: (1) high likelihood that the putative embedded target gene (pETaG) is associated with a gene cluster (e.g., BGC) ("Tier A+"), (2) moderately high likelihood that the pETaG is associated with a gene cluster (e.g., BGC) ("Tier 1"), (3) moderately low likelihood that the pETaG is associated with a gene cluster (e.g., BGC) ("Tier 2"), and (4) low likelihood that the pETaG is associated with a gene cluster (e.g., BGC) ("Tier 3"). In some embodiments, the heatmap and / or one or more subsections thereof may be used as input for the LSTM.

[0227] As shown in FIG. 8B-1 and FIG. 8B-2, the values ​​of the cells in each column of the heatmap are sequentialized into a plurality of input arrays, each corresponding to a query gene across multiple genomes, and each of the plurality of input arrays is input to an LSTM cell. Each input array also holds location information for pETaG and core synthase, represented by two scalars, with 1 and 0 corresponding to the presence and absence of a particular gene (pETaG or core synthase). The plurality of LSTM cells provide an output hierarchy. In some examples, after computing the heatmap, a vector representation of the data contained in the heatmap (e.g., one or more patterns, one or more colors, a table of values ​​indicating pETaG locations, etc.) may be provided to one or more neural networks to perform an embedding classification of the pETaGs in the heatmap. For example, in some examples, the one or more neural networks may include, for example, a long short-term memory (LSTM) model, a convolutional neural network, or other recurrent neural network (RNN), which may be suitable for processing genomic data or other text-based data represented as an array of elements. For example, in one embodiment, the LSTM model can receive as input a linear vector representation of the heatmap data (e.g., values ​​indicative of one or more patterns, one or more colors, pETaG location, etc.) and output an embeddedness classification (e.g., a probability value of how pETaG is embedded in a cluster of genes). In some examples, the LSTM model performs embeddedness classification by classifying pETaGs into one of four classes: (1) "True Positive" (e.g., there is a high likelihood that pETaG is associated with BGCs ("Tier A+")); (2) "Probable" (e.g., there is a moderately high likelihood that pETaG is associated with BGCs ("Tier 1")); (3) "Uncertain" (e.g., there is a moderately low likelihood that pETaG is associated with BGCs ("Tier 2")); (4) "True Negative" (e.g., there is a low likelihood that pETaG is associated with BGCs ("Tier 3")).

[0228] 8B-1 and 8B-2 and 8C-1 and 8C-2 show one or more example implementations of an LSTM model that outputs an embedding classification (e.g., a probability value of how pETaGs are embedded in a cluster of genes) based on an input linear vector representation of heat map data (e.g., values ​​indicative of one or more patterns, one or more colors, pETaG location, etc.), according to examples of the present disclosure. For example, in some examples, as shown by FIG. 8B-1 and 8B-2, an LSTM model can include a series of memory hierarchies, each including, e.g., a respective memory cell. In some examples, each memory cell can include, e.g., its cell state (e.g., C t-1 ~C t ), the LSTM model may include the ability to remove or add information to the cell state, carefully regulated by structures called gates. In some examples, gates in each memory cell may be provided to optionally pass information through. For example, in some examples, each memory cell may include a sigmoid neural net layer and a pointwise multiplication operation. The sigmoid layer outputs a number between "0" and "1" that describes how much of each component's data should be passed through. For example, in one embodiment, a value of "0" means "do not pass data" and a value of "1" means "pass data". In some examples, each respective memory cell may include these gates to protect and control the cell state, for example.

[0229] In some examples, during operation, each layer and memory cell of the LSTM model can begin by determining which of the heat map data to discard from the cell state. For example, in some examples, the determination may be performed by a sigmoid layer called a forget gate layer, which looks at the input data and discards each number (e.g., C t-1 ~C t), where "1" represents keeping this data entirely and "0" represents discarding this data entirely. Each layer and memory cell of the LSTM model can then decide what new information to store in the cell state. For example, in some instances, a sigmoid layer, called the input gate layer, decides which value to update, and a tan h layer creates a vector of new candidate values ​​that can be added to the state.

[0230] Then, each layer and memory cell of the LSTM model is assigned the old cell state C t-1 New cell state C t The LSTM model can then determine what to output based on each layer of the LSTM model and each cell state of the memory cell. For example, in some instances, a sigmoid layer determines what portion of the cell state to output, and then that cell state goes through a tan h layer to set the value between "-1" and "+1" and multiply the value by the output of a sigmoid gate.

[0231] As further illustrated by the prediction tables in Figures 9A and 9B according to an example of the present disclosure, Figures 8C-1 and 8C-2 show example outputs of an LSTM model that can represent the embeddability classification of pETaG into one of four classes as described above: (1) "true positive" (e.g., there is a high likelihood that pETaG is associated with BGCs ("Tier A+")), (2) "probable" (e.g., there is a moderately high likelihood that pETaG is associated with BGCs ("Tier 1")), (3) "uncertain" (e.g., there is a moderately low likelihood that pETaG is associated with BGCs ("Tier 2")); (4) "true negative" (e.g., there is a low likelihood that pETaG is associated with BGCs ("Tier 3")).

[0232] Specifically, FIG. 9A shows a table of predictive embeddedness benchmark values ​​for “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”. Similarly, FIG. 9B shows a table of final predictive embeddedness benchmark values ​​for “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”, including a positive predictive value (i.e., precision) of 48.91%, a negative predictive value of 99.49%, a sensitivity value of 91.82%, and a specificity value (i.e., recall) of 94.28%. The values ​​in the table in FIG. 9A are calculated from comparing the manually annotated hierarchies with the prediction results from the LSTM model. In the table in FIG. 9B, the results of Tier A+ and Tier 1 are combined, and the results of Tier 2 and Tier 3 are combined. Sensitivity is calculated as the true predicted positives divided by the sum of actual positives. Specificity is calculated as the true predicted negatives divided by the sum of actual negatives. Positive predictive value is calculated as the true predicted positives divided by the sum of predicted positives. Negative predictive value is calculated as the true predicted negatives divided by the sum of the predicted negatives.

[0233] Figures 10A and 10B show an example of a heatmap depicting a BGC for lovastatin production predicted by antiSMASH in Aspergillus terreus. Using the methods described herein, a smaller set of genes is identified as associated with the BGC. Figures 11A-1 to 11D-2 show the heatmaps classified into Tier A+, Tier 1, Tier 2, and Tier 3, respectively.

[0234] In some examples, based on an LSTM model that outputs an embeddedness classification (e.g., a probability value of how pETaG is embedded in a cluster of genes), the output of the LSTM model may be organized into one or more tables (as shown in FIG. 12A ) representing the four features “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”, and combined with a predetermined set of additional features (e.g., up to 27 or more features including the four features “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”).

[0235] Machine learning model for predicting ETaG likelihood FIG. 12A illustrates a data table including a combination set of features. In some examples, the combination set of features (e.g., up to 27 or more features) can be organized in a data table as shown in FIG. 12A and can be utilized to train a machine learning model (e.g., an artificial neural network (ANN), a multi-layer perceptron (MLP), a deep neural network (DNN), a convolutional neural network, etc.) to output an ETaG or pETaG probability value based on an input of the combination set of features shown in FIG. 12A. In some examples, the data can be utilized to train other types of machine learning models (e.g., Bayesian inference, decision tree-based methods such as XGBoost or random forests, etc.). In some examples, the data can be utilized to train a logistic regression model or other types of supervised models. As further shown in FIG. 12A, the data table of the combination set of features can also include annotated datasets of known ETaG or pETaG label values ​​that can be utilized, for example, as ground truth or other references for training the machine learning model.

[0236] FIG. 12B illustrates an initial training stage of a neural network trained to output ETaG or pETaG probability values, according to an example of the present disclosure. As shown, a training data set corresponding to the features included in the data table of FIG. 12A may be input to an input layer of the neural network. For structured or semi-structured input data, in one embodiment, the neural network may include a multi-layer perceptron (MLP) or other layered neural network including at least one hidden layer. For example, during training, a table of features may be input to each neuron or node. Specifically, in some examples, each neuron or node may take the table of features as input and execute one or more specified activation functions (e.g., computational functions) to generate an output based thereon. For example, in some examples, the specified activation functions (e.g., computational functions) may specifically determine the value of the output of the input neuron or node.

[0237] In some examples, each input neuron or node can be connected to a set of hidden neurons or nodes that can receive the output of the input neuron or node. In some examples, the hidden neurons or nodes can constitute a hidden layer of the neural network, e.g., each can include a weight that determines the relative strength (e.g., positive or negative) of the input of each connection to the input neuron or node. For example, in some examples, the weights of the hidden layer can affect, e.g., the effect that each input has on the hidden neuron or node, and can be iteratively adjusted for the neural network to learn over time. In some examples, as further illustrated by FIG. 12B, the neural network can be trained, e.g., based on a forward propagation technique. In some examples, the combined output of the hidden neuron or node through forward propagation and extension can include a weighted sum of the output of the hidden neuron or node and a predicted ETaG or pETaG label value (e.g., compared to a ground truth ETaG or pETaG label value).

[0238] FIG. 12C further illustrates a training stage of a neural network trained to output ETaG or pETaG probability values ​​according to an example of the present disclosure. For example, as shown in FIG. 12C, a neural network can be evaluated by utilizing a loss or cost function (e.g., supervised learning) to compare predicted ETaG or pETaG label values ​​with ground truth ETaG or pETaG label values ​​to calculate a loss (e.g., error). In some examples, the weights of hidden neurons or nodes can be iteratively adjusted to minimize the loss (e.g., calculated based on a comparison of predicted ETaG or pETaG label values ​​with ground truth ETaG or pETaG labels) to the extent that the neural network is properly and accurately trained.

[0239] FIG. 12D illustrates an inference stage of a neural network trained to output ETaG or pETaG probability values, according to an example of the present disclosure. As shown, an unknown data set of features, e.g., corresponding to one or more features included in the data table of FIG. 12A, may be input to an input layer of the neural network. For example, during inference, the table of features may be input to each neuron or node. Specifically, as described above with respect to FIG. 12B, each neuron or node may take the table of features as input and execute one or more specified activation functions (e.g., computational functions) to generate an output based thereon. For example, in some examples, the specified activation functions (e.g., computational functions) may specifically determine the value of the output of the input neuron or node. In some examples, each input neuron or node may be connected to a set of hidden neurons or nodes that may receive the output of the input neuron or node. In some examples, the hidden neurons or nodes may each include a weight that determines, for example, the relative strength (e.g., positive or negative) of the input of each connection to the input neuron or node. In some examples, as further shown in FIG. 12D, the combined output of the hidden neuron or node may include a weighted sum of the output of the hidden neuron or node and a predicted ETaG or pETaG label value (e.g., compared to a ground truth ETaG or pETaG label value). Specifically, according to examples of the present disclosure, the neural network may output an ETaG or pETaG probability value based on an unknown input data set of features. FIG. 12E shows an example data table including an unknown input data set of features and a corresponding output of ETaG or pETaG probability values ​​for each of the features, according to examples of the present disclosure. In this manner, the examples may identify and determine the likelihood that an ETaG or pETaG is associated with one or more features that represent a BGC.

[0240] Purpose The computer-based methods described herein have a variety of applications, including identification of orthologs of one or more query sequences (e.g., gene sequences of interest) in one or more query (or target) genomes, identification of BGCs, identification of resistance genes to secondary metabolites produced by BGCs in a query genome, identification and / or characterization of therapeutic targets, identification and / or characterization of targets of secondary metabolites, and drug discovery.

[0241] In some examples, the present disclosure provides a method and system for determining the boundaries of a gene cluster (e.g., a BGC) by identifying a plurality of genes associated with the gene cluster (e.g., a BGC). In some examples, the method includes: (1) identifying a plurality of query genes that co-localize with anchor genes (e.g., core synthase genes) of the BGC in a query (or reference) genome; (2) for each of the plurality of query genes, performing any one of the computer-implemented methods described in the section "Grid Representation Analysis Method" to determine the likelihood that the query gene is associated with the BGC; and (c) identifying the query genes with a specified high likelihood, which is a likelihood that is greater than a threshold value, as the plurality of genes associated with the BGC. For example, if the query gene or its ortholog has a probability of (1) exceeding any one of 30%, 40%, 50%, 60%, 70%, 80%, 90% or more for the high likelihood category, the query gene is associated with the BGC. In some examples, if the query gene or its ortholog has a combined probability of more than about 50%, 60%, 70%, 80%, 90% or more for the (1) high likelihood and (2) moderately high likelihood categories, the query gene is associated with the BGC. In some examples, if the query gene or its ortholog has a probability of more than about 30%, 40%, 50%, 60%, 70%, 80%, 90% or more for the (4) low likelihood category, the query gene is rejected as not associated with the BGC. The boundaries of the BGC (i.e., upstream and downstream limits) can be determined based on the positions of all genes associated with the BGC determined using this method.

[0242] In some examples, the present disclosure provides a method and system for identifying resistance genes to secondary metabolites produced by a BGC in a query (or reference) genome, comprising: (a) identifying a putative embedded gene that is not involved in the production of secondary metabolites by the BGC and that is co-localized (e.g., within a proximity zone of about 100 kb, 50 kb, 20 kb, or less than a user-specified distance) with an anchor gene (e.g., a core synthase gene) in the BGC in the query genome; (b) performing any one of the computer-implemented methods described in the "Grid Representation Analysis Method" section to determine the likelihood that the putative embedded gene is associated with the BGC; and (c) identifying the putative embedded gene as a resistance gene based at least in part on the likelihood that the putative embedded gene is associated with the BGC. In some examples, the putative embedded gene is identified as a resistance gene if the likelihood that the putative embedded gene is associated with the BGC exceeds a specified threshold. In some examples, the likelihood that the putative embedded gene is associated with the BGC is one of a plurality of factors used to identify the putative embedded gene as a resistance gene associated with the BGC.In some examples, the method further comprises experimentally verifying that the putative embedded gene is a resistance gene to the secondary metabolite produced by the BGC in the query genome.For example, the putative embedded gene can be expressed and contacted with the secondary metabolite to determine whether binding occurs between the product of the putative embedded gene and the secondary metabolite.

[0243] In some examples, the present disclosure provides methods and systems for the identification and / or characterization of mammalian (e.g., human) targets. For example, resistance genes (e.g., fungal resistance genes) with homologs in the human genome identified using the methods described herein provide a link between the resistance gene, the human homolog, and the secondary metabolite produced by the BGC. This link suggests that the human homolog may be a human target of the secondary metabolite, and that the secondary metabolite may interact with and / or regulate the human homolog.

[0244] In some examples, the present disclosure provides a method for identifying and / or characterizing a mammalian (e.g., human) target of a secondary metabolite of a BGC or an analog of a BGC product, comprising: (1) identifying a putative embedded target gene (pETaG) that is co-localized (e.g., within a contiguous zone of about 200 kb, 100 kb, 50 kb, 40 kb, 30 kb, 20 kb or less) with a BGC in a query genome, is homologous to a mammalian (e.g., human) gene, and does not encode an enzyme that produces a secondary metabolite of the BGC; (2) performing any one of the computer-based methods described in the "Grid Representation Analysis Method" section to determine the likelihood that the pETaG is associated with the BGC; and (3) identifying the mammalian (e.g., human) gene as a target of the secondary metabolite of the BGC based at least in part on the likelihood that the pETaG is associated with the BGC. In some examples, the mammalian (e.g., human) gene is identified as a target if the likelihood that it is associated with the BGC exceeds a threshold value. In some examples, the likelihood that pETaG is associated with the BGC is one of a plurality of factors for identifying a mammalian (e.g., human) gene as a target of a secondary metabolite of the BGC. In some examples, the method further comprises identifying a mammalian (e.g., human) homolog of pETaG in a mammalian (e.g., human) genome. In some examples, the method further comprises assaying the effect of a secondary metabolite produced by the BGC or an analog of a BGC product on a mammalian (e.g., human) target.

[0245] In some examples, the disclosure provides methods and systems for drug discovery, for example, identifying small molecule modulators of a mammalian target gene (or a reptilian target gene, an avian target gene, an amphibian target gene, or a target gene from any other organism). In some examples, the disclosure provides a method of identifying a small molecule modulator of a mammalian (e.g., human) target gene or its product, comprising: (a) identifying a homologous gene of the mammalian target gene in a fungal genome (or in an archaeal, bacterial, plant, or other genome that contains a BGC) that is co-localized (e.g., within a proximity zone of about 100 kb, 50 kb, 40 kb, 30 kb, 20 kb or less) with an anchor gene (e.g., a core synthase gene) of a BGC of the fungal genome and is not involved in the production of a secondary metabolite by the BGC; (b) performing any one of the computer-based methods described in the "Grid Representation Analysis Methods" section to determine a likelihood that the homologous gene is associated with the BGC; and (c) identifying a secondary metabolite or an analog thereof as a small molecule modulator of the mammalian target gene or its product based at least in part on the likelihood that the homologous gene is associated with the BGC. In some examples, if the likelihood that the homologous gene is associated with the BGC exceeds a threshold, the secondary metabolite or analog thereof is identified as a small molecule modulator of the mammalian (e.g., human) gene or its product. In some examples, the likelihood that the homologous gene is associated with the BGC is one of a plurality of factors used to identify the secondary metabolite or analog thereof as a small molecule modulator of the mammalian (e.g., human) gene or its product. In some examples, the method further comprises evaluating the interaction of the mammalian target gene product with a compound derived from the secondary metabolite produced by the BGC. In some embodiments, the method comprises contacting the secondary metabolite or analog thereof with a protein encoded by the mammalian target gene and detecting the activity of the protein encoded by the mammalian target gene. In some examples, the activity is binding of the protein encoded by the mammalian target gene with the secondary metabolite or analog thereof.

[0246] In some instances, the secondary metabolite is a product of an enzyme encoded by the BGC or a salt thereof, including a non-naturally occurring salt. In some instances, the secondary metabolite or analog thereof is an analog of the product of the enzyme encoded by the BGC, such as a small molecule compound having the same core structure as the secondary metabolite or a salt thereof.

[0247] In some examples, the disclosure provides a method of regulating a human target, comprising providing a secondary metabolite or an analog thereof produced by an enzyme encoded by a BGC, wherein the human target (or a nucleic acid sequence encoding the human target) is homologous to an ETaG associated with the BGC as determined using any one of the methods described herein.

[0248] In some examples, the disclosure provides a method of treating a condition, disorder, or disease associated with a human target, comprising administering to a subject susceptible to or suffering from a secondary metabolite produced by an enzyme encoded by a BGC, or an analog thereof, wherein the human target (or a nucleic acid sequence encoding the human target) is homologous to an ETaG associated with the BGC as determined using any one of the methods described herein.

[0249] In some instances, the secondary metabolite is produced by a fungus. In some instances, the secondary metabolite is acyclic. In some instances, the secondary metabolite is a polyketide. In some instances, the secondary metabolite is a terpene compound. In some instances, the secondary metabolite is a non-ribosomally synthesized peptide.

[0250] In some examples, an analog of a substance (e.g., a secondary metabolite) that shares one or more specific structural features, elements, components or moieties with a reference substance. Typically, an analog shows significant structural similarity with a reference substance, e.g., shares a core or consensus structure, but differs in a specific individual manner. In some examples, an analog is a substance that can be generated from a reference substance, e.g., by chemical manipulation of the reference substance. In some examples, an analog is a substance that can be generated by carrying out a synthetic process that is substantially similar (e.g., shares multiple steps) to that which produces the reference substance. In some examples, an analog is generated or can be generated by carrying out a synthetic process that is different from that used to generate the reference substance. In some examples, an analog of a substance is a substance that is substituted at one or more of its substitutable positions.

[0251] In some examples, the analog of the product includes the structural core of the product. In some examples, the biosynthetic product is cyclic, e.g., monocyclic, bicyclic, or polycyclic, and the structural core of the product is or includes a monocyclic, bicyclic, or polycyclic ring system. In some examples, the structural core of the product includes one ring of the bicyclic or polycyclic ring system of the product. In some examples, the product is or includes a polypeptide, and the structural core is the backbone of the polypeptide. In some examples, the product is or includes a polyketide, and the structural core is the backbone of the polyketide. In some examples, the analog is a substituted biosynthetic product that includes one or more suitable substitutents.

[0252] Identification of ETaG As noted above, in some examples, the disclosure provides methods for identifying embedded target genes ("ETaGs") or mammalian (e.g., human) target genes corresponding to ETaGs. The disclosure provides methods for identifying and / or characterizing ETaGs, databases including biosynthetic gene clusters and / or ETaG gene sequences (and optionally associated annotations), systems for identifying and / or characterizing human target genes corresponding to ETaGs, and methods for making and / or using such human target genes and / or systems containing and / or expressing them, and the like. ETaGs are described, for example, in WO201955816, the contents of which are incorporated herein by reference. The methods described herein provide improved methods for identifying ETaGs that are truly associated with BGCs and reducing the calling of false positive ETaGs identified based on co-localization and / or co-regulation with one or more biosynthetic genes in BGCs in a particular genome.

[0253] In some examples, the methods described herein are applied to identify ETaGs from fungal genomes. In some examples, ETaGs from eukaryotic fungi can have more similarity to mammalian genes than their counterparts (if any) in prokaryotes, such as certain bacteria. In some examples, fungi contain and / or contain more therapeutically relevant ETaGs than organisms that are evolutionarily more distant from humans.

[0254] In some examples, the method includes (a) identifying a putative ETaG (pETaG) sequence in a fungal genome, where (1) pETaG is co-localized (i.e., in a relative proximity zone) with an anchor gene (e.g., a core synthase gene) in a BGC, (2) pETaG is not involved in the production of secondary metabolites by the BGC, and (3) pETaG is homologous to an expressed mammalian nucleic acid sequence; (b) determining the likelihood that pETaG is associated with the BGC using any one of the computer-based methods described in the section "Grid Representation Analysis Methods"; and (c) identifying pETaG as an ETaG based on the likelihood that pETaG is associated with the BGC. For example, pETaG can be identified as an ETaG if the likelihood is above a threshold value or if the likelihood is one of multiple factors used to identify pETaG as an ETaG. In some examples, pETaG is co-regulated with at least one biosynthetic gene in the BGC. In some examples, pETaG is not co-regulated with at least one biosynthetic gene in the BGC. In some examples, the method is repeated for multiple pETaGs and used to prioritize pETaGs for experimental validation based on the likelihood that the pETaG is associated with the BGC. In some examples, ETaG is compared to mammalian, e.g., human, nucleic acid sequences to identify homologous mammalian nucleic acid sequences. In some examples, such methods can be used to identify ETaGs on a genome-wide scale, e.g., from sequences of many (e.g., hundreds, thousands, or more) genomes. The identified ETaGs can be prioritized based on the therapeutic importance of their mammalian homologs, particularly human homologs. In some examples, the biosynthetic products (i.e., secondary metabolites) or analogs thereof produced by enzymes encoded by the relevant biosynthetic gene clusters are modulators (e.g., activators, inhibitors, etc.) of human targets. In some embodiments, the biosynthetic products (i.e., secondary metabolites) or analogs thereof produced by the enzymes encoded by the relevant biosynthetic gene clusters are modulators (e.g., activators, inhibitors, etc.) of animal, bacterial, fungal or plant targets.

[0255] As will be readily understood by those skilled in the art, once the association between the biosynthetic product from the biosynthetic gene cluster, ETaG, and the human target has been established, it can be exploited in a variety of ways. For example, starting from the biosynthetic product produced by the enzymes encoded by the biosynthetic gene cluster, one can identify the ETaG located within the designated proximity zone of the biosynthetic genes of the biosynthetic gene cluster, and then identify the human target that is homologous to the ETaG. Once the human target is identified, it can be prioritized (even if previously thought to be undruggable) and the biosynthetic product can be used to develop a modulator of the human target, including any further optimization of the biosynthetic product, for medical use, for example, by preparing and assaying analogs of the product using any of a variety of methods known to those skilled in the art. Starting from the human target of therapeutic interest, one can also identify the ETaG that is homologous to the human target, and then identify the biosynthetic gene cluster that contains the biosynthetic gene within the designated proximity zone of the ETaG. Once the biosynthetic gene cluster is identified, one can characterize the biosynthetic product produced by the enzymes encoded by the biosynthetic gene cluster, and assay for the regulation of the human target or its product. The biosynthetic products can be used as lead compounds for the optimization of drug candidates using any of a variety of methods known to those of skill in the art to provide agents useful for many medical purposes, e.g., therapeutic purposes. In some embodiments, targets may be derived from other kingdoms of life, such as animals, plants, fungi, bacteria, archaea, etc.

[0256] In some instances, the present disclosure provides particular insight into targets that were previously thought to be undruggable, by providing methods for identifying their homologous ETaGs in fungi and elucidating the associated biosynthetic gene clusters. In some instances, the present disclosure significantly improves the druggability of targets that were previously thought to be undruggable, by, for example, identifying their homologous ETaGs in fungi, elucidating the associated biosynthetic gene clusters, and testing the biosynthetic products of the relevant biosynthetic gene clusters, in some cases essentially converting them into druggable targets (where they can be used directly as modulators of the human target and / or where their analogs can be used as modulators).

[0257] The ETaG is within a proximal zone relative to an anchor gene (e.g., a core synthase gene) in the BGC, is homologous to an expressed mammalian nucleic acid sequence, and is optionally co-regulated with at least one biosynthetic gene in the BGC. In some examples, the ETaG is located no more than about 100 kb, 50 kb, 40 kb, 30 kb, 20 kb, 10 kb, or less from an anchor gene (e.g., a core synthase gene) in the BGC.

[0258] In some examples, ETaG is a product of, or is homologous to a human nucleic acid sequence encoding, an existing target of therapeutic interest. In some examples, ETaG is a product of, or is homologous to a human nucleic acid sequence encoding, a novel target of therapeutic interest. In some examples, ETaG is a product of, or is homologous to a human nucleic acid sequence encoding, a target that was considered undruggable prior to the present disclosure. In some examples, ETaG is a product of, or is homologous to a human nucleic acid sequence encoding, a target that was considered undruggable by small molecules prior to the present disclosure.

[0259] In some examples, the ETaG sequence is homologous to the expressed mammalian nucleic acid sequence in that the sequence or a portion thereof is at least 20%, 30%, 40%, 50%, 60%, 70%, 80%, or 90% identical to that of the expressed mammalian nucleic acid sequence. In some examples, the ETaG sequence is homologous to the mammalian nucleic acid sequence in that the mRNA produced from ETaG or a portion thereof is homologous to that of the mammalian nucleic acid sequence. In some examples, the homologous portion is at least 50, 100, 150, 200, 500, 1000, 2000, 3000, or 5000 base pairs in length. In some examples, the homologous portion encodes a conserved protein or a portion of a conserved protein, such as a protein domain, a set of residues involved in a function (e.g., interaction with another molecule (e.g., protein, small molecule, etc.), enzymatic activity, etc.), from fungi to mammals. In some examples, the mammalian nucleic acid, such as a human nucleic acid sequence, is associated with a human disease, disorder, or condition. In some instances, such human nucleic acid sequences are existing targets for therapeutic purposes. In some instances, such human nucleic acid sequences are novel targets for therapeutic purposes. In some instances, such human nucleic acid sequences are targets previously thought to be insensitive to targeting, e.g., by small molecules.

[0260] In some examples, the ETaG sequence is homologous to a mammalian nucleic acid sequence in that a product encoded by ETaG, or a portion thereof, is homologous to a product encoded by a mammalian nucleic acid sequence. In some examples, the ETaG sequence is homologous to a mammalian nucleic acid sequence in that a protein encoded by ETaG, or a portion thereof, is homologous to a protein encoded by a mammalian nucleic acid sequence. In some examples, the ETaG sequence is homologous to a mammalian nucleic acid sequence in that a portion of a protein encoded by ETaG is homologous to a protein encoded by a mammalian nucleic acid sequence.

[0261] In some instances, the portion of the protein is a protein domain. In some instances, the protein domain is an enzyme domain. In some instances, the protein domain interacts with one or more factors, such as small molecules, lipids, carbohydrates, nucleic acids, proteins, etc.

[0262] In some instances, a portion of a protein is a functional and / or structural domain that defines the protein family that the protein belongs to. Amino acids within a particular catalytic or structural domain that defines a patent family can be selected based on predicted subfamily domain architecture and, optionally, verified by various assays, for use in homology alignment analysis.

[0263] In some instances, the portion of the protein is a set of critical residues, contiguous or non-contiguous, that are important for the function of the protein. In some instances, the function is an enzymatic activity and the portion of the protein is a set of residues necessary for the activity. In some instances, the function is an enzymatic activity and the portion of the protein is a set of residues that interact with a substrate, intermediate, or product. In some instances, the set of residues interacts with a substrate. In some instances, the set of residues interacts with an intermediate. In some instances, the set of residues interacts with a product.

[0264] In some examples, the function of the protein is the interaction with one or more factors, such as small molecules, lipids, carbohydrates, nucleic acids, proteins, etc., and a part of the protein is a set of residues required for the interaction. In some examples, the set of residues each independently contacts an interacting agent. For example, in some examples, each of the residues of the set independently contacts an interacting small molecule. In some examples, the protein is a kinase, the interacting small molecule is or includes a nucleic acid base, and each of the set of residues independently contacts the nucleic acid base, for example, via hydrogen bonds, electrostatic forces, van der Waals forces, aromatic stacking, etc. In some examples, the interacting agent is another macromolecule. In some examples, the interacting agent is a nucleic acid. In some examples, the set of residues are residues that contact an interacting nucleic acid, such as residues in a transcription factor. In some examples, the set of residues are residues that contact an interacting protein.

[0265] In some instances, the portion of the protein is or includes an essential structural element for protein effector recruitment and / or binding, for example, based on the tertiary protein structure of the human target.

[0266] Portions of proteins, such as protein domains, sets of residues that are responsible for biological functions, can be conserved from species to species, for example, from fungi to humans, in some instances, as shown in this disclosure.

[0267] In some examples, protein homology is measured based on exact identity, e.g., the same amino acid residue at a given position. In some examples, homology is measured based on amino acid residues having one or more properties, e.g., one or more identical or similar properties (e.g., polar, non-polar, hydrophobic, hydrophilic, size, acidic, basic, aromatic, etc.). Exemplary methods for assessing homology are widely known in the art and can be utilized in accordance with the present disclosure, e.g., MAFFT, MUSCLE, TCoffee, ClustalW, etc.

[0268] In some examples, the protein or portion thereof encoded by ETaG (e.g., those described in this disclosure) is at least 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, or 99%, or 100% (in the case of 100%, it is identical) homologous to that encoded by the mammalian nucleic acid sequence. In some examples, the protein or portion thereof encoded by ETaG is at least 50%, 60%, 70%, 80%, 85%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, or 99%, or 100% homologous to that encoded by the expressed mammalian nucleic acid sequence.

[0269] In some examples, ETaG is co-regulated with at least one biosynthetic gene in the biosynthetic gene cluster. In some examples, ETaG is co-regulated with two or more genes in the biosynthetic gene cluster. In some examples, ETaG is co-regulated with the biosynthetic gene cluster in that expression of ETaG increases or is turned on when a biosynthetic product produced by an enzyme encoded by the biosynthetic gene cluster (a biosynthetic product of the biosynthetic gene cluster) is produced. In some examples, ETaG is co-regulated with the biosynthetic gene cluster in that expression of ETaG increases or is turned on when the level of a biosynthetic product of the biosynthetic gene cluster increases.

[0270] In some examples, an organism comprising ETaG comprises one or more homologous genes of ETaG. In some examples, the ETaG gene sequence may optionally be more than about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, or 99% homologous to one or more gene sequences in the same genome. In some examples, the ETaG gene sequence is optionally more than about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, or 99% homologous to 2, 3, 4, 5, 6, 7, 8, 9 or more gene sequences in the same genome. In some examples, the homology is more than 10%. In some examples, the homology is more than 20%. In some examples, the homology is more than 30%. In some examples, the homology is more than 40%. In some examples, the homology is more than 50%. In some instances, the homology is greater than 60%. In some instances, the homology is greater than 70%. In some instances, the homology is greater than 80%. In some instances, the homology is greater than 90%.

[0271] In some examples, the ETaG gene sequence is about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, or 99% identical or less to any expressed gene sequence in at least 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99%, 99.1%, 99.2%, 99.3%, 99.4%, 99.5%, 99.6%, 99.7%, 99.8%, or 99.9% of the fungal nucleic acid sequences in a set that optionally are derived from different fungal strains and that include homologous biosynthetic gene clusters. In some examples, the ETaG gene sequence is about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, or 99% identical to at least 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99%, 99.1%, 99.2%, 99.3%, 99.4%, 99.5%, 99.6%, 99.7%, 99.8%, or 99.9% of the fungal gene sequences that are optionally within a proximity zone to biosynthetic genes of a homologous biosynthetic gene cluster from a different fungal strain. In some examples, the ETaG gene sequence is about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, or 99% identical to at least 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99%, 99.1%, 99.2%, 99.3%, 99.4%, 99.5%, 99.6%, 99.7%, 99.8%, or 99.9% of the fungal gene sequences that are optionally within a proximity zone to biosynthetic genes of a homologous biosynthetic gene cluster from a different fungal strain. In some examples, the ETaG gene sequence is about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%), 95%), or 99% or less identical to any expressed gene sequence in any fungal nucleic acid sequence in the set that optionally is derived from a different fungal strain and comprises a homologous biosynthetic gene cluster.In some examples, the ETaG gene sequence is optionally less than about 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 85%, 90%, 95%, or 99% identical to any expressed gene sequence that is within a proximity zone to a biosynthetic gene of a homologous biosynthetic gene cluster from a different fungal strain. In some examples, it is less than about 10% identical. In some examples, it is less than about 20% identical. In some examples, it is less than about 30% identical. In some examples, it is less than about 40%) identical. In some examples, it is less than about 50% identical. In some examples, it is less than about 60% identical. In some examples, it is less than about 70%) identical. In some examples, it is less than about 80% identical. In some examples, it is less than about 90% identical.

[0272] In some examples, the human target gene and / or its product is sensitive to regulation by the biosynthetic product of the biosynthetic gene cluster or its analog, and the human target gene has its homologous ETaG embedded in the biosynthetic gene cluster or located in a proximal zone designated for the biosynthetic gene of the cluster. In some examples, the protein encoded by the human target gene is sensitive to regulation by the biosynthetic product of the biosynthetic gene cluster or its analog, and the human target gene has its homologous ETaG embedded in the biosynthetic gene cluster or located in a proximal zone designated for the biosynthetic gene of the cluster. Thus, in some examples, the present disclosure not only provides novel human targets, but also methods and agents for regulating such human targets. In some examples, the compounds produced by the enzymes of the biosynthetic gene cluster interact with and / or regulate targets encoded by mammalian, e.g., human, nucleic acid sequences that are homologous to ETaGs associated with the biosynthetic gene cluster.

[0273] In some examples, the present disclosure provides a method for evaluating compounds using the identified ETaG and the product encoded thereby.In some examples, the present disclosure provides a method, comprising: contacting at least one test compound with the gene product encoded by the embedded target gene of fungal nucleic acid sequence; and determining that the level or activity of the gene product is changed in the presence of the test compound compared to the absence of the test compound, or determining that the level or activity of the gene product is equivalent to that observed in the presence of a reference agent that has a known effect on the level or activity.

[0274] In some examples, the disclosure provides a method for identifying and / or characterizing a mammalian, e.g., human target of a product or product analog produced by an enzyme encoded by a biosynthetic gene cluster, comprising identifying a human homolog of ETaG determined to be associated with a BGC using any one of the methods described herein, and optionally assaying the effect of the product or product analog produced by an enzyme encoded by a biosynthetic gene cluster on the target.

[0275] Further analysis may include assessing the conservation / similarity of essential structural elements for protein effector recruitment / binding, for example based on examination of the tertiary protein structure of the human target. For example, in some instances, aligned sequences were compared to PDB crystal structures. In some instances, only amino acids within specific catalytic or structural domains that define the PFAM boundaries of the ETaG / target (e.g., based on predicted subfamily domain architectures) were used for alignment analysis. ETaG sequences were directly compared to their human counterparts by aligning all ETaG and human target proteins and their phylogenetic relationships to obtain quantitative correlation data (e.g., peptide sequence similarity and / or evolutionary tree visualization) corresponding to target protein residues within 4 angstroms of the corresponding engaged protein.

[0276] Without intending to be limited by any theory, if these structural motifs are conserved in fungal ETaGs, it may indicate an increased probability that metabolites produced by ETaG-related biosynthetic gene clusters are effectors of both fungal and human target proteins, and that the produced metabolites may be drug candidates or leads for drug development against the human target. In some examples, the above analysis is used to prioritize ETaGs and their related biosynthetic gene clusters, as well as metabolites produced from the biosynthetic gene clusters, for targeting of human targets.

[0277] Computer Systems In some examples, the computer-based method, sequence, genome and / or database provided is embodied in a computer-readable medium. In some examples, the present disclosure provides a system including one or more non-transient machine-readable storage media that stores data representing the computer-based method, sequence, genome and / or database provided. Non-transient machine-readable storage media suitable for embodying the data provided include all forms of non-volatile storage, including, by way of example, semiconductor storage devices, such as EPROM, EEPROM and flash storage devices, magnetic disks, such as internal hard disks or removable disks, magneto-optical disks, and CD-ROM and DVD-ROM disks. In particular, the system provided may be particularly efficient for the provided sets and databases having the specific structures described herein.

[0278] In some examples, the present disclosure provides a computer system that can perform the methods described herein. In some examples, the present disclosure provides a computer system adapted to perform the methods provided. In some examples, the present disclosure provides a computer system adapted to query genomes and / or genome databases, for example, to identify homologs of one or more query sequences. In some examples, the present disclosure provides a computer system adapted to access one or more genome databases.

[0279] A computer system that can be used to implement all or part of the provided method may include various forms of digital computers. Examples of digital computers include, but are not limited to, laptops, desktops, workstations, personal digital assistants, servers, blade servers, mainframes, smart televisions, and other suitable computers. Mobile devices can be used to implement all or part of the provided technology. Mobile devices include, but are not limited to, tablet computing devices, personal digital assistants, mobile phones, smartphones, digital cameras, digital glasses, and other portable computing devices. The computing devices, their connections and relationships, and their functions described herein are meant to be merely examples and are not meant to be limitations on the implementation of the technology.

[0280] All or a portion of the techniques described herein and various modifications thereof may be implemented, at least in part, via a computer program product, e.g., a computer program tangibly embodied in one or more information carriers, e.g., one or more tangible machine-readable storage media, for execution by or to control the operation of a data processing apparatus, e.g., a programmable processor, computer, or multiple computers.

[0281] The computer programs for the provided techniques can be written in any type of programming language, including compiled or interpreted languages, and can be deployed in any type, including a stand-alone program, or a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program can be deployed to be executed on one computer, or on multiple computers at one site, or distributed across multiple sites and interconnected by a network.

[0282] For example, operations associated with implementing the programs and techniques may be performed by one or more programmable processors executing one or more computer programs to perform the techniques provided. All or a portion of the processes may be implemented as special purpose logic circuitry, such as an FPGA (field programmable gate array) and / or an ASIC (application specific integrated circuit).

[0283] Processors suitable for executing computer programs include, by way of example, both general-purpose and special-purpose microprocessors, and any one or more processors of any kind of digital computer. Generally, a processor receives instructions and data from a read-only or random-access memory area, or both. Elements of a computer (including a server) include one or more processors for executing instructions and one or more storage devices for storing instructions and data. Generally, a computer also includes, or is operatively coupled to receive data from, or transfer data to, one or more machine-readable storage media, such as mass storage devices for storing data, e.g., magnetic, magneto-optical, or optical disks. Non-transitory machine-readable storage media suitable for embodying computer program instructions and data include, by way of example, all forms of non-volatile storage, including semiconductor storage devices, e.g., EPROM, EEPROM, and flash storage devices, magnetic disks, e.g., internal hard disks or removable disks, magneto-optical disks, and CD-ROM and DVD-ROM disks.

[0284] Each computing device, such as a tablet computer, may include a hard drive for storing data and computer programs, and a processing device (e.g., a microprocessor) and memory (e.g., RAM) for executing the computer programs. Each computing device may include an image capture device, such as a still camera or video camera. The image capture device may be built-in or simply accessible to the computing device.

[0285] Each computing device may include a graphics system including a display screen. The display screen, such as an LCD or CRT (cathode ray tube), displays to the user images generated by the graphics system of the computing device. As is well known, a display on a computer display (e.g., a monitor) physically transforms the computer display. For example, if the computer display is LCD-based, the orientation of the liquid crystals may be changed by application of a bias voltage in a physical transformation that is visually apparent to the user. As another example, if the computer display is a CRT, the state of the phosphor screen may be changed by the influence of electrons in a physical transformation that is also visually apparent. Each display screen may be touch-sensitive, which allows a user to input information into the display screen via a virtual keyboard. In some computing devices, such as desktops or smartphones, a physical QWERTY keyboard and scroll wheel may be provided to input information into the display screen. Each computing device, and the computer programs running thereon, may also be configured to receive voice commands and perform functions in response to such commands.

[0286] example Example 1 - An exemplary workflow for the calculation of phylogenetic comparison metrics based on a custom algorithm FIG. 13 provides a schematic of an exemplary workflow for the calculation of "phylogenetic signatures" using a custom algorithm. 1) Use pETaG to search for homologs in a set of positive and negative genomes. 2) Identify the last common ancestor (LCA) and designate clades as pETaG or house copy clades (e.g., housekeeping versions of pETaG). 3) Calculate copy number differences between genes in the positive vs. negative genomes from the LCA. 4) Calculate the distance of each gene to the LCA, average the distances per clade, and calculate the ratio of pETaG / house copy clade. 5) Calculate and average the distances between pairwise combinations of genes within each pETaG or house copy clade, and calculate the ratio of pETaG / house copy clade. 6) Sum the branch lengths of pETaG and house copy clades, respectively, and calculate the branch length ratio of pETaG / house copy clade.

[0287] Example 2 - Phylogenetic signatures calculated for lovastatin ETaG Figure 14 provides a non-limiting example of the use of the above custom phylogenetic algorithm for lovastatin ETaG. 1) A phylogenetic tree is created from the set of homologous ETaG genes from the set of positive and negative genomes. 2) LCA clades are identified by traversing the tree clades by working backwards from ETaG until there is a clade that contains ETaG and at least one gene from all negative genomes. 3) The clades are separated into ETaG clades and house copy clades, and the copy number difference CND between the genes in the positive genome and the genes in the negative genome is calculated. The results show that the positive lovastatin genomes have an average of 0.86 increased copy number of HMG-CoA reductase homologs. 4) Additional phylogenetic features as above are calculated. A resulting ratio of >1 indicates that the ETaG clade is evolving faster than the house copy clade.

[0288] Example 3 - An example workflow for performing co-evolutionary evaluation Figure 15 provides a non-limiting example of a workflow for assessing coevolution when comparing coevolution between pairs of COGs using percent sequence identity. 1) COGs were identified as described above. 2) After sequence alignment and trimming, the pairwise percent sequence identity was calculated for each pair of COGs. 3) The pairwise percent sequence identity was summarized in a table. Pearson R and orthogonal regression analysis were used to explore pairwise relationships between all pairwise combinations of COGs. 4) The plot in panel 4 provides a simulated example of the coevolution results of three COGs-COG 1, COG 2, and COG 3. The results show that COG 1 has coevolved with COG 2 but not with COG 3. 5) The plot provided a true example of the coevolution results of lovastatin ETaG (panel 5, left) and house copy (panel 5, right) for the core synthase of lovastatin BGC. The results indicate that lovastatin ETaG has coevolved with the core synthase, but not with the house copy.

[0289] Example 4 - Performance data for a deep learning model trained to evaluate pETaGs for the likelihood that they are true ETaGs FIG. 16 provides non-limiting examples of the number of units per hidden layer in five different versions (i.e., versions 1, 2, 3, 4, and 5) of a deep learning model including two hidden layers trained as described above to assess the likelihood that a putative ETaG is an actual ETaG.

[0290] FIG. 17 provides a non-limiting example of performance data (test loss) for five versions (i.e., versions 1, 2, 3, 4, and 5) of the deep learning model described in FIG. 16. Model performance shown as the sum of errors after each iteration of optimization. The model that shows a decrease in test loss after each iteration is determined to be the better model.

[0291] FIG. 18 provides non-limiting examples of performance data (test specificity), i.e., values ​​of true predicted negatives / actual negatives sum for different training epochs, for five different versions (i.e., versions 1, 2, 3, 4, and 5) of the deep learning model described in FIG.

[0292] FIG. 19 provides non-limiting examples of performance data (test sensitivity), i.e., values ​​of sum of true predicted positives / actual positives for different training epochs, for five different versions (i.e., versions 1, 2, 3, 4 and 5) of the deep learning model described in FIG.

[0293] FIG. 20 provides non-limiting examples of performance data (test accuracy), i.e., true predicted positives / total predicted positives values ​​for different training epochs, for five different versions (i.e., versions 1, 2, 3, 4, and 5) of the deep learning model described in FIG.

[0294] Example 5 - Summary of exemplary target evaluation values ​​of known ETaGs including BGC (lovastatin, Monascus ruber) and negative ETaG (histone H3.2, Dendryphion sp.). Table 2 shows examples of target evaluation values ​​for a known ETaG containing a BGC (lovastatin, Monascus ruber) and a negative ETaG (histone H3.2, Dendryphion sp.). Two known lovastatin candidates are shown, which had various results in terms of the comparison heatmap. In this scenario, the empirical score and deep learning probability, along with the target evaluation metric, provide confidence that the ETaG is a true ETaG. [Table 2]

[0295] Exemplary embodiments Among the embodiments provided are the following: 1. A computer-implemented method for identifying embedded target genes (ETaGs), comprising: specifying one or more query sequences or proxies thereof; selecting one or more target genomes; performing a search of one or more target genomes using the one or more query sequences or proxies thereof to identify putative embedded target gene (pETaG) sequences that are homologs of the one or more query sequences based on a comparison of one or more homology-based metrics for the candidate pETaGs to one or more predetermined homology-based metric thresholds; Determining whether a given pETaG is an ETaG based on comparative genomics analysis of multiple genomes; 4. A computer-implemented method comprising: 2. The computer-implemented method of embodiment 1, wherein the comparative genomics analysis comprises generating a comparative genomics heat map based on the plurality of genomes. 3. The computer-implemented method of embodiment 1 or embodiment 2, wherein the plurality of genomes comprises a plurality of positive genomes and a plurality of negative genomes. 4. The computer-implemented method of any one of embodiments 1-3, wherein the comparative genomics analysis comprises determining phylogenetic characteristics, co-occurrence characteristics, co-evolution characteristics, or any combination thereof, for a given pETaG based on multiple genomes. 5. The computer-implemented method of any one of embodiments 1-4, wherein the comparative genomics analysis comprises analysis of an input dataset comprising phylogenetic features, co-occurrence features, co-evolution features, comparative genomics heat maps, data derived from comparative genomics heat maps, or any combination thereof, for pETaG using a machine learning model or an empirical algorithm to predict the probability that pETaG is an ETaG. 6. The computer-implemented method of any one of embodiments 1-5, further comprising determining that the identified pETaG is associated with a resistance mechanism based on determining the copy number of the identified pETaG. 7. The computer-implemented method of any one of embodiments 1-6, further comprising determining that pETaG is involved in the resistance mechanism based on determining a copy number difference between a positive genome containing pETaG and a negative genome not containing pETaG. 8. The computer-implemented method of any one of embodiments 1-7, wherein the one or more query sequences, or proxies thereof, comprise one or more protein sequences, one or more nucleic acid sequences, one or more Universal Protein Resource (Uniprot) identification numbers, one or more profile hidden Markov models (pHMMs), a specified set of protein sequence domains, or any combination thereof. 9. The computer-implemented method of any one of embodiments 1-8, wherein the one or more query sequences or proxies thereof are selected from a bacterial genome, an archaeal genome, a fungal genome, a plant genome, an animal genome, a human genome, or any combination thereof. 10. The computer-implemented method of any one of embodiments 1-8, wherein the one or more target genomes are selected from a bacterial genome, an archaeal genome, a fungal genome, a plant genome, an animal genome, a human genome, or any combination thereof. 11. The computer-implemented method of any one of embodiments 1-10, wherein the two or more target genomes are selected based on pairwise similarity scores, pairwise phylogenetic distances, or any combination thereof. 12. The computer-implemented method of embodiment 11, further comprising filtering the two or more selected target genomes to retain only those target genomes whose (i) pairwise similarity scores are greater than a specified pairwise similarity threshold, or (ii) whose pairwise phylogenetic distances are less than a specified phylogenetic distance threshold. 13. The computer-implemented method of embodiment 12, further comprising clustering the retained target genomes into sets using a clustering algorithm, and performing a search using one or more of the sets of clustered target genomes. 14. The computer-implemented method of embodiment 13, wherein the clustering algorithm comprises a Markov cluster algorithm. 15. The computer-implemented method of any one of embodiments 1-14, wherein the search is performed using BLAST, DIAMOND, HMMER, Exonerate, or ggsearch. 16. The computer-implemented method of any one of embodiments 1 to 15, wherein the search is limited to one or more specific regions of the one or more target genomes. 17. The computer-implemented method of embodiment 16, wherein the one or more specific regions comprise one or more biosynthetic gene clusters (BGCs). 18. The computer-implemented method of embodiment 17, wherein one or more BGCs in one or more target genomes are predicted using a BGC search algorithm. 19. The computer-implemented method of embodiment 18, wherein the BGC search algorithm includes antiSMASH, SMURF, TOUCAN, or deepBGC. 20. The computer-implemented method of embodiment 17, wherein one or more BGCs are predicted for one or more target genomes by extracting sequence regions of specified length proximal to gene sequences that match known biosynthetic core synthases determined using a sequence search tool. 21. The computer-implemented method of embodiment 20, wherein the sequence search tool comprises BLAST, DIAMOND, HMMER, Exonerate or ggsearch. 22. The computer-implemented method of embodiment 17, wherein one or more BGCs are predicted for one or more query genomes using a hidden Markov model (HMM) of known core synthases. 23. The computer-implemented method of embodiment 17, wherein one or more BGCs are predicted for one or more target genomes based on co-localization of protein sequence domains related to known core synthases. 24. The computer-implemented method of any one of embodiments 1-23, wherein the one or more homology sequence-based metrics include percent sequence identity, percent sequence coverage, E-value, bit score, HMM score, or any combination thereof. 25. The computer-implemented method of any one of embodiments 1-24, wherein the one or more predetermined homologous sequence-based metric thresholds include a sequence identity percentage threshold, a sequence coverage percentage threshold, an E-value threshold, a bit score threshold, an HMM score threshold, or any combination thereof. 26. The computer-implemented method of embodiment 25, wherein the one or more predetermined homologous sequence-based metric thresholds comprise a sequence identity percentage threshold having a value of at least 20%, at least 30%, at least 40%, at least 50%, at least 60%, at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least 98%. 27. The computer-implemented method of embodiment 25, wherein the one or more predetermined homologous sequence-based metric thresholds comprise a sequence coverage percentage threshold having a value of at least 20%, at least 30%, at least 40%, at least 50%, at least 60%, at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least 98%. 28. One or more predetermined homology-based metric thresholds are less than 10, less than 9, less than 8, less than 7, less than 6, less than 5, less than 4, less than 3, less than 2, less than 1, less than 0.01, less than 0.001, less than 1e -10 Less than 1e -20 Less than 1e -30 Less than 1e -40 Less than 1e -50 Less than 1e -60 Less than 1e -70 Less than 1e -80 Less than 1e -90 Less than or equal to 1e -100 26. The computer-implemented method of embodiment 25, comprising an E value threshold having a value less than 29. The computer-implemented method of embodiment 25, wherein the one or more predetermined homologous sequence-based metric thresholds comprise a bit score threshold having a value of at least 40, at least 50, at least 60, at least 70, at least 80, at least 90, at least 100, at least 250, at least 500, or at least 1000, or at least 5000. 30. The computer-implemented method of embodiment 25, wherein the one or more predetermined homologous sequence-based metric thresholds comprise an HMM score threshold having a value of at least 10, at least 25, at least 50, at least 100, at least 250, at least 500, at least 1000, or at least 5000. 31. Conducting a search converting one or more query sequences, including protein sequences, into nucleic acid sequences; performing a search of one or more target genomes using one or more query sequences converted to nucleic acid sequences to identify homologous nucleic acid sequences based on a comparison of one or more homology-based metrics of the candidate pETaGs to one or more predetermined homology-based metric thresholds; comparing the genomic coordinates of the homologous nucleic acid sequences with genomic coordinates corresponding to predicted protein sequences in one or more target genomes; 31. The computer-implemented method of any one of embodiments 1 to 30, comprising: 32. The computer-implemented method of embodiment 31, wherein if the homologous nucleic acid sequences overlap with nucleic acid sequences corresponding to a single predicted protein sequence and the overlap is greater than a specified nucleic acid sequence overlap threshold, the predicted protein sequence is reported as pETaG. 33. The computer-implemented method of embodiment 31, wherein if a homologous nucleic acid sequence overlaps with nucleic acid sequences corresponding to multiple predicted protein sequences and each overlap is greater than a specified nucleic acid sequence overlap threshold, only one of the predicted protein sequences is reported as pETaG. 34. The computer-implemented method of embodiment 33, wherein the predicted protein sequence reported as pETaG is the predicted protein sequence for which the homologous nucleic acid sequences and the nucleic acid sequences corresponding to the predicted protein sequence exhibit the greatest percent sequence identity, percent sequence coverage, E value or bit score value. 35. The computer-implemented method of embodiment 33, wherein the predicted protein sequence reported as pETaG is the predicted protein sequence in which the homologous nucleic acid sequence and the nucleic acid sequence corresponding to the predicted protein sequence show the longest overlapping sequence. 36. The computer-implemented method of embodiment 31, wherein if the homologous nucleic acid sequence overlaps with one or more nucleic acid sequences corresponding to one or more predicted protein sequences, but the respective overlaps are less than a specified nucleic acid sequence overlap threshold, the longest predicted protein sequence is reported as pETaG. 37. The computer-implemented method of embodiment 31, wherein if the homologous nucleic acid sequence does not overlap with the nucleic acid sequence corresponding to the predicted protein sequence, the genomic coordinates of the homologous nucleic acid sequence are reported as pETaG. 38. The computer-implemented method of any one of embodiments 32-37, wherein the specified nucleic acid sequence overlap threshold has a value of at least 20%, 30%, 40%, 50%, 60%, 70%, 75%, 80%, 85%, 90%, 95%, or 98%. 39. A comparative genomics heat map generated for a given pETaG includes a plurality of cells arranged in a grid according to a first axis and a second axis, the first axis corresponds to a plurality of different target genomes, the plurality of different target genomes including a plurality of positive genomes each having an ortholog of an anchor gene sequence of one of the known BGCs of the target genome and a plurality of negative genomes having no ortholog of the anchor gene sequence, the second axis corresponds to a plurality of target gene sequences co-localized with the anchor gene sequence of the known BGC, the putative embedded target gene (pETaG) is one of the plurality of co-localized target gene sequences, and the numerical value of each cell is (i) the presence or absence of orthologs of each co-localized query gene sequence in each target genome; (ii) sequence similarity of the orthologs to each colocalized query gene sequence; (iii) whether orthologs of each query gene sequence colocalize with orthologs of the anchor gene sequence in each genome; The computer-implemented method of any one of embodiments 2 to 38, based on 40. The computer-implemented method of embodiment 39, further comprising analyzing the comparative genomics heat map or its underlying data using a trained machine learning model, wherein the machine learning model is trained to determine the likelihood that a putative embedded gene is embedded in a gene cluster based on the numerical values ​​in a plurality of cells in the grid representation. 41. The computer-implemented method of embodiment 40, wherein the trained machine learning model comprises a long short-term memory (LSTM) model or a convolutional neural network (CNN). 42. The computer-implemented method of any one of embodiments 5-41, wherein the machine learning model used to predict the probability that pETaG is an ETaG comprises a supervised learning model. 43. The computer-implemented method of embodiment 42, wherein the supervised learning model comprises a deep learning model. 44. The computer-implemented method of embodiment 42, wherein the supervised learning model includes a decision tree model. 45. A computer-implemented method for determining the likelihood that a putative embedded target gene (pETaG) is a resistance gene for a secondary metabolite produced by a biosynthetic gene cluster (BGC) in a query genome, comprising: a) Below: i) the likelihood that pETaG is associated with a BGC based on the presence or absence of each of a plurality of orthologs of a plurality of query genes that are co-localized with the BGC in a plurality of different genomes, the plurality of genomes including a plurality of positive genomes that contain an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not contain an ortholog of an anchor gene of the BGC, the anchor gene being known to be associated with the BGC; ii) one or more phylogenetic characteristics of the last common ancestor (LCA) of pETaG homologs in a phylogenetic tree of multiple genomes; iii) one or more scores indicating co-occurrence of pETaG orthologs and anchor gene orthologs among multiple positive genomes; iv) one or more scores indicating coevolution of sequence diversity between orthologs of pETaG with respect to sequence diversity between orthologs of the anchor gene in positive genomes that contain both orthologs of pETaG and orthologs of the anchor gene; and v) one or more scores indicating the copy number of pETaG homologs in multiple positive genomes and the copy number of pETaG homologs in multiple negative genomes determining one or more parameters selected from b) determining the likelihood that pETaG is a resistance gene for a secondary metabolite produced by the BGC based on one or more parameters; 4. A computer-implemented method comprising: 46. ​​The computer-implemented method of embodiment 45, wherein pETaG co-localizes with a BGC in the query genome. 47. The computer-implemented method of embodiment 46, wherein pETaG is not involved in the production of secondary metabolites by the BGC. 48. The computer-implemented method of any one of embodiments 45 to 47, wherein the anchor gene is a core synthase gene of the BGC. 49. The computer-implemented method of any one of embodiments 45 to 48, comprising, for each of a plurality of pETaGs, determining the likelihood that the pETaG is a resistance gene for a secondary metabolite produced by a BGC in the target genome. 50.a) identifying putative BGCs in multiple genomes that have pairwise sequence similarity above a threshold; b) identifying a non-biosynthetic gene that co-localizes with an ortholog of the anchor gene in the putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in the organism of interest, and the non-biosynthetic gene is not involved in the production of a secondary metabolite by the BGC; c) for each of a plurality of query genes, identifying a non-biosynthetic gene encoding a protein having the highest sequence similarity to the protein of each of the target genes as pETaG, and identifying a genome encoding the non-biosynthetic gene as a target genome; d) for each of the plurality of query genes, determining the likelihood that each pETaG is a resistance gene for a secondary metabolite produced by each BGC in each target genome; 50. The computer-implemented method of embodiment 49, comprising: 51.a) clustering the genomes in the database into a number of clusters, each cluster containing genomes that have pairwise sequence similarity above a threshold; b) for each of the multiple clusters, i) identifying a non-biosynthetic gene that co-localizes with an ortholog of an anchor gene in a putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in an organism of interest, and the non-biosynthetic gene is not involved in the production of a secondary metabolite by the BGC; ii) for each of a plurality of query genes, identifying as a candidate pETaG a non-biosynthetic gene encoding a protein having the highest sequence similarity to the protein of each of the query genes; c) clustering the candidate pETaGs into a plurality of clusters based on sequence similarity between the pETaGs, and identifying the candidate pETaGs encoding proteins with the highest sequence similarity to the proteins of the respective query genes in each cluster as pETaGs, and each genome encoding the pETaGs as a target genome; d) for each of the plurality of query genes, determining the likelihood that each pETaG is a resistance gene for a secondary metabolite produced by each BGC in each target genome; 50. The computer-implemented method of embodiment 49, comprising: 52. The computer-implemented method of embodiment 50 or embodiment 51, wherein the threshold is a pairwise sequence similarity of at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least at least 98%. 53. The computer-implemented method of any one of embodiments 50-52, wherein each of the identified non-biosynthetic genes encodes a protein having at least about 30% sequence identity to the protein encoded by the respective query gene. 54. The computer-implemented method of any one of embodiments 50-53, wherein the plurality of query genes are all protein-coding genes in the organism of interest. 55. The computer-implemented method of any one of embodiments 50 to 54, wherein the organism of interest is a mammal. 56. The computer-implemented method of embodiment 55, wherein the organism of interest is a human. 57. The computer-implemented method of any one of embodiments 50-54, wherein the organism of interest is a reptile, bird, amphibian, animal, plant, fungus or bacterium. 58. The computer-implemented method of embodiment 55 or embodiment 56, wherein the plurality of genomes are fungal genomes. 59. The computer-implemented method of embodiment 55 or embodiment 56, wherein the plurality of genomes are bacterial genomes. 60. The computer-implemented method of embodiment 55 or embodiment 56, wherein the plurality of genomes are plant genomes. 61. The computer-implemented method of any one of embodiments 54-60, wherein each of the plurality of clusters comprises about 10 to about 100 genomes. 62. A computer-implemented method for identifying a druggable target in an organism of interest, comprising performing any one of the methods of embodiments 45 to 61, and identifying the query gene as a druggable target based on the likelihood that the respective pETaG of the query gene is a resistance gene to a secondary metabolite produced by a BGC in the target genome. 63. The computer-implemented method of embodiment 62, further comprising identifying the secondary metabolite or an analog thereof as a small molecule modulator of the query gene or a protein encoded by the query gene. 64. The computer-implemented method of embodiment 63, further comprising contacting the secondary metabolite or an analog thereof with a protein encoded by the query gene, and detecting the activity of the protein encoded by the query gene. 65. The computer-implemented method of any one of embodiments 45 to 64, wherein the number of positive genomes is equal to the number of negative genomes. 66. The computer-implemented method of embodiment 65, comprising selecting a plurality of positive genomes and a plurality of negative genomes from a database of genomes. 67. The computer-implemented method of embodiment 66, comprising clustering a database of genomes into a plurality of clusters based on sequence similarity, and selecting one positive genome per cluster to provide a plurality of positive genomes. 68. The computer-implemented method of embodiment 67, comprising selecting a negative genome having the highest sequence similarity to each positive genome in the cluster. 69. The computer-implemented method of any one of embodiments 66 to 68, wherein the average pairwise sequence identity percentage of orthologs of one or more single-copy genes in the positive genome is about 95% or less, and / or the average pairwise sequence identity percentage of orthologs of one or more single-copy genes in the negative genome is about 95% or less. 70. The computer-implemented method of any one of embodiments 45 to 69, wherein the number of positive genomes is at least 5. 71. The computer-implemented method of any one of embodiments 45-70, wherein the one or more parameters include a likelihood that pETaG is associated with the BGC based on the presence or absence of an ortholog of each of the multiple query genes in the BGC in the multiple different genomes. 72. Determining the likelihood that pETaG is associated with BGC a) receiving a grid representation including a plurality of cells arranged according to a first axis and a second axis, the first axis corresponding to a plurality of genomes, the second axis corresponding to a plurality of query genes in a BGC in a query genome, each cell comprising: i) the presence or absence of an ortholog of each query gene in each genome; and ii) sequence similarity of orthologs to each query gene; iii) whether the orthologs of each query gene colocalize with the orthologs of the anchor genes in each genome; receiving, b) inputting the grid representation into a machine learning model, where the machine learning model is trained to determine a likelihood that the pETaG is associated with the BGC based on values ​​of a plurality of cells in the grid representation, thereby providing a likelihood that the pETaG is associated with the BGC; 72. The computer-implemented method of embodiment 71, comprising: 73.a) Identifying a putative BGC containing pETaG from a library of putative BGCs derived from multiple genomes, and identifying the longest biosynthetic gene in the putative BGC as a core synthase gene; b) obtaining a plurality of positive genomes comprising an ortholog of the core synthase gene and a plurality of negative genomes not comprising an ortholog of the core synthase gene, wherein the plurality of positive genomes have pairwise sequence similarity below a threshold, and the plurality of negative genomes are selected based on sequence similarity to the plurality of positive genomes; c) creating a grid representation comprising a plurality of cells arranged according to a first axis and a second axis, the first axis corresponding to all protein-coding genes co-localized with the core synthase gene in the putative BGC in the query genome, the second axis corresponding to a plurality of positive genomes and a plurality of negative genomes, each cell comprising: i) the presence or absence of orthologs of each protein-coding gene in each genome; and ii) sequence similarity of the orthologs to their respective protein-coding genes; iii) whether orthologs of each protein-coding gene colocalize with orthologs of core synthase genes in each genome; creating a grid representation, the grid representation being calculated based on 73. The computer-implemented method of embodiment 72, further comprising: 74. The computer-implemented method of embodiment 72 or embodiment 73, wherein the machine learning model is a classification model configured to output a probability for each of a plurality of predefined likelihood categories. 75. The computer-implemented method of embodiment 74, wherein the classification model is a long short-term memory (LSTM) model. 76. The computer-implemented method of embodiment 74, wherein the classification model is a convolutional neural network (CNN). 77. The computer-implemented method of embodiment 74, wherein the classification model is an artificial neural network (ANN), a multi-layer perceptron (MLP), a deep neural network (DNN), a vision transformer model, a generative adversarial network (GAN) model, a variational autoencoder model, or a latent diffusion model. 78. The computer-implemented method of any one of embodiments 74-77, wherein the plurality of predefined likelihood categories include: (1) high likelihood, (2) somewhat high likelihood, (3) somewhat low likelihood, and (4) low likelihood. 79. The computer-implemented method of any one of embodiments 45 to 78, wherein the one or more parameters include one or more phylogenetic characteristics of the last common ancestor (LCA) of the homologs of pETaG in the phylogenetic trees of the multiple positive and negative genomes. 80. The computer-implemented method of embodiment 79, wherein the one or more phylogenetic features are selected from the group consisting of the average copy number difference (CND) between genes in the multiple positive genomes and genes in the multiple negative genomes and a value determined from the multiple positive genomes, the ratio of the average to LCA, the ratio of the standard deviation to LCA, the ratio of the average of the adjacent distances, the standard deviation of the ratio of the adjacent distances, and the sum of the clade ratios. 81. The computer-implemented method of any one of embodiments 45-80, wherein the one or more parameters include one or more scores indicating co-occurrence of orthologs of pETaG and orthologs of the anchor gene in multiple positive genomes. 82. The computer-implemented method of embodiment 81, wherein the one or more scores indicative of co-occurrence are selected from the group consisting of co-occurrence pETaG distance, co-occurrence pETaG rank, co-occurrence core distance, and co-occurrence core rank. 83. The computer-implemented method of any one of embodiments 45-82, wherein the one or more parameters include one or more scores indicative of co-evolution of sequence diversity between orthologs of pETaG with respect to sequence diversity between orthologs of the anchor gene in a positive genome that contains both an ortholog of pETaG and an ortholog of the anchor gene. 84. The computer-implemented method of embodiment 83, wherein the one or more scores indicative of coevolution are selected from the group consisting of coevolution correlation, coevolution rank, and coevolution gradient. 85. The method of any one of embodiments 79-84, wherein the one or more parameters further comprise one or more characteristics of the plurality of positive genomes and the plurality of negative genomes. 86. The computer-implemented method of embodiment 85, wherein the one or more features are selected from the group consisting of the number of positive genomes, the average pairwise genomic identity (PGI) between the positive genomes, the standard deviation of PGI between the positive genomes, the number of negative genomes, the average PGI between the negative genomes, and the standard deviation of PGI between the negative genomes. 87. The computer-implemented method of any one of embodiments 45-86, wherein determining the likelihood based on the one or more parameters comprises inputting the one or more features into a machine learning model, and the machine learning model is trained to determine the likelihood that pETaG is a resistance gene. 88. The computer-implemented method of embodiment 87, wherein the machine learning model is a deep learning model. 89. The computer-implemented method of embodiment 87, wherein the machine learning model is a decision tree model. 90. The computer-implemented method of embodiment 87, wherein the machine learning model is a Bayesian estimation model. 91. The computer-implemented method of embodiment 87, wherein the machine learning model is a logistic regression model. 92. The computer-implemented method of any one of embodiments 45 to 91, wherein whether a gene co-localizes with an anchor gene of a BGC is determined using antiSMASH. 93. The computer-implemented method of any one of embodiments 45-92, wherein whether a gene co-localizes with an anchor gene of a BGC is determined based on whether the gene is located within a close distance from the anchor gene. 94. The computer-implemented method of embodiment 93, wherein the proximity zone is about 50 kb or less. 95. The computer-implemented method of embodiment 93, wherein the proximity zone is approximately 20 kb. 96. A system comprising: one or more processors; Memory and Equipped with The memory is communicatively coupled to one or more processors and, when executed by the one or more processors, provides the system with: i) receiving as input one or more query sequences, or proxies thereof; ii) receiving as input a selection of one or more target genomes; iii) performing a search of one or more target genomes using the one or more query sequences or proxies thereof to identify putative embedded target gene (pETaG) sequences that are homologs of the one or more query sequences based on a comparison of one or more homology-based metrics for the candidate pETaGs to one or more predetermined homology-based metric thresholds; iv) determining whether a given pETaG is an actual ETaG based on comparative genomics analysis of multiple genomes related to one or more target genomes. The system is configured to store instructions to cause a 97. A system comprising: one or more processors; Memory and Equipped with The memory is communicatively coupled to the one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to perform a method for determining a likelihood that a putative embedded target gene (pETaG) is a resistance gene for a secondary metabolite produced by a biosynthetic gene cluster (BGC) in a query genome, the method comprising: a) Below: i) the likelihood that pETaG is associated with a BGC based on the presence or absence of each of a plurality of orthologs of a plurality of query genes that are co-localized with the BGC in a plurality of different genomes, the plurality of genomes including a plurality of positive genomes that contain an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not contain an ortholog of an anchor gene of the BGC, the anchor gene being known to be associated with the BGC; ii) one or more phylogenetic characteristics of the last common ancestor (LCA) of pETaG homologs in a phylogenetic tree of multiple genomes; iii) one or more scores indicating co-occurrence of pETaG orthologs and anchor gene orthologs among multiple positive genomes; iv) one or more scores indicating coevolution of sequence diversity between orthologs of pETaG with respect to sequence diversity between orthologs of the anchor gene in positive genomes that contain both orthologs of pETaG and orthologs of the anchor gene; and v) one or more scores indicating the copy number of pETaG homologs in multiple positive genomes and the copy number of pETaG homologs in multiple negative genomes determining one or more parameters selected from b) determining the likelihood that pETaG is a resistance gene for a secondary metabolite produced by the BGC based on one or more parameters; Including, the system. 98. A system comprising: one or more processors; Memory and Equipped with A system, wherein the memory is communicatively coupled to one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to perform any one of the methods of embodiments 1 to 95. 99. A non-transitory computer-readable storage medium storing one or more programs, the one or more programs including instructions that, when executed by one or more processors of an electronic device, cause the electronic device to perform any one of the methods of embodiments 1 to 95.

[0296] The above description has been described with reference to specific examples or embodiments for the purpose of explanation. However, the above exemplary description is not intended to be exhaustive or to limit the invention to the precise form disclosed. For clarity and concise description, features are described herein as parts of the same or separate variants. However, it will be understood that the scope of the present disclosure includes variants having all or a partial combination of the described features. Many modifications and variations are possible in light of the above teachings. The variants have been selected and described in order to best explain the principles of the technology and their practical application. This will enable others skilled in the art to best utilize the technology and various variants with various modifications suited to the particular use envisaged.

[0297] Although the present disclosure and examples have been fully described with reference to the accompanying drawings, it should be noted that various changes and modifications will become apparent to those skilled in the art. Such changes and modifications should be understood to be included within the scope of the present disclosure and examples as defined by the claims. Finally, the entire disclosures of the patents and publications referenced in this application are incorporated herein by reference.

Claims

1. 1. A computer-implemented method for determining the likelihood that a putative embedded target gene (pETaG) is a resistance gene for a secondary metabolite produced by a biosynthetic gene cluster (BGC) in a query genome, comprising: a) Below: i) the likelihood that the pETaG is associated with the BGC based on the presence or absence of orthologs of each of a plurality of query genes that co-localize with the BGC in a plurality of different genomes, the plurality of genomes comprising a plurality of positive genomes that contain an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not contain an ortholog of the anchor gene of the BGC, and the anchor gene is known to be associated with the BGC; ii) one or more phylogenetic characteristics of the last common ancestor (LCA) of the pETaG homologs in the phylogenetic tree of the plurality of genomes; iii) one or more scores indicating co-occurrence of the ortholog of the pETaG and the ortholog of the anchor gene among the plurality of positive genomes; iv) one or more scores indicating co-evolution of sequence diversity between orthologs of the pETaG with respect to sequence diversity between orthologs of the anchor gene in a positive genome containing both an ortholog of the pETaG and an ortholog of the anchor gene; and v) one or more scores indicating the copy number of the pETaG homologue in the plurality of positive genomes and the copy number of the pETaG homologue in the plurality of negative genomes. determining one or more parameters selected from b) determining the likelihood that the pETaG is a resistance gene to the secondary metabolite produced by BGC based on the one or more parameters; 11. A computer-implemented method comprising:

2. The computer-implemented method of claim 1 , wherein the pETaG co-localizes with the BGC in the query genome.

3. The computer-implemented method of claim 2 , wherein the pETaG is not involved in the production of the secondary metabolite by the BGC.

4. 2. The computer-implemented method of claim 1, wherein the anchor gene is a core synthase gene of the BGC.

5. 2. The computer-implemented method of claim 1, comprising determining, for each of a plurality of pETaGs, a likelihood that the pETaG is a resistance gene to a secondary metabolite produced by a BGC in the target genome.

6. a) identifying putative BGCs in multiple genomes that have pairwise sequence similarity above a threshold; b) identifying a non-biosynthetic gene that co-localizes with an ortholog of the anchor gene in the putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in the organism of interest, and the non-biosynthetic gene is not involved in the production of secondary metabolites by the BGC; c) for each of the plurality of query genes, identifying the non-biosynthetic gene encoding a protein having the highest sequence similarity to the protein of the respective target gene as the pETaG, and identifying the genome encoding the non-biosynthetic gene as the target genome; d) for each of the plurality of query genes, determining the likelihood that each of the pETaGs is a resistance gene to a secondary metabolite produced by each of the BGCs in each of the target genomes; The computer-implemented method of claim 5 , comprising:

7. a) clustering the genomes in the database into a plurality of clusters, each cluster containing genomes with pairwise sequence similarity above a threshold; b) for each of said plurality of clusters: i) identifying a non-biosynthetic gene that co-localizes with an ortholog of the anchor gene in the putative BGC, wherein the non-biosynthetic gene is homologous to any one of a plurality of query genes in an organism of interest, and the non-biosynthetic gene is not involved in the production of secondary metabolites by the BGC; ii) for each of a plurality of query genes, identifying as candidate pETaG the non-biosynthetic gene encoding a protein having the highest sequence similarity to the protein of the respective query gene; c) clustering the candidate pETaGs into a plurality of clusters based on sequence similarity among the pETaGs, and identifying the candidate pETaG encoding a protein having the highest sequence similarity to the protein of the respective query gene in each cluster as the pETaG, and each genome encoding the pETaG as the target genome; d) for each of the plurality of query genes, determining the likelihood that each of the pETaGs is a resistance gene to a secondary metabolite produced by each of the BGCs in each of the target genomes; The computer-implemented method of claim 5 , comprising:

8. 7. The computer-implemented method of claim 6, wherein the threshold is a pairwise sequence similarity of at least 70%, at least 75%, at least 80%, at least 85%, at least 90%, at least 95%, or at least 98%.

9. 7. The computer-implemented method of claim 6, wherein each of the identified non-biosynthetic genes encodes a protein having at least about 30% sequence identity to the protein encoded by the respective query gene.

10. The computer-implemented method of claim 6 , wherein the plurality of query genes is all protein-coding genes in the organism of interest.

11. The computer-implemented method of claim 6 , wherein the organism of interest is a mammal.

12. The computer-implemented method of claim 11 , wherein the organism of interest is a human.

13. The computer-implemented method of claim 6 , wherein the organism of interest is a reptile, bird, amphibian, plant, fungus, or bacterium.

14. The computer-implemented method of claim 12 , wherein the plurality of genomes are fungal genomes, bacterial genomes, or plant genomes.

15. 11. The computer-implemented method of claim 10, wherein each of the plurality of clusters comprises about 10 to about 100 genomes.

16. 10. A computer-implemented method for identifying a druggable target in an organism of interest, comprising: performing the method of claim 1; and identifying the query gene as a druggable target based on the likelihood that the respective pETaG of the query gene is a resistance gene to a secondary metabolite produced by the BGC in the target genome.

17. 17. The computer-implemented method of claim 16, further comprising identifying the secondary metabolite or analog thereof as a small molecule modulator of the query gene or the protein encoded by the query gene.

18. 18. The computer-implemented method of claim 17, further comprising contacting the secondary metabolite or an analog thereof with a protein encoded by the query gene; and detecting the activity of the protein encoded by the query gene.

19. The computer-implemented method of claim 1 , wherein the number of positive genomes is equal to the number of negative genomes.

20. 20. The computer-implemented method of claim 19, comprising selecting a plurality of positive genomes and a plurality of negative genomes from a database of genomes.

21. 21. The computer-implemented method of claim 20, comprising clustering the database of genomes into a plurality of clusters based on sequence similarity, and selecting one positive genome per cluster to provide the plurality of positive genomes.

22. 22. The computer-implemented method of claim 21, comprising selecting a negative genome that has the highest sequence similarity to each positive genome in the cluster.

23. The computer-implemented method of claim 20, wherein the average pairwise sequence identity percentage of orthologs of one or more single-copy genes in the positive genome is about 95% or less and / or the average pairwise sequence identity percentage of orthologs of one or more single-copy genes in the negative genome is about 95% or less.

24. The computer-implemented method of claim 1 , wherein the number of positive genomes is at least five.

25. 2. The computer-implemented method of claim 1, wherein the one or more parameters comprise a likelihood that the pETaG is associated with the BGC based on the presence or absence of orthologs of each of a plurality of query genes in the BGC in a plurality of different genomes.

26. 2. The computer-implemented method of claim 1, wherein the one or more parameters comprise one or more phylogenetic characteristics of the last common ancestor (LCA) of the pETaG homologs in the phylogenetic trees of the plurality of positive and negative genomes.

27. 27. The computer-implemented method of claim 26, wherein the one or more phylogenetic features are selected from the group consisting of the mean copy number difference (CND) between genes in the plurality of positive genomes and genes in the plurality of negative genomes and a value determined from the plurality of positive genomes, the mean ratio to LCA, the standard deviation ratio to LCA, the mean ratio of adjacent distances, the standard deviation of adjacent distance ratios, and the sum of clade ratios.

28. 2. The computer-implemented method of claim 1, wherein the one or more parameters comprise one or more scores indicative of co-occurrence of the ortholog of the pETaG and the ortholog of the anchor gene in the plurality of positive genomes.

29. 29. The computer-implemented method of claim 28, wherein the one or more scores indicative of co-occurrence are selected from the group consisting of co-occurrence pETaG distance, co-occurrence pETaG rank, co-occurrence core distance, and co-occurrence core rank.

30. 2. The computer-implemented method of claim 1, wherein the one or more parameters comprise one or more scores indicative of co-evolution of sequence diversity between orthologs of the pETaG with respect to sequence diversity between orthologs of the anchor gene in a positive genome containing both an ortholog of the pETaG and an ortholog of the anchor gene.

31. 31. The computer-implemented method of claim 30, wherein the one or more scores indicative of coevolution are selected from the group consisting of coevolution correlation, coevolution rank, and coevolution gradient.

32. 27. The method of claim 26, wherein the one or more parameters further comprise one or more characteristics of the plurality of positive genomes and the plurality of negative genomes.

33. 33. The computer-implemented method of claim 32, wherein the one or more features are selected from the group consisting of the number of positive genomes, the mean pairwise genome identity (PGI) between the positive genomes, the standard deviation of PGI between the positive genomes, the number of negative genomes, the mean PGI between the negative genomes, and the standard deviation of PGI between the negative genomes.

34. 2. The computer-implemented method of claim 1, wherein determining the likelihood based on the one or more parameters comprises inputting the one or more features into a machine learning model, the machine learning model being trained to determine the likelihood that the pETaG is a resistance gene.

35. 35. The computer-implemented method of claim 34, wherein the machine learning model is a deep learning model.

36. 35. The computer-implemented method of claim 34, wherein the machine learning model is a decision tree model.

37. 35. The computer-implemented method of claim 34, wherein the machine learning model is a Bayesian inference model.

38. 35. The computer-implemented method of claim 34, wherein the machine learning model is a logistic regression model.

39. 2. The computer-implemented method of claim 1, wherein whether a gene co-localizes with an anchor gene of a BGC is determined using antiSMASH.

40. 2. The computer-implemented method of claim 1, wherein whether a gene co-localizes with an anchor gene of a BGC is determined based on whether the gene is located within a close distance from the anchor gene.

41. 41. The computer-implemented method of claim 40, wherein the proximity zone is about 50 kb or less.

42. 41. The computer-implemented method of claim 40, wherein the proximity zone is approximately 20 kb.

43. 1. A system comprising: one or more processors; Memory and Equipped with The memory is communicatively coupled to the one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to perform a method for determining the likelihood that a putative embedded target gene (pETaG) is a resistance gene for a secondary metabolite produced by a biosynthetic gene cluster (BGC) in a query genome, the method comprising: a) Below: i) the likelihood that the pETaG is associated with the BGC based on the presence or absence of orthologs of each of a plurality of query genes that co-localize with the BGC in a plurality of different genomes, the plurality of genomes comprising a plurality of positive genomes that contain an ortholog of an anchor gene of the BGC and a plurality of negative genomes that do not contain an ortholog of the anchor gene of the BGC, and the anchor gene is known to be associated with the BGC; ii) one or more phylogenetic characteristics of the last common ancestor (LCA) of the pETaG homologs in the phylogenetic tree of the plurality of genomes; iii) one or more scores indicating co-occurrence of the ortholog of the pETaG and the ortholog of the anchor gene among the plurality of positive genomes; iv) one or more scores indicating co-evolution of sequence diversity between orthologs of the pETaG with respect to sequence diversity between orthologs of the anchor gene in a positive genome containing both an ortholog of the pETaG and an ortholog of the anchor gene; and v) one or more scores indicating the copy number of the pETaG homologue in the plurality of positive genomes and the copy number of the pETaG homologue in the plurality of negative genomes. determining one or more parameters selected from b) determining the likelihood that the pETaG is a resistance gene to the secondary metabolite produced by the BGC based on the one or more parameters; Including, the system.

44. 1. A system comprising: one or more processors; Memory and Equipped with 10. A system, wherein the memory is communicatively coupled to the one or more processors and configured to store instructions that, when executed by the one or more processors, cause the system to perform the method of claim 1.

45. 10. A non-transitory computer-readable storage medium storing one or more programs comprising instructions that, when executed by one or more processors of an electronic device, cause the electronic device to perform the method of claim 1.

Images