Computational methods for identifying gene networks containing functionally related genes

Abstract

We describe a computational method for identifying biosynthetic gene clusters (BGCs) independent of prior identification of core synthases. The disclosed method also facilitates linking BGCs to their potential downstream targets. The approach integrates information from co-evolution, co-occurrence, co-regulation, co-localization, and functional enrichment to group functionally related genes across gene networks, delineate BGCs, and propose putative secondary metabolite targets.

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of priority to U.S. Provisional Patent Application No. 63 / 297,565, filed January 7, 2022, the contents of which are incorporated herein by reference in their entirety.

[0002] The present disclosure relates generally to methods and systems for identifying genes associated with gene networks, such as biosynthetic gene clusters, and uses thereof, including identifying potential therapeutic targets and drug candidates. [Background technology]

[0003] Secondary metabolites identified from bacteria, fungi and plants have found use in a wide variety of applications in medicine and agriculture (see, for example, Katz, et al. (2016), "Natural Product Discovery: Past, Present, and Future", J Ind Microbiol Biotechnol. 43(2-3):155-76). These secondary metabolites are often synthesized via metabolic pathways consisting of core biosynthetic enzymes (such as polyketide synthases and nonribosomal peptide synthases) and various regulatory enzymes, which are co-localized in genomes as biosynthetic gene clusters (BGCs) (see, for example, Scherrach, et al. (2021), "Mining and Unearthing Hidden Biosynthetic Potential", Nat Commun. 12(1):3864). With the availability of a vast number of sequenced genomes, genomics database mining techniques have become crucial for the discovery of networks of functionally related genes, including BGCs. Traditional approaches to identifying BGCs have relied on the identification of core synthases using profile hidden Markov models and the subsequent inclusion of genes proximal to these core enzymes in the genome (see, e.g., Blin, et al. (2021), "antiSMASH 6.0: Improving Cluster Detection and Comparative Capabilities," Nucleic Acids Res. 49(W1):W29-W35; Scherrach, et al. (2021), ibid.). However, such approaches have several drawbacks.First, not all BGCs contain core enzymes, and many secondary metabolites can be synthesized from scaffolds derived from central metabolism (see, e.g., Wasil, et al. (2018), "Oryzines A&B, Maleidride Congeners from Aspergillus oryzae and Their Putative Biosynthesis", J. Fungi 4:96; Lim, et al. (2018), "Fungal Isocyanide Synthases and Xanthocillin Biosynthesis in Aspergillus fumigatus", mBio 9(3):e00785-18). Currently, approximately 11% of experimentally validated BGCs in the MiBiG database (see, for example, Kautsar, et al. (2020), "MIBiG 2.0: A Repository for Biosynthetic Gene Clusters of Known Function", Nucleic Acids Res. 48(D1):D454-D458) are classified as "other", i.e., BGCs without an associated core synthase. Second, these approaches rely on the colocalization of genes for core synthases and associated regulatory enzymes, which may limit the scope of metabolic interactions that may be involved in the synthesis of secondary metabolites. Finally, a major challenge in genome mining-based BGC identification is to link the identified BGCs with potential downstream targets of the synthesized secondary metabolites. One approach to making this association is the embedded resistance gene hypothesis (Yan, et al. (2020), "Recent Developments in Self-Resistance Gene Directed Natural Product Discovery", Nat Prod Rep. 37(7):879-892), in which genes encoding protein targets of secondary metabolites undergo duplication events, and paralogous copies become embedded within BGCs and undergo accelerated evolution, thereby resulting in the acquisition of mutations that abolish the effects of the secondary metabolites without affecting the primary protein function.This hypothesis was validated by genomic analysis of BGCs associated with previously discovered secondary metabolites and has also proven effective in identifying novel targets (Yan et al. (2020) ibid.), however, identification of resistance genes has been limited to those embedded within the BGCs. Summary of the Invention

[0004] To avoid the shortcomings outlined above, a computational approach for identifying networks of functionally related genes is described. In some examples, the computational approach may be used to identify BGCs, for example, independent of their association with core synthases, and may also facilitate linking BGCs to their potential downstream targets. The approach integrates information from coevolution, co-occurrence, co-regulation, co-localization, and functional enrichment to (i) group functionally related genes across gene networks, (ii) delineate specific types of gene networks (e.g., BGCs), and (iii) suggest related secondary metabolite targets if the method is used to identify BGCs. The disclosed method (and systems designed to implement the disclosed method) is based on the observation that functionally related genes coevolve (Steenwyk, et al. (2022), "An Orthologous Gene Coevolution Network Provides Insight Into Eukaryotic Cellular and Genomic Structure and Function", Sci.Adv.8, eabn0105) and are often coregulated. Furthermore, as previously discussed, genes within BGCs frequently co-localize and co-occur in related organisms. The disclosed methods provide a pipeline for the discovery and functional assignment of gene networks (e.g., BGCs).

[0005] Disclosed herein is a computer-implemented method for identifying networks of functionally related genes, the method comprising: receiving as input a selection of genomes for analysis, the selection of genomes comprising a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining a pairwise coevolution metric, a pairwise co-regulation metric, a pairwise co-occurrence metric, a pairwise co-localization metric, or any combination thereof, for the identified COGs; determining pairwise functional association scores for the identified COGs based on the determined pairwise coevolution metric, pairwise co-regulation metric, pairwise co-occurrence metric, pairwise co-localization metric, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group them into functionally related COGs; and outputting a determination that a COG cluster is a network of functionally related genes in a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes in the particular functional category.

[0006] In some embodiments, functional enrichment analysis does not require the identification of genes known to be associated with a particular functional category of gene network.

[0007] In some embodiments, the functional enrichment analysis includes testing for enrichment of genes within functional categories known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. In some embodiments, the functional categories known to be associated with BGCs include gene ontology terms or KEGG pathways known to be associated with BGCs. In some embodiments, examples of gene ontology terms known to be associated with BGCs include GO:0019748 (secondary metabolic process), GO:0044550 (secondary metabolite biosynthetic process), GO:0030639 (polyketide biosynthetic process), GO:0030638 (polyketide metabolic process), GO:0043455 (regulation of secondary metabolic process), GO:1900539 (fumonisin metabolic process), or any combination thereof. In some embodiments, KEGG pathways known to be associated with BGCs include M00778 (type II polyketide backbone biosynthesis) or M00095 (C5 isoprenoid biosynthesis, mevalonate pathway), M00937 (aflatoxin biosynthesis), M00893 (lovastatin biosynthesis), or any combination thereof.

[0008] In some embodiments, the functional enrichment analysis includes testing for enrichment of protein domain representations known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. In some embodiments, the protein domain representations known to be associated with BGCs include PFAM domain representations, Conserved Domain Database (CDD) domain representations, or TIGRFAM domain representations known to be associated with BGCs.

[0009] In some embodiments, the computer-implemented method further includes identifying a putative target of a secondary metabolite synthesized by the putative BGC by identifying a protein sequence that is not a component of a known BGC, determining a pairwise functional association score for the identified protein sequence and the putative BGC, and identifying the putative target of the secondary metabolite based on a comparison of the pairwise functional association score to a first predetermined threshold. In some embodiments, the pairwise functional association score includes a co-regulation score, and the first predetermined threshold of the co-regulation score corresponds to a p-value of 0.05 or less. In some embodiments, the pairwise functional association score includes a co-evolution score, and the first predetermined threshold of the co-evolution score corresponds to a co-evolution score value of 0.7 or more. In some embodiments, the pairwise functional association score includes a co-occurrence score, and the first predetermined threshold of the co-occurrence score corresponds to a co-occurrence score value of 0.5 or more.

[0010] In some embodiments, identification of a putative BGC does not require identification of the associated core synthase.

[0011] In some embodiments, the plurality of related genomes comprises fungal, bacterial, or plant genomes.

[0012] In some embodiments, identifying the COGs comprises using BLAST to identify orthologous genes in the multiple related genomes as bidirectional best hits, followed by clustering the identified orthologous genes. In some embodiments, identifying the COGs comprises using orthoMCL or orthoFinder to identify orthologous genes in the multiple related genomes.

[0013] In some embodiments, determining the pairwise coevolution metric for the COGs comprises: calculating the percentage identity between each pair of protein sequences in each COG of the pair of COGs to identify shared protein sequences; calculating the Pearson correlation coefficient for each pair of COGs that contains a certain minimum number of shared protein sequences to estimate the coevolution rate; filtering the COGs by removing COGs whose pairwise Pearson correlation coefficient is less than a second predetermined threshold and clustering the remaining COGs according to the estimated coevolution rate; and performing a functional enrichment analysis to remove clusters of COGs enriched for essential metabolic functional categories. In some embodiments, the second predetermined threshold corresponds to a Pearson correlation coefficient value of 0.7, 0.8, 0.9, 0.95, 0.98, or 0.99. In some embodiments, clustering the remaining COGs according to the estimated coevolution rate comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

[0014] In some embodiments, determining the pairwise co-regulation metric for the COGs comprises extracting intergenic regions within each COG, performing de novo detection of sequence motifs within the extracted intergenic regions to identify putative cis-regulatory elements or transcription factor binding sites (TFBS), comparing the putative cis-regulatory elements or TFBS identified for each COG with those identified across all other COGs to determine pairwise motif similarity scores between the COGs, filtering the COGs to exclude COGs whose pairwise motif similarity scores have a p-value equal to or less than a third predetermined threshold, and clustering the filtered COGs based on the pairwise motif similarity scores to identify co-regulated COG clusters. In some embodiments, the third predetermined threshold corresponds to a p-value of 0.05. In some embodiments, the third predetermined threshold corresponds to a p-value of 0.01. In some embodiments, clustering the remaining COGs according to the motif similarity scores comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

[0015] In some embodiments, determining the pairwise co-occurrence metric of the COGs comprises calculating a Jaccard coefficient for each pair of COGs (COG A and COG B) based on the following relationship:

number

[0016] In some embodiments, determining the pairwise colocalization metric of the COGs comprises calculating a proximity score for each pair of corresponding gene sequences in the pair of COGs (COG A and COG B) based on the following relationship:

number

[0017] In some embodiments, the pairwise functional association scores of the identified COGs are based on the addition of determined pairwise coevolution metrics, pairwise co-regulation metrics, pairwise co-occurrence metrics, pairwise co-localization metrics, or any combination thereof.

[0018] In some embodiments, clustering the identified COGs according to their pairwise functional association scores comprises the use of Markov clustering (MCL) or hierarchical clustering algorithms.

[0019] In some embodiments, the computer-implemented method further comprises determining a horizontal introgression metric based on a calculation of a codon adaptation index (CAI) or a dinucleotide signature divergence index (DSDI). In some embodiments, the horizontal introgression metric is used to further refine the clustering of co-localizing, co-occurring and / or co-evolving COGs. In some embodiments, the horizontal introgression metric is used as part of a post-processing step to search for nearby horizontally transferred genes that were missed in the upstream clustering step.

