Toxin Scanner

BioBrick discovery

For iGEM 2016 we designed a high-throughput pipeline for the identification of novel proteins directly from raw genome sequencing data. Given the specificity of our tool and the importance of biobrick discovery in iGEM, we made it publicly available for everyone to modify and use. The tool can be found in this Gitlab repository .
For the purpose of the BeeT project, we use it as a cry toxin predictor, given genomes of selected bacteria.

Toxin Specificity

For this project we need a toxin that specifically targets Varroa destructor. There are many insecticidal protein families, but the group we are most interested in are the crystal (cry) proteins. These are usually found on plasmids from Bacillus Thuringiensis and related species. ¹⁴

We know Varroa-specific insecticidal activity exists in Bacillus thuringiensis and related species, as shown in: "In vitro susceptibility of Varroa destructor and Apis mellifera to native strains of Bacillus thuringiensis." by Alquisira-Ramírez et al. ² In this paper, several isolates are described that cause a mite mortality of up to 100%. Importantly, the strains also showed no insecticidal activity against bee larvae. Because of this we started several sub-projects in parallel to maximize our chances of finding a viable V. destructor-killer. The specificity part of our project focuses on creating V. destructor-binding cry toxins and finding V. destructor-specific insecticidal proteins.

There already exists a publication about a tool called “Bt Toxin Scanner”¹. This tool does not fully support local deployment, which is needed for high-throughput analysis. Also, because of the relatively basic analysis done by the tool, we decided to develop our own tool that is fully open-source and improves upon the analysis techniques used in Bt Toxin Scanner. Our goal with this tool is to run raw sequencing files, and deliver potential cry proteins with just the click of a button.

This tool was made in preparation for results of the V. destructor isolates part of the project. We assume these V. destructor-specific insecticidal proteins, to be of the cry protein family. Or at the very least quite similar to them.
These cry proteins are a diverse group, but are known to be highly specific for individual insects and nematodes. Cry proteins are not necessarily a group of proteins that all perform the same function in the same manner. The distinction between cry and non-cry proteins is defined by a committee: Cry Protein website Based on 45% sequence similarity there are over 70 groups. This high amount of diversity makes it hard to predict when something is or isn't a cry protein.

Testing the tool

We tested the tool on a genome sample from a study with accession number PRJEB5931. ¹³ This genome was found after a co-evolution experiment, and a Bacillus thuringiensis with known nematicidal cry proteins present.
ERX463573 is the accession code of the experiment from which these raw read files came. The experiment this came from was about a Population of Bacillus thuringiensis which were coevolved with Caenorhabditis, which is a kind of nematode, as host. From the study we know to expect at the very least the following two proteins: Cry35Aa4 and Cry21Aa2.
According to the cry protein toxin list it is known that Cry35Aa4 is a binary toxin with Cry34Aa4 , and as such we may expect to find this protein, or one like it, as well.

Visualizing the result

The easiest way to visualize the results from the three separate components of the tool is to use a Venn diagram.

Figure 1: A Venn diagram showing the overlap of the output of the various methods used to predict whether or not certain genes from the genome are cry proteins or not.

As Blast had only 5 results it is easier to examine this method in detail. Two of the five were indeed proteins we expected from the paper: Cry35Aa4 and Cry21Aa2. Both of these were also picked up by the Hidden Markov Model method of cry protein detection as shown in figure 1, but not by the RandomForest method.
Next we found a protein that matched very well with the Cry34Aa group, which is a complimentary protein to Cry35Aa4.
The two other proteins which were found: Cry38Aa1, which has no known insecticidal target, but is highly similar to proteins that do: Cry15Aa1, Cry23Aa1, and Cry33Aa1. ¹⁷
and "Gene_5518" which was 85.87% identical to: Cry14Ab1. This protein is only mentioned in a patent by Sampson et al. 2012.¹⁸ This protein is quite interesting because it might be a completely new protein or a variation of Cry14Ab1 specific to nematodes.

Tool overview

We use a combination of existing tools to come to the prediction of novel cry proteins. The entire pipeline consists of four scripts in total, one of which is entirely dedicated to analysis of the Random Forest model and not further used in the main program. The others are there to group the Machine Learning specific functions, the functions that handle known cry proteins, and the functions that handle raw sequence data. Figure 2 gives a graphical representation of the pipeline, though some modules have been left out for the sake of readability. Here, we shall go through each part of the process in a step-by-step manner.

Figure 2: A graphical representation of the pipeline showing the various methods and tools used. All the pink diamonds are clickable and will take you to the respective tool's homepage. The known cry proteins box will take you to the cry protein database.

Software description

Click here for a highly detailed overview of how the pipeline works.

