<jupyter_start><jupyter_text>**Run the application on BRCA1 dataset**<jupyter_code>#import dataset
def create_dataset (file_path):
    data=pd.read_csv(file_path, header=None)
    data.columns = ['Data_{}'.format(int(i)+1) for i in data.columns]
    data.columns = [*data.columns[:-1], 'labels']
    return data
data = create_dataset('/home/k.benrayana/PathogenicAaMutations/data/BRCA1_CTerDom.csv') #BRCA1 dataset

data

#prepare the sequences and labels

d = np.array(data)
sample_len = d.shape[0] # number of protein sequences in the whole dataset
time_len = d.shape[1]-1 #  length of the protein seq in the whole  dataset
#Hydrophilic coding vocabulary for hydrophobicity representation
vocab1 = {'A': 1.8, 'R':-4.5,'N':-3.5,'D':-3.5,'C':2.5,'Q':-3.5,'E':-3.5,'G':-0.4,'H' :-3.2,'I':4.5,
         'L':3.8,'K':-3.9,'M':1.9,'F':2.8,'P':-1.6,'S':-0.8,'T':-0.7,'W':-0.9,'Y':-1.3,'V':4.2,'X':-3.5,'*':-100}
#Mutation type Benign:0; Pathogenic:1
tags = {'0':0, '1':1}
#prepare data splits for embeddings + hydrophbicity index
data=np.array(data)
data_len = data.shape[0]
data_seqlen = data.shape[1] -1

sequences = data[:, 0:data_seqlen]
labels = data[:, data_seqlen]

#replace each Aa of sequences with hydrophobicity index
Aa_hydro = [[vocab1[word] for word in sentence] for sentence in sequences]
Aa_hydro = np.array(Aa_hydro).reshape(data_len,data_seqlen)
Aa_hydro
#Prepare labels reshape (data_len, 1)
labels = labels.reshape(data_len,1)
#compute embeddings of all sequences
transformer = BioTransformers(backend="protbert_bfd", device="cuda")
embeddings= transformer.compute_embeddings(sequences)["mean"]
labels
Aa_hydro
embeddings
#embedding representation  with hydrophibicity encoder
sequences_hydro = np.concatenate((embeddings,Aa_hydro),axis=1)

sequences_hydro

#split data 70% training set 15% validation set 15 % test or evaluation set
over = RandomOverSampler()  # balance training set

def split_dataset(sequences, labels):
    x_train, x_test, y_train, y_test = train_test_split(sequences, np.array(labels), test_size=0.15, random_state=1)
    x_train, x_val, y_train, y_val = train_test_split(x_train, y_train, test_size=0.176, random_state=2) 
    y_train=y_train.astype(float).reshape(y_train.size,1)
    y_test=y_test.astype(float).reshape(y_test.size,1)
    y_val=y_val.astype(float).reshape(y_val.size,1)
    x_train, y_train = over.fit_resample(x_train, y_train)
    return x_train, y_train, x_val, y_val, x_test, y_test
#reshape final representation to 3D as an input for CNN/BiLSTM classifier 
def model_input(x_train, x_test, x_val):

    x_train=np.reshape(x_train, (x_train.shape[0], 1, x_train.shape[1]))
    x_test=np.reshape(x_test, (x_test.shape[0], 1, x_test.shape[1]))
    x_val=np.reshape(x_val, (x_val.shape[0], 1, x_val.shape[1]))
    return x_train, x_test, x_val
x_train, y_train, x_test, y_test, x_val, y_val=split_dataset(sequences_hydro, labels) #split dataset to train/validation/test
x_train, x_test, x_val=model_input(x_train, x_test, x_val)

#verify if the train set is balanced
unique, counts = np.unique(y_train, return_counts=True)
dict(zip(unique, counts))