[0020] In some embodiments, the computer-implemented method further comprises evaluating the gene identified as belonging to the putative BGC to determine whether it is a resistance gene. In some embodiments, the resistance gene is an embedded target gene (ETaG) or a non-embedded target gene (NETaG). In some embodiments, the computer-implemented method further comprises performing an in vitro assay to test a secondary metabolite produced by the putative BGC for activity against a resistance gene homolog identified in the target genome or a protein encoded thereby. In some embodiments, the computer-implemented method further comprises performing an in vivo assay to test a secondary metabolite produced by the putative BGC for activity against a resistance gene homolog identified in the target genome or a protein encoded thereby. In some embodiments, the target genome comprises a mammalian genome, a human genome, an avian genome, a reptile genome, an amphibian genome, a plant genome, a fungal genome, a bacterial genome, or a viral genome.

[0021] The invention also includes a system that includes one or more processors, and that is communicatively coupled to the one or more processors, which, when executed by the one or more processors, causes the system to receive as input a selection of genomes for analysis, where the selection of genomes includes a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining a pairwise coevolution metric, a pairwise coregulation metric, a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof, for the identified COGs; and determining a pairwise coevolution metric, a pairwise coregulation metric, a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof, for the determined COGs. and a memory configured to store instructions that cause: determining pairwise functional association scores for the identified COGs based on a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group into functionally related COGs; and outputting a determination that the COG cluster is a network of functionally related genes in a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes in the particular functional category.

[0022] In some embodiments, the functional enrichment analysis does not require the identification of genes known to be associated with a particular functional category of gene network. In some embodiments, the functional enrichment analysis includes testing for enrichment of genes in functional categories known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. In some embodiments, the functional categories known to be associated with BGCs include gene ontology terms or KEGG pathways known to be associated with BGCs. In some embodiments, the functional enrichment analysis includes testing for enrichment of protein domain representations known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. In some embodiments, the protein domain representations known to be associated with BGCs include PFAM domain representations, Conserved Domain Database (CDD) domain representations, or TIGRFAM domain representations known to be associated with BGCs.

[0023] In some embodiments, the system further comprises instructions for identifying a putative target of a secondary metabolite synthesized by the putative BGC by identifying a protein sequence that is not a component of a known BGC, determining a pairwise functional association score for the identified protein sequence and the putative BGC, and identifying the putative target of the secondary metabolite based on a comparison of the pairwise functional association score to a first predetermined threshold. In some embodiments, the pairwise functional association score comprises a co-regulation score, and the first predetermined threshold of the co-regulation score corresponds to a p-value of 0.05 or less. In some embodiments, the pairwise functional association score comprises a co-evolution score, and the first predetermined threshold of the co-evolution score corresponds to a co-evolution score value of 0.7 or more. In some embodiments, the pairwise functional association score comprises a co-occurrence score, and the first predetermined threshold of the co-occurrence score corresponds to a co-occurrence score value of 0.5 or more.

[0024] In some embodiments, identification of a putative BGC does not require identification of the associated core synthase.

[0025] A non-transitory computer-readable medium storing one or more programs, the one or more programs which, when executed by one or more processors of the system, provide the system with: A non-transitory computer readable medium is disclosed, the medium comprising instructions to cause: receiving as an input a selection of genomes for analysis, the selection of genomes including a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining a pairwise coevolution metric, a pairwise co-regulation metric, a pairwise co-occurrence metric, a pairwise co-localization metric, or any combination thereof, for the identified COGs; determining pairwise functional association scores for the identified COGs based on the determined pairwise coevolution metric, pairwise co-regulation metric, pairwise co-occurrence metric, pairwise co-localization metric, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group into functionally related COGs; and outputting a determination that a COG cluster is a network of functionally related genes in a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes in a particular functional category.

[0026] In some embodiments, the functional enrichment analysis does not require the identification of genes known to be associated with a particular functional category of gene network. In some embodiments, the functional enrichment analysis includes testing for enrichment of genes in functional categories known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. In some embodiments, the functional categories known to be associated with BGCs include gene ontology terms or KEGG pathways known to be associated with BGCs. In some embodiments, the functional enrichment analysis includes testing for enrichment of protein domain representations known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. In some embodiments, the protein domain representations known to be associated with BGCs include PFAM domain representations, Conserved Domain Database (CDD) domain representations, or TIGRFAM domain representations known to be associated with BGCs.

[0027] In some embodiments, the non-transitory computer readable medium further comprises instructions for identifying a putative target of a secondary metabolite synthesized by the putative BGC by identifying a protein sequence that is not a component of a known BGC, determining a pairwise functional association score for the identified protein sequence and the putative BGC, and identifying the putative target of the secondary metabolite based on a comparison of the pairwise functional association score to a first predetermined threshold. In some embodiments, the pairwise functional association score comprises a co-regulation score, and the first predetermined threshold of the co-regulation score corresponds to a p-value of 0.05 or less. In some embodiments, the pairwise functional association score comprises a co-evolution score, and the first predetermined threshold of the co-evolution score corresponds to a co-evolution score value of 0.7 or more. In some embodiments, the pairwise functional association score comprises a co-occurrence score, and the first predetermined threshold of the co-occurrence score corresponds to a co-occurrence score value of 0.5 or more.

[0028] In some embodiments, identification of a putative BGC does not require identification of the associated core synthase.

[0029] 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.

[0030] 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.

[0031] 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]

[0032] [Figure 1] A non-limiting example of a process flow chart for identifying networks of functionally related genes according to one or more examples of the present disclosure is presented.

[0033] [Diagram 2] 1 presents a non-limiting schematic diagram of a computing device according to one or more examples of the present disclosure.

[0034] [Diagram 3]1 presents a non-limiting schematic of a pipeline for BGC discovery and function assignment.

[0035] [Figure 4] Non-limiting examples of BGC identification pipeline performance metrics for known and predicted BGCs are presented. The figure shows recall (proportion of genes known to be involved in the biosynthesis of the identified target molecule), precision (proportion of predicted genes that are true BGC genes for the target molecule) and recall (core synthase clusters) (proportion of core synthase clusters predicted by antiSMASH for identification).

[0036] [Diagram 5] We present non-limiting examples of data on the contribution of different components of the BGC identification pipeline to prediction performance. This figure shows the performance of the BGC identification pipeline when colocalization (Coloc only), colocalization and coregulation (Coloc+Coreg), colocalization and coevolution (Coloc+Coevo), or colocalization, coregulation and coevolution are all applied to predict BGCs (Coloc+Coreg+Coevo). Recall and precision are described elsewhere herein. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0037] A computational approach for identifying networks of functionally related genes is described. In some examples, the computational approach may be used to identify BGCs, for example, independent of their association with core synthases, and may also facilitate linking BGCs to their potential downstream targets. The approach integrates information from coevolution, co-occurrence, co-regulation, co-localization, and functional enrichment to (i) group functionally related genes across gene networks, (ii) delineate specific types of gene networks (e.g., BGCs), and (iii) suggest associated secondary metabolite targets if the method is used to identify BGCs. The disclosed method (and systems designed to implement the disclosed method) is based on the observation that functionally related genes coevolve (Steenwyk et al. (2022), supra) and are often co-regulated. Furthermore, as previously discussed, genes within gene networks such as BGCs frequently co-localize and co-occur in related organisms. Details regarding the individual components of the approach and how they can be combined into a pipeline, for example, for BGC discovery and functional assignment, are described below.

[0038] In some examples, a disclosed method (e.g., a computer-implemented method) for identifying a network of functionally related genes may include receiving as input a selection of genomes for analysis, where the selection of genomes includes a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining a pairwise coevolution metric, a pairwise co-regulation metric, a pairwise co-occurrence metric, a pairwise co-localization metric, or any combination thereof, for the identified COGs; determining pairwise functional association scores for the identified COGs based on the determined pairwise coevolution metric, pairwise co-regulation metric, pairwise co-occurrence metric, pairwise co-localization metric, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group them into functionally related COGs; and outputting a determination that the COG cluster is a network of functionally related genes in a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes in the particular functional category.

[0039] In some examples, the disclosed method for identifying putative biosynthetic gene clusters (BGCs) may include receiving as input a selection of genomes for analysis, where the selection of genomes includes a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining a pairwise coevolution metric, a pairwise coregulation metric, a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof, for the identified COGs; determining a pairwise functional association score for the identified COGs based on the determined pairwise coevolution metric, a pairwise coregulation metric, a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group into functionally related COGs; and outputting a determination that the COG cluster is a putative BGC based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes within a particular functional category and / or protein domain known to be associated with BGCs.

[0040] In some examples, the disclosed methods further include identifying a putative target of a secondary metabolite synthesized by the putative BGC by identifying a protein sequence that is not a component of a known BGC, determining a pairwise functional association score for the identified protein sequence and the putative BGC, and identifying the putative target of the secondary metabolite based on a comparison of the pairwise functional association score to a first predetermined threshold.

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

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

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

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

[0045] As used herein, "secondary metabolite" refers to a small organic molecule or 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.

[0046] The term "biosynthetic gene cluster" or "BGC" is used interchangeably herein and refers to a locally clustered group of one or more genes that together 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, for example, 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," in: 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, BGCs may also contain other genes, e.g., genes encoding products not involved in the biosynthesis of secondary metabolites, interspersed among the biosynthetic genes. These genes are referred to herein as "associated" with a BGC if their products are functionally related to the secondary metabolites of the BGC. Some genes, e.g., genes not involved in the biosynthesis of secondary metabolites produced by a BGC, are referred to herein as "embedded" in a BGC if their products are functionally related to the secondary metabolites of the BGC and are located in physical proximity to the biosynthetic genes of the cluster.Some genes, e.g., genes not involved in the biosynthesis of secondary metabolites produced by the BGC, are referred to herein as "non-embedded" if their products are functionally related to the secondary metabolites of the BGC but are not located in physical proximity to the biosynthetic genes of the BGC. "Anchor genes" refer to biosynthetic genes (e.g., core synthases) involved in the biosynthesis of secondary metabolites produced by the BGC that are known to be co-localized and functionally related (i.e., associated) with the BGC.

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

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

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

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

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

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

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

[0054] 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.

[0055] Computational methods for identifying gene networks containing functionally related genes Coevolution of functionally related genes: Genes involved in the same biological process or pathway generally face similar evolutionary pressures and are therefore likely to evolve at similar rates. Previous studies analyzing the overall coevolution of genes in budding yeast have shown that functionally related genes exhibit highly correlated evolutionary rates and that this property is not necessarily linked to the physical location of the genes in the genome (see, e.g., Steenwyk et al. (2022), ibid.). Networks derived from clustering of coevolving genes exhibited a high degree of functional consistency, with many essential genes involved in important biological processes grouped together. This association of coevolving genes mirrored observations from gene interaction networks obtained from studying the effects on growth of single and double gene knockout mutants (Costanzo, et al. (2016), "A Global Genetic Interaction Network Maps a Wiring Diagram of Cellular Function", Science 353(6306):aaf1420). Thus, coevolution can be used as an approach to group functionally related genes.

