Methods and systems for identifying genes associated with biosynthetic gene clusters
Patent Information
- Application Number
- JP2024527067
- Authority / Receiving Office
- JP · JP
- Patent Type
- Applications
- Current Assignee / Owner
- Priority Date
- 2021-11-10
- Filing Date
- 2022-11-04
- Publication Date
- 2025-10-27
AI Technical Summary
Current methods struggle to accurately define the genomic boundaries of biosynthetic gene clusters (BGCs) and identify non-biosynthetic genes associated with them, which are crucial for understanding the function and therapeutic potential of secondary metabolites produced by microorganisms.
A computer-implemented method using grid representations and machine learning models, such as LSTM, to determine the likelihood of gene association with BGCs by analyzing colocalization and sequence similarity across diverse genomes, enhancing the accuracy of identifying embedded genes and resistance genes.
Improves the precision of defining BGC boundaries and identifies functionally related genes, enabling the identification of resistance genes and potential therapeutic targets, thereby facilitating drug discovery and understanding microbial interactions.
Smart Images

Figure 00000000_0000_ABST
Abstract
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] Disclosed herein are exemplary methods, systems, and non-transitory storage media for identifying genes associated with (or embedded in) a gene cluster (e.g., a biosynthetic gene cluster (BGC) involved in the synthesis of a primary or secondary metabolite). The disclosed methods and systems can be used to determine the boundaries of a gene cluster in a genome (e.g., the boundaries of a BGC), to identify embedded genes of a cluster, to identify resistance genes to secondary metabolites produced by a BGC, or to aid in the identification of a small molecule modulator of a target gene (i.e., a gene of interest) or a protein encoded by a target gene.
[0007] One aspect of the present application is a computer-implemented method for determining a likelihood that a putative embedded gene is associated with a gene cluster (e.g., BGC), wherein the putative embedded gene co-localizes with an anchor gene (e.g., core synthase gene) known to be associated with a gene cluster (e.g., BGC) in a query genome, comprising: a) receiving a grid representation (e.g., heat map) including a plurality of cells arranged according to a first axis and a second axis, wherein the first axis corresponds to a plurality of different genomes, and the second axis corresponds to a plurality of query genes that co-localize with an anchor gene (e.g., core synthase gene) of the gene cluster (e.g., BGC) in the query genome, wherein the putative embedded gene is one of a plurality of query gene orthologs, and each cell represents (i) a respective query gene in each genome; (i) receiving a grid representation of a putative embedded gene based on the presence or absence of an ortholog (e.g., bidirectional best hit or "BBH") of the child, (ii) sequence similarity of the ortholog (e.g., BBH) to the respective query gene, and (iii) whether the ortholog (e.g., BBH) of the respective query gene is co-localized with an ortholog (e.g., BBH) of an anchor gene (e.g., core synthase gene) in the respective genome; and (b) inputting the grid representation into a machine learning model, the machine learning model being trained to determine a likelihood that the putative embedded gene is embedded in a gene cluster (e.g., BGC) based on values of a plurality of cells in the grid representation, thereby providing a likelihood that the putative embedded gene is associated with the gene cluster (e.g., BGC). In some embodiments, the anchor gene is a core synthase gene (i.e., the longest biosynthetic gene) of the gene cluster. In some embodiments, the ortholog of the gene is a BBH of the gene. In some embodiments, the grid representation is a data matrix (e.g., a table). In some embodiments, the grid representation is a heat map. In some embodiments, the grid representation is a subset of a larger grid representation (eg, a subsection of a larger data matrix or heatmap).In some embodiments, the grid representation and / or one or more subsections thereof may be used as input for a machine learning model.
[0008] In some embodiments according to any one of the computer-based methods above, the method further comprises creating a grid representation. In some embodiments, the method includes a) identifying a putative gene cluster (e.g., a putative BGC) comprising a putative embedded gene in the query genome from the plurality of genomes, the putative gene cluster comprising an anchor gene (e.g., a core synthase gene) known to be associated with the gene cluster, the anchor gene co-localizing with the putative embedded gene, and b) identifying a plurality of positive genomes comprising an ortholog of the anchor gene (e.g., BBH) and a plurality of negative genomes not comprising an ortholog of the anchor gene (e.g., BBH), the plurality of positive genomes having a pairwise sequence similarity below a threshold, and the plurality of negative genomes being ranked based on sequence similarity or phylogenetic distance to the plurality of positive genomes. and c) creating a grid representation comprising a plurality of cells arranged according to a first axis and a second axis, where the first axis corresponds to all protein-coding genes that co-localize with the anchor gene in the putative gene cluster in the query 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 (e.g., BBH) of the respective protein-coding gene in the respective genome, (2) the sequence similarity of the ortholog (e.g., BBH) to the respective protein-coding gene, and (3) whether the ortholog of the respective protein-coding gene co-localizes with the ortholog of the anchor gene in the respective genome.
[0009] In some embodiments according to any one of the computer-based methods above, 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., a 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.
[0010] In some embodiments according to any one of the computer-based methods above, the machine learning model is a regression model configured to output a probability that a putative embedded gene is associated with a gene cluster (e.g., a BGC). In some embodiments, the regression model is a logistic regression model.
[0011] In some embodiments according to any one of the computer-based methods above, the method further comprises displaying the grid representation and the likelihood.
[0012] In some embodiments of any one of the computer-based methods above, the grid representation is ordered along a first axis and / or a second axis. In some embodiments, the grid representation is hierarchically clustered, for example along the first axis or the second axis. In some embodiments, the grid representation is ordered along the first axis based on phylogenetic relationships between the multiple genomes (e.g., based on a phylogenetic tree of the multiple genomes). In some embodiments, the grid representation is ordered along the first axis based on pairwise sequence similarity between the multiple genomes. In some embodiments, the grid representation is ordered along the second axis based on the positions of the multiple query genes in the query genome. In some embodiments, the grid representation is ordered along the second axis based on functional annotations of the multiple query genes in the query genome.
[0013] In some embodiments according to any one of the computer-based methods above, the plurality of genomes includes a plurality of positive genomes each having an ortholog of the anchor gene (e.g., a core synthase gene) and a plurality of negative genomes each having no ortholog of the anchor gene (e.g., a core synthase gene). In some embodiments, the number of positive genomes is equal to the number of negative genomes. In some embodiments, the plurality of positive genomes are selected from a plurality of genome clusters based on sequence similarity of the genomes in the database, and no two positive genomes in the grid representation belong to the same genome cluster. In some embodiments, each negative genome is selected by identifying a genome in the database that has the highest sequence similarity or shortest phylogenetic distance to the positive genome, but does not have an ortholog of the anchor gene. In some embodiments, the average pairwise sequence identity percentage of orthologs of one or more single copy genes in the positive genome is about 99.5% or less (e.g., about any one of 99%, 98%, 95%, 90%, 85%, 80%, 75%, 70%, 65%, 60%, 55%, or 50% 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 99.5% or less (e.g., about any one of 99%, 98%, 95%, 90%, 85%, 80%, 75%, 70%, 65%, 60%, 55%, or 50% or less).
[0014] In some embodiments of any one of the computer-based methods above, the first axis corresponds to at least 2, 4, 8, 16, 20, 30, 40, 50, 75, 100, 150, 200, 250 or more genomes. In some embodiments, the first axis corresponds to about 50 genomes.
[0015] In some embodiments according to any one of the computer-based methods above, the plurality of genomes are fungal genomes. In some embodiments, the plurality of genomes are plant kingdom (green algae and / or plant) genomes. In some embodiments, the plurality of genomes are bacterial genomes. In some embodiments, the plurality of genomes are archaeal genomes. In some embodiments, the plurality of genomes are protozoan genomes. In some embodiments, the plurality of genomes are Chromista (e.g., brown algae, diatoms, cryptophytes, etc.) genomes. In some embodiments, the plurality of genomes are animal genomes.
[0016] In some embodiments of any one of the above computer-based methods, whether a gene co-localizes with an anchor gene (e.g., a core synthase gene) of a gene cluster (e.g., a BGC) is determined using antiSMASH, SMURF, TOUCAN, or deepBGC. In some embodiments, whether a gene co-localizes with an anchor gene (e.g., a core synthase gene) of a gene cluster (e.g., a BGC) is determined using antiSMASH. In some embodiments, whether a gene co-localizes with an anchor gene (e.g., a core synthase gene) of a gene cluster (e.g., a BGC) is determined based on whether the gene is located within a proximal zone upstream or downstream of the anchor gene (e.g., a core synthase gene). In some embodiments, the proximity zone is no more than about any one of 200kb, 100kb, 90kb, 80kb, 70kb, 60kb, 50kb, 45kb, 40kb, 35kb, 30kb, 25kb, 20kb, 15kb, 10kb, or 5kb. In some embodiments, the proximity zone is at least about any one of 5kb, 10kb, 15kb, 20kb, 25kb, 30kb, 35kb, 40kb, 45kb, 50kb, 60kb, 70kb, 80kb, 90kb, 100kb, or more. In some embodiments, the proximity zone is any one of about 5 kb to 20 kb, 5 kb to 50 kb, 5 kb to 100 kb, 5 kb to 200 kb, 20 kb to 50 kb, 20 kb to 100 kb, 20 kb to 200 kb, 50 kb to 100 kb, 50 kb to 200 kb, 10 kb to 50 kb, or 10 kb to 100 kb. In some embodiments, the proximity zone is about 50 kb. In some embodiments, the proximity zone is about 20 kb.
[0017] Another aspect of the present application provides a method for identifying a resistance gene to a secondary metabolite produced by a BGC in a query genome, the method comprising: (a) identifying a putative embedded gene that is not involved in the production of a secondary metabolite by the BGC that co-localizes (e.g., within a proximity zone of 50 kb, 20 kb, or any user-specified distance or less) with an anchor gene (e.g., a core synthase gene) in a BGC in the query genome; (b) performing a method according to any one of the computer-based methods described herein to determine a 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 embedded gene is associated with the BGC.
[0018] Another aspect of the application provides a method of identifying a small molecule modulator of a mammalian target gene (i.e., a mammalian gene of interest), comprising: (a) identifying a homologous gene of the mammalian target gene in the fungal genome that is co-localized 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 a method according to any one of the computer-based methods described herein to determine the likelihood that the homologous gene is associated with the BGC; and (c) identifying a secondary metabolite or analog thereof as a small molecule modulator of the mammalian target gene based at least in part on the likelihood that the homologous gene is associated with the BGC. In some embodiments, the homologous gene encodes a protein having at least about 20%, 25%, 30%, 35%, 40%, 45%, 50%, 55%, 60%, 65%, 70%, 85%, 90%, 95% or more sequence identity or homology to a protein encoded by the mammalian target gene. In some embodiments, the homologous gene encodes a protein having at least about 30% sequence identity or homology to the protein encoded by the mammalian target gene. In some embodiments, the method further comprises contacting a secondary metabolite or analog thereof with the protein encoded by the mammalian target gene and detecting an activity of the protein encoded by the mammalian target gene (e.g., binding to the secondary metabolite or analog thereof).
[0019] Another aspect of the present application provides a method for identifying a plurality of genes associated with a BGC, the method comprising: (a) identifying a plurality of query genes that co-localize with an anchor gene (e.g., a core synthase gene) of the BGC in a query genome; (b) for each of the plurality of query genes, determining a likelihood that each query gene is associated with the BGC using a method according to any one of the computer-based methods described above; and (c) identifying query genes having a high likelihood, above a threshold, of being associated with the BGC as the plurality of genes associated with the BGC.
[0020] One aspect of the present application includes the steps of: a) identifying a putative biosynthetic gene cluster (BGC) comprising a putative embedded gene in a query genome from a plurality of genomes, wherein the putative BGC comprises an anchor gene known to be associated with the BGC, and the anchor gene co-localizes with the putative embedded gene; b) identifying a plurality of positive genomes comprising an ortholog of the anchor gene and a plurality of negative genomes not comprising an ortholog of the anchor 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 or phylogenetic distance to the plurality of positive genomes; and c) identifying a first axis and a second axis and a putative embedded gene in a query genome from a plurality of genomes, wherein the putative BGC comprises an anchor gene known to be associated with the BGC, and the anchor gene co-localizes with the putative embedded gene; and creating a grid representation including 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 anchor genes in a putative BGC in a query genome, and the second axis corresponding to a plurality of positive genomes and a plurality of negative genomes, each cell being based on (1) the presence or absence of an ortholog of each 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 each protein-coding gene co-localizes with an ortholog of the anchor gene in the respective genome. In some embodiments, the anchor gene is a core synthase gene (e.g., the longest biosynthetic gene) of the BGC. In some embodiments, the ortholog of the gene is a bidirectional best hit (BBH) of the gene. In some embodiments, the method further comprises hierarchically clustering the grid representation. In some embodiments, the method further includes ordering the grid representation along a first axis, e.g., based on phylogenetic relationships between the multiple genomes (e.g., a phylogenetic tree) or based on pairwise sequence similarity between the multiple genomes (e.g., a cladogram derived from pairwise genome comparison). In some embodiments, the grid representation is a data matrix (e.g., a table). In some embodiments, the grid representation is a heat map. In some embodiments, the method further includes displaying the grid representation.In some embodiments, the grid representation is a subset of a larger grid representation (e.g., a subsection of a larger data matrix or heatmap). In some embodiments, the grid representation and / or one or more subsections thereof can be used as input for a machine learning model.
[0021] Also 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 techniques or methods, or one or more computer-based steps of the methods, described herein.
[0022] 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.
[0023] 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.
[0024] 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]
[0025] [Figure 1] FIG. 1 shows exemplary putative biosynthetic gene clusters (BGCs) predicted by antiSMASH.
[0026] [Diagram 2] FIG. 1 illustrates an exemplary method for generating a grid representation (e.g., a heat map) of genomic data.
[0027] [Diagram 3] 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 gene cluster (e.g., a BGC) in a query genome is associated with the BGC, according to some examples.
[0028] [Figure 4A-1] FIG. 1 illustrates an exemplary heat map. [Figure 4A-2] FIG. 1 illustrates an exemplary heat map. [Figure 4A-3] FIG. 1 illustrates an exemplary heat map.
[0029] [Figure 4B-1] FIG. 1 illustrates an exemplary long short-term memory (LSTM) model used to classify an input heatmap into one of multiple likelihood categories (e.g., four likelihood categories or hierarchies). FIG. 2 illustrates an exemplary LSTM model with a series of memory hierarchies. [Figure 4B-2] FIG. 1 illustrates an exemplary long short-term memory (LSTM) model used to classify an input heatmap into one of multiple likelihood categories (e.g., four likelihood categories or hierarchies). FIG. 2 illustrates an exemplary LSTM model with a series of memory hierarchies. [Figure 4C-1] FIG. 1 illustrates an exemplary long short-term memory (LSTM) model used to classify an input heatmap into one of multiple likelihood categories (e.g., four likelihood categories or a hierarchy). FIG. 2 illustrates an exemplary diagram of the output of the LSTM model. [Figure 4C-2]FIG. 1 illustrates an exemplary long short-term memory (LSTM) model used to classify an input heatmap into one of multiple likelihood categories (e.g., four likelihood categories or a hierarchy). FIG. 2 illustrates an exemplary diagram of the output of the LSTM model.
[0030] [Figure 5A] 13 is a table comparing manual and machine learning-based classification of heatmaps for “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3” categories, respectively.
[0031] [Figure 5B] FIG. 1 shows a table comparing manual and machine learning based classification of heat maps for “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”, including positive predictive value, negative predictive value, sensitivity value, and specificity value.
[0032] [Figure 6A] FIG. 13 shows an exemplary heatmap of the BGC of lovastatin, identifying the true boundaries of the lovastatin BGC compared to the predicted lovastatin BGC by antiSMASH. [Figure 6B] FIG. 13 shows an exemplary heatmap of the BGC of lovastatin, identifying the true boundaries of the lovastatin BGC compared to the predicted lovastatin BGC by antiSMASH.
[0033] [Figure 7A-1] 1 shows an example heatmap that has been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap classified as "Tier A+". [Figure 7A-2] 1 shows an example heatmap that has been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap classified as "Tier A+". [Figure 7B-1] 1 shows an example heatmap that has been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap that is classified as "Tier 1". [Figure 7B-2] 1 shows an example heatmap that has been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap that is classified as "Tier 1". [Figure 7C-1] 1 shows an example heatmap that has been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap that is classified as "Tier 2". [Figure 7C-2] 1 shows an example heatmap that has been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap that is classified as "Tier 2". [Figure 7D-1] 1 shows an example heatmap of samples that have been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap of samples that have been classified as "Tier 3". [Figure 7D-2] 1 shows an example heatmap of samples that have been manually reviewed and classified as different likelihood categories. FIG. 2 shows an example heatmap of samples that have been classified as "Tier 3".
[0034] [Figure 8A] FIG. 1 illustrates a data table of a set of features (e.g., up to 26 or more features) that may be organized into a data table and utilized to train a neural network.
[0035] [Figure 8B] FIG. 1 illustrates the initial training stage of a neural network trained to output a probability value that a putative embedded gene (e.g., pETaG) is associated with a BGC (i.e., an "embedded gene probability value" (e.g., ETaG) probability value").
[0036] [Figure 8C] FIG. 13 illustrates a further training stage of a neural network trained to output probability values (e.g., ETaG probability values) for embedded genes.
[0037] [Figure 8D]FIG. 1 illustrates the inference stage of a neural network trained to output ETaG probability values.
[0038] [Figure 8E] 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. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
[0039] The present disclosure provides methods for identifying genes associated with gene clusters (e.g., biosynthetic gene clusters (BGCs)) by comparative genomics analysis using grid representations (e.g., heat maps). The grid representations allow for visual and machine learning-based assessment of co-occurrence and co-localization between genes found in close proximity to the gene cluster (BGC), its biosynthetic or non-biosynthetic genes, or other genes of interest.
[0040] 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.
[0041] 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, e.g., bidirectional best hits (BBHs), of genes that co-localize with anchor genes (e.g., core synthase genes) known to be associated with gene clusters (e.g., BGCs) across multiple diverse genomes to determine the likelihood that an ortholog of a query gene (e.g., a gene in a reference or query genome) that co-localizes with an anchor gene (e.g., core synthase gene) of a gene cluster (e.g., BGC) in a query genome is associated with the gene cluster (BGC) in the targeted or queried 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.
[0042] The methods 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 BGCs in different genomes by identifying genes associated with the BGC.
[0043] Furthermore, the method can be used to identify resistance genes embedded in the BGC 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 the BGC 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.
[0044] definition 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.
[0045] The terms "biosynthetic gene cluster" or "BGC" are used interchangeably herein and refer to a locally clustered group of one or more genes that together encode a biosynthetic pathway for the production of a secondary metabolite. Exemplary BGCs include, but are not limited to, biosynthetic gene clusters for synthesizing nonribosomal 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. BGCs contain genes encoding 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 not involved in the biosynthesis of secondary metabolites that are 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. "Anchor genes" refer to biosynthetic or non-biosynthetic genes that are known to be co-localized and functionally related (i.e., associated) with the BGC.
[0046] The term "co-localized" refers to the presence of two or more genes in close proximity to one another, e.g., genes that are separated in the genome by about 200 kb or less, about 100 kb or less, about 50 kb or less, about 40 kb or less, about 30 kb or less, about 20 kb or less, about 10 kb or less, about 5 kb or less, or less.
[0047] The term "homolog" refers to a gene that is part of a set of genes whose gene sequences (i.e., nucleic acid sequences) and / or the sequences of their protein products are inherited from a common origin. Homologs can arise through speciation events, or through gene duplication events, or through horizontal gene transfer events. Homologs can be identified by phylogenetic methods, or through the identification of common functional domains in aligned nucleic acid or protein sequences, or through sequence comparison.
[0048] The term "ortholog" refers to two or more genes that are predicted to have evolved from a common ancestral gene by speciation. 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 (or BBH) of the second gene, and the second gene is the bidirectional best hit (or BBH) of the first gene. Identifying BBH is a commonly used method to infer orthology.
[0049] "Percent sequence identity" or "percent sequence homology" with respect to the protein sequences described herein is defined as the percentage of amino acid residues in a candidate polypeptide sequence that are identical or homologous to the amino acid residues in the polypeptide to which it is compared, after aligning the sequences and considering any conservative substitutions as part of the sequence identity. Homology between different amino acid residues is determined based on an alternative matrix, such as BLOSUM (BLOcks SUbstitution Matrix). Alignment for determining percent amino acid sequence identity can be achieved in various ways known to those skilled in the art, for example, using publicly available computer software such as BLAST, BLAST-2, ALIGN or Megalign (DNASTAR) software. Those skilled in the art can determine the appropriate parameters for measuring alignment, including any algorithms required to achieve maximum alignment over the entire length of the sequences being compared.
[0050] As used herein, "sequence similarity" between two genes means similarity in either the nucleic acid (e.g., DNA or mRNA) sequences encoded by the genes or the amino acid sequences of the gene products.
[0051] In the following description, it is to be understood that the singular forms "a", "an" and "the" used in the following description are intended to include the plural forms unless the context clearly indicates otherwise. It is also to be understood that the term "and / or" as used herein refers to and encompasses any and all possible combinations of one or more of the associated listed items. Furthermore, it is to be understood that the terms "includes", "including", "comprises" and / or "comprising", as used herein, 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.
[0052] 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, or hardware, and when embodied in software, may be downloaded to reside and operate on different platforms used by various operating systems. As will be apparent from the following description, unless otherwise stated, descriptions utilizing terms such as "processing," "calculating," "computing," "determining," "displaying," "generating," and the like throughout the description are 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.
[0053] The section headings used herein are for organizational purposes only and are not to be construed as limiting the subject matter described. The description is presented to enable one of ordinary skill in the art to make and use the invention and is provided in the context of a patent application and its requirements.
[0054] Grid Representation Analysis Method The systems and methods described herein relate to the identification of 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 the BGC across diverse genomes.
[0055] FIG. 1 shows exemplary putative BGC regions predicted by antiSMASH. BGCs contain a set of genes encoding enzymes, including signature biosynthetic enzymes, in a biosynthetic pathway for the production of secondary metabolites, the longest of which is referred to herein as the "core biosynthetic protein." In a BGC, non-biosynthetic genes that co-localize with the biosynthetic genes may be homologs of human proteins that include therapeutic targets of interest. Such non-biosynthetic genes may be functionally related to the secondary metabolites produced by the BGC or may encode functionally unrelated products. A non-biosynthetic gene in a BGC that is a homolog of a human protein of interest is a putative embedded target gene (pETaG). The methods described herein leverage comparative genomics to determine the likelihood that a gene of interest is associated with a BGC.
[0056] FIG. 2 shows an exemplary method 200 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 pETaG is associated with a gene cluster (e.g., a BGC) in the genome.
[0057] FIG. 3 illustrates an exemplary method 300 for determining the likelihood that a putative embedded gene is associated with a gene cluster (e.g., BGC). Process 200 and process 300 are performed, for example, using one or more electronic devices implementing a software platform. In some examples, process 200 and / or process 300 are performed using a client-server system, with blocks of process 200 and / or process 300 being split in any manner between a server and one or more client devices. In some examples, process 200 and / or process 300 are performed using only one client device or only multiple client devices. In process 200 and / or 300, 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 200 and / or process 300. Thus, the operations shown (and described in more detail below) are exemplary in nature and, therefore, should not be considered limiting.
[0058] Grid Representation In block 302 of FIG. 3, an exemplary system (e.g., 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.
[0059] 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 may 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 embodiments, 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, may be stored in a separate table and used as input for the machine learning-based methods described herein. In some embodiments, the grid representation may be a physical representation of the underlying data matrix, such as a heat map, that facilitates visualization of the data.
[0060] 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.
[0061] In some embodiments, the orthologue of a query gene in a given genome can 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, 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.
[0062] 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 the target genome. In some embodiments, 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) of a BGC. In some embodiments, 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 embodiments, the proximal zone is 5 kb or less upstream or downstream of a gene. In some embodiments, the proximal zone is no more than 10 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 15 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 20 kb upstream or downstream of the gene.In some embodiments, the proximal zone is no more than 25 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 30 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 35 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 40 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 45 kb upstream or downstream of the gene. In some embodiments, the proximal zone is no more than 50 kb upstream or downstream of the gene.
[0063] 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.
[0064] For example, the grid representation is genome Q1 to Q nThe first axis (e.g., Y axis) corresponds to the gene G1 to G m and 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.
[0065] 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, genome phylogeny, or the presence or absence of orthologs corresponding to all query genes in the grid representation. In some embodiments, 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.
[0066] In some embodiments, 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 embodiments, pETaG is homologous to an expressed mammalian nucleic acid sequence. In some embodiments, the mammalian nucleic acid sequence is an expressed mammalian nucleic acid sequence. In some embodiments, the mammalian nucleic acid sequence is a mammalian gene. In some embodiments, the mammalian nucleic acid sequence is an expressed mammalian gene. In some embodiments, the mammalian nucleic acid is a human nucleic acid sequence. In some embodiments, the human nucleic acid sequence is an expressed human nucleic acid sequence. In some embodiments, the human nucleic acid sequence is a human gene. In some embodiments, the human nucleic acid sequence is an expressed human gene.
[0067] An example of a genome heat map is shown in Figure 4A-1 to Figure 4A-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 embodiments, the BGC is identified by antiSMASH. In some embodiments, a gene in a 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 designated positive genomes. Half of the genomes do not contain the BBH of the core synthase gene and are designated 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 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 genomes.
[0068] Each genome shown in the grid representation may correspond to an assembled genome, or multiple genome fragments obtained from genome sequencing. In some embodiments, 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.
[0069] The methods described herein are suitable for any genome that contains a BGC. Bacterial, plant, and fungal genomes are known to encode biosynthetic gene clusters. In some embodiments, the query (or reference) genome and the multiple queried (or target) genomes used to generate the grid representation belong to the same kingdom. In some embodiments, the query genome and the multiple queried 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 more related to mammalian genomes than to bacterial or plant genomes. Thus, fungal genomes may be preferred for identifying ETaGs that correspond to human target genes (i.e., human genes of interest) for secondary metabolites produced by BGCs carrying ETaGs.
[0070] At least two genomes are required to construct the grid representation. In some embodiments, the first axis of the grid representation corresponds to at least about any one of 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 150, 200, 250 or more genomes. In some embodiments, the first axis of the grid representation corresponds to at least 20 genomes. In some embodiments, 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 each other in terms of their sequence similarity and / or phylogenetic relationship to generate the grid representation to balance between the performance of the method (e.g., accuracy of prediction) and computational resources.
[0071] 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 embodiments, the positive genome and the negative genome are selected from a database of genomes. In some embodiments, 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 embodiments, 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.
[0072] 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).
[0073] 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.
[0074] 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.
[0075] The grid representation can be generated, for example, using the method illustrated in FIG.
[0076] 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.
[0077] 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.
[0078] 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 embodiments, 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.
[0079] 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.
[0080] 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.
[0081] Gene clusters (e.g., BGCs) protein or gene orthologous group clusters (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 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 may contain, for example, protein or gene IDs from the same COG in each row of a table.
[0082] In block 202 of FIG. 2, 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 202. 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.
[0083] 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.
[0084] 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.
[0085] In block 204 of Figure 2, 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 204. 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.
[0086] 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.
[0087] 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.
[0088] 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.
[0089] 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).
[0090] In block 206 of FIG. 2, 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.
[0091] 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 (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 (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.
[0092] 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.
[0093] 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.
[0094] A final table storing the values of 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.
[0095] The generated grid representation and / or a subset of the generated grid representation (e.g., a heatmap, or a data matrix) can be input into a machine learning model, such as an LSTM, to provide the likelihood that a putative embedded gene (e.g., pETaG) is associated with a BGC.
[0096] Machine learning model for predicting "embeddedness" In block 304 of FIG. 3, 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., a BGC) based on the values of a number of cells in the grid representation.
[0097] 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 embodiments, 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, e.g., Figures 7A-1 and 7A-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, e.g., Figures 7B-1 and 7B-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, e.g., Figures 7C-1 and 7C-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 7D-1 and Figure 7D-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%.
[0098] In block 306 of FIG. 3, 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.
[0099] In block 308 of FIG. 3, 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).
[0100] 4A-1-4A-3 show heatmaps depicting the degree of "embeddedness" (i.e., association with BGCs) of pETaG in a calculated heatmap according to an embodiment 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 the BBH of pETaG and are referred to herein as "positive genomes." Half of the genomes do not contain the BBH of pETaG 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.
[0101] 4B-1-4C-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.
[0102] As shown in FIG. 4B-1 and FIG. 4B-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 embodiments, after computing the heatmap, a vector representation of the data contained in the heatmap (e.g., a table of values indicating one or more patterns, one or more colors, pETaG locations, etc.) may be provided to one or more neural networks to perform an embedding classification of pETaGs in the heatmap. For example, in some embodiments, the one or more neural networks may include, for example, a long short-term memory (LSTM) model, a convolutional neural network (CNN), or other recurrent neural network (RNN) that 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 embodiments, 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")).
[0103] 4B-1-4B-2 and 4C-1-4C-2 illustrate 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 embodiments of the present disclosure. For example, in some embodiments, as illustrated by FIGS. 4B-1-4B-2, an LSTM model can include a series of memory hierarchies, each including, e.g., a respective memory cell. In some embodiments, each memory cell can include, e.g., a memory hierarchies that store 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 embodiments, gates in each memory cell may be provided to optionally pass information through. For example, in some embodiments, 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 embodiments, each respective memory cell may include these gates, for example, to protect and control the cell state.
[0104] In some embodiments, 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 embodiments, the determination may be performed by a sigmoid layer called a forget gate layer, which looks at the input data and determines which of the cell states (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 embodiments, a sigmoid layer, called the input gate layer, decides which values to update, and a tan h layer creates a vector of new candidate values that can be added to the state.
[0105] 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 embodiments, 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 a value between "-1" and "+1" and multiply the value by the output of a sigmoid gate.
[0106] As further illustrated by the prediction tables in Figures 5A and 5B according to an embodiment of the present disclosure, Figures 4C-1 to 4C-2 show examples of output 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")), and (4) "true negative" (e.g., there is a low likelihood that pETaG is associated with BGCs ("Tier 3")).
[0107] Specifically, FIG. 5A shows a table of predictive embeddedness benchmark values for “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”. Similarly, FIG. 5B 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. 5A are calculated from comparing the manually annotated hierarchies with the prediction results from the LSTM model. In the table in FIG. 5B, 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.
[0108] Figures 6A and 6B show an example of 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 7A-1 to 7D-2 show heatmaps classified into Tier A+, Tier 1, Tier 2, and Tier 3, respectively.
[0109] In some embodiments, 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. 8A ) representing the four features “Tier A+”, “Tier 1”, “Tier 2”, and “Tier 3”, and combined with a predefined 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”).
[0110] Machine learning model for predicting ETaG likelihood Specifically, FIG. 8A illustrates a data table including a combination set of features. In some embodiments, the combination set of features (e.g., up to 27 or more features) can be organized in a data table as shown in FIG. 8A and can be utilized to train a 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.) to output an ETaG or pETaG probability value based on an input of the combination set of features shown in FIG. 8A. 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 can be utilized to train a logistic regression model or other types of supervised models. As further illustrated in FIG. 8A, 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 neural network.
[0111] FIG. 8B illustrates an initial training stage of a neural network trained to output ETaG or pETaG probability values, according to an embodiment of the present disclosure. As shown, a training data set corresponding to the features included in the data table of FIG. 8A 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 embodiments, 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 embodiments, the specified activation functions (e.g., computational functions) may specifically determine the value of the output of the input neuron or node.
[0112] In some embodiments, 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 embodiments, the hidden neurons or nodes can constitute a hidden layer of the neural network, and can 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. For example, in some embodiments, the weights of the hidden layer can affect, for example, 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 embodiments, as further illustrated by FIG. 8B, the neural network may be trained, for example, based on a forward propagation technique. In some embodiments, 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).
[0113] FIG. 8C further illustrates a training stage of a neural network trained to output ETaG or pETaG probability values according to an embodiment of the present disclosure. For example, as shown in FIG. 8C, 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 embodiments, 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.
[0114] FIG. 8D illustrates an inference stage of a neural network trained to output ETaG or pETaG probability values, according to an embodiment 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. 8A, 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. 8B, 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 embodiments, the specified activation functions (e.g., computational functions) may specifically determine the value of the output of the input neuron or node. In some embodiments, 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 embodiments, 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 embodiments, as further shown in FIG. 8D, 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 an embodiment 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. 8E illustrates an exemplary 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 an embodiment of the present disclosure. In this manner, the embodiment may identify and determine the likelihood that an ETaG or pETaG is associated with one or more features that represent a BGC.
[0115] Purpose The computer-based methods described herein in the "Grid Representation Analysis Methods" section have a variety of applications.
[0116] In some embodiments, 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 embodiments, the method includes: (1) identifying a plurality of query genes that are co-localized 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 or its ortholog is associated with the corresponding BGC in the plurality of genomes; and (c) identifying the query gene or its ortholog with a specified high likelihood, which is a likelihood above a threshold, associated with the BGC as a 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 embodiments, if the query gene or its orthologue has a combined probability of (1) high likelihood and (2) more than about 50%, 60%, 70%, 80%, 90% or more for the category of moderately high likelihood, the query gene is associated with the BGC. In some embodiments, if the query gene or its orthologue has a probability of (4) more than about 30%, 40%, 50%, 60%, 70%, 80%, 90% or more for the category of low likelihood, the query gene is rejected as not associated with the BGC. The boundaries (i.e., upstream and downstream limits) of the BGC can be determined based on the positions of all genes associated with the BGC determined using this method.
[0117] In some embodiments, 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 putative embedded genes that are not involved in the production of secondary metabolites by the BGC and that are 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 a 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 embodiments, 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 embodiments, 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 embodiments, 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.
[0118] In some embodiments, 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.
[0119] In some embodiments, 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 embodiments, a 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 embodiments, the likelihood that pETaG is associated with the BGC is one of multiple factors used to identify a mammalian (e.g., human) gene as a target of a secondary metabolite of the BGC. In some embodiments, the method further comprises identifying a mammalian (e.g., human) homolog of pETaG in the mammalian (e.g., human) genome. In some embodiments, the method further comprises assaying the effect of a secondary metabolite produced by the BGC or an analog of a BGC product on the mammalian (e.g., human) target.
[0120] In some embodiments, 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 embodiments, 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 containing a BGC) that is co-localized (e.g., within a contiguous 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 embodiments, 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 embodiments, 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 embodiments, 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 embodiments, the activity is the binding of the protein encoded by the mammalian target gene with the secondary metabolite or analog thereof.
[0121] In some embodiments, the secondary metabolite is a product of an enzyme encoded by a BGC or a salt thereof, including a non-naturally occurring salt. In some embodiments, the secondary metabolite or analog thereof is an analog of the product of an enzyme encoded by a BGC, such as a small molecule compound having the same core structure as the secondary metabolite, or a salt thereof.
[0122] In some embodiments, 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.
[0123] In some embodiments, 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.
[0124] In some embodiments, the secondary metabolite is produced by a fungus. In some embodiments, the secondary metabolite is acyclic. In some embodiments, the secondary metabolite is a polyketide. In some embodiments, the secondary metabolite is a terpene compound. In some embodiments, the secondary metabolite is a non-ribosomally synthesized peptide.
[0125] In some embodiments, an analog of a substance (e.g., a secondary metabolite) 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 embodiments, an analog is a substance that can be generated from a reference substance, e.g., by chemical manipulation of the reference substance. In some embodiments, 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 embodiments, 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 embodiments, an analog of a substance is a substance that is substituted at one or more of its substitutable positions.
[0126] In some embodiments, the analog of the product comprises the structural core of the product. In some embodiments, the biosynthetic product is cyclic, e.g., monocyclic, bicyclic, or polycyclic, and the structural core of the product is or comprises a monocyclic, bicyclic, or polycyclic ring system. In some embodiments, the structural core of the product comprises one ring of the bicyclic or polycyclic ring system of the product. In some embodiments, the product is or comprises a polypeptide, and the structural core is the backbone of the polypeptide. In some embodiments, the product is or comprises a polyketide, and the structural core is the backbone of the polyketide. In some embodiments, the analog is a substituted biosynthetic product that comprises one or more suitable substitutents.
[0127] Identification of ETaG In some embodiments, the present disclosure provides methods for identifying embedded target genes ("ETaGs") or mammalian (e.g., human) target genes (i.e., (human) genes of interest) corresponding to ETaGs. 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.
[0128] In some embodiments, the methods described herein are applied to identify ETaGs from fungal genomes. In some embodiments, ETaGs from eukaryotic fungi can have more similarity to mammalian genes than their counterparts (if any) in prokaryotes, such as certain bacteria. In some embodiments, fungi contain and / or contain more therapeutically relevant ETaGs than organisms that are evolutionarily more distant from humans.
[0129] In some embodiments, 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 "Grid Representation Analysis Methods" section; 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 embodiments, pETaG is co-regulated with at least one biosynthetic gene in the BGC. In some embodiments, pETaG is not co-regulated with at least one biosynthetic gene in the BGC. In some embodiments, 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 embodiments, ETaG is compared to mammalian, e.g., human, nucleic acid sequences to identify homologous mammalian nucleic acid sequences. In some embodiments, 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 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 the human target.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, archaeal, fungal or plant targets.
[0130] 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.
[0131] In some embodiments, 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 embodiments, the present disclosure significantly improves the druggability of targets that were previously thought to be undruggable, for example, by identifying their homologous ETaGs in fungi, elucidating the associated biosynthetic gene clusters, and testing the biosynthetic products of the relevant biosynthetic gene clusters, essentially converting them, in some cases, 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).
[0132] 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 embodiments, the ETaG is located no more than about any one of 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.
[0133] In some embodiments, ETaG is a product of an existing target of therapeutic interest or is homologous to a human nucleic acid sequence that encodes it. In some embodiments, ETaG is a product of a novel target of therapeutic interest or is homologous to a human nucleic acid sequence that encodes it. In some embodiments, ETaG is a product of a target that was considered undruggable prior to this disclosure or is homologous to a human nucleic acid sequence that encodes it. In some embodiments, ETaG is a product of a target that was considered undruggable by small molecules prior to this disclosure or is homologous to a human nucleic acid sequence that encodes it.
[0134] In some embodiments, 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 embodiments, 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 embodiments, the homologous portion is at least 50, 100, 150, 200, 500, 1000, 2000, 3000, or 5000 base pairs in length. In some embodiments, 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 embodiments, the mammalian nucleic acid, e.g., a human nucleic acid sequence, is associated with a human disease, disorder, or condition. In some embodiments, such human nucleic acid sequences are existing targets for therapeutic purposes. In some embodiments, such human nucleic acid sequences are novel targets for therapeutic purposes. In some embodiments, such human nucleic acid sequences are targets previously thought to be insensitive to targeting, e.g., by small molecules.
[0135] In some embodiments, 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 the mammalian nucleic acid sequence. In some embodiments, 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 product encoded by the mammalian nucleic acid sequence. In some embodiments, 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 the mammalian nucleic acid sequence.
[0136] In some embodiments, the portion of the protein is a protein domain. In some embodiments, the protein domain is an enzyme domain. In some embodiments, the protein domain interacts with one or more factors, such as small molecules, lipids, carbohydrates, nucleic acids, proteins, etc.
[0137] In some embodiments, a portion of a protein is a functional and / or structural domain that defines the protein family to which the protein belongs. 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.
[0138] In some embodiments, 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 embodiments, the function is an enzymatic activity and the portion of the protein is a set of residues required for the activity. In some embodiments, 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 embodiments, the set of residues interacts with a substrate. In some embodiments, the set of residues interacts with an intermediate. In some embodiments, the set of residues interacts with a product.
[0139] In some embodiments, the function of the protein is interaction with one or more factors, e.g., 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 embodiments, each of the set of residues independently contacts an interacting agent. For example, in some embodiments, each of the residues of the set independently contacts an interacting small molecule. In some embodiments, 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, e.g., via hydrogen bonds, electrostatic forces, van der Waals forces, aromatic stacking, etc. In some embodiments, the interacting agent is another macromolecule. In some embodiments, the interacting agent is a nucleic acid. In some embodiments, the set of residues are residues that contact an interacting nucleic acid, e.g., residues in a transcription factor. In some embodiments, the set of residues are residues that contact an interacting protein.
[0140] In some embodiments, 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 a human target.
[0141] 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, in some embodiments, from fungi to humans as shown in this disclosure.
[0142] In some embodiments, protein homology is measured based on exact identity, e.g., the same amino acid residue at a given position. In some embodiments, 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.).
[0143] In some embodiments, the protein or portion thereof encoded by ETaG (e.g., as 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 embodiments, 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.
[0144] In some embodiments, ETaG is co-regulated with at least one biosynthetic gene in the biosynthetic gene cluster. In some embodiments, ETaG is co-regulated with two or more genes in the biosynthetic gene cluster. In some embodiments, ETaG is co-regulated with the biosynthetic gene cluster in that the expression of ETaG increases or is turned on when a biosynthetic product (biosynthetic product of the biosynthetic gene cluster) produced by an enzyme encoded by the biosynthetic gene cluster is produced. In some embodiments, ETaG is co-regulated with the biosynthetic gene cluster in that the expression of ETaG increases or is turned on when the level of the biosynthetic product of the biosynthetic gene cluster increases.
[0145] In some embodiments, the organism comprising ETaG comprises one or more homologous genes of ETaG. In some embodiments, 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 embodiments, 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 embodiments, the homology is more than 10%. In some embodiments, the homology is more than 20%. In some embodiments, the homology is more than 30%. In some embodiments, the homology is more than 40%. In some embodiments, the homology is greater than 50%. In some embodiments, the homology is greater than 60%. In some embodiments, the homology is greater than 70%. In some embodiments, the homology is greater than 80%. In some embodiments, the homology is greater than 90%.
[0146] In some embodiments, 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 comprise homologous biosynthetic gene clusters. In some embodiments, 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 within a proximity zone to biosynthetic genes of a homologous biosynthetic gene cluster from a different fungal strain. In some embodiments, 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 within a proximity zone to biosynthetic genes of a homologous biosynthetic gene cluster from a different fungal strain. In some embodiments, 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 embodiments, the ETaG gene sequence is optionally no more 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 biosynthetic genes of a homologous biosynthetic gene cluster from a different fungal strain. In some embodiments, it is no more than about 10%) identical. In some embodiments, it is no more than about 20% identical. In some embodiments, it is no more than about 30% identical. In some embodiments, it is no more than about 40%) identical. In some embodiments, it is no more than about 50% identical. In some embodiments, it is no more than about 60% identical. In some embodiments, it is no more than about 70%) identical. In some embodiments, it is no more than about 80% identical. In some embodiments, it is no more than about 90% identical.
[0147] In some embodiments, the human target gene and / or its product is sensitive to regulation by the biosynthetic product of the biosynthetic gene cluster or its analog, the human target gene having its homologous ETaG embedded in the biosynthetic gene cluster or located in a designated proximal zone relative to the biosynthetic gene of the cluster. In some embodiments, the protein encoded by the human target gene is sensitive to regulation by the biosynthetic product of the biosynthetic gene cluster or its analog, the human target gene having its homologous ETaG embedded in the biosynthetic gene cluster or located in a designated proximal zone relative to the biosynthetic gene of the cluster. Thus, in some embodiments, the present disclosure not only provides novel human targets, but also methods and agents for regulating such human targets. In some embodiments, 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.
[0148] In some embodiments, the present disclosure provides a method for evaluating compounds using identified ETaG and the product encoded thereby.In some embodiments, 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.
[0149] In some embodiments, 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 the enzyme encoded by the biosynthetic gene cluster on the target.
[0150] Further analysis may include evaluating 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 embodiments, the aligned sequences were compared to PDB crystal structures. In some embodiments, 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 architecture) 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 yielded 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.
[0151] 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 embodiments, 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.
[0152] Computer Systems In some embodiments, the computer-based method, sequence, genome and / or database provided is embodied in a computer-readable medium. In some embodiments, 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.
[0153] In some embodiments, the present disclosure provides a computer system capable of executing the provided methods described herein. In some embodiments, the present disclosure provides a computer system adapted to execute the provided methods. In some embodiments, the present disclosure provides a computer system adapted to query the provided genome and / or database. In some embodiments, the present disclosure provides a computer system adapted to access the provided database.
[0154] Computer systems that can be used to implement all or part of the provided technology can include digital computers of various forms. 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.
[0155] 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.
[0156] 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.
[0157] 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).
[0158] 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.
[0159] 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.
[0160] 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.
[0161] Exemplary embodiments Among the embodiments provided are the following: 1.a) identifying a putative gene cluster comprising a putative embedded gene in a query genome from a plurality of genomes, the putative gene cluster comprising an anchor gene known to be associated with the gene cluster, the anchor gene co-localizing with the putative embedded gene; b) identifying a plurality of positive genomes that contain an ortholog of the anchor gene and a plurality of negative genomes that do not contain an ortholog of the anchor 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 or phylogenetic distance to the plurality of positive genomes; c) creating a grid representation including 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 anchor genes 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: (1) the presence or absence of orthologs of each protein-coding gene in each genome; (2) sequence similarity of orthologs to their respective protein-coding genes; (3) whether orthologs of each protein-coding gene colocalize with orthologs of anchor genes in each genome; Creating a grid representation based on 4. A computer-implemented method comprising: 2. A computer-implemented method for determining the likelihood that a putative embedded gene is associated with a gene cluster, the putative embedded gene co-localizing with an anchor gene known to be associated with the gene cluster in a query genome; 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 different genomes, the plurality of genomes including a plurality of positive genomes each having an ortholog of an anchor gene and a plurality of negative genomes not having an ortholog of the anchor gene, the second axis corresponding to a plurality of query gene orthologs co-localized with an anchor gene of a BGC in a query genome, the putative embedded gene being one of the plurality of query genes, each cell including: (i) the presence or absence of an ortholog of each query gene in each genome; (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 or a subsection thereof into a machine learning model, the machine learning model being trained to determine a likelihood that a putative embedded gene is embedded in a gene cluster based on values of a plurality of cells in the grid representation, thereby providing a likelihood that the putative embedded gene is associated with the gene cluster; A method comprising: 3. The computer-implemented method of embodiment 2, further comprising generating a grid representation. 4. Generating a grid representation a) identifying a putative gene cluster comprising a putative embedded gene in a query genome from a plurality of genomes, the putative gene cluster comprising an anchor gene known to be associated with the gene cluster, the anchor gene co-localizing with the putative embedded gene; b) identifying a plurality of positive genomes that contain an ortholog of the anchor gene and a plurality of negative genomes that do not contain an ortholog of the anchor 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 or phylogenetic distance 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 anchor genes in a putative gene cluster in the query genome, the second axis corresponding to a plurality of positive genomes and a plurality of negative genomes, each cell comprising: (1) the presence or absence of orthologs of each protein-coding gene in each genome; (2) sequence similarity of orthologs to their respective protein-coding genes; (3) whether orthologs of each protein-coding gene colocalize with orthologs of anchor genes in each genome; Creating a grid representation based on 4. The computer-implemented method of embodiment 3, comprising: 5. The computer-implemented method of any one of embodiments 2-4, wherein the machine learning model is a classification model configured to output a probability for each of a plurality of predefined likelihood categories. 6. The computer-implemented method of embodiment 5, wherein the classification model is a long short-term memory (LSTM) model or a convolutional neural network (CNN) model. 7. The computer-implemented method of any one of embodiments 5 or 6, wherein the plurality of predefined likelihood categories include: (1) high likelihood, (2) somewhat high likelihood, (3) somewhat low likelihood, and (4) low likelihood. 8. The computer-implemented method of any one of embodiments 1-7, wherein the grid representation is a heatmap representation. 9. The computer-implemented method of any one of embodiments 2-8, further comprising displaying the grid representation and the likelihood. 10. The computer-implemented method of any one of embodiments 1-9, wherein the grid representation is hierarchically clustered. 11. The computer-implemented method of any one of embodiments 1 to 10, wherein the number of positive genomes is equal to the number of negative genomes. 12. The computer-implemented method of embodiment 10 or 11, wherein a plurality of positive genomes are selected from a plurality of genome clusters based on sequence similarity of the genomes in the database, and no two positive genomes in the grid representation belong to the same genome cluster. 13. The computer-implemented method of embodiment 12, wherein each negative genome is selected by identifying a genome in the database that has the highest sequence similarity or shortest phylogenetic distance to the positive genome, but does not have an ortholog of the anchor gene. 14. The computer-implemented method of any one of embodiments 1 to 13, wherein 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. 15. The computer-implemented method of any one of embodiments 1-14, wherein the first axis corresponds to at least 20 genomes. 16. The computer-implemented method of embodiment 15, wherein the first axis corresponds to approximately 50 genomes. 17. The computer-implemented method of any one of embodiments 1-16, wherein the plurality of genomes are fungal genomes. 18. The computer-implemented method of any one of embodiments 1-16, wherein the plurality of genomes are plant genomes. 19. The computer-implemented method of any one of embodiments 1-16, wherein the plurality of genomes are bacterial genomes. 20. The computer-implemented method of any one of embodiments 1 to 19, wherein whether a gene co-localizes with an anchor gene of a gene cluster is determined using antiSMASH. 21. The computer-implemented method of any one of embodiments 1 to 20, wherein whether a gene co-localizes with an anchor gene of a gene cluster is determined based on whether the gene is located within a proximity zone upstream or downstream of the anchor gene. 22. The computer-implemented method of embodiment 21, wherein the proximity zone is 50 kb or less. 23. The computer-implemented method of embodiment 22, wherein the proximity zone is approximately 20 kb. 24. The computer-implemented method of any one of embodiments 1-23, wherein the gene cluster is a biosynthetic gene cluster (BGC). 25. A method for identifying resistance genes to secondary metabolites produced by a BGC in a query genome, comprising: (a) Identifying putative embedded genes that are not involved in the production of secondary metabolites by BGCs that colocalize with anchor genes in BGCs in the query genome; and (b) performing the method of any one of embodiments 2 to 24 to determine the likelihood that a putative embedded gene is associated with a 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; and 4. A computer-implemented method comprising: 26. A computer-implemented method for identifying a small molecule modulator of a target gene, comprising: (a) Identifying homologous genes of target genes in fungal genomes that are co-localized with anchor genes of BGCs in fungal genomes and are not involved in the production of secondary metabolites by BGCs; (b) performing the method of any one of embodiments 2 to 24 to determine the likelihood that a homologous gene is associated with a BGC; and (c) identifying a secondary metabolite, or an analog thereof, as a small molecule modulator of the target gene based at least in part on the likelihood that the homologous gene is associated with the BGC; and 4. A computer-implemented method comprising: 27. The computer-implemented method of embodiment 26, wherein the homologous gene encodes a protein having at least about 30% sequence identity to the protein encoded by the target gene. 28. The computer-implemented method of embodiment 26 or 27, further comprising contacting the secondary metabolite or an analog thereof with the protein encoded by the target gene, and detecting the activity of the protein encoded by the target gene. 29. The computer-implemented method of any one of embodiments 26-28, wherein the target gene is a mammalian gene. 30. The computer-implemented method of embodiment 29, wherein the mammalian gene is a human gene. 31. The computer-implemented method of any one of embodiments 26-28, wherein the target gene is a reptilian gene, an avian gene, or an amphibian gene. 32. The computer-implemented method of any one of embodiments 26-28, wherein the target gene is a bacterial gene. 33. The computer-implemented method of any one of embodiments 26-28, wherein the target gene is a fungal gene. 34. The computer-implemented method of any one of embodiments 26-28, wherein the target gene is a plant gene. 35. A computer-implemented method for identifying a plurality of genes associated with BGC, comprising: (a) identifying multiple query genes that co-localize with anchor genes of a BGC in a query genome; (b) for each of a plurality of query genes, determining a likelihood that each query gene is associated with a BGC using a method according to any one of embodiments 2 to 24; (c) identifying a plurality of query genes having a high likelihood of being associated with the BGC, the high likelihood being higher than a threshold, as the plurality of genes associated with the BGC; 4. A computer-implemented method comprising: 36. The computer-implemented method of any one of embodiments 2 to 35, wherein the query gene is a gene homologous to a mammalian protein. 37. The computer-implemented method of embodiment 36, wherein the mammalian protein is a human protein. 38. The computer-implemented method of any one of embodiments 2-35, wherein the query gene is a gene homologous to a reptilian protein, an avian protein, or an amphibian protein. 39. The computer-implemented method of any one of embodiments 2 to 35, wherein the query gene is a gene corresponding to a bacterial protein. 40. The computer-implemented method of any one of embodiments 2 to 35, wherein the query gene is a gene corresponding to a fungal protein. 41. The computer-implemented method of any one of embodiments 2 to 35, wherein the query gene is a gene corresponding to a plant protein. 42. The computer-implemented method of any one of embodiments 1 to 41, wherein the anchor gene is a core synthase gene of the BGC. 43. 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 42. 44. 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 42.
[0162] 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.
[0163] 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 gene is associated with a gene cluster, wherein the putative embedded gene co-localizes with an anchor gene known to be associated with the gene cluster in a query genome; a) receiving a grid representation including a plurality of cells arranged according to a first axis and a second axis, wherein the first axis corresponds to a plurality of different genomes, the plurality of genomes including a plurality of positive genomes each having an ortholog of an anchor gene and a plurality of negative genomes not having an ortholog of the anchor gene, the second axis corresponds to a plurality of query gene orthologs co-localized with the anchor gene of a BGC in the query genome, the putative embedded gene is one of the plurality of query genes, and each cell (i) the presence or absence of an ortholog of each of the query genes in each of the genomes; and (ii) the sequence similarity of the orthologs to the respective query genes; (iii) whether the orthologs of the respective query genes are co-localized with the orthologs of the anchor genes in the respective genomes; receiving, based on b) inputting the grid representation or a subsection thereof into a machine learning model, wherein the machine learning model is trained to determine a likelihood that the putative embedded gene is embedded in the gene cluster based on the values of the cells in the grid representation, thereby providing the likelihood that the putative embedded gene is associated with the gene cluster; A method comprising:
2. The computer-implemented method of claim 1 , further comprising generating the grid representation.
3. generating the grid representation, a) identifying a putative gene cluster comprising a putative embedded gene in a query genome from a plurality of genomes, wherein the putative gene cluster comprises an anchor gene known to be associated with the gene cluster, and the anchor gene co-localizes with the putative embedded gene; b) identifying a plurality of positive genomes that contain an ortholog of the anchor gene and a plurality of negative genomes that do not contain an ortholog of the anchor 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 or phylogenetic distance 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, wherein the first axis corresponds to all protein-coding genes that co-localize with the anchor genes in the putative gene cluster in the query genome, and the second axis corresponds to the plurality of positive genomes and the plurality of negative genomes, and each cell comprises: (1) the presence or absence of orthologs of each protein-coding gene in each genome; (2) the sequence similarity of the ortholog to the respective protein-coding gene; (3) whether the orthologs of the respective protein-coding genes co-localize with the orthologs of the anchor genes in the respective genomes; and creating a grid representation based on The computer-implemented method of claim 2 , comprising:
4. The computer-implemented method of claim 1 , wherein the machine learning model is a classification model configured to output a probability for each of a plurality of predefined likelihood categories.
5. 5. The computer-implemented method of claim 4, wherein the classification model is a long short-term memory (LSTM) model or a convolutional neural network (CNN) model.
6. 5. The computer-implemented method of claim 4, wherein the plurality of predefined likelihood categories include: (1) high likelihood, (2) somewhat high likelihood, (3) somewhat low likelihood, and (4) low likelihood.
7. The computer-implemented method of claim 1 , wherein the grid representation is a heatmap representation.
8. The computer-implemented method of claim 1 , further comprising displaying the grid representation and the likelihood.
9. The computer-implemented method of claim 1 , wherein the grid representation is hierarchically clustered.
10. The computer-implemented method of claim 1 , wherein the number of positive genomes is equal to the number of negative genomes.
11. 10. The computer-implemented method of claim 9, wherein the plurality of positive genomes are selected from a plurality of genome clusters based on sequence similarity of genomes in a database, and no two positive genomes in the grid representation belong to the same genome cluster.
12. 12. The computer-implemented method of claim 11, wherein each negative genome is selected by identifying a genome in the database that has the highest sequence similarity or shortest phylogenetic distance to a positive genome but does not have an ortholog of the anchor gene.
13. The computer-implemented method of claim 3, 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.
14. The computer-implemented method of claim 1 , wherein the first axis corresponds to at least 20 genomes.
15. 15. The computer-implemented method of claim 14, wherein the first axis corresponds to approximately 50 genomes.
16. The computer-implemented method of claim 1 , wherein the plurality of genomes are fungal genomes.
17. The computer-implemented method of claim 1 , wherein the plurality of genomes are plant genomes.
18. The computer-implemented method of claim 1 , wherein the plurality of genomes are bacterial genomes.
19. 2. The computer-implemented method of claim 1, wherein whether a gene co-localizes with an anchor gene of a gene cluster is determined using antiSMASH.
20. The computer-implemented method of claim 1, wherein whether the putative embedded gene co-localizes with an anchor gene of a gene cluster is determined based on whether the gene is located within a proximity zone upstream or downstream of the anchor gene.
21. 21. The computer-implemented method of claim 20, wherein the proximity zone is 50 kb or less.
22. 22. The computer-implemented method of claim 21, wherein the proximity zone is approximately 20 kb.
23. The computer-implemented method of claim 1 , wherein the gene cluster is a biosynthetic gene cluster (BGC).
24. 1. A computer-implemented method for identifying resistance genes to secondary metabolites produced by BGCs in a query genome, comprising: (a) identifying putative embedded genes that co-localize with anchor genes in the BGC in the query genome, wherein the putative embedded genes are not involved in the production of the secondary metabolite by the BGC; (b) performing the method of claim 1 to determine the likelihood that the putative embedded gene is associated with the BGC; (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; and 11. A computer-implemented method comprising:
25. 1. A computer-implemented method for identifying small molecule modulators of a target gene, comprising: (a) identifying a homologous gene of the target gene in the fungal genome that is co-localized with an anchor gene of a BGC in the fungal genome and is not involved in the production of secondary metabolites by the BGC; (b) performing the method of claim 1 to determine the likelihood that the homologous gene is associated with the BGC; (c) identifying the secondary metabolite or analog thereof as a small molecule modulator of the target gene based at least in part on the likelihood that the homologous gene is associated with the BGC; and 11. A computer-implemented method comprising:
26. 26. The computer-implemented method of claim 25, wherein the homologous gene encodes a protein having at least about 30% sequence identity to the protein encoded by the target gene.
27. 26. The computer-implemented method of claim 25, further comprising contacting the secondary metabolite or analog thereof with a protein encoded by the target gene and detecting the activity of the protein encoded by the target gene.
28. 26. The computer-implemented method of claim 25, wherein the target gene is a mammalian gene.
29. 29. The computer-implemented method of claim 28, wherein the mammalian gene is a human gene.
30. 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.
31. 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.