Objectives

Our aim includes the development of computational model and software that will enable us to design plasmids with required copy number via the manipulation plasmid’s origin of replication (Ori).

By using our model, we can understand how copy number is encoded in the origin of replication; we can use this knowledge to predict the copy number from Ori sequence and design new Oris. This will allow synthetic biologists to tune the expression of genes of interest for various objectives such as alternative proteins for food tech, biosensors, metabolic engineering, and DNA vaccines, among others via optimizing protein production with minimal metabolic burden on the cell.

A gradient of plasmid copy number can be achieved by controlling the interaction between priming RNA and inhibitory RNA. To maximize copy number, RNAp should have a good affinity to the DNA strand, and RNAi should be repressed. The amount of RNAp and RNAi in the host cell is controlled by their promoter sequences, located in the plasmid’s ORI (see more details in the Project Description page).

Model Features

We extracted thousands of sequence-related features for RNAp, RNAi promoters including promoter binding energy, promoter strength, PSSM, nucleotide counts, one hot encoding, motifs, features based on the secondary structure of RNA (see more details in the RNA secondary structure features section) and more. Below we briefly review some of these features.

Motifs

Motifs files were downloaded from (1) DPINTERACT database [3] which is a comprehensive library of binding sites for E. coli DNA-binding proteins containing a set of 68 motifs, between 10 and 50 in width (average width 24.5), and (2) SwissRegulon E.coli database [4] which is a database containing genome-wide annotations of regulatory motifs, promoters and TF binding sites (TFBSs) in promoter regions across model organisms, contains a set of 97 motifs, between 6 and 30 in width (average width 20.1).

The score and P-value of each motif were calculated with FIMO module, and sequences with no matching motif were replaced with 0 score, and p-value of 1.

PSSM

The promoter sequences of priming RNA and inhibitory RNA corresponding to the top and bottom 20% copy number values were used to compute a position-specific scoring matrix (PSSM), a commonly used representation of motifs (patterns) in biological sequences, for differentiating between high and low copy number plasmids. The PSSM score of each sequence is calculated based on the 20% high copy number matrix.
The PSSM results are illustrated in the figures:

Promoter strength

The promoter strength was calculated using an energy matrix [KBT] taken from an article by Brewster et al. [5]. The matrix covers base pairs [-41:-1] where 1 denotes the transcription start site. Each row corresponds to a given nucleotide and each column corresponds to a given position. The promoter strength was calculated for the entire sequence and for the two edited zones of the priming and inhibitory RNAs.

One hot encoding features

each promoter sequence was split into 144 (=36*4) columns so that each column represents the index and nucleotide combination, and the value is set to 1 if it matches the current sequence and 0 otherwise.

Nucleotide counts features

We used nucleotide counts from an article by Muhammod et al. [6] including:

Table 2: Nucleotide counts features and their description.

Feature Name	Number of features	Feature Description
zCurve	3 features	x_axis =(ΣA+ΣG)-(ΣC+ΣT) y_axis =(ΣA+ΣC)-(ΣG+ΣT) z_axis =(ΣA+ΣT)-(ΣG+ΣC)
gcContent	1 feature	(ΣG+ΣC)/(ΣA+ΣC+ΣG+ΣT) * 100%
ATGC ratio	1 feature	(ΣA+ΣT)/(ΣG+ΣC)
Cumulative Skew	2 features	GC Skew=(ΣG-ΣC)/(ΣG+ΣC) AT Skew=(ΣA-ΣT)/(ΣA+ΣT)
Pseudo KNC	when –kGap = 2, 20 features	Features will be numbers of A, C, G, T, AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, and TT of the whole sequence of DNA respectively.
monoMonoKGap	when –kGap = 2, 32 features	Features will be numbers of A_A, A_C, A_G, A_T, C_A, C_C, C_G, C_T, G_A, G_C, G_G, G_T, T_A, T_C, T_G, T_T, A__A, A__C, A__G, A__T, C__A, C__C, C__G, C__T, G__A, G__C, G__G, G__T, T__A, T__C, T__G, T__T of the whole sequence of DNA respectively.
monoDiKGap	when –kGap = 2, 128 features	Feature structure will be X_XX, and X__XX of the whole sequence of DNA respectively.
monoTriKGap	when –kGap = 2, 512 features	Feature structure will be X_XXX, and X__XXX of the whole sequence of DNA respectively.
diMonoKGap	when –kGap = 2, 128 features	Feature structure will be XX_X, and XX__X of the whole sequence of DNA respectively.
diDiKGap	when –kGap = 2, 512 features	Feature structure will be XX_XX, and XX__XX of the whole sequence of DNA respectively.
diTriKGap	when –kGap = 2, 2048 features	Feature structure will be XX_XXX, and XX__XXX of the whole sequence of DNA respectively.
triMonoKGap	when –kGap = 2, 512 features	Feature structure will be XXX_X, and XXX__X of the whole sequence of DNA respectively.
triDiKGap	when –kGap = 2, 2048 features	Feature structure will be XXX_XX, and XXX__XX of the whole sequence of DNA respectively.

