Feature_Generation_for_Protein_peptide_Sequences.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMIQF34ci5ZYiCI8SHXLdyT",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/shilpasy/Building-and-Interpreting-Random-Forest-based-Cell-Penetrating-Peptide-Prediction-Model/blob/master/Feature_Generation_for_Protein_peptide_Sequences.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Simple feature calculations for protein and peptide sequences to use in AI/ML related projects. These can be a basic set of features on which you can develop further depending on the problem of interest."
      ],
      "metadata": {
        "id": "nZCYYJPSqsAb"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#!pip install biopython"
      ],
      "metadata": {
        "id": "Y1XSa-E6vQ2u"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "import random\n",
        "import math\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import itertools\n",
        "from collections import namedtuple\n",
        "from Bio.Seq import Seq\n",
        "from Bio.SeqUtils.ProtParam import ProteinAnalysis # you can install this library using pip\n",
        "from matplotlib import pyplot as plt\n",
        "import seaborn as sns\n",
        "%matplotlib inline\n",
        "import warnings\n",
        "#warnings.filterwarnings('ignore')\n",
        "#pd.set_option('display.max_rows', 1000)"
      ],
      "metadata": {
        "id": "qsCgwUF7vHvS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Your input file i.e. training data which has sequences and target-labels is \"train\" here and the way in which you call the function is all_features_df = generate_features(train). This dataframe will have the calculated features. Features calculated here: amino acid frequencies, dipeptide frequencies, tripeptide frequencies, common biochemical properties"
      ],
      "metadata": {
        "id": "ztCvrvycrqqc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def calculate_aa_frequencies(sequence, peptides):\n",
        "    aa_frequencies = [sequence.count(aa) / len(sequence) for aa in peptides]\n",
        "    return aa_frequencies\n",
        "\n",
        "def generate_features(train):\n",
        "    dipeptide = [''.join(x) for x in itertools.product('ACDEFGHIKLMNPQRSTVWY', repeat=2)]\n",
        "    # tripeptide = [''.join(y) for y in itertools.product('ACDEFGHIKLMNPQRSTVWY', repeat=3)]\n",
        "    biochem_features = [\"net_charge\", \"IsoelectricPoint\", \"Hbdonors\", \"Hbacceptors\", \"HbDonorAcceptorDiff\", \"hydropathy\", \"MolWt\", \"SS_B\", \"SS_T\", \"SS_H\"]\n",
        "\n",
        "    features = []\n",
        "\n",
        "    for i, row in train.iterrows():\n",
        "        sequence = row['protein_sequence']\n",
        "        labelclass = row['tm']\n",
        "\n",
        "        aa_frequencies = calculate_aa_frequencies(sequence, 'ACDEFGHIKLMNPQRSTVWY')\n",
        "        #dipeptide_count = calculate_aa_frequencies(sequence, dipeptide)\n",
        "       # tripeptide_count = calculate_aa_frequencies(sequence, tripeptide)\n",
        "\n",
        "        # Biochemical properties features\n",
        "        eachSeq = sequence.replace('B', 'N').replace('Z', 'Q').replace('J', 'I').replace('X', 'G').replace('m', '')\n",
        "        positive_aa_R = eachSeq.count('R')\n",
        "        positive_aa_K = eachSeq.count('K')\n",
        "        negative_aa_D = eachSeq.count('D')\n",
        "        negative_aa_E = eachSeq.count('E')\n",
        "        aromaticnum_W = eachSeq.count('W')\n",
        "        aromaticnum_Y = eachSeq.count('Y')\n",
        "        aromaticnum_F = eachSeq.count('F')\n",
        "        count_aa_H = eachSeq.count('H'); count_aa_N = eachSeq.count('N'); count_aa_Q = eachSeq.count('Q'); count_aa_S = eachSeq.count('S'); count_aa_T = eachSeq.count('T');\n",
        "        net_charge = (positive_aa_R + positive_aa_K - negative_aa_D - negative_aa_E) / len(eachSeq)\n",
        "        PI = ProteinAnalysis(eachSeq).isoelectric_point()\n",
        "        Hbdonors = (positive_aa_R * 5 + count_aa_Q * 2 + count_aa_S + count_aa_T + count_aa_H * 2 +\n",
        "                    positive_aa_K * 3 + aromaticnum_W + aromaticnum_Y + count_aa_N * 2) / len(eachSeq)\n",
        "        Hbacceptors = (count_aa_N * 2 + negative_aa_D * 4 + count_aa_Q * 2 + negative_aa_E * 4 +\n",
        "                       count_aa_H * 2 + count_aa_S * 2 + count_aa_T * 2 + aromaticnum_Y) / len(eachSeq)\n",
        "        Hbdonaccdiff = (Hbdonors - Hbacceptors) / len(eachSeq)\n",
        "        SS = ProteinAnalysis(eachSeq).secondary_structure_fraction()  # [Sheet, Turn, Helix]\n",
        "\n",
        "        hydropathy = ProteinAnalysis(eachSeq).gravy()\n",
        "        MW = ProteinAnalysis(eachSeq).molecular_weight()\n",
        "        MW = MW / len(eachSeq)\n",
        "\n",
        "        # Combine all features\n",
        "        #all_features = aa_frequencies + dipeptide_count + [ .... ]\n",
        "        all_features = aa_frequencies + [net_charge, PI, Hbdonors, Hbacceptors, Hbdonaccdiff, hydropathy, MW, SS[0], SS[1], SS[2]]\n",
        "        all_features.append(labelclass)\n",
        "        features.append(all_features)\n",
        "    #columns = ['aa_' + aa for aa in 'ACDEFGHIKLMNPQRSTVWY'] + ['dipeptide_' + dp for dp in dipeptide] + [f\"biochem_{feature}\" for feature in biochem_features] + ['LABEL'] # # ['tripeptide_' + tp for tp in tripeptide] +\n",
        "    columns = ['aa_' + aa for aa in 'ACDEFGHIKLMNPQRSTVWY'] + [f\"biochem_{feature}\" for feature in biochem_features] + ['LABEL']\n",
        "\n",
        "    features_df = pd.DataFrame(features, columns=columns)\n",
        "\n",
        "    return features_df\n",
        "\n",
        "###### Generate all features ########\n",
        "all_features_df = generate_features(train)"
      ],
      "metadata": {
        "id": "PL9cBfYgs6Cg"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}