[0056] The disclosed methodology for evaluating co-evolved genes may include the following. i. A set of diverse but related genomes, e.g. from fungal, bacterial or plant species, is selected to facilitate the capture of true evolutionary signals. ii. Orthologs shared between these species are then identified by bidirectional best BLAST hits and / or approaches to group them into clusters of orthologous genes (COGs) using software tools such as orthoMCL and orthoFinder (e.g., Lee, et al. (2003), "OrthoMCL: Identification of Ortholog Groups for Eukaryotic Genomes", Genome Res. 13(9):2178-89; Emms, et al. (2019), "OrthoFinder: Phylogenetic Orthology Inference for Comparative Genomics", Genome Biol. 20(1):238). iii.To estimate the co-evolutionary rates between COGs, the percentage of identity between all pairs of proteins in each COG is calculated. For each COG, the Pearson correlation coefficient is then calculated between it and all other identified COGs with the minimum number of shared pairs (i.e., equivalent genome pairs), generating a correlation matrix. This provides an estimate of how similar the rates of evolution are between each pair of COGs. iv. To generate a network of coevolving COGs, Markov Clustering (MCL) or hierarchical clustering can be used to group COGs that show the highest degree of coevolution. To ensure that only strong coevolutionary links are considered, Pearson's correlation coefficients below a specified threshold (e.g., 0.75) can be filtered out. This approach groups COGs into clusters of coevolving genes. v. A functional enrichment analysis performed on the network (e.g., a statistical analysis to identify functional categories (e.g., based on Gene Ontology terms) of genes that are over-represented in the network of co-evolving COGs, where functions are inferred from the set of annotated genomes selected for analysis) can then reveal a high degree of functional coherence within the clusters, with essential functions grouped into primary clusters and genes encoding auxiliary functions grouped into smaller clusters. vi. In some instances, clusters of COGs enriched for essential functions may be filtered from downstream analyses as these COGs are unlikely to be involved in secondary metabolite biosynthesis.

[0057] Co-regulation of functionally related genes: Functionally related genes are often co-regulated, which can serve as an additional layer of information in the analysis. This can be achieved by identifying signatures of common regulation (i.e., shared putative cis-regulatory elements or transcription factor binding sites (TFBS)) within and between COGs, as these binding sites are often conserved between related species.

[0058] The disclosed methodology for identifying co-regulated genes can include the following. i. To identify co-regulated genes, first extract the intergenic regions of genes (identified as above) within each COG. ii. De novo motif detection can be performed on these intergenic regions using motif detection software such as MEME (see, e.g., Bailey, et al. (2015), "The MEME Suite", Nucleic Acids Res. 43(W1):W39-49) or HOMER (see, e.g., Heinz, et al. (2010), "Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities", Mol Cell. 38(4):576-589). iii. The putative TFBSs identified for each COG by this analysis can then be compared to the TFBSs identified across all other COGs to generate a pairwise similarity matrix between COGs based on the similarity of their identified motifs. iv. To identify networks of co-regulated COGs, these pairwise motif similarity scores can then be used to cluster the COGs (e.g., via MCL or hierarchical clustering) to generate co-regulated COG clusters. To ensure that only highly significant motif similarities are considered when clustering co-regulated COGs, a p-value cutoff (e.g., p-value of 0.01 or greater) can be utilized to filter the COGs by motif similarity scores. v. Functional enrichment analysis performed on gene networks can be used to identify COG clusters in which a high degree of functional coherence exists.

[0059] Co-occurrence of COGs: Because genes of a gene network (e.g., a BGC) need to function together (e.g., to produce a target secondary metabolite), they may be expected to show correlated patterns in their occurrence across genomes, such that all genomes share a core set of genes that can, for example, produce a secondary metabolite.

[0060] The disclosed methodology for assessing co-occurring genes may include the following. i. To utilize co-occurrence in identifying gene networks, a co-occurrence score can be calculated for each pair of COGs by calculating the Jaccard coefficient of the set of genomes in each COG. For example, to generate a co-occurrence score between COG A and COB A, the following Jaccard coefficient is calculated: Co-occurrence score = (|A∩B|) / (|A∪B|) where A∩B is the number of genomes shared between COG A and COG B, and A∪B is the set of all genomes present in A or B. ii. It provides a measure of pairwise co-occurrence between COGs that ranges from 0 (for COGs that occur in completely different genomes) to 1 (for COGs that occur in the exact same genome). Thus, a pairwise co-occurrence matrix can be calculated between all COGs.

[0061] COG colocalization: Networks of functionally related genes (e.g., BGCs) are often composed of genes that are located close to each other in the genome. This information can also be integrated to facilitate grouping of functionally related genes.

[0062] The disclosed methodology for incorporating colocalization information may include the following. To leverage the colocalization information, a proximity score can be generated that captures how closely orthologous genes are located to each other across species. Non-limiting examples of such proximity scores can be given by: Proximity score = 1 / (1 + number of interval genes) In the formula, adjacent genes receive a score of 1, more distant genes receive a score less than 1 (the more distant the gene, the smaller the score), and genes in separate contigs receive a value of 0. ii. Thus, in comparing two COGs, one can calculate the proximity score of each gene in COG A with its corresponding gene from the same genome in COG B, and take the average of all the computable pairs to get the average proximity score of COG A to COG B. Thus, one can calculate a pairwise colocalization matrix between all COGs. iii. Because this measure of colocalization is calculated across species, it also captures synteny between genes. It can therefore be used to cluster COGs, grouping SYNTENOUS COGs together.

[0063] Signatures of horizontal gene transfer: BGCs can be exchanged between microbial species through the process of horizontal gene transfer, thereby driving increased diversity of natural products. Because organisms have intrinsic genomic characteristics, in some instances, these laterally transferred genes may still have the genomic signature of the original donor, especially if they were recently acquired or derived from very distantly related organisms or the process of improvement is slow. Thus, genomic characteristics such as codon usage, GC content and / or dinucleotide ratios may differ significantly between, for example, the horizontally transferred BGC and the rest of the host genome, thereby providing another metric that can be used to delineate gene networks such as BGCs.

[0064] For example, the codon adaptation index (CAI) of a gene sequence can be calculated relative to the rest of the genome (which serves as a reference) by first calculating the relative synonymous codon usage (RSCU) across the genome as previously described (see, e.g., Sharp, et al. (1987), "The Codon Adaptation Index--A Measure of Directional Synonymous Codon Usage Bias, and Its Potential Applications", Nucleic Acids Res. 15(3):1281-95). The CAI is calculated by comparing the codon usage of the gene to the RSCU table. CAI, a value between 0 and 1, serves as a measure of how similar the codon usage of a gene or locus is to the reference set. Genes with higher CIA values ​​have more similar codon usage compared to the rest of the genome. CIA can be calculated for all genes in the genome, and those with lower than expected values ​​can be considered candidates for horizontally transferred genes.

[0065] Co-localized, co-occurring, co-regulated and / or co-evolved gene clusters can then be further refined by CAI-based clustering. Alternatively, CAI can be used as part of post-processing (described below) to search for nearby horizontally transferred genes that were missed in the upstream clustering step.

[0066] Dinucleotide relative abundance (DRA) is the ratio of observed dinucleotide frequency to expected frequency derived from single nucleotide frequencies under independent assumptions and can be calculated using the following formula:

number

[0067] In the formula, p xy is the DRA of dinucleotide xy, and f xy is the dinucleotide frequency, and f x、f y is the single nucleotide frequency. The calculated DRA value for the target genome locus, e.g., the whole BGC or a specific gene, can be compared with the overall / background DRA value from the whole genome to generate a dinucleotide signature difference index (DSDI), e.g., by taking the Euclidean distance between the two vectors. In some examples, the DSDI can be generated by taking the chi-square / DRA divergence, delta distance, or quadratic discriminant between the two vectors (see, e.g., Baran, et al. (2008), "Detecting Horizontally Transferred and Essential Genes Based on Dinucleotide Relative Abundance", DNA Research 15:267-276). The DSDI can be calculated for all genes in the genome, and those with values ​​higher than the expected value can be considered candidates for horizontally transferred genes.

[0068] Co-localizing, co-occurring, co-regulated and / or co-evolving COGs can then be further refined by DSDI-based clustering. Alternatively, DSDI can be used as part of post-processing (described below) to search for nearby horizontally transferred genes that were missed in the upstream clustering step. In some instances, CAI and DSDI can be combined for this analysis.