dG total and TX rate

ΔGtotal is the difference in free energy between an unbound promoter and a promoter, RNA polymerase and the subunit of RNA polymerase σ70 factor complex with a stable transcriptional bubble, used to predict a promoter’s TX rate compared to a reference promoter sequence with calculated ΔGtotal,ref and measured TXref. ΔGtotal,ref was set to zero arbitrarily so that stronger and more favourable interaction free energies have more negative values compared with the reference promoter, which has a low transcription rate [7].

Entropy

Entropy measures the variability of nucleotides at a given position. It’s a way to assess how conserved or variable a position is among the sequences. In our case the data is promoter sequences and it’s relevant for the positions where researchers inserted mutations.
First based on our sequences (all with the same length) we get the probability of each nucleotide to appear in each index (per each index we calculate the ratio: \(\dfrac{\text{count of nucleotide}N}{\text{number of sequences}}\) so that we get per each index 4 ratios, one for each nucleotide). We use these probabilities to calculate the entropy per sequence as follows:

$$S = -\frac{\sum_{i=0}^{n} p_i \log_2(p_i)}{2}$$

Where \(p_i\) is the probability we calculated before of the nucleotide at index i. We divide the value by 2 for normalization purposes [8].

De Novo Motif Detection

We wanted to search for de novo motifs that appear at higher rates in groups of promoters that have high copy number against a reference group of promoters with low copy number, and vice versa, and use those motifs’ scores as features for a different set of promoters to avoid overfitting.

To do so, we have used HOMER (Hypergeometric Optimization of Motif EnRichment), which is a suite of tools for Motif Discovery and next-gen sequencing analysis. It is a collection of command-line programs. HOMER was primarily written as a de novo motif discovery algorithm and is well suited for finding 8-20 bp motifs in large-scale genomics data [9].

Model selection

We had 3 groups of candidate models - Neural network, Linear Regression models, Random Forest and Ensemble Learning models. We wanted to choose the model that answered several criteria: model performance, model complexity (we want the model to be as explainable as possible) and preventing overfitting.

Neural Networks

A Neural Network is a machine-learning technique that is composed of layers of interconnected processing units, with input, hidden, and output layers. Each unit is connected to the units in adjacent layers, with each connection having a weight value. The weights are estimated with a backpropagation process, where they are adjusted in order to minimize the average classification or regression error. Inputs are multiplied by these weights and summed at each unit, and then passed through an activation function, such as a sigmoid, tanh, or rectified linear unit (ReLU), before being passed to the next layer [10].

Neural Networks excel in performance for complex tasks, but they are highly complex and can easily overfit if not regulated. They are challenging to interpret due to their "black-box" nature.

Linear Regression models

Linear regression is widely used in predictive analysis methods when prediction results can be quantified and modeled using multiple input variables. It is a method of identifying linear correlations between dependent and independent variables through analyzing and modeling data. Therefore, this method is capable of modeling relationships between dependent variables and independent variables by analyzing and learning from existing training data [11].