Genome assembly

Raw genome data directly from a Next Generation Sequencing device comes in the form of files containing many pieces of DNA called ‘reads’. These reads are usually obtained through cutting up the genome and sequencing the small bits at random after amplifying them many times. The idea is that many reads will overlap and can be assembled like a big puzzle. Our pipeline uses a de novo assembler, meaning that it can assemble the reads from scratch, without a reference genome to map them to. This is done by a 'graph-based approach' The nodes of the graph represent whole reads or parts of each read called k-mers, while the edges represent the overlap between them.. ³ Once the genome is assembled into contigs A contig is a closed loop in the graph, the longest stretches of DNA possible from all of the reads. we can start doing gene prediction. In this pipeline we used the IDBA-UD The Iterative De Bruijn Graph De Novo Assembler for highly Uneven sequencing Depth (IDBA-UD) is a tool for assembling genomes in a de novo manner . ¹⁵

Gene prediction

The pipeline does gene prediction through the use of Hidden Markov Models (HMM) ⁴ that have been trained to distinguish between prokaryotic coding and non-coding regions. These regions are then translated using BioPython's ⁶ translate function. These are the proteins we compare to our existing Cry Protein database ⁵ In this pipeline we used the GeneMarkS utilizes a non-supervised training procedure and can be used for a newly sequenced prokaryotic genome with no prior knowledge of any protein or rRNA genes. program for this purpose. ¹⁵

Machine Learning

As cry proteins are a group that is hard to define biologically, other than having specific insecticidal activity, it is somewhat of a challenge to find them through means other than direct sequence comparisons. This is why we use three methods to find potential cry proteins:

A direct sequence comparison with a stringent cut-off (E valueExpect value is the number of hits one can "expect" to see by chance when searching a database of a particular size. 0.1 or lower): BLAST The Basic Local Alignment Search Tool (BLAST) is an algorithm for comparing sequence information, such as amino-acid (protein) or nucleic-acid (DNA/RNA) sequences. ⁸
Using a database of only cry proteins, this tool will allow us to rapidly identify other proteins with highly similar sequences.
A probabilistic approach based on the primary subdivision between cry proteins (45% sequence similarity): HMMER.
The HMMER toolkit contains several programs that the pipeline uses. ⁹ The program first builds many different profiles based on the multiple sequence alignment of the primary cry protein groups. This is handled by Clustal Omega. ¹⁰
- HMMbuild: Builds a profile HMM from an input multiple alignment.
- HMMpress: Formats an HMM database into a binary format for hmmscan.
- HMMscan: Searches a protein sequence against a protein profile HMM database.
All of these methods are still based on the same Hidden Markov Model principle ⁴.
A machine learning approach based on protein features using many weak classifiers to give a majority vote: "RandomForest".⁷ This part is explained further in the protein classification section.

Protein Classification

The protein classification part of the machine learning is done by the pipeline in two steps: First the features of a protein are predicted using Biopython ⁶ , and then the pre-trained model predicts a class based on those features.

Features

After some experimentation with the features chosen for modelling, we came up with the following list which we think will capture essential structural elements and allow for good predictions.

Molecular weight of the protein.
Aromaticity, ie the relative frequency of aromatic amino acids (Phenylalanine, Tryptophan and Tyrosine)
Instability index, ie the relative frequency of instability causing dipeptides, for more information see Guruprasad et al 1990.¹¹
Isoelectric point, which is the pH at which a particular molecule carries no net electrical charge.
Alpha Helix fraction, the relative frequency of amino acids that tend to be in found in the secondary structure known as alpha helices.
Turn fraction, the relative frequency of amino acids that tend to be in found in the secondary structure known as turns.
Beta Sheet fraction, the relative frequency of amino acids that tend to be in found in the secondary structure known as beta sheets.

Scikit-Learn

We made the classification model using scikit-learn or sklearn.¹² The class used for this model is called the RandomForestClassifier which works by fitting a number of decision trees, in our case 2000, on various sub-samples of the dataset and averaging these weak classifiers to improve accuracy and control for. overfitting Overfitting is what happens when a model follows its training data so well it can not predict anything new with any accuracy because too much of the noise from the training data has been modelled, rather than the underlying relationship. We assess the predictive accuracy of our model by doing a stratified 10-fold crossvalidation technique from the StratifiedKFold class. These 10-fold cross-validation runs were assessed with the roc-curves method. A ROC curve is a plot of the true positive rate (TPR) The proportion of correctly predicted instances of the testing data. against the false positive rate (FPR) The proportion of wrongly predicted instances of the testing data. ¹⁹ These kinds of plots are extremely useful in determining whether or not a Machine Learning model such as RandomForest is well trained.