#Classifier CNN/BiLSTM
input_shape = x_train.shape[2] # shape of input: 1024 (protbert_bfd embeddings) + size of protein sequence (BRCA1: 225/PTEN: 403)
timesteps = x_train.shape[0] # number of protein sequences
model = Sequential()
model.add(InputLayer(input_shape=(timesteps, input_shape)))
model.add(Conv1D(filters=128, kernel_size=1, activation='relu')) # expected input data shape: (batch_size, timesteps, input_shape)
model.add(MaxPooling1D(pool_size=3, padding="same"))
model.add(Bidirectional(LSTM(100, return_sequences=True)))
model.add(Dense(1, activation='sigmoid'))
bce_loss = tf.keras.losses.BinaryCrossentropy() #BCE loss function
accuracy = tf.keras.metrics.BinaryAccuracy() # Binary Accuracy (expressed in the [0,1] interval)
optimizer = tf.optimizers.Adam(learning_rate = 0.001) #Adam optimizer
model.compile(optimizer, bce_loss, accuracy)  
model.fit(x_train, y_train,validation_data=(x_val, y_val), epochs=300, batch_size=32) #Training of the model

#architecture of the model

model.summary()
#architecture of the model
tf.keras.utils.plot_model(
    model, to_file='app_BRCA1.png', show_shapes=True, 
    show_layer_names=True, rankdir='TB', expand_nested=False, dpi=96
)<jupyter_output><empty_output><jupyter_text>model evaluationy_pred = model(x_test).squeeze().detach().numpy()model_evaluation_accuracy(y_test, y_pred)<jupyter_code>#model evaluation on test set
y_pred = model.predict(x_test)   
prediction=pd.DataFrame(y_pred.flatten()).rename(columns={0:'probability'}) #convert predictionss to dataframe
y_pred=prediction['probability']
#labels
y_test=y_test.astype(float).reshape(y_test.size,1) # convert numpy array elements from object to float for labels

prediction['targets']=pd.DataFrame(y_test) #convert labels from numpy array to DataFrame and convert elements from object to float
y_pred=prediction['probability']

y_true=prediction['targets']
#Plot Precision-Recall

# calculate model precision-recall curve
precision, recall, _ = precision_recall_curve(y_true, y_pred)

plt.figure()
plt.plot(precision, recall, color='green', label='protbert_bfd_CNN_BiLSTM')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.title('Precision-Recall Curve ')
plt.legend()
plt.show()
plt.savefig('PR_curve.png')
# Plot ROC_AUC

fpr, tpr, threshold = roc_curve(y_true, y_pred)
#roc_auc = auc(fpr, tpr)
plt.figure()
plt.plot(fpr, tpr, color='orange', label='ROC')
plt.plot([0, 1], [0, 1], color='darkblue', linestyle='--')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve')
plt.legend()
plt.show()
plt.savefig('ROC_curve.png')
auc = roc_auc_score(y_true, y_pred)
print('AUC Score: %.2f' % auc)
# Plot confusion matrix
model_evaluation_accuracy(y_true, y_pred)

confusion_matrix_plot(y_true, (y_pred> 0.5).astype(int), ["0", "1"])
matrix_index = ["0","1"]
from sklearn.metrics import classification_report 
y_pred= np.around(y_pred)
print(classification_report(y_true, y_pred,target_names=matrix_index))

#Compare Predictions and targets distribution

fig, axes = plt.subplots(1, 2)
fig.suptitle('Probability and targets distribution for BRCA1 dataset', fontsize=16, fontweight='bold')
fig.set_figheight(7)
fig.set_figwidth(12)
axes[0].set_title('Targets distribution')
axes[1].set_title('Prediction distributions')


sns.set(font_scale=1.2)
sns.set_style("whitegrid")
plt.xlim(0, 1)
plt.ylim(0, 0.8)

plt.xlabel('Probability')
plt.ylabel('Count')

sns.histplot(ax=axes[1],data=prediction, x="probability", hue="targets", multiple="stack", stat='probability')
plt.legend(['Pathogenic', 'Benign'],loc='best', fontsize=8)

plt.tight_layout()
plt.figsize=(10,5)
sns.countplot(ax=axes[0],x=y_true)
y_true.value_counts()<jupyter_output><empty_output><jupyter_text>**Run the application on PTEN dataset**<jupyter_code>def create_dataset (file_path):
    data=pd.read_csv(file_path, header=None)
    data.columns = ['Data_{}'.format(int(i)+1) for i in data.columns]
    data.columns = [*data.columns[:-1], 'labels']
    return data