Linear regression models have low to moderate performance but are straightforward and easy to understand. They have low complexity, low risk of overfitting, and high interpretability via coefficient analysis.

Lasso

Lasso (Least Absolute Shrinkage and Selection Operator) is a linear regression model that performs feature selection and regularization, making it a powerful tool for predicting outcomes in high-dimensional datasets.
The lasso model uses L1 regularization, which adds a penalty term equal to the absolute value of the coefficients of the regression parameters in order to shrink the coefficients of less important features to zero. This results in a more parsimonious model that only includes the most relevant features in the final prediction while also improving the model's performance by reducing overfitting.

Ridge regression

Ridge regression is a method of estimating the coefficients of multiple-regression models in scenarios where the independent variables are highly correlated [12]. It is especially useful for mitigating multicollinearity in linear regression, which commonly occurs in models with large numbers of parameters [13]. Generally, the method provides improved efficiency in parameter estimation problems in exchange for a tolerable amount of bias [14].

Elastic Net

The elastic net is a regularization and variable selection method. The elastic net is particularly useful when the number of predictors is much bigger than the number of observations [15]. The elastic net is a regularization method that linearly combines L1 and L2 penalties, which are used in LASSO and ridge regression methods respectively [16].

Random Forest model

The random forest algorithm works by constructing a multitude of decision trees on different subsets of the training data. Each decision tree is constructed using a random subset of the features, and at each split in the tree, a random subset of the features is considered. This process of feature randomness helps to reduce overfitting and improve the generalization performance of the model.

Once the individual decision trees have been constructed, the random forest algorithm makes predictions by averaging the predictions of the individual trees. In the case of regression tasks, the average of the individual predictions is used to determine the final prediction [17].

Random forests offer a good balance between performance and complexity. They are moderately complex and less prone to overfitting. Their interpretability is moderate through feature importance scores.

Ensemble Learning models

Ensemble Learning is a technique that uses the “Wisdom of the crown” principle and builds a strong predictive model by combining the predictions of multiple weak models, typically decision trees, in a sequential manner.
There are two main classes of ensemble learning methods:

• Bagging method builds models in parallel using a random subset of data (sampling with replacement).
• Boosting method builds “weak” models in sequence using the whole data, with each model improving on the previous model’s error.

Gradient Boosting Trees, CatBoost, and XGBoost are all variations of gradient boosting algorithms. All of them build decision trees, which are trained sequentially and use the error from the previous tree to adjust its learning and eventually minimize the loss function. Gradient boosting is primarily used to reduce the bias error of the model. Based on the bias-variance tradeoff, it is a greedy algorithm that can overfit a training dataset quickly. However, this overfitting can be controlled by shrinkage (learning rate), tree constraint (pruning), and regularization, which are implemented in CatBoost and XGBoost models [18].

Ensemble learning combines models for high performance. They generally have moderate complexity, a moderate risk of overfitting, and moderate interpretability through feature importance.

Gradient Boosting Trees

In each iteration, a tree is built using all the features, based on the residual errors of the previous tree. Each tree in the ensemble is shrunk after it is multiplied by the learning rate (ƞ) which ranges between 0 to 1.
The final prediction is given by the formula:
y(pred) = y1 + ( ƞ* r1) + (ƞ * r2) + ....... + (ƞ * rN)
where y1 is the prediction of the first tree, r1 is the prediction of the second tree that was trained based on the residual of the first one, and so on. Below is an illustration of ensemble models. This process is illustrated in Figure 8.

XGBoost

Based on gradient boosting trees and additionally contains regularization term and tree pruning to prevent overfitting.

Cat-Boost

Based on gradient boosting trees. Adapts the learning rate during training and incorporates built-in cross-validation and other techniques to mitigate overfitting.

Table 3: A summary of the models’ properties.

	Neural network	Linear regression models	Ensemble learning models	Random Forest
Model performance	High	Low to moderate	Moderate to high	Moderate to high
Model complexity	High	Low	Moderate to High	Moderate
Model overfitting	High	Low	Moderate to low	Moderate
Model interpretability	Low	High	Moderate	Moderate

