<jupyter_start><jupyter_text>Regression on RBS sequences - <NAME> - Predict TIRThis notebook shows the regression of RBS sequences.There are several parts of the implementation. - Reading the dataset- Embedding    - Label    - One-hot embedding    - K-mer     - PMW (to be added)    - Unsupervised (to be added)- Kernels    - Spectrum    - Mixed Spectrum    - Weighted Degree    - Weighted Degree with Shifting- Evaluate    - RMSE score    - Plot<jupyter_code># direct to proper path
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import numpy as np

import pandas as pd
import matplotlib.pyplot as plt
import itertools
from collections import defaultdict

from codes.embedding import Embedding
from codes.environment import Rewards_env
from codes.ucb import GPUCB
from codes.evaluations import evaluate, plot_eva
from codes.regression import Regression
from codes.kernels_pairwise import spectrum_kernel_pw, sum_spectrum_kernel_pw, sum_onehot_spectrum_kernel_pw, mixed_spectrum_kernel_pw, WD_kernel_pw, sum_onehot_WD_kernel_pw, WD_shift_kernel_pw

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import PairwiseKernel, DotProduct, RBF
from sklearn.kernel_ridge import KernelRidge

from ipywidgets import IntProgress
from IPython.display import display
import warnings
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Reading the Dataset<jupyter_code># Data downloaded from https://github.com/synbiochem/opt-mva
# Paper https://pubs.acs.org/doi/abs/10.1021/acssynbio.8b00398

# A: whole RBS sequence (len: 29)
# B: extracted 20-base RBS seqeunce (A[7:27], len: 20), training features
# C: TIR labels
# D: the design part RBS (B[7:13], len: 6) 

Path = '../data/RBS_seqs.csv'

df = pd.read_csv(Path)
df.columns = ['A', 'B', 'C']
df['D'] = df['B'].str[7:13]
df.head()
# drop the exact same record (same B and C)

df = df.drop_duplicates(subset = ['B', 'C'])
df.shape
Log_flag = False # indicates whether take log label
Norm_method = 'minmax' # indicates how to normalize label (one of 'mean', 'minmax', None)

def normalize(df):
    # take log FC -- possiblely provide Gaussain distribution?
    if Log_flag:
        df['C'] = np.log(df['C'])
    if Norm_method == 'mean':
        # mean normalization
        df['C'] = (df['C'] - df['C'].mean())/df['C'].std()
    elif Norm_method == 'minmax':
        # min-max normalization 
        df['C'] = 2 * (df['C'] - df['C'].min())/(df['C'].max() - df['C'].min()) -1
    else:
        assert Norm_method == None
        
    return df
# data: num_data * 2, columns: [20-base RBS seq (B), TIR labels (C)]
data = np.asarray(normalize(df)[['B', 'C']])
data.shape
plt.hist(data[:,1])<jupyter_output><empty_output><jupyter_text>Regression- Kernel Ridge Regression- Gaussian Process Regression<jupyter_code># indicates whether cross validation (KFOLD)
cross_val_flag = False 

# indicates whether plot predict_label vs. true label
plot_flag = True 

# string kernel list 
# kernels = [spectrum_kernel_pw, mixed_spectrum_kernel_pw, WD_kernel_pw, WD_shift_kernel_pw]
kernels = [sum_onehot_spectrum_kernel_pw, sum_onehot_WD_kernel_pw]

# name dictionaries
regression_name = {KernelRidge: 'KR',
                  GaussianProcessRegressor: 'GPR'}
embedding_name = ['onehot', 'kmer', 'label']
kernel_name = {
               spectrum_kernel_pw: 'spec',
               sum_spectrum_kernel_pw: 'sspec',
               sum_onehot_spectrum_kernel_pw: 'sospec',
               mixed_spectrum_kernel_pw: 'mspec',
               WD_kernel_pw: 'WD',
               sum_onehot_WD_kernel_pw: 'sWD',
               WD_shift_kernel_pw: 'WDshift'}

train_scores = {}
test_scores = {}
def run_regression(model, kernel, embedding):
    if model == KernelRidge:
        reg = Regression(model(kernel = kernel), data, embedding)
    elif model == GaussianProcessRegressor:
        reg = Regression(model(kernel = PairwiseKernel(metric = kernel)), data, embedding)
    reg.train()
    train_score, test_score = reg.evaluate(cross_val_flag = cross_val_flag, plot_flag = plot_flag)
    key = regression_name[model] + '_' + kernel_name[kernel] + '_' + embedding
    train_scores[key] = train_score
    test_scores[key] = test_score
    print()

