<jupyter_start><jupyter_text>**Predicting the activity of antimicrobial peptides with machine learning models****Author:** [<NAME>](https://sayalaruano.github.io/)This notebook is part of my capstone project for the [Machine Learning Zoomcamp](https://datatalks.club/courses/2021-winter-ml-zoomcamp.html).[Antimicrobial peptides](https://en.wikipedia.org/wiki/Antimicrobial_peptides) (AMPs) are small bioactive drugs, commonly with fewer than 50 amino acids, which have appeared as promising compounds to control infectious disease caused by multi-drug resistant bacteria or superbugs. These superbugs are not treatable with the available drugs because of the development of some mechanisms to avoid the action of these compounds, which is known as antimicrobial resistance (AMR). According to the World Health Organization, AMR is one of the [top ten global public health threats facing humanity in this century](https://www.who.int/news-room/fact-sheets/detail/antimicrobial-resistance), so it is important to search for AMPs that combat these superbugs and prevent AMR.However, the search for AMPs to combat superbugs by experimental methods is unfeasible because of the huge number of these compounds. So, it is required to explore other methods for handling this problem, and machine learning (ML) models could be great candidates for this task. Thus, in this project I created some ML binary classifiers to predict the activity of antimicrobial peptides.You can find more details about my capstone project in this [GitHub repository](https://github.com/sayalaruano/CapstoneProject-MLZoomCamp).For this work, I took as a reference the [notebook](https://github.com/dataprofessor/peptide-ml) and [video](https://www.youtube.com/watch?v=0NrFIGLwW0Q&feature=youtu.be) from [Dataprofessor](https://github.com/dataprofessor) about this topic. Also, the datasets, some ideas, and references to compare the performance of the best model were obtained from this [article](https://pubs.acs.org/doi/10.1021/acs.jcim.1c00251).**Datasets:** https://biocom-ampdiscover.cicese.mx/dataset**Note:** `mkdir`, `wget`, `unzip`, `make`, `cd`, `grep`, and `cat` are Unix based OS commands, so if you are in Windows, maybe these terminal tools won't work.  **Installation of libraries**  **[LazyPredict](https://github.com/shankarpandala/lazypredict)**Library to run many sklearn's binary classifiers with default parameters<jupyter_code>!pip install lazypredict<jupyter_output><empty_output><jupyter_text>**[Pfeature](https://github.com/raghavagps/Pfeature)**Library to calculate molecular descriptors (features) for the AMPs<jupyter_code># Download and unzip the files to install Pfeature
! wget https://github.com/raghavagps/Pfeature/raw/master/PyLib/Pfeature.zip
! unzip Pfeature.zip
# Move to the Pfeature folder - % is required if you are working in colab
% cd Pfeature
# Run the installation script
! python setup.py install
# Move to the content directory - % is required if you are working in colab
% cd ../<jupyter_output>UsageError: Line magic function `%` not found.<jupyter_text>**[CD-HIT](https://github.com/weizhongli/cdhit)**Program to remove the redundant sequences of the AMPs<jupyter_code># # Download and unzip the files to install CD-HIT
!wget https://github.com/weizhongli/cdhit/releases/download/V4.8.1/cd-hit-v4.8.1-2019-0228.tar.gz
!tar -xzvf cd-hit-v4.8.1-2019-0228.tar.gz
# Move to the Pfeature folder
% cd cd-hit-v4.8.1-2019-0228/
# Run the installation script
!make
# Move to the content directory
% cd ../<jupyter_output>/content<jupyter_text>**Imports, data obtention, and pre-processing**<jupyter_code># Data manipulation and EDA
import pandas as pd
import numpy as np

# Plots
import matplotlib.pyplot as plt
import seaborn as sns

# Features calculation
from Pfeature.pfeature import aac_wp
from Pfeature.pfeature import atc_wp
from Pfeature.pfeature import pcp_wp
from Pfeature.pfeature import rri_wp
from Pfeature.pfeature import ddr_wp
from Pfeature.pfeature import sep_wp
from Pfeature.pfeature import ser_wp
from Pfeature.pfeature import spc_wp
from Pfeature.pfeature import ctc_wp
from Pfeature.pfeature import ctd_wp

# Machine learning classifers
import lazypredict
import pickle
from lazypredict.Supervised import LazyClassifier
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, f1_score, matthews_corrcoef
from sklearn.metrics import make_scorer
# Download data 
# Active AMPs
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/TR_starPep_active_AMPs.fasta
# Non-active AMPs
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/TR_starPep_nonactive_AMPs.fasta
# Print the content of one of the files
!cat TR_starPep_active_AMPs.fasta<jupyter_output><empty_output><jupyter_text>**Remove redundant sequences using CD-HIT**<jupyter_code># Active AMPs
!cd-hit-v4.8.1-2019-0228/cd-hit -i TR_starPep_active_AMPs.fasta -o TR_starPep_active_AMPs.txt -c 0.99
# Non-active AMPs
!cd-hit-v4.8.1-2019-0228/cd-hit -i TR_starPep_nonactive_AMPs.fasta -o TR_starPep_nonactive_AMPs.txt -c 0.99
# Number of active AMPs (before processing)
! grep -c ">" TR_starPep_active_AMPs.fasta
# Number of non-redundant and active AMPs (after processing)
! grep -c ">" TR_starPep_active_AMPs.txt
# Number of non-active AMPs (before processing)
! grep -c ">" TR_starPep_nonactive_AMPs.fasta
# Number of non-redundant and non-active AMPs (after processing)
!grep -c ">" TR_starPep_nonactive_AMPs.txt<jupyter_output>7666<jupyter_text>**Feature matrices preparation**For this work, I will consider ten features that require only input and output files as parameters. These features are summarized in the table below.Feature class | Description | Function---|---|---AAC | Amino acid composition | aac_wpABC | Atom and bond composition | atc_wp, btc_wpPCP | Physico-chemical properties | pcp_wpRRI | Repetitive Residue Information | rri_wpDDR | Distance distribution of residues |ddr_wpSEP | Shannon entropy | sep_wpSER | Shannon entropy of residue level | ser_wpSPC | Shannon entropy of physicochemical property | spc_wpCTC | Conjoint Triad Calculation | ctc_wpCTD | Composition enhanced transition distribution | ctd_wpTo know the details of these features, you can read the [Pfeature Manual](https://webs.iiitd.edu.in/raghava/pfeature/Pfeature_Manual.pdf)  **Calculation of molecular descriptors (features) with Pfeature**You can run this section by yourself, but due to the size of the dataset it can take a while, so I recommend you to use output files provided [here](https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/PF_features.zip), and move to the next section to load this data.<jupyter_code># Create a list with names of functions to obtain the features
PF_list = [aac_wp,
           atc_wp,
           pcp_wp,
           rri_wp,
           ddr_wp,
           sep_wp,
           ser_wp, 
           spc_wp,
           ctc_wp, 
           ctd_wp]<jupyter_output><empty_output><jupyter_text>**Function for calculating the different features**<jupyter_code>def PF_feature_calc(input, feature):
  temp = input.rstrip('txt')
  output = temp + feature.__name__ + '.csv'
  df_out = feature(input, output)
  df_in = pd.read_csv(output)
  return df_in<jupyter_output><empty_output><jupyter_text>**Function to create feature matrices**This function calculates features for both positive and negative classes, combines the two classes, and merges the results into a single dataframe.<jupyter_code># Create a directory to store the cv files of the fature matrices
!mkdir PF_features
# Define variables with names of the active and non-active AMPs files
pos = 'TR_starPep_active_AMPs.txt'
neg = 'TR_starPep_nonactive_AMPs.txt'

def act_non_act_concat(po, ne, feature_name):
  # Calculate feature
  po_feature = PF_feature_calc(po, feature_name)
  ne_feature = PF_feature_calc(ne, feature_name)
  # Create class labels
  po_class = pd.Series(['positive' for i in range(len(po_feature))])
  ne_class = pd.Series(['negative' for i in range(len(ne_feature))])
  # Combine po and ne
  po_ne_class = pd.concat([po_class, ne_class], axis=0)
  po_ne_class.name = 'class'
  po_ne_feature = pd.concat([po_feature, ne_feature], axis=0)
  # Combine feature and class
  df = pd.concat([po_ne_feature, po_ne_class], axis=1)
  # Export th csv file of the feature matrix
  #out_name = feature_name.__name__ + ".csv"
  #df.to_csv("PF_features/" + out_name)
  return df
# Define a dictionary to store all feature matrices
PF_features_dfs = {}

# Loop to create feature matrices with all features
for i in PF_list: 
  temp = act_non_act_concat(pos, neg, i)
  PF_features_dfs[i.__name__] = temp
!mkdir PF_features<jupyter_output><empty_output><jupyter_text>**Load molecular descriptors (features)**<jupyter_code># Download and unzip calculated molecular descriptor files
! mkdir PF_features
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/PF_features.zip --directory PF_features
! unzip PF_features/PF_features.zip -d PF_features
# Create a list with names of the feature's files
PF_files_list = [
           "aac_wp",
           "atc_wp",
           "pcp_wp",
           "rri_wp",
           "ddr_wp",
           "sep_wp",
           "ser_wp", 
           "spc_wp",
           "ctc_wp", 
           "ctd_wp"]
# Load target variable 
temp = pd.read_csv("PF_features/aac_wp.csv", index_col=0)

y = temp['class'].copy()

# Encoding the target class label
y = y.map({"positive": 1, "negative": 0})
# Load original molecular descriptor files and storage them in a dictionary
PF_features_dfs = {}

for i in PF_files_list:
  path = "PF_features/"
  feature_file_path = ''.join([path,i])
  
  feature_name_csv = ''.join([feature_file_path,'.csv'])
  print(feature_name_csv)
  
  PF_features_dfs[str(i)] = pd.read_csv(feature_name_csv, index_col=0)
  PF_features_dfs[str(i)] = PF_features_dfs[str(i)] .drop('class', axis=1)<jupyter_output><empty_output><jupyter_text>**Exploratory Data Analysis**  **Target variable**<jupyter_code># Search missing values 
PF_features_dfs["aac_wp"]["class"].isnull().sum()
# Bar plot of active and inactive AMPs

# Data
active = (PF_features_dfs["aac_wp"]["class"] == "positive").sum()
inactive = (PF_features_dfs["aac_wp"]["class"] == "negative").sum()

x = ['Active', 'Inactive']
y = [active, inactive]

# Setup plot
fig, ax = plt.subplots()

# Make bar plot
p = ax.bar(x, y, color = ['#F8766D', '#00BFC4'], ec = 'black')

ax.set_title('AMPs training dataset activity', fontsize=14, fontweight='bold', pad=15)
#ax.set_xticklabels(x, fontweight='bold')

ax.set_ylim(0,10000)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

# Label above bar
for index, data in enumerate(y):
    plt.text(x=index-0.1 , y =data+1000 , s=f"{data}" , fontdict=dict(fontsize=14))

fig.set_size_inches(5,4.5)
plt.show()<jupyter_output><empty_output><jupyter_text>The previous plot shows that this dataset is balanced regarding their classes  **Amino acid composition  feature matrix**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["aac_wp"]
# General measurements of the feature matrix 
PF_features_dfs["aac_wp"].info()
# General statistics of the feature matrix for all AMPs
feat_stat = PF_features_dfs["aac_wp"].describe().T

# Bar plot of mean aminoacid composition for all the AMPs
mean_values = feat_stat["mean"]

mean_values.plot.bar(figsize=(8,4), ec='black')

plt.title('Mean aminoacid composition for all the AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for active AMPs
feat_stat_act = PF_features_dfs["aac_wp"][PF_features_dfs["aac_wp"]["class"] 
                                          == "positive"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_act = feat_stat_act["mean"]

mean_values_act.plot.bar(figsize=(8,4), color = 'green', ec='black')

plt.title('Mean aminoacid composition for active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for non-active AMPs
feat_stat_nact = PF_features_dfs["aac_wp"][PF_features_dfs["aac_wp"]["class"] 
                                          == "negative"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_nact = feat_stat_nact["mean"]

mean_values_nact.plot.bar(figsize=(8,4), color = 'red', ec='black')

plt.title('Mean aminoacid composition for non-active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()<jupyter_output><empty_output><jupyter_text>We observed that the distribution of the mean aminoacid composition from active and non-active AMPs is different.  **Atom and bond composition**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["atc_wp"]
# General measurements of the feature matrix 
PF_features_dfs["atc_wp"].info()
# General statistics of the feature matrix for all AMPs
feat_stat = PF_features_dfs["atc_wp"].describe().T

# Bar plot of mean aminoacid composition for all the AMPs
mean_values = feat_stat["mean"]

mean_values.plot.bar(figsize=(8,4), ec='black')

plt.title('Mean atom composition for all the AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for active AMPs
feat_stat_act = PF_features_dfs["atc_wp"][PF_features_dfs["atc_wp"]["class"] 
                                          == "positive"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_act = feat_stat_act["mean"]

mean_values_act.plot.bar(figsize=(8,4), color = 'green', ec='black')

plt.title('Mean atom composition for active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for non-active AMPs
feat_stat_nact = PF_features_dfs["atc_wp"][PF_features_dfs["atc_wp"]["class"] 
                                          == "negative"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_nact = feat_stat_nact["mean"]

mean_values_nact.plot.bar(figsize=(8,4), color = 'red', ec='black')

plt.title('Mean atom composition for non-active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()<jupyter_output><empty_output><jupyter_text>There are not apparent differences in the distribution of mean atom composition from active and non-active AMPs  **Physico-chemical properties**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["pcp_wp"]
# General measurements of the feature matrix 
PF_features_dfs["pcp_wp"].info()
# General statistics of the feature matrix for all AMPs
feat_stat = PF_features_dfs["pcp_wp"].describe().T

# Bar plot of mean aminoacid composition for all the AMPs
mean_values = feat_stat["mean"]

mean_values.plot.bar(figsize=(8,4), ec='black')

plt.title('Mean physico-chemical properties for all the AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for active AMPs
feat_stat_act = PF_features_dfs["pcp_wp"][PF_features_dfs["pcp_wp"]["class"] 
                                          == "positive"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_act = feat_stat_act["mean"]

mean_values_act.plot.bar(figsize=(8,4), color = 'green', ec='black')

plt.title('Mean physico-chemical properties for active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for non-active AMPs
feat_stat_nact = PF_features_dfs["pcp_wp"][PF_features_dfs["pcp_wp"]["class"] 
                                          == "negative"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_nact = feat_stat_nact["mean"]

mean_values_nact.plot.bar(figsize=(8,4), color = 'red', ec='black')

plt.title('Mean physico-chemical properties for non-active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()<jupyter_output><empty_output><jupyter_text>There are not apparent differences in the distribution of mean physocochemical properties from active and non-active AMPs  **Repetitive Residue Information**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["rri_wp"]
# General measurements of the feature matrix 
PF_features_dfs["rri_wp"].info()
# General statistics of the feature matrix for all AMPs
feat_stat = PF_features_dfs["rri_wp"].describe().T

# Bar plot of mean aminoacid composition for all the AMPs
mean_values = feat_stat["mean"]

mean_values.plot.bar(figsize=(8,4), ec='black')

plt.title('Mean repetitive residues for all the AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for active AMPs
feat_stat_act = PF_features_dfs["rri_wp"][PF_features_dfs["rri_wp"]["class"] 
                                          == "positive"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_act = feat_stat_act["mean"]

mean_values_act.plot.bar(figsize=(8,4), color = 'green', ec='black')

plt.title('Mean epetitive residues for active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for non-active AMPs
feat_stat_nact = PF_features_dfs["rri_wp"][PF_features_dfs["rri_wp"]["class"] 
                                          == "negative"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_nact = feat_stat_nact["mean"]

mean_values_nact.plot.bar(figsize=(8,4), color = 'red', ec='black')

plt.title('Mean epetitive residues for non-active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()<jupyter_output><empty_output><jupyter_text>We observed that the distribution of the mean repetitive residue information from active and non-active AMPs is different.  **Distance distribution of residues**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["ddr_wp"]
# General measurements of the feature matrix 
PF_features_dfs["ddr_wp"].info()
# General statistics of the feature matrix for all AMPs
feat_stat = PF_features_dfs["ddr_wp"].describe().T

# Bar plot of mean aminoacid composition for all the AMPs
mean_values = feat_stat["mean"]

mean_values.plot.bar(figsize=(8,4), ec='black')

plt.title('Mean distance distribution of residues for all the AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for active AMPs
feat_stat_act = PF_features_dfs["ddr_wp"][PF_features_dfs["ddr_wp"]["class"] 
                                          == "positive"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_act = feat_stat_act["mean"]

mean_values_act.plot.bar(figsize=(8,4), color = 'green', ec='black')

plt.title('Mean distance distribution of residues for active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# General statistics of the feature matrix for non-active AMPs
feat_stat_nact = PF_features_dfs["ddr_wp"][PF_features_dfs["ddr_wp"]["class"] 
                                          == "negative"].describe().T

# Bar plot of mean aminoacid composition for active AMPs
mean_values_nact = feat_stat_nact["mean"]

mean_values_nact.plot.bar(figsize=(8,4), color = 'red', ec='black')

plt.title('Mean distance distribution of residues for non-active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()<jupyter_output><empty_output><jupyter_text>We observed that the distribution of the mean distnace distribution of residues from active and non-active AMPs is different.  **Shannon entropy**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["sep_wp"]
# General measurements of the feature matrix 
PF_features_dfs["sep_wp"].info()
# Histogram of shanon entropy for all the AMPs
SEP = PF_features_dfs["sep_wp"]["SEP"]

SEP.hist(figsize=(8,4), ec='black')

plt.title('Histogram of shanon entropy values for all the AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
PF_features_dfs["sep_wp"]["SEP"][PF_features_dfs["sep_wp"]["class"] 
                                          == "positive"]
# Histogram of shanon entropy for active AMPs
SEP_act = PF_features_dfs["sep_wp"]["SEP"][PF_features_dfs["sep_wp"]["class"] 
                                          == "positive"]

SEP_act.hist(figsize=(8,4),color='green', ec='black')

plt.title('Histogram of shanon entropy values for active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
# Histogram of shanon entropy for non-active AMPs
SEP_nact = PF_features_dfs["sep_wp"]["SEP"][PF_features_dfs["sep_wp"]["class"] 
                                          == "negative"]

SEP_nact.hist(figsize=(8,4),color='red', ec='black')

plt.title('Histogram of shanon entropy values for non-active AMPs', fontsize=14, fontweight='black', pad=15)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()<jupyter_output><empty_output><jupyter_text>There are not apparent differences in the distribution of shanon entropy from active and non-active AMPs  **Shannon entropy of residue level**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["ser_wp"]
# General measurements of the feature matrix 
PF_features_dfs["ser_wp"].info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 14609 entries, 0 to 7665
Data columns (total 21 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   SER_A   14609 non-null  float64
 1   SER_C   14609 non-null  float64
 2   SER_D   14609 non-null  float64
 3   SER_E   14609 non-null  float64
 4   SER_F   14609 non-null  float64
 5   SER_G   14609 non-null  float64
 6   SER_H   14609 non-null  float64
 7   SER_I   14609 non-null  float64
 8   SER_K   14609 non-null  float64
 9   SER_L   14609 non-null  float64
 10  SER_M   14609 non-null  float64
 11  SER_N   14609 non-null  float64
 12  SER_P   14609 non-null  float64
 13  SER_Q   14609 non-null  float64
 14  SER_R   14609 non-null  float64
 15  SER_S   14609 non-null  float64
 16  SER_T   14609 non-null  float64
 17  SER_V   14609 non-null  float64
 18  SER_W   14609 non-null  float64
 19  SER_Y   14609 non-null  float64
 20  class   14609 non-null  object 
dtypes: float64(20), object(1)
memory usa[...]<jupyter_text>**Shannon entropy of physicochemical property**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["spc_wp"]
# General measurements of the feature matrix 
PF_features_dfs["spc_wp"].info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 14609 entries, 0 to 7665
Data columns (total 26 columns):
 #   Column     Non-Null Count  Dtype  
---  ------     --------------  -----  
 0   SEP_PC     14609 non-null  float64
 1   SEP_NC     14609 non-null  float64
 2   SEP_NE     14609 non-null  float64
 3   SEP_PO     14609 non-null  float64
 4   SEP_NP     14609 non-null  float64
 5   SEP_AL     14609 non-null  float64
 6   SEP_CY     14609 non-null  float64
 7   SEP_AR     14609 non-null  float64
 8   SEP_AC     14609 non-null  float64
 9   SEP_BS     14609 non-null  float64
 10  SEP_NE_pH  14609 non-null  float64
 11  SEP_HB     14609 non-null  float64
 12  SEP_HL     14609 non-null  float64
 13  SEP_NT     14609 non-null  float64
 14  SEP_HX     14609 non-null  float64
 15  SEP_SC     14609 non-null  float64
 16  SEP_SS_HE  14609 non-null  float64
 17  SEP_SS_ST  14609 non-null  float64
 18  SEP_SS_CO  14609 non-null  float64
 19  SEP_SA_BU  14609 non-null  float64
 20  SEP_SA_[...]<jupyter_text>**Conjoint Triad Calculation**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["ctc_wp"]<jupyter_output><empty_output><jupyter_text>**Composition enhanced transition distribution**<jupyter_code># Take a look at the feature matrix
PF_features_dfs["ctd_wp"]<jupyter_output><empty_output><jupyter_text>---  **Machine learning models**First, I tested more than 30 ML binary classifiers using the [LazyPredict](https://github.com/shankarpandala/lazypredict) Python library. I chose the best models according some performance metrics such as accuracy, ROC AUC, F1 score, and Matthews Correlation Coefficient. Then, I fine-tuned the hyperparameters of the best models using sklearn's class [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.htmlsklearn.model_selection.GridSearchCV). Finally, considering the results of hyperparameter tuning and performance metrics, I obtained the best ML model to predict AMPs activity.<jupyter_code># Run ML classifiers for all the feature matrices with LazyPredict 
Results = {}

for i, j in PF_features_dfs.items():
    
    # Data splittung
    X_train, X_test, y_train, y_test = train_test_split(j, y, test_size=0.2, 
                                                        random_state =42, 
                                                        stratify=y)

    # Define and build the Lazyclassifier
    clf = LazyClassifier(verbose=0,ignore_warnings=True, custom_metric=matthews_corrcoef)
    
    models,predictions = clf.fit(X_train, X_test, y_train, y_test)

    # Store the results in the dictionary
    Results[i] = models
# Extract the best 3 models for all the feature matrices and put them on a dataset
results_all_fm = pd.DataFrame()

for i,j in Results.items():
    temp = Results[i].head(3)

    temp["Dataset"] = i

    results_all_fm = results_all_fm.append(temp)
# Sort the dataframe by the Matthews correlation coefficient and identify
# the best models
results_all_fm = results_all_fm.sort_values("matthews_corrcoef", ascending=False)

# Save the table of results
results_all_fm.to_csv("results_best_models_all_feat_mat.csv")
results_all_fm
# Plot of Accuracy
plt.figure(figsize=(5, 10))
sns.set_theme(style="whitegrid")
ax = sns.barplot(y=results_all_fm.index+"_"+ results_all_fm.Dataset, x="Accuracy", 
                 data=results_all_fm)
ax.set(xlim=(0, 1))
# Plot of MCC
plt.figure(figsize=(5, 10))
sns.set_theme(style="whitegrid")
ax = sns.barplot(y=results_all_fm.index+"_"+ results_all_fm.Dataset, 
                 x="matthews_corrcoef", data=results_all_fm)
ax.set(xlim=(0, 1))<jupyter_output><empty_output><jupyter_text>---  **Hyperparameter tuning of the best ML models**According to the previous analysis, the best models were **ExtraTreesClassifier** and **RandomForestClassifier** with ser_awp and aac_wp feature matrices.I performed feature tuning with sklearn class [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.htmlsklearn.model_selection.GridSearchCV). This class allows to specify a grid with features for tuning, and also the number of folds if we want to use kfold cross validation for this process. So, I obtained the best set of features for the best 4 models in this way.<jupyter_code># Dictionary with the best 2 feature matrices
Best2featmat = {}
Best2featmat["ser_wp"] = PF_features_dfs["ser_wp"]
Best2featmat["aac_wp"] = PF_features_dfs["aac_wp"]
# Dataframe to save performance metrics of the best models
results_bestmodels = pd.DataFrame(columns=['Dataset', 'Model', 'ROC_AUC_cv',
                                'Accuracy_cv', 'Precision_cv', 
                                 'Recall_cv', 'F1score_cv', 'MCC_cv'])<jupyter_output><empty_output><jupyter_text>**RandomForest**<jupyter_code># A different model for each feature matrix
for i, j in Best2featmat.items():

    # Data splitting 
    X_train, X_test, y_train, y_test = train_test_split(j, y, test_size=0.2, 
                                                        random_state=42,
                                                        stratify=y)

    # Set the grid of parameters of GridSearchCV
    tuned_parameters_rf = [
        {"max_depth": [10, 20, 30, 50], "n_estimators": [10, 30, 50, 100, 200], 
          "max_features": ["sqrt", "log2", 0.2, 0.5], "n_jobs":[-1], "random_state": [10]}
    ]

    print(f'# Tuning hyper-parameters Random Forest {i}')
    print()

    # Scoring metrics to be calculated
    scoring = {"AUC": "roc_auc", "Accuracy": "accuracy", 
               "Precision": "precision", "Recall": "recall", 
               "F1score": "f1", "MCC": make_scorer(matthews_corrcoef)}

    # Create GridSearchCV object with its parameters 
    rf = GridSearchCV(RandomForestClassifier(), tuned_parameters_rf, scoring=scoring, 
                      cv=3, refit="AUC", return_train_score=True)
    
    # Train models 
    rf.fit(X_train, y_train)

    print("Best parameters set found on development set:")
    print()
    print(rf.best_params_)
    print(rf.best_score_)

    # Obtain performance metrics of the best model  
    roc_auc_cv = max(rf.cv_results_["mean_test_AUC"])
    print("roc_auc_cv: ", roc_auc_cv)

    acc_cv = max(rf.cv_results_["mean_test_Accuracy"])
    print("acc_cv: ", acc_cv)

    prec_cv = max(rf.cv_results_["mean_test_Precision"])
    print("prec_cv: ", prec_cv)

    recall_cv = max(rf.cv_results_["mean_test_Recall"])
    print("recall_cv: ", recall_cv)

    f1score_cv = max(rf.cv_results_["mean_test_F1score"])
    print("f1score_cv: ", f1score_cv)

    mcc_cv = max(rf.cv_results_["mean_test_MCC"])
    print("mcc_cv: ", mcc_cv)

    # Add information to the results dataframe
    results_bestmodels = results_bestmodels.append({'Dataset': i, 'Model': 'RandomForest', 
                                                    'ROC_AUC_cv': roc_auc_cv,
                                                    'Accuracy_cv': acc_cv, 
                                                    'Precision_cv':prec_cv, 
                                                    'Recall_cv': recall_cv, 
                                                    'F1score_cv': f1score_cv,
                                                    'MCC_cv': mcc_cv}, ignore_index=True)<jupyter_output># Tuning hyper-parameters Random Forest ser_wp

Best parameters set found on development set:

{'max_depth': 50, 'max_features': 'sqrt', 'n_estimators': 200, 'n_jobs': -1, 'random_state': 10}
0.9420116388059961
roc_auc_cv:  0.9420116388059961
acc_cv:  0.8832032502752348
prec_cv:  0.9156911972881921
recall_cv:  0.8327317283012822
f1score_cv:  0.8710845947819701
mcc_cv:  0.7679060782911419
# Tuning hyper-parameters Random Forest aac_wp

Best parameters set found on development set:

{'max_depth': 20, 'max_features': 'sqrt', 'n_estimators': 200, 'n_jobs': -1, 'random_state': 10}
0.9423206051122271
roc_auc_cv:  0.9423206051122271
acc_cv:  0.8832888300784013
prec_cv:  0.9160226062327795
recall_cv:  0.8321905074174293
f1score_cv:  0.8713093123197407
mcc_cv:  0.7679423909102394<jupyter_text>**ExtraTreesClassifier**<jupyter_code># A different model for each feature matrix
for i, j in Best2featmat.items():

    # Data splitting 
    X_train, X_test, y_train, y_test = train_test_split(j, y, test_size=0.2, 
                                                        random_state=42,
                                                        stratify=y)

    # Set the grid of parameters of GridSearchCV
    tuned_parameters_rf = [
        {"max_depth": [10, 20, 30, 50], "n_estimators": [10, 30, 50, 100, 200], 
          "max_features": ["auto", "sqrt", "log2", 0.2], "n_jobs":[-1], "random_state": [10]}
    ]

    print(f'# Tuning hyper-parameters ExtraTreesClassifier {i}')
    print()

    # Scoring metrics to be calculated
    scoring = {"AUC": "roc_auc", "Accuracy": "accuracy", 
               "Precision": "precision", "Recall": "recall", 
               "F1score": "f1", "MCC": make_scorer(matthews_corrcoef)}

    # Create GridSearchCV object with its parameters 
    etc = GridSearchCV(ExtraTreesClassifier(), tuned_parameters_rf, scoring=scoring, 
                      cv=3, refit="AUC", return_train_score=True)
    
    # Train models 
    etc.fit(X_train, y_train)

    print("Best parameters set found on development set:")
    print()
    print(etc.best_params_)
    print(etc.best_score_)

    # Obtain performance metrics of the best model  
    roc_auc_cv = max(etc.cv_results_["mean_test_AUC"])
    print("roc_auc_cv: ", roc_auc_cv)

    acc_cv = max(etc.cv_results_["mean_test_Accuracy"])
    print("acc_cv: ", acc_cv)

    prec_cv = max(etc.cv_results_["mean_test_Precision"])
    print("prec_cv: ", prec_cv)

    recall_cv = max(etc.cv_results_["mean_test_Recall"])
    print("recall_cv: ", recall_cv)

    f1score_cv = max(etc.cv_results_["mean_test_F1score"])
    print("f1score_cv: ", f1score_cv)

    mcc_cv = max(etc.cv_results_["mean_test_MCC"])
    print("mcc_cv: ", mcc_cv)

    # Add information to the results dataframe
    results_bestmodels = results_bestmodels.append({'Dataset': i, 'Model': 'ExtraTreesClassifier', 
                                                    'ROC_AUC_cv': roc_auc_cv,
                                                    'Accuracy_cv': acc_cv, 
                                                    'Precision_cv':prec_cv, 
                                                    'Recall_cv': recall_cv, 
                                                    'F1score_cv': f1score_cv,
                                                    'MCC_cv': mcc_cv}, ignore_index=True)
results_bestmodels
results_bestmodels.to_csv("results_bestMLclassifiers.csv")<jupyter_output><empty_output><jupyter_text>**Evaluation of the best model on test and external datasets**According to the summary tables presented above, the best model was **ExtraTreesClassifier** with `max_depth` of 50 and `n_estimators` of 200 as parameters, and low variance `Amino Acid Composition (aac_wp)` molecular descriptor as feature matrix. Thus, we trained this model again on the entire training dataset (without creating the validation partition) and evaluated it on the **test and external datasets**. AMPs can have more than one activity, including antibacterial, antifungal, antiparasitic, antiviral, among others. Test partition has active AMPs with a single activity, while the external partition has AMPs with more than one activity and it represents a more real scenario of virtual screening with much more non active AMPs than the active ones. You can find more details about these datasets in this [article](https://pubs.acs.org/doi/10.1021/acs.jcim.1c00251).<jupyter_code># Training dataset without validation partition
X = Best2featmat["aac_wp"]

# Model parameters
max_depth = 50
n_estimators = 200

# Implementation of the best model 
etc_def = ExtraTreesClassifier(max_depth=max_depth, n_estimators=n_estimators, n_jobs=-1, 
                            random_state=10).fit(X, y)
# Download test and external datasets
# Active AMPs from test dataset
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/TS_starPep_active_AMPs.fasta

# Non-active AMPs from test dataset
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/TS_starPep_nonactive_AMPs.fasta

# Active AMPs from external dataset
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/EX_starPep_active_AMPs.fasta

# Non-active AMPs from external dataset
! wget https://github.com/sayalaruano/CapstoneProject-MLZoomCamp/raw/main/Data/EX_starPep_nonactive_AMPs.fasta

# Create feature matrices of test and external datasets 
# For this, we will use the functions of the "Calculation of molecular descriptors 
# (features) with Pfeature" section

# Define variables with names of the active and non-active AMPs files from test 
# and external datasets
pos_test = 'TS_starPep_active_AMPs.fasta'
neg_test = 'TS_starPep_nonactive_AMPs.fasta'

pos_ext = 'EX_starPep_active_AMPs.fasta'
neg_ext = 'EX_starPep_nonactive_AMPs.fasta'

# Create test and external datasets
test = act_non_act_concat(pos_test, neg_test, aac_wp)
external = act_non_act_concat(pos_ext, neg_ext, aac_wp)

# Export feature matrices fro test and external datasets
test.to_csv("aac_wp_test.csv")
external.to_csv("aac_wp_external.csv")
# Obtain test and external y, and encoding them
y_test = test['class'].copy()
y_test = y_test.map({"positive": 1, "negative": 0})

y_ext = external['class'].copy()
y_ext = y_ext.map({"positive": 1, "negative": 0})

# Obtain test and external feature matrices
X_test = test.drop('class', axis=1)
X_ext = external.drop('class', axis=1)
# Predicting values with the best model on test and external datasets
y_test_pred = etc_def.predict(X_test)
y_ext_pred = etc_def.predict(X_ext)
# Dataset to store performance metrics
results_bestMLmodel = pd.DataFrame()
# Calculate performance metrics
roc_auc_ext = roc_auc_score(y_ext, y_ext_pred)
    
roc_auc_test = roc_auc_score(y_test, y_test_pred)

acc_ext = accuracy_score(y_ext, y_ext_pred)

acc_test = accuracy_score(y_test, y_test_pred)

prec_ext = precision_score(y_ext, y_ext_pred)

prec_test = precision_score(y_test, y_test_pred)

recall_ext = recall_score(y_ext, y_ext_pred)

recall_test = recall_score(y_test, y_test_pred)

f1score_ext = f1_score(y_ext, y_ext_pred)

f1score_test = f1_score(y_test, y_test_pred)

mcc_ext = matthews_corrcoef(y_ext, y_ext_pred)

mcc_test = matthews_corrcoef(y_test, y_test_pred)

# Add columns to the dataset
results_bestMLmodel["Performance_metric"] = ['ROC_AUC','Accuracy', 'Precision', 
                                             'Recall', 'F1score', 'MCC']

results_bestMLmodel["Test_dataset"] = [roc_auc_test, acc_test, prec_test, 
                                       recall_test, f1score_test, mcc_test]

results_bestMLmodel["External_dataset"] = [roc_auc_ext, acc_ext, prec_ext, 
                                       recall_ext, f1score_ext, mcc_ext]
# Show the evaluation results of the best model
results_bestMLmodel
results_bestMLmodel.to_csv("evaluation_bestMLclassifier.csv")<jupyter_output><empty_output><jupyter_text>**Feature importance of the best model**<jupyter_code># Retrieve feature importance from the RF model
importance = pd.Series(etc_def.feature_importances_, name = 'Gini')

# Retrieve feature names
feature_names = pd.Series(X.columns, name = 'Feature')
# Combine feature names and Gini values into a Dataframe
feat_import = pd.concat([feature_names, importance], axis=1, names=['Feature', 'Gini'])

# Sort the dataframe by Gini in descending order
feat_import_sorted = feat_import.sort_values('Gini', ascending=False)
# Plot of feature importance
plt.figure(figsize=(5, 10))
sns.set_theme(style="whitegrid")
ax = sns.barplot(x = 'Gini', y = 'Feature', data = feat_import_sorted)
plt.xlabel("Feature Importance")<jupyter_output><empty_output><jupyter_text>**Save the best model into a file**<jupyter_code># Define the output file's name
output_file = f'ExtraTreesClassifier_maxdepth{max_depth}_nestimators{n_estimators}'
output_file
# Export the binary file of the best model
with open(output_file, "wb") as f_out:
    pickle.dump(etc_def, f_out)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10501211387038159, 'fraction_numerical': 0.027710478497880074, 'mean_word_length': 4.192635303367842, 'pattern_counts': {'":': 28, '<': 106, '<?xml version=': 0, '>': 110, 'https://': 28, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 73, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '647372', 'license': 'MIT', 'n_tokens_mistral': 13962, 'n_tokens_neox': 12820, 'n_words': 3400, 'path': 'EDA_Binary_classifiers_AMPs_activity.ipynb', 'repo_name': 'sayalaruano/CapstoneProject-MLZoomCamp', 'size': '39639'}