“The best approach […] is to randomly divide the dataset into three parts: a training set, a validation set, and a test set. The training set is used to fit the models; the validation set is used to estimate prediction error for model selection; the test set is used for assessment of the generalization error of the final chosen model.” [19]

After splitting our data, in order to find the best model for our data, we performed a hyperparameter tuning for each model using “Optuna” [20] (see more details in the hyperparameters tuning section), and evaluated the R-squared and MAE of the train and validation sets. This approach allowed us to examine the model performance and model overfitting.

These are the results for the RNAp model:

These are the results for the RNAi model:

In light of these results, we decided to proceed to the next steps, feature selection and hyperparameter optimization, with the following models: XGBoost, Cat-Boost, Light GBM, and Random Forest.

Feature selection

In our research workflow, feature selection has been a crucial step in enhancing the effectiveness of our machine learning models. Initially, we calculated the mutual information (MI) between each feature and the copy number, selecting those with MI values greater than the mean MI as potential candidates. Subsequently, we assessed the MI between the remaining features, identifying pairs of features exhibiting strong correlations within the top percentile. For each pair of correlated features, we chose the one displaying higher MI with the copy number, ensuring the retention of the most informative attributes.

Finally, to further refine our feature set and optimize model interpretability, we applied the eBoruta approach.

Boruta is an algorithm for selecting relevant features in a dataset, primarily designed for use with Random Forests [22]. It falls under the category of "wrapper" methods, as it employs a machine learning model to filter out features that are pertinent to the model's learning objective.

Boruta-SHAP is a feature selection algorithm that combines the Boruta algorithm with SHAP values to select important features.

Boruta is a tree-based feature selection algorithm that works by creating random shadow features (noise) and comparing the importance of the original features to the importance of the shadow features. Features that are more important than the shadow features are kept, while features that are not more important than the shadow features are removed.

SHAP [23] values are a way to measure the individual impact of each feature on the prediction of a model. SHAP values are calculated by considering all possible combinations of features and their impact on the prediction.

Using e-boruta with SHAP values works by using SHAP values to rank the importance of the features. The features with the highest SHAP values are then selected using the Boruta algorithm.

eBoruta-SHAP has several advantages over other feature selection algorithms, including:

• It is robust to noise and outliers.
• It can handle high-dimensional data.
• It can select features that are important for both linear and non-linear models.
• It can select features that are important for both classification and regression tasks.

By integrating this feature selection process, we obtained a comprehensive understanding of feature importance. This final selection of features, chosen based on their impact on model predictions, allowed us to train a more efficient and interpretable machine learning model, thereby advancing the quality and relevance of our research outcomes.

Features that survive this feature selection steps are considered the most important for the model's performance and are retained for training the final machine learning model.

Hyperparameters tuning using Optuna

Optuna is a powerful and versatile Python library for hyperparameter tuning, which has gained significant attention in the field of machine learning and optimization. Optuna automates the search for optimal hyperparameter configurations using Bayesian optimization, where each trial employs a different set of hyperparameter values and evaluates the R-squared score on the validation set [20]. This way, the best set of hyperparameters for our data and for each model is selected.

Here is a glimpse of our hyperparameter optimization study.

RNA secondary structure features

The model based on inhibitory RNA promoters produced less successful results than the model based on priming RNA promoters. Therefore, we decided to add more features that could help improve our model. The inhibitory RNA (RNAi or RNA I) and its promoter are complementary to a part of the priming RNA (RNAp or RNA II) (figure 21). Thus, a mutation in the promoter of RNAi will cause a mutation in RNAp itself and could alter its secondary structure. As a result, we planned to generate features that describe RNAp secondary structure and apply them to the model based on the inhibitory RNA promoters.