[0069] Then, clustering of functionally related genes: To cluster functionally related genes, these four calculated measures for pairwise co-evolution (Pearson's correlation coefficient), pairwise co-regulation (motif similarity score), pairwise co-localization (average proximity score) and pairwise co-occurrence (co-occurrence score), or combinations or subsets thereof, can be combined into a unified score of functional association. For example, in some examples, a unified score of functional association can be calculated based on the additive scores of any two or more of the pairwise metrics for co-evolution, co-regulation, co-localization, co-occurrence, and optionally horizontal introgression, where the individual metrics can be weighted equally in some examples, or weighted differently, for example according to the rank order of their predictive values, in some examples. This pairwise functional association score can then be used to cluster the COGs, for example using either MCL or hierarchical clustering, to group functionally related COGs.

[0070] Functional enrichment analysis to identify potential gene networks: Following the method outlined above, which leads to grouping all functionally related COGs into clusters, functional gene networks, functional enrichment analysis can be performed among the identified COG clusters likely to be involved in, for example, the synthesis of secondary metabolites, to identify, for example, biosynthetic gene clusters, to test for enrichment of specific functional categories (e.g., using gene ontology terms and / or KEGG pathways) and / or specific protein domains (e.g., PFAM domains) known to be associated with BGCs, e.g., methyltransferases, monooxygenases, etc. COG clusters enriched for specific functional categories and / or protein domains can then be prepared as putative gene networks, e.g., putative BGCs.

[0071] Post-processing: To ensure that all components of a gene network, such as BGCs, are captured, post-processing steps can be performed to capture proximal genes that may be unique to a particular clade or species. Post-processing steps can include, but are not limited to, converting COG-level information (across species) into genome-level information for each genome selected for analysis. For each genome, each identified cluster of functionally related genes can be enriched with embedded genes that are found between the cluster genes of the particular genome, but not found at the COG level. Similarly, clusters can also be enriched with proximal genes that have similar DNA signatures of horizontal introgression as other genes in the identified cluster of functionally related genes. Conversely, genes found to be very far from the core set of functionally related genes in a particular genome can be considered for exclusion from the cluster.

[0072] Association of downstream targets: To associate potential targets of secondary metabolites with the identified BGCs, proteins that are typically not BGC components but have high functional association scores can be readily identified (with or without considering colocalization) to rank potential targets of secondary metabolites.

[0073] FIG. 1 presents a non-limiting example of a flow chart of a process 100 for identifying networks of functionally related genes. The process 100 can be performed as a computer-implemented method using software executing on one or more processors of one or more electronic devices, computers, or computer platforms, for example. In some examples, the process 100 is performed using a client-server system, and blocks of the process 100 are divided in any manner between a server and a client device. In other examples, blocks of the process 100 are divided between a server and multiple client devices. Thus, although parts of the process 100 are described herein as being performed by a particular device of a client-server system, it will be understood that the process 100 is not so limited. In other examples, the process 100 is performed using only one client device or only multiple client devices. In the process 100, some blocks are optionally combined, the order of some blocks is optionally changed, and some blocks are optionally omitted. In some examples, additional steps can be performed in combination with the process 100. Accordingly, the operations illustrated (and described in more detail below) are exemplary in nature and, therefore, should not be considered as limiting.

[0074] In step 102 of Figure 1, a selection of genomes to analyze is received as input, where the selection of genomes includes multiple related genomes. In some examples, for example, the selection of genomes may be input by a user of a system configured to perform the methods (e.g., computer-implemented methods) described herein.

[0075] In some examples, the multiple related genomes may include genomes of organisms known or suspected to contain a gene network of interest (e.g., a network of functionally related genes or gene products). Examples of gene networks include, but are not limited to, biochemical pathways, cellular pathways, and signal transduction pathways involved in immune response, such as gene regulatory networks, primary and secondary metabolic pathways, hormone signal transduction pathways, and JAK-STAT pathways. In some examples, the multiple related genomes may include, for example, mammalian genomes, human genomes, avian genomes, reptile genomes, amphibian genomes, plant genomes, fungal genomes, bacterial genomes, or viral genomes.

[0076] In some examples, the gene network of interest may include a BGC that produces a secondary metabolite, and the multiple related genomes may include, for example, a fungal genome, a bacterial genome, or a plant genome.

[0077] In some examples, the multiple related genomes may be input in any of a variety of formats or representations known to those of skill in the art, including, but not limited to, nucleotide sequences, amino acid sequences, or CDD, Gene3D, PANTHER, Pfam, ProSitePatterns, ProSiteProfiles, SUPERFAMILY, SMART, TIGRFAM, SFLD, Hamap, Coils, PRINTS, PIRSR, AntiFam, MobiDBLite, or PIRSF representations of protein domains encoded by genes in the genomes of the multiple related genomes, or any combination thereof. In some examples, the sequence of protein domain representations of the genome is generated by a process that includes searching the protein sequence of each gene in the genome and identifying protein domains within the protein sequence by sequence alignment against, for example, CDD, Gene3D, PANTHER, Pfam, ProSitePatterns, ProSiteProfiles, SUPERFAMILY, SMART, TIGRFAM, SFLD, Hamap, Coils, PRINTS, PIRSR, AntiFam, MobiDBLite, or PIRSF representation databases, respectively.

[0078] In some examples, the representation of the genome may further include associated gene ontology (GO) terms, identification of any known resistance genes present in the genome, identification of additional regulatory elements such as promoters, enhancers or silencers present in the genome, or identification of additional epigenetic elements such as histone folding, DNA methylation or acetylation present in the genome.

[0079] In some examples, the multiple related genomes received as input may include 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or more than 100 genomes, or any number within this range.

[0080] In step 104 of Figure 1, clusters of orthologous genes (COGs) can be identified in the multiple related genomes, as described elsewhere herein. In some examples, for example, identifying COGs can include identifying orthologous genes in the multiple related genomes as bidirectional best hits (BBHs) using BLAST, followed by clustering the identified orthologous genes. In some examples, identifying COGs may include identifying orthologous genes in multiple related genomes using software tools such as orthoMCL (Li, et al. (2003), "OrthoMCL: Identification of Ortholog Groups for Eukaryotic Genomes", Genome Research 13:2178-2189), or orthoFinder (Emms, et al. (2015), "OrthoFinder: Solving Fundamental Biases in Whole Genome Comparisons Dramatically Improves Orthogroup Inference Accuracy", Genome Biology 16:157).

[0081] In step 106 of FIG. 1, pairwise co-occurrence metrics, pairwise co-localization metrics, pairwise co-evolution metrics, pairwise co-regulation metrics, or any combination thereof, may be determined for the identified COGs, as described elsewhere herein.

[0082] In some examples, determining pairwise coevolution metrics for COGs may include calculating the percentage identity between each pair of protein sequences in each COG of a pair of COGs to identify shared protein sequences, calculating a Pearson's correlation coefficient for each pair of COGs that contains a certain minimum number of shared protein sequences to estimate a coevolution rate, filtering the COGs by removing COGs whose pairwise Pearson's correlation coefficient is below a predetermined threshold and clustering the remaining COGs according to their estimated coevolution rates, and performing a functional enrichment analysis to remove clusters of COGs enriched for essential metabolic functional categories.

[0083] In some examples, the specified minimum number of shared protein sequences used to identify COGs for which a Pearson correlation coefficient can be calculated can be 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or more than 100 shared protein sequences.

[0084] In some examples, the predetermined threshold value of the Pearson correlation coefficient used to exclude COGs from clustering may correspond to a Pearson correlation coefficient value of 0.7, 0.8, 0.9, 0.95, 0.98, or 0.99, or any value within this range.

[0085] In some examples, clustering the remaining COGs according to their estimated coevolutionary rates may include the use of Markov clustering (MCL) or hierarchical clustering algorithms.

[0086] In some examples, determining the pairwise co-regulation metrics for the COGs may include extracting intergenic regions within each COG, performing de novo detection of sequence motifs within the extracted intergenic regions to identify putative cis-regulatory elements or transcription factor binding sites (TFBS), comparing the identified putative cis-regulatory elements or TFBS for each COG with those identified across all other COGs to determine a pairwise motif similarity score between the COGs, filtering the COGs to exclude COGs whose pairwise motif similarity scores have a p-value below a predetermined threshold, and clustering the filtered COGs based on the pairwise motif similarity scores to identify co-regulated COG clusters.

[0087] In some examples, the predetermined threshold pairwise motif similarity score p-value may correspond to a p-value of 0.001, 0.01, 0.02, 0.03, 0.04, 0.05, or any p-value within this range.

[0088] In some examples, clustering the remaining COGs according to motif similarity scores may include the use of Markov clustering (MCL) or hierarchical clustering algorithms.

[0089] In some examples, determining the pairwise co-occurrence metric of the COGs includes calculating a Jaccard coefficient for each pair of COGs (COG A and COG B) based on the following relationship:

number

[0090] In some examples, the predetermined threshold for the pairwise co-occurrence score may correspond to a pairwise co-occurrence score value of 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, or 0.9, or any value within this range.

[0091] In some examples, clustering the remaining COGs according to co-occurrence scores may include the use of Markov clustering (MCL) or hierarchical clustering algorithms.

[0092] In some examples, determining the pairwise colocalization metric of the COGs includes calculating a proximity score for each pair of corresponding gene sequences in the pair of COGs (COG A and COG B) based on the following relationship:

number

[0093] In some examples, the clustering of COGs by averaged proximity scores may include the use of a Markovian clustering (MCL) or hierarchical clustering algorithm. In some examples, the clustering is performed using a Markovian clustering (MCL) algorithm and an MCL inflation parameter value ranging from 1.5 to 3.0. In some examples, the MCL inflation parameter value may be 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, or any value within this range.

[0094] In step 108 of FIG. 1, pairwise functional association scores for the identified COGs are determined based on the pairwise co-occurrence metrics, co-localization metrics, co-evolution metrics, co-regulation metrics, or combinations thereof determined in step 106.

[0095] In some examples, the pairwise functional association score of the identified COGs includes an algebraic function of the determined pairwise coevolution metric, pairwise coregulation metric, pairwise co-occurrence metric, pairwise colocalization metric, or any combination thereof. In some examples, for example, the pairwise functional association score of the identified COGs is based on the addition of the determined pairwise coevolution metric, pairwise coregulation metric, pairwise co-occurrence metric, pairwise colocalization metric, or any combination thereof. Such pairwise functional association scores can have values ​​ranging from 0 to 1, with higher values ​​indicating higher functional association. In some examples, for example, a pairwise functional association score value of 0.5, 0.6, 0.7, or 0.8 or greater can be used as a cutoff threshold applied prior to clustering. In some examples, as shown in FIG. 3, the pairwise coevolution metric, pairwise coregulation metric, pairwise co-occurrence metric, pairwise colocalization metric can be applied sequentially as an alternative to combining them in a single pairwise functional association score.

[0096] In step 110 of Figure 1, the identified COGs are clustered according to their pairwise functional association scores to group functionally related COGs. In some examples, clustering the identified COGs according to their pairwise functional association scores may include the use of Markov clustering (MCL) or hierarchical clustering algorithms.

[0097] In some examples, the method may further include determining a horizontal introgression metric, for example, based on calculation of a codon adaptation index (CAI) or dinucleotide signature dissimilarity index (DSDI). The CAI value may range, for example, from 0 to 1 (e.g., values ​​of 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0, or any value within this range), with lower values ​​indicating horizontal introgression. In some examples, the horizontal introgression metric may be used to further refine the clustering of co-localized, co-occurring and / or co-evolving COGs.

[0098] In some examples, the horizontal introgression metric can be used as part of a post-processing step to search for nearby horizontally transferred genes that were missed in the upstream clustering step.

[0099] In step 112 of Figure 1, a determination that a COG cluster is a network of functionally related genes in a particular functional category is output based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes in a particular functional category. In some examples, the identification of a putative gene network does not require the identification of associated genes known to be associated with a gene network of a given functional category.

[0100] In some examples, functional enrichment analysis may include testing for enrichment of genes within functional categories known to be associated with, for example, biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. Functional categories known to be associated with BGCs may include, for example, gene ontology terms or KEGG pathways known to be associated with BGCs.

[0101] Examples of gene ontology terms known to be associated with BGCs include, but are not limited to, GO:0019748 (secondary metabolic process), GO:0044550 (secondary metabolite biosynthetic process), GO:0030639 (polyketide biosynthetic process), GO:0030638 (polyketide metabolic process), GO:0043455 (regulation of secondary metabolic process), GO:1900539 (fumonisin metabolic process), or any combination thereof.

[0102] Examples of KEGG pathways known to be associated with BGCs include, but are not limited to, M00778 (type II polyketide backbone biosynthesis), M00095 (C5 isoprenoid biosynthesis, mevalonate pathway), M00937 (aflatoxin biosynthesis), M00893 (lovastatin biosynthesis), or any combination thereof.

[0103] In some examples, functional enrichment analysis may include testing for enrichment of protein domain representation known to be associated with, for example, biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs.

[0104] Examples of protein domains known to be associated with BGCs include, but are not limited to, PFAM domains, Conserved Domain Database (CDD) domains, or TIGRFAM domains known to be associated with BGCs. Examples of protein domains associated with BGCs include, but are not limited to, PF00550 (phosphopanteine ​​binding site), PF00501 (AMP-binding enzyme), PF07690 (Major Facilitator Superfamily), PF00067 (Cytochrome P450), PF00698 (acyltransferase domain), PF08242 (methyltransferase domain), or any combination thereof.

[0105] In some examples, the method may further include identifying a putative target of a secondary metabolite synthesized by the putative BGC by identifying a protein sequence that is not a component of a known BGC, determining a pairwise functional association score of the identified protein sequence and the putative BGC, and identifying the putative target of the secondary metabolite based on a comparison of the pairwise functional association score to a predetermined threshold. In some examples, identification of a putative BGC does not require identification of an associated core synthase.

[0106] In some examples, the pairwise functional association score may include a co-regulation score, where the pre-defined threshold for the co-regulation score corresponds to a p-value of, for example, 0.05, 0.04, 0.03, 0.02, 0.01, or 0.001 or less.

[0107] In some examples, the pairwise functional association score may include a coevolution score, and the first predetermined threshold of the coevolution score corresponds to a coevolution score value of 0.7, 0.8, 0.9, or 0.95 or greater.

[0108] In some examples, the pairwise functional association scores can include a co-occurrence score, and the first predetermined threshold of the co-occurrence score corresponds to a co-occurrence score of 0.5, 0.6, 0.7, 0.8, 0.9, or 0.95 or greater. How to use

[0109] The disclosed computer-based methods for identifying gene networks, and in particular for identifying BGCs, have a variety of applications, including, for example, performing further evaluation of genes predicted to be part of a BGC to (i) identify homologs or orthologs of one or more target sequences (e.g., gene sequences) of interest in one or more target genomes, (ii) identify resistance genes to secondary metabolites produced by the BGC in the target genome, (iii) predict the function of a secondary metabolite produced by the BGC, and / or (iv) identify BGCs that encode biosynthetic enzymes for producing a secondary metabolite having an activity of interest (e.g., a therapeutic activity of interest), etc.

[0110] Methods for evaluating genes embedded within or associated with BGCs to identify resistance genes (e.g., "embedded target genes" (ETaGs) or "non-embedded target genes" (NETaGs)) are described in International Patent Application Nos. PCT / US2022 / 049016, PCT / US2022 / 049040, PCT / US2022 / 079965 and PCT / US2022 / 080447, the contents of each of which are incorporated herein in their entirety. In some examples, for example, a method for identifying resistance genes (e.g., embedded target genes (ETaGs) and / or non-embedded target genes (NETaGs)) includes receiving a selection of at least one target sequence of interest, receiving a selection of target genomes from a genomics database, where the selection of target genomes includes a plurality of target genomes from organisms known to produce or likely to produce secondary metabolites, performing a search to identify homologs of the at least one target sequence in the plurality of target genomes, generating a phylogenetic tree based on the identified homologs of the at least one target sequence, classifying genomes of the plurality of target genomes as positive genomes or negative genomes based on the phylogenetic tree, where a positive genome is a genome that belongs to a clade in which multiple copies of at least one target sequence homolog are present, and a negative genome is a genome that belongs to a clade in which multiple copies of at least one target sequence homolog are present. a genome belonging to a clade in which a single copy of at least one target sequence homolog is present; a target sequence homolog present in multiple copies in the positive genome is a putative resistance gene (e.g., pETaG or pNETaG); and based at least in part on the classification of the positive and negative genomes, determining whether or not the positive and negative genomes have at least one of the following: i) one or more scores indicating co-occurrence of at least one target sequence homolog (putative resistance gene (e.g., pETaG or pNETaG)) with one or more genes associated with a biosynthetic gene cluster (BGC); ii) one or more scores indicating co-evolution of at least one target sequence homolog (putative resistance gene (e.g., pETaG or pNETaG)) with one or more genes associated with a BGC; iii) one or more scores indicating co-evolution of at least one target sequence homolog (putative resistance gene (e.g., pETaG or pNETaG)) with one or more genes associated with a BGC;and iv) one or more scores indicating co-regulation with one or more genes associated with the BGC, and iv) one or more scores indicating co-expression of at least one target sequence homolog (putative resistance gene (e.g., pETaG or pNETaG)) with one or more genes associated with the BGC, and determining the likelihood that the putative resistance gene (e.g., pETaG or pNETaG) is a resistance gene (e.g., an embedded target gene (ETaG) or a non-embedded target gene (NETaG)) based on the at least one genomic parameter.

