In Silico Analysis

Unsupervised Machine Learning

The goal to reduce the variability of the RBS BBa_B0034 via redesigning the spacer region would mean that our team should locate spacer mutants for which we could have quantitative results regarding their context dependency. That is no easy task and would normally require a large amount of experimental measurements with different spacer variants and genes of interest in order to create a large dataset. Of course, during the iGEM short time window, such a procedure would be impossible in parallel with the rest of the project experiments. As a result, our team decided to use an already tested-in-lab dataset and the different features of its input sequences in order to determine the context dependency of our novel RBS sequences.

The dataset was provided in the research paper by Duan et. al.[5]. In particular, they characterized a library of RBS mutants in the context of GFP and RFP fluorescent proteins. In this way, they created a labeled dataset of RBS sequences in which RBS variants with similar values of GFP and RFP fluorescence(< 20% difference) would be characterized as low context dependent (LCD) and the others would be labeled as high context dependent(HCD).

Simultaneously they also performed a feature analysis in order to determine the exact factors that play a more significant role in the context dependency of their algorithm and put an emphasis on the nucleotide sequence not only of the entire RBS sequence,but also of the RBS - Start Codon spacer region in particular, as well as characteristics of the secondary structure, such as MFE, and the binding affinity between RBS and 16S rRNA.After correlation analysis, they concluded to a 27-long vector of RBS characteristics which showed important statistical correlation regarding protein translation rates (via the Mann-Whitney Test [6]).

In the above workflow, they created a dataset that was unbalanced, meaning that the more context independent sequences were significantly fewer than the ones showing higher context independency (1:8 ratio). That type of dataset is a prime candidate for unsupervised learning models, which are capable of detecting outlier candidates. That is due to the fact that the given training model is already unbalanced, meaning that it is tuned to only detect normal cases. In the situation where the model is presented with a novel anomaly case (an LCD RBS in our case), it will create an anomaly score leading to the detection of it.

In terms of our workflow, we used the already characterized dataset provided by Duan et. al [5] and trained it for many different unsupervised learning algorithms via the pycaret library [7] . The only main difference is the binding affinity metric which we noticed that it couldn’t be replicated by any metric we tested on the Vienna RNA package. So, instead we used the pf_fold() function and in particular the ensemble energy variable it gives as an output which is proportional to the binding affinity metric used in the paper[8]. This was done not only for our own testing sequences but also for the original dataset ones.

The anomaly detection module of pycaret ranks trained and tested sequences based on an anomaly metric from 0 to 1, with the anomaly LCD values being closer to an anomaly score of 1. The 4 learning models used for our analysis via the pycaret library are briefly presented below.

1. Principal Component Analysis (PCA)

PCA is a popular technique for analyzing large datasets containing a high number of dimensions/features per observation, increasing the interpretability of data while preserving the maximum amount of information, and enabling the visualization of multidimensional data [9]. Formally, PCA is a statistical technique for reducing the dimensionality of a dataset.

In terms of using it in unsupervised learning, PCA can detect outlier inputs via calculating the euclidean distance between the novel input vector and the reduced space of the dataset provided by PCA. A high distance means a high reconstruction error and thus an anomaly score closer to one.

2. One Class Support Vector Machine (One-class SVM)

One-class SVM is a special case of SVMs, which are normally supervised learning models with associated learning algorithms that analyze data for classification and regression analysis. In brief, SVM maps training examples to points in space so as to maximize the width of the gap between the various categories. New examples are then mapped into that same space and predicted to belong to a category based on which side of the gap they fall. The kernel separating the different labels could be either linear or non-linear, as shown below.

As for one-class SVM, it is a variance of the original SVM algorithm which has as a target to find the hyperplane is as close to the original dataset points as possible [12]. As a result, when a novel point is presented that is away from the hyperplane, then one-class SVM detects it as an anomaly

3. K-nearest neighbors (knn)

The knn algorithm is a supervised learning classifier, which uses proximity to make classifications or predictions about the grouping of an individual data point The algorithm computes the euclidean distance of a novel point in order to compute its k-nearest neighbors and then use the neighbors labels in order to determine the new points’ label.

4. Isolation Forest Algorithm (IForest)

Isolation forest or IForest is an anomaly detection algorithm that uses isolation (how far a data point is to the rest of the data), rather than modeling the normal points Isolation forests are based on decision trees; randomly sub-sampled data is processed in a tree structure based on randomly selected features. he samples that travel deeper into the tree are less likely to be anomalies as they require more cuts to isolate them. Similarly, the samples which end up in shorter branches indicate anomalies as it was easier for the tree to separate them from other observations [14].

After selecting the algorithms described above, then we trained the initial dataset for each one of those algorithms and saved the models that were produced on the .pkl files. Simultaneously we created our own muted spacer sequences on the RBS_BBa0034 variants-100 with a spacer length 5 and 100 more with a spacer length 6- and tested the 4 different trained models in order to rank the novel RBS variants in terms of “anomaly”. Then kept the top 20 variants from each model, i.e. the 20 sequences with the highest anomaly scores and then kept the variables that were present on all 4 “top 20” variant lists.

The resulted sequences from the in silico analysis are presented on table 1

Due to a lack of time for testing all 5 sequences, the wet lab team initially arbitrarily decided on testing 3 variants. Since all variants have the same AC content and one T or C in the spacer region, the T containing variants were selected on the grounds of the rough estimation that an AT hydrogen bond with a neighboring sequence would offer less stability than a stronger GC bond. Hence, we were led to choose RBS variants with index values 2,3,5 on table 1 for experimental validation of our in silico analysis. In the end, due to time and wet lab engineering constraints, only the RBS variant with index value 5 (V59) was tested in comparison with the original BBa_B0034 RBS.

For a better view of the generated scripts and analysis, please refer to our improve-a-part repository on gitlab