The replication of the ColE1 plasmid family, and pUC19 as part of it, is a highly regulated process, coordinated by the interplay of two pivotal RNA molecules: RNAi, transcribed from the antisense promoter P1, and RNAp, known as the pre-primer RNA, transcribed from the sense promoter P2 [24]. This regulation hinges on a remarkable phenomenon — the persistence of RNAp on the DNA template — an essential prerequisite for initiating replication.
This persistence is made possible through a specific interaction: a G-rich loop within the RNAp transcript, situated 265 nucleotides upstream of the origin, forms a crucial partnership with a C-rich stretch on the DNA template strand, positioned just 20 nucleotides upstream of the origin [25]. A balance between two distinct RNAp conformations plays a pivotal role in determining the percentage of transcripts capable of forming the persistent hybrid necessary for DNA synthesis initiation [25].

RNAp beta (𝞫) sequences can be paired either with a complementary upstream alpha (𝞪) sequence or with an equivalent downstream gamma (𝞬) sequence. After approximately 120 nucleotides have been transcribed by RNA polymerase, the nascent RNAp transcript folds into a secondary structure composed of three stem-loop domains I, II and III. However, this conformation is only temporary, as the polymerization of several more nucleotides disrupts stem loop III and creates an extended stem loop, IV, that is stabilized by hydrogen bonding between alpha and beta sequences. As transcription progresses, two structures are possible. Changing from alpha-beta to beta-gamma pairing perturbs the RNAp conformation, determining whether the transcript is competent for RNA-DNA hybridization and primer formation. RNA molecules in the alpha-beta conformation form the persistent hybrid necessary for DNA synthesis initiation, while transcripts in the beta-gamma conformation do not [25].

RNAi conformation consists of three stem loops and an unpaired 5' region nine nucleotides long. RNAi binds to RNAp by complementary base pairing. The interaction between RNAi and RNAp affects the folding of RNAp, such that beta-gamma pairing is favored. In vitro, RNAi affects primer formation only when added within a short period of time after initiation of RNAp transcription [25].