[0111] In some examples, determining the likelihood that the putative resistance gene is a resistance gene may comprise comparing at least one determined genomic parameter to at least one predefined threshold.

[0112] In some examples, the selection of at least one target sequence of interest may be provided as an input by a user of a system configured to execute a computer-implemented method. In some examples, for example, the at least one target sequence of interest may include a sequence of a gene identified as belonging to a BGC by any of the methods described elsewhere herein.

[0113] In some examples, at least one target sequence of interest may comprise an amino acid sequence, a nucleotide sequence, or any combination thereof.In some examples, at least one target sequence of interest may comprise a peptide sequence or a portion thereof, a protein sequence or a portion thereof, a protein domain sequence or a portion thereof, a gene sequence or a portion thereof, or any combination thereof.In some examples, at least one target sequence of interest may comprise a mammalian sequence, a human sequence, a plant sequence, a fungal sequence, a bacterial sequence, an archaeal sequence, a viral sequence, or any combination thereof.

[0114] In some examples, at least one target sequence of interest may include a primary target sequence and one or more related sequences. In some examples, the one or more related sequences may include a sequence that is functionally related to the primary target sequence. In some examples, the one or more related sequences may include a sequence that is pathway related to the primary target sequence.

[0115] In some examples, the selection of the target genome may be submitted as an input by a user of a system configured to execute the computer-implemented method. In some examples, the multiple target genomes may include plant genomes, fungal genomes, bacterial genomes, or any combination thereof. In some examples, the genomics database may include a public genomics database. In some examples, the genomics database includes a proprietary genomics database.

[0116] In some examples, the search for identifying at least one target sequence homolog (e.g., a homolog of a gene sequence identified as belonging to a BGC) can include identifying the homolog based on a probabilistic sequence alignment model. In some examples, the probabilistic sequence alignment model is a profile hidden Markov model (pHMM). In some examples, the homolog is identified based on comparing the probabilistic sequence alignment model score with a predetermined threshold value.

[0117] In some examples, the search for identifying homologs of at least one target sequence may include identifying homologs based on sequence alignment using a local sequence alignment search tool, calculating sequence homology metrics based on the alignment, and comparing the calculated sequence homology metrics with a predetermined threshold value.In some examples, the local sequence alignment search tool includes BLAST, DIAMOND, HMMER, Exonerate, or ggsearch.In some examples, the predetermined threshold value may include a threshold value of sequence identity percentage, sequence coverage percentage, E value, or bit score value.

[0118] In some examples, the search for identifying homologs of at least one target sequence may include identifying homologs based on the use of gene and / or protein domain annotation tools. In some examples, the gene and / or protein domain annotation tools include InterProScan or EggNOG.

[0119] In some examples, generating a phylogenetic tree based on the identified homologs of at least one target sequence may include aligning the homolog sequences using an alignment software tool, trimming the aligned homolog sequences using a sequence trimming software tool, and constructing a phylogenetic tree using a phylogenetic tree construction software tool. In some examples, the alignment software tool includes MAFFT, MUSCLE, or ClustalW. In some examples, the sequence trimming software tool includes trimAI, GBlocks, or ClipKIT. In some examples, the phylogenetic tree construction software tool includes FastTree, IQ-TREE, RAxML, MEGA, MrBayes, BEAST, or PAUP. In some examples, the construction of the phylogenetic tree may be based on a maximum likelihood algorithm, a maximum parsimony algorithm, a neighbor-joining algorithm, a distance matrix algorithm, or a Bayesian estimation algorithm.

[0120] In some examples, the score or scores indicating co-occurrence may be determined based on the identification of a positive correlation between the presence of multiple copies of the putative resistance gene and the presence of one or more genes of the BGC in the positive genome. In some examples, identifying a positive correlation between the presence of multiple copies of the putative resistance gene and the presence of one or more genes of the BGC in the positive genome may include using a clustering algorithm to cluster aligned protein sequences, aligned nucleotide sequences, aligned protein domain sequences, or aligned pHMMs for a group of BGCs to identify a BGC community in the multiple target genomes. In some examples, identifying a positive correlation between the presence of multiple copies of the putative resistance gene and the presence of one or more genes of the BGC in the positive genome may include using a phylogenetic analysis of protein sequences or protein domains of a group of BGCs to identify a BGC community in the multiple target genomes. In some examples, identifying a positive correlation between the presence of multiple copies of the putative resistance gene and the presence of one or more genes of the BGC in the positive genome may include selecting genomes with a specific taxonomy to identify a BGC community in the multiple target genomes.

[0121] In some examples, the one or more scores indicating the co-evolution of the putative resistance gene and the one or more genes associated with the BGC may be determined based on a co-evolution correlation score, a co-evolution rank score, a co-evolution gradient score, or any combination thereof. In some examples, the co-evolution correlation score may be based on a correlation between the pairwise sequence identity percentage of the cluster of orthologous groups (COGs) for the putative resistance gene and the pairwise sequence identity percentage of the cluster of orthologous groups (COGs) for one of the one or more genes associated with the BGC. In some examples, the co-evolution rank score may be based on a ranking of correlation coefficients of the COGs including one of the one or more genes associated with the BGC in ascending order relative to the COG including the putative resistance gene. In some examples, in the event of a tie in distance scores, the rank of all the tied COGs may be set equal to the lowest rank in the group. In some examples, the co-evolution gradient score may be based on an orthogonal regression of the pairwise sequence identity percentage of the COGs for the putative resistance gene and the pairwise sequence identity percentage of the COGs for one of the one or more genes associated with the BGC. In some instances, only COGs arising from unique positive genomes with more than three genes remaining after removing the corresponding genes from the negative genome are used to evaluate the coevolutionary correlation score, coevolutionary rank score or coevolutionary gradient score.

[0122] In some examples, the score or scores indicative of co-regulation may be based on DNA motif detection from intergenic sequences of one or more genes associated with the BGC and putative resistance genes.

[0123] In some examples, the score or scores indicative of co-expression may be based on differential expression and / or clustering analysis of global transcriptome data.

[0124] In some examples, the one or more genes associated with a biosynthetic gene cluster (BGC) can include anchor genes, core synthase genes, biosynthetic genes, genes not involved in the biosynthesis of a secondary metabolite produced by the BGC, or any combination thereof.

[0125] In some examples, the putative resistance gene can be a putative embedded target gene (pETaG) or a putative non-embedded target gene (pNETaG).

[0126] In some examples, the resistance gene can be an embedded target gene (ETaG) or a non-embedded target gene (NETaG).

