05. Clustering BGCs into GCFs
In recent times (mid 2021), BiG-SCAPE (Munoz, Selem-Mjoica, Mullowney, et al. 2020) has been the popular tool for clustering of individual biosynthetic gene clusters (BGCs) into gene cluster families (GCFs) for lineage specific exploration and mapping of biosynthetic potential; however, it is designed to cluster extremely diverse BGCs from across the entire domain of bacteria and works at finding homology on the level of protein domains, specifically those noted to be of relevance to BGCs.
If only a single lineage is of interest however, the precision of clustering can be greatly improved by using whole protein sequences with a standard well established homolog identification tool. In lsaBGC-Cluster.py, and subsequent programs, we utilize homology relationships between proteins from antiSMASH (Blin et al. 2019) identified by OrthoFinder (Emms and Kelly 2019), with the goal to group together orthologous BGCs across isolate genomes into a single GCF. The simplified scope of the clustering (within a single lineage) in combination with the ability to utilize full protein sequences enables us to avoid under- and over-clustering BGCs to form GCFs, which we aim to ideally represent a set of orthologous BGCs where paralagous instances are split into separate GCFs where appropriate. Importantly, this definition of GCFs differs from how the BiG-SCAPE authors define it as sets of BGCs from "across a range of organisms that are linked to a highly similar natural product chemotype".
While BiG-SCAPE incorporates three distinct metrics to assess grouping BGCs into the same GCF, all three metrics are based on protein domains and sequences between homologous BGCs can be left un-utilized to help guide clustering. Working in the context of across the domain of bacteria, using protein domains makes a lot of sense and can help avoid incorporating noise from the sequences in between the domains; however, if working in the confines of a single lineage (e.g. species or genus), these stretches of sequence between domains can be extremely informative. By identifying homology across full protein sequences, we are further able to employ a very simple clustering approach whereby we only assess how many homolog groups BGCs share in common to determine orthologous GCFs. The clustering is controlled by two relatively simple to understand parameters, the Jaccard similarity index (how many homolog groups are shared by a pair of BGCs normalized by the number of homolog groups present in either BGC; default value = 50; range of accepted values = 0 to 100) and an inflation parameter (is a parameter for the Markov graph-based clustering algorithm - basically the larger this parameter the more granular the ; default value = 1.4; range of accepted values = 0.8 to 5).
By keeping the clustering focused on homology and ignoring syntenic similarity, we avoid the issue of overclustering and splitting orthologous GCFs into discrete categories, where issues can arise. For instance, when looking at ectoine BGCs across M. luteus genomes through a phylogenomic perspective we noticed that two very similar BGCs which were identified in phylogenetically similar isolates were assigned to different GCF categories by BiG-SCAPE:
Are you someone who likes to use a single cutoff for different taxa in an effort to standardize comparing results across them? Or are you someone who likes to tailor their parameters to each specific taxonomy of interest? Whichever one you are, we have you covered with sensible default parameters that should generally work well as well as an option to automatically test different combinations of parameters. We strongly advocate for this latter option however, as the scope of the lineage might influence which clustering parameters are best suited for the analyses of interest for you (e.g. if working in the context of just a single species or sub-species, perhaps a more course clustering configuration could lead to reduced singleton BGCs left unassigend to a GCF while still preserving the assumption of orthology between BGC instances within such clusters).
Users can simply run lsaBGC-Cluster.py with the --run_parameter_tests flag the first time around to generate a PDF report (8x11 for easy printing, for those in the US at least) to examine which combina ion of the Jaccard Similarity index and the MCL Inflation parameter best control the various tradeoffs between performing granular vs. coarse clustering.
Here is a description of the figures the report will feature:
The plot on the title page shows the number of GCFs which resulted (generally want to maximize as this will lead to more clean delineation of orthologous BGCs from paralogous BGCs) vs. the number of singleton BGCs which remain un-clustered(want to minimize to be inclusive of orthologous variants of BGCs which have significantly diverged/altered) when using different values of the Jaccard Index and MCL Inflation parameters. Each point corresponds to results from clustering the BGCs into GCFs at a specific Jaccard Index and MCL Inflation combination. A jitter is added to prevent overlapping points, so interpret the plot carefully.
The following 2nd page in the report contains four heatmaps, colored in different shades of the primary colors. Each cell in the heatmap corresponds to statistic assessing the quality of the GCF clustering performed using some specific MCL Inflation value (x-axis) and a specific Jaccard Similarity cutoff (y-axis).
- Blue Heatmap (upper left) - The shading represents the ratio of the number of BGCs which are singletons relative to the number of BGCs which were grouped into GCFs. The larger this ratio, the darker the shading.
- Green Heatmap (upper right) - The shading represents the number of GCFs formed from clustering with specific parameters. The larger this count, the darker the shading.
- Yellow Heatmap (lower left) - The shading represents the proportion of GCFs formed from clustering with specific parameters which contain BGCs with distinct annotations from AntiSMASH (e.g. an ectoine BGC and a terpene BGC being grouped together). BGCs which are not annotated are accounted for and considered as belonging to a distinct annotation category "unknown". The larger this proportion, the darker the shading.
- Red Heatmap (lower right) - The shading represents the proportion of GCFs formed from clustering with specific parameters which contain BGCs with distinct annotations from AntiSMASH (e.g. an ectoine BGC and a terpene BGC being grouped together). BGCs which are not annotated are ignored and not treated as conflicting evidence for multiple annotation categories being present in a single GCF. The larger this proportion, the darker the shading.
The 3rd page of the report features another set of overview plots, this time in the form of adjacent grids of boxplots. Each grid corresponds to a specific Jaccard Similarity cutoff values, whereas the individual boxplots correspond to different MCL Inflation parameter values. The boxplots depict a range of statistics across each GCF resulting from the respective clustering.
- First Row - Sizes of the single copy core of GCFs. This is the number of genes which are found in a single copy across all samples with the GCF.
- Second Row - Number of samples which have each GCF.
- Third Row - Average number of genes observed in each GCF.
- Fourth Row - Standard deviation for number of genes observed in each GCF.
As will often be apparent from these plots, the Jaccard Similarity cutoff parameter is by far more critical in impacting GCF clustering results than the MCL Inflation parameter which can be used more for fine tuning clustering results.
The 4th page of the report features yet another set of overview plots, this time they are alluvian plots showcasing how GCFs get collapsed or split up between two "parameter-adjacent" clustering attempts. By "parameter-adjacent", we mean that as one changes the MCL Inflation parameter you can see the influence this will have on the individual GCFs. For lineages with tons of GCFs, this visualization will not be very informative and I might incorporate a flag to filter its production if the case later on.
The 5th page and onwards, the bulk of the report, shows the individual GCF level details for each clustering attempt with a parameter combination. These pages comprise of 4 subplots.
The top row features two distinct plots. The left plot is a histogram of sample counts for each GCF (how many samples each GCF has), colored according to whether the GCF has a single copy core gene set. BGCs which are singletons are colored in grey as the notion of possessing a single copy core gene set doesn't apply to them. The right plot is a simple histogram of the standard deviation of gene counts per GCF. If GCFs have many small and large BGCs, this might imply the clustering had some issues and might be under-clustering (clustering too coarsely).
The middle row plot is a scatter plot of pie charts. Each pie chart represents a GCF and the coloring of the pie chart represents the annotation class assigned to members of that GCF, BGCs which are unannotated are colored in grey, the rest of the annotation categories are arbitrarily assigned a color. The point of this figure is not to draw biological insight, but simply to assess whether the GCFs are comprising of single annotation categories as would be expected.
The bottom row plot is a bar chart which shows a subset of the largest GCFs, ranked in order. Sometimes, for lineages with lots of BGCs and GCFs, the pie charts overlap in the middle plot and there are many annotations observed, so it becomes difficult to parse the scatter pie plot. This figure aims to alleviate such struggles.
An example report page for the species of M. luteus
An example report page for the genus of Streptomyces
Once again, the purpose of these plots is to help guide which values to select for the two parameters controlling clustering for your specific lineage of interest. However, it might be that during exploration of the individual GCFs, one of them really fascinates you, in such case an iterative process of using lsaBGC-See.py to visualize how individual BGC instances are added to a GCF or removed when using different clustering criterion could be potentially insightful. Additionally, the
usage: lsaBGC-Cluster.py [-h] -b BGC_LISTINGS -m ORTHOFINDER_MATRIX -o OUTPUT_DIRECTORY [-c CORES] [-i MCL_INFLATION] [-j JACCARD_CUTOFF] [-s] [-r]
Program: lsaBGC-Cluster.py
Author: Rauf Salamzade
Affiliation: Kalan Lab, UW Madison, Department of Medical Microbiology and Immunology
This program will cluster BGCs found by AntiSMASH using MCL based on similarity exhibited in ortholog group presence/
absence data. Clustering uses MCL.
optional arguments:
-h, --help show this help message and exit
-b BGC_LISTINGS, --bgc_listings BGC_LISTINGS
BGC listing file. Tab delimited: 1st column contains path to AntiSMASH BGC Genbank and 2nd column contains sample name.
-m ORTHOFINDER_MATRIX, --orthofinder_matrix ORTHOFINDER_MATRIX
OrthoFinder matrix.
-o OUTPUT_DIRECTORY, --output_directory OUTPUT_DIRECTORY
Output directory.
-c CORES, --cores CORES
Number of cores to use for MCL step.
-i MCL_INFLATION, --mcl_inflation MCL_INFLATION
Inflation parameter to be used for MCL.
-j JACCARD_CUTOFF, --jaccard_cutoff JACCARD_CUTOFF
Cutoff for Jaccard similarity of homolog groups shared between two BGCs.
-s, --split_by_annotation
Partition BGCs into groups based on annotation first.
-r, --run_parameter_tests
Run tests for selecting the best inflation parameter and jaccard for MCL analysis and exit.