We used the sequence of RNAi and RNAp as it was found in pUC 19 plasmid sequence, that was obtained from Addgene (Plasmid #50005) [26] and the mutations of RNAi promoter as described in the article by Rouches et al [2]. The variants were edited at 9 positions (-33,-30), (-10,-7), and +1 in RNAi promoter which corresponds to positions 112, (118,121), (142,145) in RNAp sequence. It should be noted that the promoter of RNAi that is complementary to RNAp as it was found in PUC 19 is included in the data which contains the RNAi promoter variants and the suitable copy number.

We used the python package ViennaRNA [27] which includes the tools RNAfold and RNAeval. RNAfold calculates the minimum free energy RNA’s secondary structure, MFE secondary structure, and the centroid secondary structure, which is the secondary structure with minimal base pair distance to all other secondary structures in the Boltzmann ensemble. RNAfold also calculates the probabilities of each base to bond with another base along the RNA sequence; these probabilities are presented as colors on the RNA secondary structure image. RNAeval calculates the free energy of a given secondary structure, by each base pair bond’s free energy. Our aim is to generate features that describe RNAp secondary structure and characterize the alpha, beta and gamma sequences, based on the mutations in RNAi promoter.

We calculated the RNA secondary structure of RNAp as it was found in PUC 19 plasmid using RNAfold and we got 2 secondary structures as discussed. We calculate the RNA secondary structure of RNAp at several lengths, 108 bases as the length of RNAi, 200 bases to capture the alpha beta structure and 553 bases as the full length of RNAp.

Below are some examples of the secondary structure of RNAp and RNAi for illustration and verification of the structures described in the literature.

The MFE secondary structure of RNAp as found in PUC 19 at full length, is shown on the right (figure 27), and the centroid secondary structure on the left (figure 26).

The MFE secondary structure of RNAp as found in PUC 19 at a length of 108 bases with annotations of stem loops 1,2 and 3, is shown on the right (figure 29). In addition, the centroid secondary structure is shown on the left (figure 28).

The MFE secondary structure of RNAp as found in PUC 19 at a length of 200 bases with annotations of alpha sequence, beta sequence, -10 hexamer mutations, -35 hexamer mutations and the +1 mutation, is shown on top (figure 30). In addition, the centroid secondary structure is shown below (figure 31).

It can be seen through comparison with the sequence of the ColE1 ori previously discussed that in both structures at a length of 200 bases the conformation of alpha-beta pairing and stem loop IV can be found. This structure is associated with primer formation and corresponds with the high copy number of the RNAi promoter that is complementary to this RNAp.

The MFE secondary structure of RNAi as found in PUC 19 with annotations of stem loops 1,2 and 3, is shown on the right (figure 33). In addition, the centroid secondary structure is shown on the left (figure 32).

The MFE secondary structure of RNAp corresponding to the RNAi promoter with the highest copy number at full length, is shown on the right (figure 35), and the centroid secondary structure on the left (figure 34).

The MFE secondary structure of RNAp corresponding to the RNAi promoter with the highest copy number at a length of 200 bases with annotations of alpha sequence, beta sequence, -10 hexamer mutations, -35 hexamer mutations and the +1 mutation, is shown on top (figure 36). In addition, the centroid secondary structure is shown below (figure 37).

We developed many features describing RNAp secondary structure based on RNAfold and RNAeval calculations. The features are calculated for a range of lengths of RNAp in order to describe its structure during transcription:

Table 4: Features describing RNAp secondary structure.

Feature Group Name	Number of Features	Feature Description
Base pairs match ratio	22 features	We calculate the following ratio:\(\dfrac{\text{number of paired bases in specific area}}{\text{number of paired bases in pUC19 sequence in this area}}\). We get the paired bases using RNAFold's output. Specific features examples: alpha_beta_match_seq_end_130 – it means we used the first 130 bases of the RNA sequences and checked alpha-beta area. alpha_beta_extended_match_seq_end_400 – it means we used the first 400 bases of the RNA sequences and checked extended alpha-beta area. Note: we used varied lengths of RNA because during its formation it changes shapes, and it can be interesting to see the difference in structures in different stages.
Base pairing probabilities	22,877 features	RNAFold can output in an additional file the probability of each nucleotide to pair with a nearby nucleotide. We parsed the file and took each such probability as a feature. Specific feature example: seq_end_400_375_385 – it means we used the first 400 bases of the RNA sequences and checked the probability of the nucleotide in index 375 to pair with nucleotide in index 385.
MFE and Centroid comparison	1,730 features	RNAFold can predict the secondary structure of a given RNA sequence using ‘MFE’ and ‘centroid’. We compared how similar the two predictions are. RNAFold’s documentation states that the more similar they are the prediction is more reliable, so we thought it could be a good feature as well. The features are whether base pairs in different indexes exist in both predicted structures. Specific feature example: seq_end_180 mfe == centroid 118_141 – it means that for each RNA sequence we used the first 180 bases as an input for RNAFold and got two predicted secondary structures: MFE and Centorid, for each base pair we check if it appears in both forms, in this case we check if the pair combined of nucleotides at indexes 118 and 141 appears in both forms.
Stem loops MFE	20 features	Another output we can extract is the minimal free energy (MFE) of specific areas in the secondary structure (using RNAFold and then RNAeval – another tool in ViennaRNA package). Specific feature example: mfe_seq_end_300_sl_IV – it means we used the first 300 bases of the RNA sequences and calculated MFE for stem loop IV.
Base pairing of specific areas	6 features	We checked how many bases are paired in specific areas in the predicted secondary structure of the RNA sequences. Specific feature example: base_pairing_mfe_alpha_range – it means we checked in the MFE predicted structure how many bases are paired in the alpha area.
RNA form representation	36 features	We tried to represent the form using some basic characteristics of the secondary structure output. We calculated: Number of unpaired bases Number of ‘arcs’ – a term we use for continuous sequence of unpaired bases that form an open ended circular shape. Number of ‘interior loops’ – a term we use for continuous sequences of unpaired bases that form an interior loop shape (similar to the tip of RNA hairpin) Specific features examples: unpaired_cnt_seq_end_82 - counts number of bases that don’t pair with other bases in the first 82 bases of the RNA sequences. arc_cnt_seq_end_102 – counts how many open ended circles we find in the first 102 bases of the RNA sequences. interior_loop_cnt_seq_end_172 – counts how many interior loops we find in the first 172 bases of the RNA sequences.
RNA form in Window	28 features	These features include: mfe, unpaired bases count, ‘arcs’ count and ‘interior loops’ count but in local windows of size 70. Since the bonds of RNA’s secondary structure are usually local it makes sense to check these features in small windows. Specific features examples: mfe_window_130_200 – calculates MFE of RNA between index 130 to 200. unpaired_cnt_window_120_190 – counts number of bases that don’t pair with other bases in the secondary structure of RNA between index 120 to 190. arc_cnt_window_110_180 – counts how many open ended circles we find in the window (here between index 110 to 180). interior_loop_cnt_window_100_170 – counts how many interior loops we find in the window (here between 100 to 170).

ColE1-like Plasmids Bioinformatics Analysis

Motivation

Based on our discussion with Alagene, we were encouraged to expand our tool to different plasmids and hosts that the scientific community can use. We have gathered ColE1-like plasmids and analyzed them against E. coli puc19 plasmid to detect evolutionarily conserved locations.

We have created a bioinformatics data analysis pipeline of ColE1-like plasmids to gain new insights about the central modulators of the plasmid copy number - RNAp, RNAi, and their corresponding. Specifically, as explained below, we clustered RNAp promoters based on their sequence similarity to find distinct promoter groups. We believe that each cluster represents a novel family of promoters that can regulate RNAp expression. We plan to test the effect of each such family in the lab. This information will improve our model in two dimensions:

1. Adding additional degrees of freedom to the modifications we can do to the RNAp promoter.
2. Finding promoters that can work in various additional hosts.

We have created a code pipeline that includes the following components (more details appear in the next sub-sections):

• Download ColE-like plasmids fasta files.
• Develop a null model for evaluating the quality of the alignment score of a sequence to our RNAi or RNAp sequence.
• Given a plasmid, we extract RNAi and RNAp above our null model threshold.
• We extract the promoters of RNAi and RNAp 40bp upstream of RNAi and RNAp.
• For each pair of promoters, we calculated a distance measure based on the chimera algorithm.
• We cluster the promoters using calculated chimera scores using different clustering algorithms.
• We extract interesting consensus sequences from different clusters using Multiple Sequence Alignment.
• We plan to test the consensus sequences in the lab.

Data

We have gathered ColE1-like plasmids from the paper “Evolution of ColE1-like plasmids across γ-Proteobacteria: From bacteriocin production to antimicrobial resistance” [28]. The initial dataset included more than 1,000 plasmids from 62 host species that share RNAi and RNAp replication system (see figures 51, 52 and table 5).

RNAI And RNAP Detection

Given a plasmid, we find RNAi and RNAp by local alignment to the puc19 sequences:

RNAI (forward) = AACAGTATTTGGTATCTGCGCTCTGCTGAAGCCAGTTACCTTCGGAAAAAGAGTTGGTAGCTCTTGATCCGGCAAACAAACCACCGCTGGTAGCGGTGGTTTTTTTGTTTGCAAGCAGCAGATTACGCGCAGAAAAAAAGGATCTCA

RNAP_seq (forward) = GCAAACAAAAAAACCACCGCTACCAGCGGTGGTTTGTTTGCCGGATCAAGAGCTACCAACTCTTTTTCCGAAGGTAACTGGCTTCAGCAGAGCGCAGATACCAAATACTGTTCTTCTAGTGTAGCCGTAGTTAGGCCACCACTTCAAGAACTCTGTAGCACCGCCTACATACCTCGCTCTGCTAATCCTGTTACCAGTGGCTGCTGCCAGTGGCGATAAGTCGTGTCTTACCGGGTTGGACTCAAGACGATAGTTACCGGATAAGGCGCAGCGGTCGGGCTGAACGGGGGGTTCGTGCACACAGCCCAGCTTGGAGCGAACGACCTACACCGAACTGAGATACCTACAGCGTGAGCTATGAGAAAGCGCCACGCTTCCCGAAGGGAGAAAGGCGGACAGGTATCCGGTAAGCGGCAGGGTCGGAACAGGAGAGCGCACGAGGGAGCTTCCAGGGGGAAACGCCTGGTATCTTTATAGTCCTGTCGGGTTTCGCCACCTCTGACTTGAGCGTCGATTTTTGTGATGCTCGTCAGGGGGGCGGAGCCTATGGAAA

We ensure that in our data, RNAi and RNAp are on different strands on each plasmid and that their alignment score is higher than our null model.

Null Model

To create a null model, we randomly sampled 1000 locations in the E. coli genome and performed alignments against RNAi and RNAP above. Our null model thresholds (the minimum thresholds for RNAP and RNAi) were established based on their maximum alignment scores in the random sampling.

As seen in the figures below, the minimum threshold for RNAi was set to 374, corresponding to the maximum alignment score between the forward and reversed RNAi sequences, while the threshold for RNAp was established at 1117.

Promoters Extraction

After we found RNAi and RNAp on different strands above the threshold, we extracted 40bp upstream to find their corresponding promoters.

Then, we created the following table with the promoters:

ChimeraARS Score

ChimeraARS (Chimera Average Repetitive Substring) [29] is a measure of the tendency of a genomic sequence to include long sub-sequences that appear in a reference sequence or a reference set of genomic sequences. It can replace classical alignment scores when the homology is not very high.

For each extracted promoter, we have calculated its corresponding chimeraARS score using the following algorithm:

1. For position i in the gene, look for the longest nucleotide substring which begins at position i and appears in the reference.
2. The substring’s length is this position’s score.
3. The sequence’s score (including start and stop codons) is the average position score.

The chimera score calculation used puc19 promoters as a reference. The advantages of chimera are that it is a fast and robust way to find similarities to a reference sequence.

Clustering

To achieve robust insights and characterize the generated promoters, we divide them into distinct groups and analyze each group individually.

We computed chimera scores for each promotor by comparing them to the other promotors and then employed them as features in our clustering algorithm. We have done this for both RNAi and RNAp.

Our clustering algorithm leverages a dataset comprising 800 features representing chimera scores, a crucial aspect in understanding the characteristics of promoters. To ensure the selection of the most effective model, we harnessed the power of seven state-of-the-art clustering algorithms as hyperparameters, employing a meticulous grid search parameter tuning algorithm. Throughout the evaluation process, each model underwent further optimization of its own hyperparameters through a dedicated grid search algorithm. The performance metric employed to gauge the quality of clustering results consisted of two key components: a Silhouette score greater than 0.4 combined with the minimum of Davies-Bouldin score divided by Calinski-Harabasz score. The algorithm systematically evaluated various clustering approaches, including DBSCAN, HDBSCAN, K-Means, DPCA (as introduced in [30], [31]), Affinity Propagation, Gaussian Mixture, and Spectral Clustering, to identify the most suitable clustering technique for the given dataset. To plot the clusters, a PCA algorithm was used. [32]

The figure above reveals distinct clusters formed by all four algorithms. Each color group in a cluster plot represents a group of observations that the algorithm clustered together based on the used metric. The top-performing algorithms were HDB and SP, achieving scores of 0.32 and 0.52, respectively, with cluster numbers of 17 and 46.

Clusters Analysis

The clustering algorithm resulted in different promoter clusters, and on each, we ran Multiple Sequence Alignments using Clustalo [33] against the pUC19 RNAp promoter and found the consensus sequence for each cluster. As stated above, one of the goals of this analysis was an attempt to find genomic evidence of ColE1-like plasmids in other hosts, therefore we decided to focus our efforts on clusters with substantial representation of hosts other than E.coli. One of the clusters had sequences belonging exclusively to the Yersiniaceae family, a member of the order Enterobacteriales.

Part of the clusters of our final RNAp promoters clustering are shown below:

Each color corresponds to a different base, as shown in the legend. Each cluster comprises a series of RNAp promoter sequences with minor changes on top of each other after multiple sequence alignment.

Though we haven’t performed any experiments yet, due to the lack of time, our findings do suggest that a ColE1 system can be found in other organisms.