[0127] In some examples, a method for predicting a function of a secondary metabolite includes receiving a selection of at least one target sequence of interest, where the at least one target sequence of interest corresponds to a gene sequence associated with a biosynthetic gene cluster (BGC) known to produce a secondary metabolite; receiving a selection of target genomes from a genomics database, where the selection of target genomes includes a plurality of target genomes from organisms known to produce a secondary metabolite; performing a search to identify homologs of the at least one target sequence in the plurality of target genomes; generating a phylogenetic tree based on the identified homologs of the at least one target sequence; classifying genomes of the plurality of target genomes as positive genomes or negative genomes based on the phylogenetic tree, where a positive genome is a genome that belongs to a clade in which multiple copies of at least one target sequence homolog are present and a negative genome is a genome that belongs to a clade in which a single copy of at least one target sequence homolog is present; and determining, based at least in part on the classification of positive and negative genomes, at least one genomic parameter selected from the following: i) one or more scores indicative of co-occurrence of at least one target sequence homolog (putative resistance gene) and one or more genes associated with the BGC; ii) one or more scores indicative of co-evolution of at least one target sequence homolog (putative resistance gene) and one or more genes associated with the BGC; iii) one or more scores indicative of co-regulation of at least one target sequence homolog (putative resistance gene) and one or more genes associated with the BGC; and iv) one or more scores indicative of co-expression of at least one target sequence homolog (putative resistance gene) and one or more genes associated with the BGC; and determining, based on the at least one genomic parameter, a likelihood that the putative resistance gene is a resistance gene encoding a protein target acted upon by a secondary metabolite.

[0128] In some examples, a method for identifying biosynthetic gene clusters (BGCs) encoding biosynthetic enzymes for producing a secondary metabolite having an activity of interest includes receiving a selection of at least one target sequence of interest, wherein the at least one target sequence of interest comprises a sequence encoding a therapeutic target of interest; receiving a selection of target genomes from a genomics database, wherein the selection comprises a plurality of target genomes from organisms known to produce secondary metabolites; performing a search to identify homologs of the at least one target sequence in the plurality of target genomes; generating a phylogenetic tree based on the identified homologs of the at least one target sequence; classifying genomes of the plurality of target genomes as positive genomes or negative genomes based on the phylogenetic tree, wherein a positive genome is a genome belonging to a clade in which multiple copies of at least one target sequence homolog are present, a negative genome is a genome belonging to a clade in which a single copy of at least one target sequence homolog is present, and a target sequence homolog present in multiple copies in a positive genome is a putative genome. and determining, based at least in part on the classification of the positive and negative genomes, at least one genomic parameter selected from the following: i) one or more scores indicative of co-occurrence of at least one target sequence homolog (putative resistance gene) and one or more genes associated with a biosynthetic gene cluster (BGC); ii) one or more scores indicative of co-evolution of at least one target sequence homolog (putative resistance) and one or more genes associated with a BGC; iii) one or more scores indicative of co-regulation of at least one target sequence homolog (putative resistance gene) and one or more genes associated with a BGC; and iv) one or more scores indicative of co-expression of at least one target sequence homolog (putative resistance gene) and one or more genes associated with a BGC; and determining, based on the at least one genomic parameter, a likelihood that the putative resistance gene is an actual resistance gene associated with a BGC that produces a secondary metabolite that acts on the protein product encoded by the resistance gene.

[0129] In some examples, the methods of the disclosure may further include performing an in vitro assay, e.g., an assay to detect or measure the activity (e.g., receptor binding activity, enzyme activating activity, enzyme inhibitory activity, etc.) of the secondary metabolite (or analog thereof) against a mammalian (e.g., human) protein encoded by a mammalian (e.g., human) gene that is homologous to ETaG or NETaG identified in an organism that includes a biosynthetic gene cluster (BGC) that produces the secondary metabolite. In some examples, the methods may further include performing an in vitro assay to detect or measure the activity (e.g., receptor binding activity, enzyme activating activity, enzyme inhibitory activity, etc.) of the secondary metabolite (or analog thereof) against a protein (e.g., reptile, avian, amphibian, plant, fungal, bacterial, or viral protein) encoded by a reptile, avian, amphibian, plant, fungal, bacterial, or viral gene that is homologous to ETaG or NETag identified in an organism that includes a biosynthetic gene cluster (BGC) that produces the secondary metabolite.

[0130] In some examples, the methods of the disclosure may further include performing an in vivo assay, e.g., an assay to detect or measure the activity (e.g., receptor binding activity, enzyme activating activity, enzyme inhibitory activity, intracellular signaling pathway activity, disease response, etc.) of the secondary metabolite (or analog thereof) against a mammalian (e.g., human) protein encoded by a mammalian (e.g., human) gene that is homologous to ETaG or NETaG identified in an organism that includes a biosynthetic gene cluster (BGC) that produces the secondary metabolite. In some examples, the methods may further include performing an in vivo assay to detect or measure the activity (e.g., receptor binding activity, enzyme activating activity, enzyme inhibitory activity, intracellular signaling pathway activity, disease response, etc.) of the secondary metabolite (or analog thereof) against a protein (e.g., reptile, avian, amphibian, plant, fungal, bacterial, or viral protein) encoded by a reptile, avian, amphibian, plant, fungal, bacterial, or viral gene that is homologous to ETaG or NETaG identified in an organism that includes a biosynthetic gene cluster (BGC) that produces the secondary metabolite.

[0131] In some examples, the methods of the disclosure may be used to identify and / or characterize, for example, a mammalian (e.g., human) target of a secondary metabolite (or analog thereof) produced by a BGC. In some examples, the methods of the disclosure may be used to identify and / or characterize a reptile, bird, amphibian, plant, fungal, bacterial, viral target of a secondary metabolite (or analog thereof) produced by a BGC, or a target from any other organism.

[0132] In some examples, the methods of the present disclosure may be used, for example, in drug discovery efforts to identify small molecule modulators of mammalian (e.g., human) target genes. In some examples, the methods of the present disclosure may be used to identify small molecule modulators of reptile target genes, avian target genes, amphibian target genes, plant target genes, fungal target genes, bacterial target genes, viral target genes, or target genes from any other organism.

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

[0134] In some examples, the disclosure provides a method of regulating a human target (or a target from another organism), 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 or NETaG associated with the BGC as determined using any one of the methods described herein.

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

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

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

[0138] In some examples, the analog of the product includes the structural core of the product. In some examples, the biosynthetic product is cyclic, e.g., monocyclic, bicyclic, or polycyclic, and the structural core of the product is or includes a monocyclic, bicyclic, or polycyclic ring system. In some examples, the structural core of the product includes one ring of the bicyclic or polycyclic ring system of the product. In some examples, the product is or includes a polypeptide, and the structural core is the backbone of the polypeptide. In some examples, the product is or includes a polyketide, and the structural core is the backbone of the polyketide. In some examples, the analog is a substituted biosynthetic product that includes one or more suitable substitutents. A system for identifying gene networks containing functionally related genes

[0139] Also disclosed herein are systems designed to implement any of the disclosed methods for identifying gene networks (e.g., BGCs) that comprise a set of functionally related genes. The systems may, for example, include one or more processors, and communicatively coupled to the one or more processors, and when executed by the one or more processors, the systems may include: receiving as input a selection of genomes for analysis, the selection of genomes comprising a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining a pairwise coevolution metric, a pairwise coregulation metric, a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof, for the identified COGs; and determining a pairwise coevolution metric, a pairwise coregulation metric, a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof, for the determined COGs. and a memory unit configured to store instructions to cause: determining pairwise functional association scores for the identified COGs based on a pairwise co-occurrence metric, a pairwise colocalization metric, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group into functionally related COGs; and outputting a determination that the COG cluster is a network of functionally related genes in a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes in the particular functional category.

[0140] Computer Processors and Systems FIG. 2 illustrates an example of a computing device according to one or more examples of the present disclosure. The device 200 may be a host computer connected to a network. The device 200 may be a client computer or a server. As shown in FIG. 2, the device 200 may be any suitable type of microprocessor-based device, such as a personal computer, a workstation, a server, or a handheld computing device (portable electronic device) such as a phone or tablet. The device may include, for example, one or more of a processor 210, an input device 220, an output device 230, a storage 240, and a communication device 260. The input device 220 and the output device 230 may generally correspond to those described above, and they may be connectable or integrated with the computer.

[0141] The input device 220 may be any suitable device for conveying input, such as a touch screen, a keyboard or keypad, a mouse, or a voice recognition device. The output device 230 may be any suitable device for conveying output, such as a touch screen, a tactile device, or a speaker.

[0142] Storage 240 may be any suitable device providing storage, such as RAM, cache, electrical, magnetic, or optical memory, including a hard drive, or removable storage disk. Communications device 260 may include any suitable device capable of sending and receiving signals over a network, such as a network interface chip or device. The components of the computer may be connected in any suitable manner, such as via a physical bus 270 or wirelessly.

[0143] The software 250 that may be stored in the memory / storage 240 and executed by the processor 210 may include, for example, programming that embodies functions of the present disclosure (e.g., as embodied in the devices described above).

[0144] The software 250 may also be stored and / or propagated in any non-transitory computer-readable storage medium for use by or in connection with an instruction execution system, apparatus, or device, such as those described above, that can fetch instructions associated with the software from and execute the instructions. In the context of this disclosure, a computer-readable storage medium may be any medium, such as storage 240, that can contain or store programming for use by or in connection with an instruction execution system, apparatus, or device.

[0145] The software 250 may also be propagated in any transport medium for use by or in connection with an instruction execution system, apparatus, or device, such as those described above, that can fetch instructions associated with the software from the instruction execution system, apparatus, or device and execute the instructions. In the context of this disclosure, a propagation medium may be any medium that can communicate, propagate, or transmit programming for use by or in connection with an instruction execution system, apparatus, or device. Propagation readable media may include, but are not limited to, electronic, magnetic, optical, electromagnetic, or infrared wired or wireless propagation media.

[0146] The device 200 may be connected to a network, which may be any suitable type of interconnected communication system. The network may implement any suitable communication protocol and may be protected by any suitable security protocol. The network may include any suitable configuration of network links capable of implementing transmission and reception of network signals, such as wireless network connections, T1 or T3 lines, cable networks, DSL, or telephone lines.

[0147] Device 200 may implement any operating system suitable for operating on a network. Software 350 may be written in any suitable programming language, such as C, C++, Java, or Python. In various embodiments, application software embodying functionality of the present disclosure may be deployed in different configurations, such as, for example, as a web-based application or web service, in a client / server configuration, or via a web browser.