for model in [KernelRidge, GaussianProcessRegressor]:
    for kernel in kernels:
        run_regression(model, kernel, 'label')<jupyter_output>Model:  KernelRidge(alpha=1, coef0=1, degree=3, gamma=None,
            kernel=<function sum_onehot_spectrum_kernel_pw at 0x7fea4b9b3400>,
            kernel_params=None)
Train RMSE:  0.16841515144307403
Test RMSE:  0.27223719302209237

Model:  KernelRidge(alpha=1, coef0=1, degree=3, gamma=None,
            kernel=<function sum_onehot_WD_kernel_pw at 0x7fea4b9b3598>,
            kernel_params=None)
Train RMSE:  0.16327937702573053
Test RMSE:  0.23693472392610324

Model:  GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
                         kernel=PairwiseKernel(gamma=1.0, metric=<function sum_onehot_spectrum_kernel_pw at 0x7fea4b9b3400>),
                         n_restarts_optimizer=0, normalize_y=False,
                         optimizer='fmin_l_bfgs_b', random_state=None)
Train RMSE:  0.15807836844196377
Test RMSE:  0.24392759138634879

Model:  GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
                         kernel=PairwiseKernel(gamma=1.0, metric=<funct[...]<jupyter_text>Kernel Ridge Regression<jupyter_code># Baseline: onehot embedding with DotProduct kernel

reg = Regression(KernelRidge(kernel = DotProduct()), data, 'onehot')
reg.train()
train_score, test_score = reg.evaluate(cross_val_flag = cross_val_flag)
train_scores['KR_onehot_DP'] = train_score
test_scores['KR_onehot_DP'] = test_score
# kmer embedding with DotProduct kernel (l = 3)

reg = Regression(KernelRidge(kernel = DotProduct()), data, 'kmer')
reg.train()
train_score, test_score = reg.evaluate(cross_val_flag = cross_val_flag)
train_scores['KR_kmer_DP'] = train_score
test_scores['KR_kmer_DP'] = test_score<jupyter_output>Model:  KernelRidge(alpha=1, coef0=1, degree=3, gamma=None,
            kernel=DotProduct(sigma_0=1), kernel_params=None)
Train RMSE:  0.2088366335310831
Test RMSE:  0.24754401163846318<jupyter_text>Gaussian Process Regression<jupyter_code># # baseline: onehot embedding with DotProduct Kernel 
reg = Regression(GaussianProcessRegressor(kernel = DotProduct()), data, 'onehot')
reg.train()
train_score, test_score = reg.evaluate(cross_val_flag = cross_val_flag)
train_scores['GPR_onehot_DP'] = train_score
test_scores['GPR_onehot_DP'] = test_score
# kmer embedding (l = 3)
reg = Regression(GaussianProcessRegressor(kernel = DotProduct()), data, 'kmer')
reg.train()
train_score, test_score = reg.evaluate(cross_val_flag = cross_val_flag)
train_scores['GPR_kmer_DP'] = train_score
test_scores['GPR_kmer_DP'] = test_score
sorted_train_scores = sorted(train_scores.items(), key=lambda kv: kv[1])

import collections

sorted_train_scores = collections.OrderedDict(sorted_train_scores)

plt.bar(range(len(sorted_train_scores)), list(sorted_train_scores.values()), align='center')
plt.xticks(range(len(sorted_train_scores)), list(sorted_train_scores.keys()), rotation = 70)
plt.ylabel('RMSE')
plt.title('Train Scores')
sorted_test_scores = sorted(test_scores.items(), key=lambda kv: kv[1])

import collections

sorted_test_scores = collections.OrderedDict(sorted_test_scores)

plt.bar(range(len(sorted_test_scores)), list(sorted_test_scores.values()), align='center')
plt.xticks(range(len(sorted_test_scores)), list(sorted_test_scores.keys()), rotation = 70)
plt.ylabel('RMSE')
plt.title('Test Scores')
import pickle

#with open('sum_kernels.pickle', 'wb') as handle:
#    pickle.dump(sorted_test_scores, handle, protocol=pickle.HIGHEST_PROTOCOL)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1137549306527548, 'fraction_numerical': 0.029774780506425753, 'mean_word_length': 4.839524517087667, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 23, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '24527', 'license': 'Apache-2.0', 'n_tokens_mistral': 2904, 'n_tokens_neox': 2647, 'n_words': 581, 'path': 'notebook/others/regression/Regression_RBS_TIR_[-1,1].ipynb', 'repo_name': 'mholowko/Solaris', 'size': '7869'}