data = create_dataset('/home/k.benrayana/PathogenicAaMutations/data/Pten_dataset.csv') #PTEN dataset
data
#prepare the sequences and labels

d = np.array(data)
sample_len = d.shape[0] # number of protein sequences in the whole dataset
time_len = d.shape[1]-1 #  length of the protein seq in the whole  dataset
#Hydrophilic coding vocabulary for hydrophobicity representation
vocab1 = {'A': 1.8, 'R':-4.5,'N':-3.5,'D':-3.5,'C':2.5,'Q':-3.5,'E':-3.5,'G':-0.4,'H' :-3.2,'I':4.5,
         'L':3.8,'K':-3.9,'M':1.9,'F':2.8,'P':-1.6,'S':-0.8,'T':-0.7,'W':-0.9,'Y':-1.3,'V':4.2,'X':-3.5,'*':-100}
#Mutation type Benign:0; Pathogenic:1
tags = {'0':0, '1':1}
#prepare data splits for embeddings + hydrophbicity index
data=np.array(data)
data_len = data.shape[0]
data_seqlen = data.shape[1] -1

sequences = data[:, 0:data_seqlen]
labels = data[:, data_seqlen]

#replace each Aa of sequences with hydrophobicity index
Aa_hydro = [[vocab1[word] for word in sentence] for sentence in sequences]
Aa_hydro = np.array(Aa_hydro).reshape(data_len,data_seqlen)
Aa_hydro
#Prepare labels reshape (data_len, 1)
labels = labels.reshape(data_len,1)
#compute embeddings of all sequences
transformer = BioTransformers(backend="protbert_bfd", device="cuda")
embeddings= transformer.compute_embeddings(sequences)["mean"]
labels
Aa_hydro
embeddings
#embedding representation  with hydrophibicity encoder
sequences_hydro = np.concatenate((embeddings,Aa_hydro),axis=1)
sequences_hydro
#split data 70% training set 15% validation set 15 % test or evaluation set
over = RandomOverSampler()  # balance training set
def split_dataset(sequences, labels):
    x_train, x_test, y_train, y_test = train_test_split(sequences, np.array(labels), test_size=0.15, random_state=1)
    x_train, x_val, y_train, y_val = train_test_split(x_train, y_train, test_size=0.176, random_state=2) 
    y_train=y_train.astype(float).reshape(y_train.size,1)
    y_test=y_test.astype(float).reshape(y_test.size,1)
    y_val=y_val.astype(float).reshape(y_val.size,1)
    x_train, y_train = over.fit_resample(x_train, y_train)
    return x_train, y_train, x_val, y_val, x_test, y_test
#reshape final representation as an input for CNN/BiLSTM classifier i
def model_input(x_train, x_test, x_val):

    x_train=np.reshape(x_train, (x_train.shape[0], 1, x_train.shape[1]))
    x_test=np.reshape(x_test, (x_test.shape[0], 1, x_test.shape[1]))
    x_val=np.reshape(x_val, (x_val.shape[0], 1, x_val.shape[1]))
    return x_train, x_test, x_val
x_train, y_train, x_test, y_test, x_val, y_val=split_dataset(sequences_hydro, labels) #split dataset to train/validation/test

 x_train, x_test, x_val=model_input(x_train, x_test, x_val)      #reshape final representation to 3D as an input for CNN/BiLSTM classifier 

#verify if the train set is balanced
unique, counts = np.unique(y_train, return_counts=True)
dict(zip(unique, counts))