[0148] example Example 1 - Benchmarking known and predicted BGC Figure 3 presents a non-limiting schematic of a pipeline for BGC discovery and function assignment, based on the methods described elsewhere herein. The process 300 begins with the selection of a set of related genomes in step 302, followed by identification of clusters of orthologous genes (COGs) in step 304, co-occurrence analysis in step 306, co-localization analysis in step 308, co-regulation analysis in step 310, co-evolution analysis in step 312, and identification of a candidate set of functionally related loci (i.e., candidate gene networks) in step 314. In step 316, statistical analysis (and optionally other post-processing steps) is performed to identify specific types of gene networks, e.g., biosynthetic gene clusters (BGCs) in step 318.

[0149] To evaluate the above approach, closely related fungal genomes known to have BGCs for atpenins, cladosporins, citreoviridins, cyclosporins, restricticins or xanthocilins were identified. For each BGC type, 20 genomes showing approximately 90% genome identity were selected. In addition, six very distantly related genomes (approximately 60% identical to the target set) were also selected as outgroups. The best BLAST hits of the complements of proteins from all 26 genomes were calculated, and the resulting pairwise identity percentages were used as input for MCL clustering to group orthologous proteins (or related genes) in COGs. COGs containing proteins from all 20 target genomes and the six outgroup genomes were considered to be part of the core genomes of these species and therefore unlikely to play a role in secondary metabolism. These COGs were therefore excluded from downstream analyses to speed up the process.

[0150] Colocalization analysis was then performed by calculating the average proximity score between all pairs of COGs, as described elsewhere herein. This pairwise matrix was then used to cluster the COGs using MCL into clusters of colocalized (syntenous) COGs. These represent loci whose genetic composition is conserved among all or a subset of the 20 target genomes considered for the analysis of each target. These syntenous COGs then served as the starting point for calculating the other metrics used to delineate BGCs.

[0151] To assess co-regulation, motif similarity scores were calculated in MEME against a background distribution generated using all intergenic sequences from all genes across target species (as described elsewhere herein) by first extracting the intergenic regions of all genes in the syntenous COG and using this as input for de novo motif detection. If a significant motif was detected and at least one of the genes in the COG significantly matched the motif in its intergenic region, we retained the COG as a member of the syntenous COG, otherwise we dropped the COG. This generated syntenous co-regulated COGs.

[0152] To assess co-evolution (as described elsewhere herein), for each syntenous COG, all protein sequences in the COG were aligned using MAFFT and trimmed to remove gaps using trimAl. Pairwise percentage identity was then calculated between all pairs of genes in the COG. To estimate co-evolution, the Pearson correlation coefficient was calculated from the percentage identity between all pairs of COGs in the synthetic COG for the matched genome pairs. The generated pairwise matrix was filtered for correlation coefficients ≦0.9 and then used to cluster the COGs using MCL. Only the first-order clusters of co-evolving COGs were retained.

[0153] To determine which of the identified loci were candidate BGCs, we first converted the syntenous COG information into gene-level clusters for each of the 20 target species used in the analysis. We then used a hypergeometric test to evaluate each gene-level cluster for enrichment of PFAM domains known to be present within known BGCs. Clusters with p-values ​​≦0.01 were considered candidate BGCs.

[0154] To evaluate the performance of the above pipeline, we compute three metrics: i. Recall of target cluster genes (i.e., the proportion of known cluster genes identified)

number

number

number

[0155] The above BGC identification pipeline achieved a high degree of recall of the target cluster genes (i.e., the proportion of genes known to be involved in the biosynthesis of the captured target molecule, e.g., cyclosporine). The recall for the six evaluated target BGCs ranged between 0.88 and 1, with an average of 0.95 ± 0.04 (Figure 4). Furthermore, the accuracy of the approach (i.e., the proportion of predicted genes that are true target genes known to be involved in the biosynthesis of the target molecule) was also determined to be very high, ranging from 0.71 to 0.91 with an average accuracy of 0.83 ± 0.09 (Figure 4). Thus, overall, the proposed pipeline performed well in identifying specific BGC genes with an overall F-score of about 0.89.

[0156] To assess the performance of the pipeline more holistically, its predictions were compared to the complement of BGCs predicted by antiSMASH (Blin et al. (2021), supra) (i.e., core synthases containing BGCs). Based on this comparison, we observed that our BGC detection approach was able to detect approximately 70% of the antiSMASH predicted clusters on average across the 120 genomes selected for our analysis (Figure 4). Overall, these results show that our approach is able to detect clusters with their appropriate boundaries (only one or two additional genes or missing genes) while still performing well overall.

[0157] To evaluate how each component of our pipeline improves prediction performance, we evaluated the use of colocalization alone or in combination with coregulation and coevolution (Figure 5). Here, we see that the recall of the target cluster genes is generally not affected by the application of coregulation or coevolution, except in the case of cyclosporine, where the recall falls from 0.88 to 0.73 (Figure 5). This may be expected since the recall value was already very high. On the other hand, the accuracy is generally improved by applying coregulation and / or coevolution information for all targets except atpenine, which was already very accurate using only colocalization information. These results indicate that, depending on the cluster, the incorporation of orthogonal sources of genomic data can substantially improve the ability to accurately predict target cluster genes.

[0158] Finally, one of the key advantages of our proposed pipeline is that it can detect BGCs independent of core synthases, unlike many of the other currently available BGC detection algorithms such as antiSMASH. To verify this, we included the BGC of xanthocilline, which lacks a core synthase, in our benchmark test set. Given that this BGC does not have a canonical core synthase, it will not be picked up by algorithms such as antiSMASH. However, as shown in Figure 4, our proposed pipeline is able to identify this cluster with similar accuracy (0.71) and recall (1.0) as a core synthase-containing cluster. This indicates that our pipeline is a more comprehensive tool for BGC detection than the current state-of-the-art.

[0159] Exemplary embodiments Among the embodiments provided are the following: 1. A computer-implemented method for identifying networks of functionally related genes, comprising: receiving as input a selection of genomes for analysis, the selection of genomes including a plurality of related genomes; Identifying clusters of orthologous genes (COGs) in a plurality of related genomes; determining, for the identified COGs, a pairwise coevolution metric, a pairwise coregulation metric, a pairwise cooccurrence metric, a pairwise colocalization metric, or any combination thereof; determining pairwise functional association scores for the identified COGs based on the determined pairwise coevolution metric, pairwise coregulation metric, pairwise co-occurrence metric, pairwise colocalization metric, or any combination thereof; Clustering the identified COGs according to their pairwise functional association scores to group them into functionally related COGs; and outputting a determination that the COG cluster is a network of functionally related genes within a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes within a particular functional category; A method comprising: 2. The computer-implemented method of embodiment 1, wherein the functional enrichment analysis does not require identification of genes known to be associated with gene networks of a particular functional category. 3. The computer-implemented method of embodiment 1 or embodiment 2, wherein the functional enrichment analysis comprises testing for enrichment of genes within functional categories known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. 4. The computer-implemented method of embodiment 3, wherein the functional categories known to be associated with BGCs include gene ontology terms or KEGG pathways known to be associated with BGCs. 5. The computer-implemented method of embodiment 4, wherein the gene ontology terms known to be associated with the BGC include GO:0019748 (secondary metabolic process), GO:0044550 (secondary metabolite biosynthetic process), GO:0030639 (polyketide biosynthetic process), GO:0030638 (polyketide metabolic process), GO:0043455 (regulation of secondary metabolic process), GO:1900539 (fumonisin metabolic process), or any combination thereof. 6. The computer-implemented method of embodiment 4, wherein the KEGG pathways known to be associated with BGCs include M00778 (type II polyketide backbone biosynthesis) or M00095 (C5 isoprenoid biosynthesis, mevalonate pathway), M00937 (aflatoxin biosynthesis), M00893 (lovastatin biosynthesis), or any combination thereof. 7. The computer-implemented method of any one of embodiments 1 to 6, wherein the functional enrichment analysis comprises testing for enrichment of protein domain representation known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs. 8. The computer-implemented method of embodiment 7, wherein the protein domain representation known to be associated with a BGC comprises a PFAM domain representation, a Conserved Domain Database (CDD) domain representation, or a TIGRFAM domain representation known to be associated with a BGC. 9. Identify putative targets of secondary metabolites synthesized by putative BGCs Identifying protein sequences that are not components of known BGCs; Determining pairwise functional association scores of the identified protein sequences and putative BGCs; and Identifying putative targets of the secondary metabolites based on a comparison of the pairwise functional association scores to a first predetermined threshold. 9. The computer-implemented method of any one of claims 2 to 8, further comprising identifying the 10. The computer-implemented method of embodiment 9, wherein the pairwise functional association score comprises a co-regulation score, and the first predetermined threshold of the co-regulation score corresponds to a p-value of 0.05 or less. 11. The computer-implemented method of embodiment 9, wherein the pairwise functional association scores include a coevolution score, and the first predetermined threshold of the coevolution score corresponds to a coevolution score value of 0.7 or greater. 12. The computer-implemented method of embodiment 9, wherein the pairwise functional association scores include co-occurrence scores, and the first predetermined threshold of the co-occurrence score corresponds to a co-occurrence score of 0.5 or greater. 13. The computer-implemented method of any one of embodiments 2 to 12, wherein identification of a putative BGC does not require identification of an associated core synthase. 14. The computer-implemented method of any one of embodiments 1 to 13, wherein the plurality of related genomes comprises fungal, bacterial or plant genomes. 15. The computer-implemented method of any one of embodiments 1 to 14, wherein identifying COGs comprises using BLAST to identify orthologous genes in multiple related genomes as bidirectional best hits, followed by clustering the identified orthologous genes. 16. The computer-implemented method of any one of embodiments 1 to 15, wherein identifying COGs comprises using orthoMCL or orthoFinder to identify orthologous genes in a plurality of related genomes. 17. Determining pairwise coevolution metrics for COGs calculating the percentage of identity between each pair of protein sequences within each COG of the COG pairs to identify shared protein sequences; calculating the Pearson correlation coefficient for each pair of COGs that contain a certain minimum number of shared protein sequences to estimate the coevolution rate; filtering the COGs by removing COGs whose pairwise Pearson correlation coefficient is below a second predetermined threshold and clustering the remaining COGs according to their estimated coevolutionary rates; and performing functional enrichment analysis to filter out clusters of COGs enriched for essential metabolic function categories; 17. The computer-implemented method of any one of embodiments 1 to 16, comprising: 18. The computer-implemented method of embodiment 17, wherein the second predetermined threshold corresponds to a Pearson correlation coefficient value of 0.7, 0.8, 0.9, 0.95, 0.98, or 0.99. 19. The computer-implemented method of embodiment 17 or embodiment 18, wherein clustering the remaining COGs according to the estimated coevolutionary rates comprises using a Markov clustering (MCL) or hierarchical clustering algorithm. 20. Determining pairwise co-regulation metrics for COGs Extracting intergenic regions within each COG; performing de novo detection of sequence motifs within the extracted intergenic regions to identify putative cis-regulatory elements or transcription factor binding sites (TFBS); comparing the putative cis-regulatory elements or TFBSs identified for each COG with those identified across all other COGs to determine pairwise motif similarity scores between the COGs; filtering the COGs to exclude COGs whose pairwise motif similarity scores have a p-value below a third predetermined threshold; and Clustering the filtered COGs based on pairwise motif similarity scores to identify co-regulated COG clusters 20. The computer-implemented method of any one of embodiments 1 to 19, comprising: 21. The computer-implemented method of embodiment 20, wherein the third predetermined threshold corresponds to a p-value of 0.05. 22. The computer-implemented method of embodiment 20, wherein the third predetermined threshold corresponds to a p-value of 0.01. 23. A computer-implemented method according to any one of embodiments 20 to 22, wherein clustering the remaining COGs according to motif similarity scores comprises using a Markov clustering (MCL) or hierarchical clustering algorithm. 24. Determining pairwise co-occurrence metrics for COGs includes calculating the Jaccard coefficient for each pair of COGs (COG A and COG B) based on the following relationship:

