Peptide_Prediction_BigQuery_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo: Peptide Prediction with BigQuery ML\n",
    "\n",
    "### Supporting research: Covid19 and beyond for vaccine candidates \n",
    "\n",
    "BigQuery ML enables users to create and execute machine learning models in BigQuery using SQL queries. The goal is to democratize machine learning by enabling SQL practitioners to build models using their existing tools and to increase development speed by eliminating the need for data movement.\n",
    "In this tutorial, you use the sample Covid19 dataset for BigQuery.\n",
    "Comments & Feedback @jigmehta\n",
    "## Objectives\n",
    "In this tutorial, you use:\n",
    "+ BigQuery `ML.CREATE` to create a classification model using the `CREATE MODEL` statement\n",
    "+ The `ML.EVALUATE` function to evaluate the ML model\n",
    "+ Use `ML.TRANSFORM`feature engineering functions to improve model performance\n",
    "+ The `ML.PREDICT` function to make predictions using the ML model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jupyter magics are notebook-specific shortcuts that allow you to run commands with minimal syntax. Jupyter notebooks come with many [built-in commands](https://ipython.readthedocs.io/en/stable/interactive/magics.html). The BigQuery client library, `google-cloud-bigquery`, provides a cell magic, `%%bigquery`. The `%%bigquery` magic runs a SQL query and returns the results as a pandas `DataFrame`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Explore Peptidic Epitope Data\n",
    "For our coronavirus subset of data, let's examine data first:\n",
    "+ Following query shows information about antigen.\n",
    "+ Observe specific antigen that might of more interest to our study for coronavirus.\n",
    "+ Can you find it?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read GCP project id from env.\n",
    "shell_output=!gcloud config list --format 'value(core.project)' 2>/dev/null\n",
    "GCP_PROJECT_ID=shell_output[0]\n",
    "print(\"GCP project ID:\" + GCP_PROJECT_ID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define project__id variable\n",
    "project_id = 'covid-19-271622'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT * \n",
    "FROM `corona.antigen`\n",
    "LIMIT 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually inspect number of epitopes by antigen name in our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery epitopes --project $project_id\n",
    "SELECT\n",
    " Antigen_Name,\n",
    " sum(__Epitopes) AS epitope\n",
    "FROM `corona.antigen`\n",
    "GROUP BY Antigen_Name\n",
    "ORDER BY epitope DESC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epitopes.plot(kind='bar', x='Antigen_Name', y='epitope');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Following query shows epitope detail for those linear peptides whose parent protein\n",
    "# came from Spike glycoprotein, which is our interests to investigate for possible antigen\n",
    "# to protect against coronavirus\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT Epitope_ID, Object_Type, Description, Starting_Position,\n",
    "            Ending_Position, Antigen_Name,Parent_Protein, Organism_Name\n",
    "FROM `corona.epitope` \n",
    "WHERE Parent_Protein = 'Spike glycoprotein'\n",
    "LIMIT 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the resulting query, we have a set of peptides generated from parent antigen protein. Now, we need to identify how these peptides bind with HLA molecule.\n",
    "\n",
    "Following query shows HLA-peptide binding affinity information. The goal is to leverage machine learning models to predict binding affinity for many different and new peptides with HLA MHC Class-I molecules so that testing for vaccine candidates can be prioritized to accelerate solutions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT Reference_ID, Description, Allele_Name,MHC_allele_class, Assay_Group,\n",
    "       Qualitative_Measure, Quantitative_measurement \n",
    "FROM `corona.mhc_ligand`\n",
    "WHERE Epitope_ID = 12967 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our dataset there are multiple of alleles and, for each allele we have epitope binding affinity measure. Lete inspect how epitopes and allele bindings are distributed in the dataset. This is one way to short list possible allele to consider for peptide binding research.\n",
    "\n",
    "We will use dataframe to hold query result to build our visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from google.cloud import bigquery\n",
    "client = bigquery.Client()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql = \"\"\"\n",
    "SELECT Qualitative_Measure, \n",
    "       Allele_Name, \n",
    "       count(1) as count \n",
    "FROM `corona.mhc_ligand`\n",
    "GROUP BY  1,2\n",
    "ORDER BY 3 DESC\n",
    "\"\"\"\n",
    "df = client.query(sql).to_dataframe()\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pivot_table = df.pivot(index='Allele_Name', columns='Qualitative_Measure', values='count')\n",
    "pivot_table.plot(kind='bar', stacked=True, figsize=(15, 7));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on data insepction, it seems like for coronavirus data set, HLA-A-01-01, HLA-A-02-01, HLA-A-02-02, HLA-A-02-03 and HLA-A-020-06 alleles are most likely candidates for MHC-I molecules.\n",
    "\n",
    "It is important to narrow down focus peptides and HLA binding to speed up testing on most probable vaccine candidate. Lets identify set of peptides that binds well with HLA, you can further narrow it to specific HLA as well. Following qury execution suggest that peptides with length of 9 and 10 are having better binding affinity. We should focus our research on those for now. Of course, more factors, features as well as 3D structure of aplha/beta components of HLA+Peptide bindings are useful for more complex modeling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql = \"\"\"\n",
    "SELECT Qualitative_Measure, \n",
    "       length(Description) as Peptide_Length, \n",
    "       count(1) as count \n",
    "FROM `corona.mhc_ligand`\n",
    "GROUP BY  1,2\n",
    "ORDER BY 3 DESC\n",
    "\"\"\"\n",
    "df = client.query(sql).to_dataframe()\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pivot_table = df.pivot(index='Peptide_Length', columns='Qualitative_Measure', values='count')\n",
    "pivot_table.plot(kind='bar', stacked=False, figsize=(15, 7));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets build feature table that translate Peotide Amino Acid sqquence into positioning columns. Explore amino acid properties table, it is prepared with allocation of acid sequnece number as well as additional properites that you can leverage for cross feature engineering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT *\n",
    "FROM `hla_peptide_generic.amino_acid`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will examine pre-built feature table called mhc_qual_feature which has column presentation of 9 and 10 letter (MER) peptide sequence which are most suitable for MHC-1 binding. Table also contains qualitative measure and quantitative measure.\n",
    "\n",
    "+ Peptides are transformed into column A1 to A10 based on amino acid position value in a chain\n",
    "+ result_score is quantitative measure: lower the number, higher binding affinity\n",
    "+ result columns is binary classification based on 5 class result_quality score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT *\n",
    "FROM `corona.mhc_qual_feature`\n",
    "LIMIT 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Build First Classification Model\n",
    "Building ML models with BigQuery is as simple as writing SQL statements; makes ML modeling accessible to even SQL developers and analysts. We will create a simple classification model to predict for a given peptide if there is strong binding affinity with certain HLA Allele.\n",
    "\n",
    "Following statement creates a classification model using logistic regression. We are selecting feature columns of Allele and amino acid positioning within a peptide to classify if a peptide is a good candidate for vaccine testing.\n",
    "\n",
    "+ Note: we are filtering data for peptides with length of 9 or 10 mers only. Also, since we run multiple samples, we are randomizing samples by 80% of data for learning.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "CREATE OR REPLACE MODEL `corona.Classification_model_bq1`\n",
    "OPTIONS\n",
    "(\n",
    "model_type='logistic_reg',\n",
    "input_label_cols=['result_quality']\n",
    ")\n",
    "AS\n",
    "SELECT\n",
    " result_quality, Allele, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10\n",
    " FROM\n",
    "  `corona.mhc_qual_feature`\n",
    " WHERE length(peptide) IN (9,10)\n",
    " AND rand() < 0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The query takes several minutes to complete. After the first iteration is\n",
    "complete, your model (`Classification_model_bq1`) appears in the navigation panel of the\n",
    "BigQuery web UI. Because the query uses a `CREATE MODEL` statement to create a\n",
    "table, you do not see query results. The output is an empty `DataFrame`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get training statistics\n",
    "\n",
    "To see the results of the model training, you can use the\n",
    "[`ML.TRAINING_INFO`](https://cloud.google.com/bigquery/docs/reference/standard-sql/bigqueryml-syntax-train)\n",
    "function, or you can view the statistics in the BigQuery web UI.\n",
    "In this tutorial, you use the `ML.TRAINING_INFO` function.\n",
    "\n",
    "A machine learning algorithm builds a model by examining many examples and\n",
    "attempting to find a model that minimizes loss. This process is called empirical\n",
    "risk minimization.\n",
    "\n",
    "Loss is the penalty for a bad prediction &mdash; a number indicating\n",
    "how bad the model's prediction was on a single example. If the model's\n",
    "prediction is perfect, the loss is zero; otherwise, the loss is greater. The\n",
    "goal of training a model is to find a set of weights that have low\n",
    "loss, on average, across all examples.\n",
    "\n",
    "To see the model training statistics that were generated when you ran the\n",
    "`CREATE MODEL` query:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT\n",
    "  *\n",
    "FROM\n",
    "  ML.TRAINING_INFO(MODEL `corona.Classification_model_bq1`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery loss_curve --project $project_id\n",
    "SELECT\n",
    "  iteration, loss, eval_loss\n",
    "FROM\n",
    "  ML.TRAINING_INFO(MODEL `corona.Classification_model_bq1`)\n",
    "ORDER BY iteration ASC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_curve.plot(x='iteration');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `loss` column represents the loss metric calculated after the given iteration\n",
    "on the training dataset. Since you performed a logistic regression, this column\n",
    "is the [log loss](https://en.wikipedia.org/wiki/Cross_entropy#Cross-entropy_error_function_and_logistic_regression).\n",
    "The `eval_loss` column is the same loss metric calculated on\n",
    "the holdout dataset (data that is held back from training to validate the model).\n",
    "\n",
    "For more details on the `ML.TRAINING_INFO` function, see the\n",
    "[BigQuery ML syntax reference](https://cloud.google.com/bigquery/docs/reference/standard-sql/bigqueryml-syntax-train)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate your model\n",
    "\n",
    "After creating your model, you evaluate the performance of the classifier using\n",
    "the [`ML.EVALUATE`](https://cloud.google.com/bigquery/docs/reference/standard-sql/bigqueryml-syntax-evaluate)\n",
    "function. You can also use the [`ML.ROC_CURVE`](https://cloud.google.com/bigquery/docs/reference/standard-sql/bigqueryml-syntax-roc)\n",
    "function for logistic regression specific metrics.\n",
    "\n",
    "A classifier is one of a set of enumerated target values for a label. For\n",
    "example, in this tutorial you are using a binary classification model that\n",
    "detects transactions. The two classes are the values in the `label` column:\n",
    "`0` (no transactions) and not `1` (transaction made).\n",
    "\n",
    "To run the `ML.EVALUATE` query that evaluates the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT\n",
    "  *\n",
    "FROM ML.EVALUATE(MODEL `corona.Classification_model_bq1`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT roc_auc,\n",
    "       CASE WHEN roc_auc > .85 THEN 'good'\n",
    "            WHEN roc_auc > .7 THEN 'fair'\n",
    "            WHEN roc_auc > .5 THEN 'not great'\n",
    "            ELSE 'poor' END AS model_quality\n",
    "FROM ML.EVALUATE(MODEL `corona.Classification_model_bq1`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Improve Model Performance with Feature Engineering\n",
    "BigQuery offers many [transform / preprocessing](https://cloud.google.com/bigquery-ml/docs/reference/standard-sql/bigqueryml-preprocessing-functions) functions for feature engineering on data. Advantage of transform functions is that once you build a model with preprocessing as part of model definition, prediction data does not need to be prepared as the model will apply transformation for the input. Lets see one example of a transform feature and rebuild our model to check if we get better model performance.\n",
    "\n",
    "Following statement will create another classification model with preprocessing of the result score to normalize its deviation with respect to min-max value of an attribute.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "CREATE OR REPLACE MODEL `corona.Classification_model_bq2`\n",
    "TRANSFORM (result_quality, Allele, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, \n",
    "           ML.MIN_MAX_SCALER(result_score) OVER() AS rs\n",
    ")\n",
    "OPTIONS\n",
    "(\n",
    "model_type='logistic_reg',\n",
    "input_label_cols=['result_quality']\n",
    ")\n",
    "AS\n",
    "SELECT result_quality, Allele, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, result_score\n",
    "FROM `corona.mhc_qual_feature`\n",
    "WHERE length(peptide) IN (9,10)\n",
    "AND rand() < 0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The query takes several minutes to complete. After the first iteration is\n",
    "complete, your model (`Classification_model_bq2`) appears in the navigation panel of the\n",
    "BigQuery web UI. The output is an empty `DataFrame`.\n",
    "\n",
    "Check model performance after feature engineering, run the `ML.EVALUATE` query that evaluates the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery loss_curve2 --project $project_id\n",
    "SELECT\n",
    "  iteration, loss, eval_loss\n",
    "FROM\n",
    "  ML.TRAINING_INFO(MODEL `corona.Classification_model_bq2`)\n",
    "ORDER BY iteration ASC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_curve2.plot(x='iteration');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT roc_auc,\n",
    "       CASE WHEN roc_auc > .85 THEN 'good'\n",
    "            WHEN roc_auc > .7 THEN 'fair'\n",
    "            WHEN roc_auc > .5 THEN 'not great'\n",
    "            ELSE 'poor' END AS model_quality\n",
    "FROM ML.EVALUATE(MODEL `corona.Classification_model_bq2`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see model performance has improved!\n",
    "\n",
    "With BigQuery you can take advantage of an already available highly powerful computer data processing and analysis platform to build a machine learning model, without moving your data! You can learn more about BQML here. For our data set, you can build a DNN_Regression model to predict the result_score for a peptide HLA binding. Try that as a practice!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Run Prediction on BQML Model\n",
    "Now that you have evaluated your model, the next step is to use it to predict\n",
    "outcomes. \n",
    "\n",
    "To run the query that uses the model to predict the number of transactions:\n",
    "Following example demonstrate leveraging BQ model for prediction. Optionally, you can export model and publish it on to Google AI Platform for serving prediction.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bigquery --project $project_id\n",
    "SELECT\n",
    "  predicted_result_quality, predicted_result_quality_probs, result_quality as original_result, result_score\n",
    "FROM ML.PREDICT(MODEL `corona.Classification_model_bq2`, (\n",
    "  SELECT result_quality, Allele, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, result_score\n",
    "  FROM `corona.mhc_qual_feature`\n",
    "  WHERE length(peptide) IN (9,10)\n",
    "  AND rand() < 0.0009))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result shows predicted quality class with confidence. You can compare that with original result. Next step is to operationalize ML pipeline so that you can efficiently perform data updates and model updates. Check out AI Pipeline example for peptide prediction to learn more!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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
}