Venn-Diagrams

All three separate methods give a different output. In order to make an easy visual comparison, we make use of a web tool that can easily calculate and draw custom Venn-diagrams. This tool was made by the university of Gent and can be found here.

References

1. Ye, W., Zhu, L., Liu, Y., Crickmore, N., Peng, D., Ruan, L., & Sun, M. (2012). Mining new crystal protein genes from Bacillus thuringiensis on the basis of mixed plasmid-enriched genome sequencing and a computational pipeline. Applied and environmental microbiology, 78(14), 4795-4801. ↩

2. Alquisira-Ramírez, E. V., Paredes-Gonzalez, J. R., Hernández-Velázquez, V. M., Ramírez-Trujillo, J. A., & Peña-Chora, G. (2014). In vitro susceptibility of Varroa destructor and Apis mellifera to native strains of Bacillus thuringiensis. Apidologie, 45(6), 707-718. ↩

3. Compeau, P. E., Pevzner, P. A., & Tesler, G. (2011). How to apply de Bruijn graphs to genome assembly. Nature biotechnology, 29(11), 987-991. ↩

4. Rabiner, L., & Juang, B. (1986). An introduction to hidden Markov models. ieee assp magazine, 3(1), 4-16. ↩

5. Crickmore, N., Baum, J., Bravo, A., Lereclus, D., Narva, K., Sampson, K., Schnepf, E., Sun, M. and Zeigler, D.R. " Bacillus thuringiensis toxin nomenclature" (2016) http://www.btnomenclature.info/ ↩

6. Cock PA, Antao T, Chang JT, Bradman BA, Cox CJ, Dalke A, Friedberg I, Hamelryck T, Kauff F, Wilczynski B and de Hoon MJL (2009) Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25, 1422-1423 ↩

7. Liaw, A., & Wiener, M. (2002). Classification and regression by randomForest. R news, 2(3), 18-22. ↩

8. Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Basic local alignment search tool. Journal of molecular biology, 215(3), 403-410. ↩

9. Eddy, S. R. (1998). Profile hidden Markov models. Bioinformatics, 14(9), 755-763. ↩

10. Sievers, F., & Higgins, D. G. (2014). Clustal Omega, accurate alignment of very large numbers of sequences. Multiple sequence alignment methods, 105-116. ↩

11. Guruprasad, K., Reddy, B. B., & Pandit, M. W. (1990). Correlation between stability of a protein and its dipeptide composition: a novel approach for predicting in vivo stability of a protein from its primary sequence. Protein engineering, 4(2), 155-161. ↩

12. Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., ... & Vanderplas, J. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12(Oct), 2825-2830. ↩

13. Masri, L., Branca, A., Sheppard, A. E., Papkou, A., Laehnemann, D., Guenther, P. S., ... & Brzuszkiewicz, E. (2015). Host–pathogen coevolution: the selective advantage of Bacillus thuringiensis virulence and its cry toxin genes. PLoS Biol, 13(6), e1002169. ↩

14.de Maagd, R. A., Bravo, A., & Crickmore, N. (2001). How Bacillus thuringiensis has evolved specific toxins to colonize the insect world. TRENDS in Genetics, 17(4), 193-199. ↩

15. Peng, Y., Leung, H. C., Yiu, S. M., & Chin, F. Y. (2012). IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics, 28(11), 1420-1428. ↩

16. Besemer, J., Lomsadze, A., & Borodovsky, M. (2001). GeneMarkS: a self-training method for prediction of gene starts in microbial genomes. Implications for finding sequence motifs in regulatory regions. Nucleic acids research, 29(12), 2607-2618. ↩

17. Baum, J. A., Chu, C. R., Rupar, M., Brown, G. R., Donovan, W. P., Huesing, J. E., ... & Vaughn, T. (2004). Binary toxins from Bacillus thuringiensis active against the western corn rootworm, Diabrotica virgifera virgifera LeConte. Applied and environmental microbiology, 70(8), 4889-4898. ↩

18. Sampson, K. S., Tomso, D. J., & Dumitru, R. V. (2012). U.S. Patent No. 8,318,900. Washington, DC: U.S. Patent and Trademark Office. ↩

19. Bradley, A. P. (1997). The use of the area under the ROC curve in the evaluation of machine learning algorithms. Pattern recognition, 30(7), 1145-1159. ↩

Team:Wageningen UR/Software

Toxin Scanner

Toxin Scanner

BioBrick discovery

Toxin Specificity

Testing the tool

Visualizing the result

Tool overview

Software description

Genome assembly

Gene prediction

Machine Learning

Protein Classification

Features

Scikit-Learn

Venn-Diagrams

References