number

number

[0160] From the above, it should be understood that although specific implementations of the disclosed method and system have been illustrated and described, various modifications can be made thereto and are contemplated herein. Also, the present invention is not intended to be limited by the specific examples provided herein. Although the present invention has been described with reference to the above specification, the description and illustration of the preferred embodiments herein are not meant to be construed in a limiting sense. Furthermore, it should be understood that all aspects of the present invention are not limited to the specific depictions, configurations or relative proportions described herein which depend upon various conditions and variables. Various changes in form and details of the embodiments of the present invention will be apparent to those skilled in the art. Therefore, it is contemplated that the present invention will cover any such modifications, variations, and equivalents.

Claims

1. 1. A computer-implemented method for identifying networks of functionally related genes, comprising: receiving as input a selection of genomes for analysis, the selection of genomes including a plurality of related genomes; identifying clusters of orthologous genes (COGs) in the plurality of related genomes; determining pairwise coevolution metrics, pairwise co-regulation metrics, pairwise co-occurrence metrics, pairwise co-localization metrics, or any combination thereof, for the identified COGs; determining pairwise functional association scores for the identified COGs based on the determined pairwise coevolution metrics, pairwise co-regulation metrics, pairwise co-occurrence metrics, pairwise co-localization metrics, or any combination thereof; clustering the identified COGs according to their pairwise functional association scores to group functionally related COGs; and outputting a determination that a COG cluster is a network of functionally related genes within a particular functional category based on a functional enrichment analysis performed on at least one COG cluster to identify COG clusters that are enriched for genes within said particular functional category; 11. A computer-implemented method comprising:

2. The computer-implemented method of claim 1 , wherein the functional enrichment analysis does not require identification of genes known to be associated with gene networks of a particular functional category.

3. 2. The computer-implemented method of claim 1, wherein the functional enrichment analysis comprises testing for enrichment of genes within functional categories known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs.

4. The computer-implemented method of claim 3 , wherein the functional categories known to be associated with BGC include gene ontology terms or KEGG pathways known to be associated with BGC.

5. 5. The computer-implemented method of claim 4, wherein gene ontology terms known to be associated with the BGC include GO:0019748 (secondary metabolic process), GO:0044550 (secondary metabolite biosynthetic process), GO:0030639 (polyketide biosynthetic process), GO:0030638 (polyketide metabolic process), GO:0043455 (regulation of secondary metabolic process), GO:1900539 (fumonisin metabolic process), or any combination thereof.

6. The computer-implemented method of claim 4, wherein the KEGG pathways known to be associated with the BGC include M00778 (type II polyketide backbone biosynthesis) or M00095 (C5 isoprenoid biosynthesis, mevalonate pathway), M00937 (aflatoxin biosynthesis), M00893 (lovastatin biosynthesis), or any combination thereof.

7. 2. The computer-implemented method of claim 1, wherein the functional enrichment analysis comprises testing for enrichment of protein domain representation known to be associated with biosynthetic gene clusters (BGCs), thereby identifying those COG clusters as putative BGCs.

8. 8. The computer-implemented method of claim 7, wherein the protein domain representation known to be associated with a BGC comprises a PFAM domain representation, a Conserved Domain Database (CDD) domain representation, or a TIGRFAM domain representation known to be associated with a BGC.

9. The putative targets of secondary metabolites synthesized by the putative BGCs were: Identifying protein sequences that are not components of known BGCs; determining a pairwise functional association score for the identified protein sequence and the putative BGC; and identifying putative targets of the secondary metabolites based on a comparison of the pairwise functional association scores to a first predetermined threshold. The computer-implemented method of claim 2 , further comprising identifying the target by:

10. 10. The computer-implemented method of claim 9, wherein the pairwise functional association scores comprise co-regulation scores, and the first predetermined threshold of the co-regulation scores corresponds to a p-value of less than or equal to 0.

05.

11. 10. The computer-implemented method of claim 9, wherein the pairwise functional association scores comprise a coevolution score, and wherein the first predetermined threshold of the coevolution score corresponds to a coevolution score value of 0.7 or greater.

12. 10. The computer-implemented method of claim 9, wherein the pairwise functional association scores comprise co-occurrence scores, and the first predetermined threshold for the co-occurrence scores corresponds to a co-occurrence score of 0.5 or greater.

13. 3. The computer-implemented method of claim 2, wherein identification of a putative BGC does not require identification of an associated core synthase.

14. The computer-implemented method of claim 1 , wherein the plurality of related genomes comprises fungal, bacterial, or plant genomes.

15. 2. The computer-implemented method of claim 1, wherein identifying COGs comprises using BLAST to identify orthologous genes in the plurality of related genomes as bidirectional best hits, followed by clustering the identified orthologous genes.

16. 2. The computer-implemented method of claim 1, wherein identifying COGs comprises using orthoMCL or orthoFinder to identify orthologous genes in the plurality of related genomes.

17. Determining pairwise coevolution metrics for the COGs calculating the percentage identity between each pair of protein sequences within each COG of the COG pairs to identify shared protein sequences; calculating the Pearson correlation coefficient for each pair of COGs that contain a certain minimum number of shared protein sequences to estimate the coevolution rate; filtering the COGs by excluding COGs whose pairwise Pearson correlation coefficient is less than a second predetermined threshold and clustering the remaining COGs according to the estimated coevolutionary rates; and performing functional enrichment analysis to filter out clusters of COGs enriched for essential metabolic function categories; The computer-implemented method of claim 1 , comprising:

18. 18. The computer-implemented method of claim 17, wherein the second predetermined threshold corresponds to a Pearson correlation coefficient value of 0.7, 0.8, 0.9, 0.95, 0.98, or 0.

99.

19. 18. The computer-implemented method of claim 17, wherein clustering the remaining COGs according to the estimated coevolutionary rates comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

20. Determining pairwise co-regulation metrics for COGs Extracting intergenic regions within each COG; performing de novo detection of sequence motifs within said extracted intergenic regions to identify putative cis-regulatory elements or transcription factor binding sites (TFBS); comparing the putative cis-regulatory elements or TFBSs identified for each COG with those identified across all other COGs to determine pairwise motif similarity scores between COGs; filtering the COGs to exclude COGs whose pairwise motif similarity scores have a p-value below a third predetermined threshold; and Clustering the filtered COGs based on the pairwise motif similarity scores to identify co-regulated COG clusters.

20. The computer-implemented method of claim 1, comprising:

21. 21. The computer-implemented method of claim 20, wherein the third predetermined threshold corresponds to a p-value of 0.

05.

22. 21. The computer-implemented method of claim 20, wherein the third predetermined threshold corresponds to a p-value of 0.

01.

23. 21. The computer-implemented method of claim 20, wherein clustering the remaining COGs according to the motif similarity scores comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

24. Determining pairwise co-occurrence metrics of the COGs Calculating the Jaccard coefficient for each pair of COGs (COG A and COG B) based on the following relationship: [Equation 1] filtering the COGs to exclude COGs whose pairwise co-occurrence scores are below a fourth predetermined threshold; and clustering the filtered COGs based on the pairwise co-occurrence scores to identify co-occurring COG clusters. The computer-implemented method of claim 1 , comprising:

25. 25. The computer-implemented method of claim 24, wherein the fourth predetermined threshold corresponds to a pairwise co-occurrence score value of 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, or 0.

9.

26. 26. The computer-implemented method of claim 24 or claim 25, wherein clustering the remaining COGs according to the co-occurrence scores comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

27. Determining pairwise colocalization metrics of COGs Calculating a proximity score for each pair of corresponding gene sequences in a pair of COGs (COG A and COG B) based on the following relationship: [Equation 2] determining an averaging colocalization metric for each pair of COGs; and clustering the COGs based on the averaged proximity scores to identify co-localized COG clusters; The computer-implemented method of claim 1 , comprising:

28. 28. The computer-implemented method of claim 27, wherein the clustering of COGs by averaged proximity scores comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

29. 29. The computer-implemented method of claim 27 or claim 28, wherein the clustering is performed using a Markov Clustering (MCL) algorithm and an MCL inflation value in the range of 1.5 to 5.

0.

30. 2. The computer-implemented method of claim 1, wherein the pairwise functional association scores of the identified COGs are based on adding the determined pairwise coevolution metrics, pairwise co-regulation metrics, pairwise co-occurrence metrics, pairwise co-localization metrics, or any combination thereof.

31. 2. The computer-implemented method of claim 1, wherein clustering the identified COGs according to their pairwise functional association scores comprises using a Markov clustering (MCL) or hierarchical clustering algorithm.

32. 10. The computer-implemented method of claim 1, further comprising determining horizontal introgression metrics based on calculating a codon adaptation index (CAI) or a dinucleotide signature divergence index (DSDI).

33. 33. The computer-implemented method of claim 32, wherein the horizontal introgression metric is used to further refine the clustering of co-localizing, co-occurring and / or co-evolving COGs.

34. 33. The computer-implemented method of Claim 32, wherein the horizontal introgression metric is used as part of a post-processing step to search for nearby horizontally transferred genes that were missed in an upstream clustering step.

35. The computer-implemented method of claim 3 , further comprising evaluating a gene identified as belonging to the putative BGC to determine whether it is a resistance gene.

36. 36. The computer-implemented method of claim 35, wherein the resistance gene is an embedded target gene (ETaG) or a non-embedded target gene (NETaG).

37. 36. The computer-implemented method of claim 35, further comprising performing an in vitro assay to test secondary metabolites produced by the putative BGC for activity against a resistance gene homolog or protein encoded thereby identified in the target genome.

38. 36. The computer-implemented method of claim 35, further comprising performing an in vivo assay to test secondary metabolites produced by the putative BGC for activity against a resistance gene homolog or protein encoded thereby identified in the target genome.

39. 39. The computer-implemented method of claim 37 or claim 38, wherein the target genome comprises a mammalian genome, a human genome, an avian genome, a reptile genome, an amphibian genome, a plant genome, a fungal genome, a bacterial genome, or a viral genome.

Images