3-Regressor-Performance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regressor Performance Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook guides through the analysis of an existing random forest regressor. The performance evaluation is based on the R^2 score from sklearn. The correlation of measured and predicted expression values is plotted. The feature importance from the random forest regression represent the contributions of each nucleotide-position to the prediction. They are extracted and visualized with a Logo-plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System initiation\n",
    "\n",
    "Loading all necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import logomaker as lm\n",
    "import matplotlib.pyplot as plt\n",
    "import joblib\n",
    "import pickle\n",
    "from math import sqrt\n",
    "from ExpressionExpert_Functions import Data_Src_Load, split_train_test, list_onehot, Insert_row_, my_CrossValScore, make_DataDir\n",
    "from sklearn.model_selection import cross_val_score, GroupShuffleSplit\n",
    "from sklearn.metrics import r2_score, mean_squared_error, make_scorer\n",
    "%matplotlib inline\n",
    "my_r2_score = make_scorer(r2_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variable setting\n",
    "\n",
    "We load the naming conventions from 'config.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Name_Dict = dict()\n",
    "with open('config.txt') as Conf:\n",
    "    myline = Conf.read().splitlines()\n",
    "    for line in myline:\n",
    "        if not line.startswith('#'):\n",
    "            (key, val) = line.split(':', 1)\n",
    "            Name_Dict[str(key.strip())] = val.strip()\n",
    "        \n",
    "\n",
    "Data_File = Name_Dict['Data_File']\n",
    "# extract the filename for naming of newly generated files\n",
    "File_Base = Name_Dict['File_Base']\n",
    "# the generated files will be stored in a subfolder with custom name\n",
    "Data_Folder = Name_Dict['Data_Folder']\n",
    "# column name of expression values\n",
    "Y_Col_Name = eval(Name_Dict['Y_Col_Name'])\n",
    "# Extracting entropy cutoff for removal of non-informative positions\n",
    "Entropy_cutoff = float(Name_Dict['Entropy_cutoff'])\n",
    "# figure file type\n",
    "Fig_Type = Name_Dict['Figure_Type']\n",
    "# Figure font size\n",
    "FigFontSize = Name_Dict['Figure_Font_Size']\n",
    "make_DataDir(Name_Dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data manipulation\n",
    "\n",
    "For the machine learning the data is first separated into training and test sets. The training set is used to generate a standard scaler for expression standardization to zero mean and unit variance. On each position the entropy is calculated to assess how much nucleotide diversity has been sampled on each position. If at any position the entropy is zero, i.e. only one nucleotide is present in all samples, this position is removed because it is non-informative for further analysis (Position entropy analysis). \n",
    "\n",
    "### Loading train and test set used for training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TrainTest_File = os.path.join(Data_Folder, '{}_{}_TrainTest-Data.pkl'.format(Name_Dict['ML_Date'], File_Base))\n",
    "TrainTest_Data = pickle.load(open(TrainTest_File,'rb'))\n",
    "SeqTrain, SeqTest = TrainTest_Data['Train'], TrainTest_Data['Test']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regressor analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random Forest Regressor file identifier for date\n",
    "ML_Date = Name_Dict['ML_Date']\n",
    "ML_Type = Name_Dict['ML_Regressor']\n",
    "\n",
    "# Number of independent promoter library measurements\n",
    "Measure_Numb = int(Name_Dict['Library_Expression'])\n",
    "ML_Best = dict()\n",
    "Expr_Scaler = dict()\n",
    "Y_train = np.empty((SeqTrain.shape[0], Measure_Numb))\n",
    "Y_train_pred = np.empty((SeqTrain.shape[0], Measure_Numb))\n",
    "Y_test = np.empty((SeqTest.shape[0], Measure_Numb))\n",
    "Y_test_pred = np.empty((SeqTest.shape[0], Measure_Numb))\n",
    "CVsplit = 25\n",
    "scores = np.empty((CVsplit, Measure_Numb))\n",
    "\n",
    "r2_train = np.empty(Measure_Numb)\n",
    "r2_test = np.empty(Measure_Numb)\n",
    "rmse_train = np.empty(Measure_Numb)\n",
    "rmse_test = np.empty(Measure_Numb)\n",
    "for Meas_Idx in range(Measure_Numb): \n",
    "    # loading correct ML regressor file and parameters for data preparation\n",
    "    Regressor_File = os.path.join(Data_Folder, '{}_{}_{}_{}-Regressor.pkl'.format(ML_Date, File_Base, Y_Col_Name[Meas_Idx].replace(' ','-'), ML_Type))\n",
    "    Parameter_File = os.path.join(Data_Folder, '{}_{}_{}_{}-Params.pkl'.format(ML_Date, File_Base, Y_Col_Name[Meas_Idx].replace(' ','-'), ML_Type))\n",
    "\n",
    "    try:\n",
    "        ML_DictName = '{}_Regressor'.format(Y_Col_Name[Meas_Idx])\n",
    "        ML_Best[ML_DictName] = joblib.load(Regressor_File)\n",
    "        # I assume the parameters have been generated in the same run as the regressor itself and is located in the same directory following the default naming scheme\n",
    "        Data_Prep_Params = pickle.load(open(Parameter_File,'rb'))\n",
    "        # extracting the positions that were removed because of insufficient information content\n",
    "        Positions_removed = Data_Prep_Params['Positions_removed']\n",
    "        # if the data was standardized we load the corresponding function\n",
    "        if eval(Name_Dict['Data_Standard']):\n",
    "            # extracting standard scaler from existing random forest regressor\n",
    "            # The standard scaler default name is the name of the expression measurement column with suffix: '_Scaler'\n",
    "            Scaler_DictName = '{}_Scaler'.format(Y_Col_Name[Meas_Idx])\n",
    "            Expr_Scaler[Scaler_DictName] = Data_Prep_Params[Scaler_DictName]\n",
    "    except FileNotFoundError:\n",
    "        print('Regressor file not found. Check parameter \"ML_Date\" in \"config.txt\"')\n",
    "\n",
    "    X_tmp = list_onehot(list(np.delete(np.array(list(SeqTrain['Sequence_label-encrypted'])),Positions_removed, axis=1)))\n",
    "    X_train = np.array(X_tmp).reshape(len(SeqTrain.index),-1)\n",
    "    AddFeat = eval(Name_Dict['Add_Feat'])\n",
    "    # adding the additional feature, here GC-content\n",
    "    X_train = np.append(X_train,SeqTrain[AddFeat].values, axis=1)\n",
    "    # activity prediction of training set with best random forest estimator\n",
    "    Y_train_pred[:,Meas_Idx] = ML_Best[ML_DictName].predict(X_train)\n",
    "    # correcting the prediction for standardized data\n",
    "    if eval(Name_Dict['Data_Standard']):\n",
    "        Y_train_pred[:,Meas_Idx] = Expr_Scaler[Scaler_DictName].inverse_transform(Y_train_pred[:,Meas_Idx])\n",
    "\n",
    "    # Test set prediction\n",
    "    # removing sequence positions that were missing in the feature vector for ml\n",
    "    # getting one-hot encodings from the original train and test data\n",
    "    X_tmp = list_onehot(list(np.delete(np.array(list(SeqTest['Sequence_label-encrypted'])),Positions_removed, axis=1)))\n",
    "    X_test = np.array(X_tmp).reshape(len(SeqTest.index),-1)\n",
    "    # adding the additional feature, here GC-content\n",
    "    X_test = np.append(X_test,SeqTest[AddFeat].values, axis=1)\n",
    "    # activity prediction of training set with best random forest estimator\n",
    "    Y_test_pred[:,Meas_Idx] = ML_Best[ML_DictName].predict(X_test)\n",
    "    # correcting the prediction for standardized data\n",
    "    if eval(Name_Dict['Data_Standard']):\n",
    "        Y_test_pred[:,Meas_Idx] = Expr_Scaler[Scaler_DictName].inverse_transform(Y_test_pred[:,Meas_Idx])\n",
    "\n",
    "    # corresponding observations scaled\n",
    "#     Scaler_DictName = '{}_Scaler'.format(Y_Col_Name[Meas_Idx])\n",
    "    Y_train[:, Meas_Idx] = SeqTrain[Y_Col_Name[Meas_Idx]].values\n",
    "    Y_test[:, Meas_Idx] = SeqTest[Y_Col_Name[Meas_Idx]].values\n",
    "\n",
    "    r2_train[Meas_Idx] = r2_score(Y_train[:, Meas_Idx], Y_train_pred[:, Meas_Idx])\n",
    "    r2_test[Meas_Idx] = r2_score(Y_test[:, Meas_Idx], Y_test_pred[:, Meas_Idx])\n",
    "    rmse_train[Meas_Idx] = sqrt(mean_squared_error(Y_train[:, Meas_Idx], Y_train_pred[:, Meas_Idx]))\n",
    "    rmse_test[Meas_Idx] = sqrt(mean_squared_error(Y_test[:, Meas_Idx], Y_test_pred[:, Meas_Idx]))\n",
    "    \n",
    "    # cross-validation scoring\n",
    "    cv = GroupShuffleSplit(n_splits=CVsplit, test_size=.1, random_state=42)\n",
    "    # if applicable correcting target variable according to the standardization\n",
    "    if eval(Name_Dict['Data_Standard']):\n",
    "        Y_train2 = np.ravel(Expr_Scaler[Scaler_DictName].transform(SeqTrain[Y_Col_Name[Meas_Idx]].values.reshape(-1, 1)))\n",
    "    else:\n",
    "        Y_train2 = Y_train[:, Meas_Idx]\n",
    "    groups = SeqTrain['Sequence_letter-encrypted']\n",
    "    scores[:,Meas_Idx] = cross_val_score(ML_Best[ML_DictName], X_train, Y_train2, groups=groups, scoring=my_r2_score, cv=cv)# , groups=groups, scoring=my_r2_score\n",
    "\n",
    "    print('Cross validation for ML-Type: {}'.format(ML_Type))\n",
    "    print('R2 Statistic: {:.2f} (+/-{:.2f})'.format(scores[:,Meas_Idx].mean(), scores[:,Meas_Idx].std()*2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance visualization\n",
    "### Calculation of model predictions\n",
    "\n",
    "Plot of predicted to measured expression strength for training and test data sets and R$^2$ correlation coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of independent promoter library measurements\n",
    "Expr_Unit = Name_Dict['Expression_Unit']\n",
    "\n",
    "for Meas_Idx in range(Measure_Numb): \n",
    "#     fig, axs = plt.subplots(nrows=1, ncols=1)\n",
    "\n",
    "    plt.scatter(Y_train[:, Meas_Idx], Y_train_pred[:, Meas_Idx], marker='+')\n",
    "    plt.scatter(Y_test[:, Meas_Idx], Y_test_pred[:, Meas_Idx], marker='o', c='r')\n",
    "    plt.title('{} of {}'.format(ML_Type, Y_Col_Name[Meas_Idx]), fontsize=18)\n",
    "    plt.xlabel('measured {}'.format(Expr_Unit), fontsize=FigFontSize)\n",
    "    plt.ylabel('predicted {}'.format(Expr_Unit), fontsize=FigFontSize)\n",
    "    plt.legend(['Training set, R$^2$={:.2f}'.format(r2_train[Meas_Idx]),'Test set, R$^2$={:.2f}'.format(r2_test[Meas_Idx])], loc='upper left', fontsize=FigFontSize)\n",
    "\n",
    "    # saving the figure\n",
    "    Fig_ID = Name_Dict['CorrPlot_File']\n",
    "    CorrPlot_File = os.path.join(Data_Folder, '{}_{}_{}_{}_{}.{}'.format(time.strftime('%Y%m%d'), File_Base, ML_Type, Fig_ID, Y_Col_Name[Meas_Idx].replace(' ','-'), Fig_Type))\n",
    "    plt.savefig(CorrPlot_File, bbox_inches='tight', format=Fig_Type)\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature importance for random forest\n",
    "### Importance of GC-content for prediction "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if ML_Type == 'RFR' or ML_Type == 'GBR':\n",
    "    Nr_EngFeat = len(eval(Name_Dict['Add_Feat']))\n",
    "    Nucleotide_Importance = dict()\n",
    "    AddFeat_Importance = dict()\n",
    "    GC_Importance = np.empty((Measure_Numb,1))\n",
    "\n",
    "    for Meas_Idx in range(Measure_Numb): \n",
    "        ML_DictName = '{}_Regressor'.format(Y_Col_Name[Meas_Idx])\n",
    "        #number of engineered additional features\n",
    "    #     FI_DictName = '{}_FI'.format(Y_Col_Name[Meas_Idx])\n",
    "    #     AddFI_DictName = '{}_FI'.format(Y_Col_Name[Meas_Idx])\n",
    "        Nucleotide_Importance[Meas_Idx] = np.array(ML_Best[ML_DictName].feature_importances_[0:-Nr_EngFeat]).reshape(-1,4)\n",
    "        AddFeat_Importance = ML_Best[ML_DictName].feature_importances_[-Nr_EngFeat:]\n",
    "    #     print(Nucleotide_Importance)\n",
    "        MeasName = '{}_FI'.format(Y_Col_Name[Meas_Idx])\n",
    "    #     Feat_Nucl = [len(i) for i in Nucleotide_Importance[Measure_Numb]]\n",
    "    #     Feat_All = Feat_Nucl[Meas_Idx]*4 + len(eval(Name_Dict['Add_Feat']))\n",
    "        FeatList_All = np.append(np.hstack(Nucleotide_Importance[Meas_Idx]),AddFeat_Importance)\n",
    "        GC_Importance[Meas_Idx] = len(FeatList_All) - np.arange(len(FeatList_All))[np.argsort(FeatList_All)==len(FeatList_All)-1]\n",
    "        print('Importance of GC-content in {}: {}'.format(Y_Col_Name[Meas_Idx], GC_Importance[Meas_Idx]))\n",
    "else:\n",
    "    print('GC-content importance available only for decision tree methods')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequence logo of nucleotide importance\n",
    "\n",
    "The feature importance of the random forest regressor, i.e. the y-axis in the Logo-plot, is normalized to sum over all nucleotide-positions to one.\n",
    "\n",
    "The logos are generated with [Logomaker](https://logomaker.readthedocs.io/en/latest/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if ML_Type == 'RFR' or ML_Type == 'GBR':\n",
    "    # Number of independent promoter library measurements\n",
    "    for Meas_Idx in range(Measure_Numb): \n",
    "        MeasName = '{}_FI'.format(Y_Col_Name[Meas_Idx])\n",
    "        PWM_tmp = pd.DataFrame(Nucleotide_Importance[Meas_Idx], columns=['A','C','G','T'])\n",
    "        PWM_best = Insert_row_(Positions_removed, PWM_tmp, np.zeros([len(Positions_removed),4]))\n",
    "        nn_logo = lm.Logo(PWM_best)\n",
    "        nn_logo.ax.set_xlabel('Position', fontsize=FigFontSize)\n",
    "        nn_logo.ax.set_ylabel(r'$\\frac{Importance_i}{\\sum Importance}$', fontsize=FigFontSize)\n",
    "        nn_logo.ax.set_title('{} of {}'.format(ML_Type, Y_Col_Name[Meas_Idx]), fontsize=FigFontSize)\n",
    "\n",
    "        # saving the figure\n",
    "        Fig_ID = Name_Dict['LogoPlot_File']\n",
    "        LogoPlot_File = os.path.join(Data_Folder, '{}_{}_{}_{}_{}.{}'.format(time.strftime('%Y%m%d'), File_Base, ML_Type, Fig_ID, Y_Col_Name[Meas_Idx].replace(' ','-'), Fig_Type))\n",
    "        plt.savefig(LogoPlot_File, bbox_inches='tight', format=Fig_Type)\n",
    "\n",
    "        plt.show()\n",
    "else:\n",
    "    print('Nucleotide importance available only for decision tree methods')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}