Team:Wageningen UR/Software

MENU ▤
Wageningen UR iGEM 2016

 

Introduction

iGEM and sequencing

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

Background

For this project we needed a toxin specific to Varroa destructor. There are many insecticidal protein families, but the one we were most interested in are the "crystal" (or 'cry') proteins. These are usually found on plasmids found in Bacillus Thuringiensis and related species. 14 We know varroa-specific insecticidal activity exist in Bacillus thuringiensis and related species, due to a paper called: "In vitro susceptibility of Varroa destructor and Apis mellifera to native strains of Bacillus thuringiensis." by Alquisira-Ramírez et al. 2 In this paper several isolates are described that cause a mite mortality of up to 100% mite mortality. 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 focusses on creating V. destructor-binding cry toxins and finding V. destructor-specific insecticidal proteins.
This tool was made in preparation for results of the latter, finding V. destructor-specific insecticidal proteins, which we assume to be of the cry protein family. 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 not a cry protein 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.

Current developments

There already exists publication about a tool called “Bt Toxin Scanner”1, which inspired us to make our own tool. Due to the tool not being fully supported for local deployment, which is needed for high-throughput analysis, and the relatively basic analysis done by the tool, we decided to develop our own tool that is fully open-source and improves upon the 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.

Methods

Overview

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 1 gives a graphical representation of the pipeline, though some modules have been left out for readability's sake.

Figure 1: A graphical representation of the pipeline showing the various methods and tools used.

Stepwise description

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 pieces will overlap and so the reads 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 to map them to. It does this using a 'graph based approach.'3 Once the genome is assembled into 'contigs' we can start doing gene prediction.

Gene prediction

The pipeline does gene prediction through using Hidden Markov Models (HMM) 4 that have been trained to distinguish between prokaryotic coding and non-coding regions. These regions we then translate using BioPython's 6 translate function. These are the proteins we compare to our existing Cry Protein database 5

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:

  1. A direct sequence comparison with a stringent cut-off: "Blast".
  2. A probabilistic approach based on the primary subdivision between cry proteins (45% sequence similarity): "HMM". 4
  3. A machine learning approach based on protein features using many weak classifiers to give a majority vote: "RandomForest".7

Output

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

Tools Used

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. This means that it can take raw sequencing data and use it to create longer stretches of overlapping DNA called contigs.

Gene Mark

GeneMark is gene prediction software, making use of HMM. This will allow us to analyze novel genomic sequences.

Blast

The Basic Local Alignment Search Tool (BLAST) is an algorithm for comparing sequence information, such as the amino-acid (protein) or nucleic-acid (DNA/RNA) sequences. 8 Using a database of only cry proteins, this tool will allow us to rapidly identify other proteins with highly similar sequences.

HMMER

The HMMER toolkit contains several programs that the pipeline uses. 9 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. 10

  • 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.

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 6 , 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 some 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.11
  • 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.12 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 to improve accuracy and control overfitting. We assess the predictive accuracy of our model by doing a stratified 10-fold crossvalidation technique from the StratifiedKFold class. These 10 crossvalidation runs were assessed with the ROC-curves method.

Results / 'validation'

PRJEB5931

We tested the pipeline on a genome sample from a study with accession number: PRJEB5931. 13 This genome was found after a co-evolution experiment and a bt with known nematicidal cry proteins present.

Key result

The easiest way to visualize what the three seperate components of the pipeline have done is to compare their results in a venn diagram (with accompanying lists).

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

Discussion

It seems to work, but the definition of cry protein is vague, so the ML could be further optimized based on a more rigorous examination of what is and what is not a cry protein.

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. ↩

@. 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. ↩

@. ↩

@. ↩

@. ↩

@. ↩

@. ↩ Stuff I'm not sure what to do with:

Side features

has overwrite protection can be multi threaded creates its own directory structure Has tools to do detailed analysis of the random forest/machine learning model created(plotting decision boundaries, interpreting specific predictions) The tool can also be modified to search and compare other sets of genes against raw genome sequences. This would be done by providing different input and modifying some of the internal code handling filenames.

Exploring cry protein features

Creating sub divisions based on clustering, PCA and kmeans clustering show that the only feature where to be made on the simple features. Other than the proteins clustering by protein length/molecular weight, there do not seem to be clear clusters to be made based on the basic features we chose for identifying cry toxins.