#Classifier CNN/BiLSTM
input_shape = x_train.shape[2]   # shape of input: 1024 (protbert_bfd embeddings) + size of protein sequence (BRCA1: 225/PTEN: 403)
timesteps = x_train.shape[0] # number of protein sequences
model = Sequential()
model.add(InputLayer(input_shape=(timesteps, input_shape)))
model.add(Conv1D(filters=128, kernel_size=1, activation='relu')) # expected input data shape: (batch_size, timesteps, input_shape)
model.add(MaxPooling1D(pool_size=3, padding="same"))
model.add(Bidirectional(LSTM(100, return_sequences=True)))
model.add(Dense(1, activation='sigmoid'))
bce_loss = tf.keras.losses.BinaryCrossentropy() #BCE loss function
accuracy = tf.keras.metrics.BinaryAccuracy() # Binary Accuracy (expressed in the [0,1] interval)
optimizer = tf.optimizers.Adam(learning_rate = 0.001) #Adam optimizer
model.compile(optimizer, bce_loss, accuracy)  
model.fit(x_train, y_train,validation_data=(x_val, y_val), epochs=300, batch_size=32) #Training of the model
#architecture of the model

model.summary()
#architecture of the model
tf.keras.utils.plot_model(
    model, to_file='app_BRCA1.png', show_shapes=True, 
    show_layer_names=True, rankdir='TB', expand_nested=False, dpi=96
)
#model evaluation on test set
y_pred = model.predict(x_test)   
prediction=pd.DataFrame(y_pred.flatten()).rename(columns={0:'probability'}) #convert predictionss to dataframe
y_pred=prediction['probability']
#labels
y_test=y_test.astype(float).reshape(y_test.size,1) # convert numpy array elements from object to float for labels

prediction['targets']=pd.DataFrame(y_test) #convert labels from numpy array to DataFrame and convert elements from object to float
y_pred=prediction['probability']

y_true=prediction['targets']
#Plot Precision-Recall

# calculate model precision-recall curve
precision, recall, _ = precision_recall_curve(y_true, y_pred)

plt.figure()
plt.plot(precision, recall, color='green', label='protbert_bfd_CNN_BiLSTM')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.title('Precision-Recall Curve ')
plt.legend()
plt.show()
plt.savefig('PR_curve.png')
# Plot ROC_AUC

fpr, tpr, threshold = roc_curve(y_true, y_pred)
#roc_auc = auc(fpr, tpr)
plt.figure()
plt.plot(fpr, tpr, color='orange', label='ROC')
plt.plot([0, 1], [0, 1], color='darkblue', linestyle='--')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve')
plt.legend()
plt.show()
plt.savefig('ROC_curve.png')
auc = roc_auc_score(y_true, y_pred)
print('AUC Score: %.2f' % auc)
# Plot confusion matrix
model_evaluation_accuracy(y_true, y_pred)

confusion_matrix_plot(y_true, (y_pred> 0.5).astype(int), ["0", "1"])
matrix_index = ["0","1"]
from sklearn.metrics import classification_report 
y_pred= np.around(y_pred)
print(classification_report(y_true, y_pred,target_names=matrix_index))

#compare Predictions and targets distribution

fig, axes = plt.subplots(1, 2)
fig.suptitle('Probability and targets distribution for PTEN dataset', fontsize=16, fontweight='bold')
fig.set_figheight(7)
fig.set_figwidth(12)
axes[0].set_title('Targets distribution')
axes[1].set_title('Prediction distributions')


sns.set(font_scale=1.2)
sns.set_style("whitegrid")
plt.xlim(0, 1)
plt.ylim(0, 0.8)

plt.xlabel('Probability')
plt.ylabel('Count')

sns.histplot(ax=axes[1],data=prediction, x="probability", hue="targets", multiple="stack", stat='probability')
plt.legend(['Pathogenic', 'Benign'],loc='best', fontsize=8)

plt.tight_layout()
plt.figsize=(10,5)
sns.countplot(ax=axes[0],x=y_true)
y_true.value_counts()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14672270146722702, 'fraction_numerical': 0.0247732802477328, 'mean_word_length': 6.8449971081550025, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '464125', 'license': 'Apache-2.0', 'n_tokens_mistral': 5333, 'n_tokens_neox': 5007, 'n_words': 1040, 'path': 'Application_notebook.ipynb', 'repo_name': 'KhalilBR94/PathogenicAaMutation', 'size': '13563'}