04_Cider_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7b1c073a",
   "metadata": {},
   "source": [
    "# Cider analysis on the full proteome\n",
    "\n",
    "<a rel=\"license\" href=\"https://creativecommons.org/licenses/by-sa/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://licensebuttons.net/l/by-sa/4.0/88x31.png\" title='This work is licensed under a Creative Commons Attribution 4.0 International License.' align=\"right\"/></a>\n",
    "\n",
    "Author: Dr Antonia Mey   \n",
    "Email: antonia.mey@ed.ac.uk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "51210d5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import glob \n",
    "import os\n",
    "import urllib\n",
    "from collections import Counter\n",
    "import urllib\n",
    "import json\n",
    "\n",
    "# import the relevant code\n",
    "from localcider.sequenceParameters import SequenceParameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0ecb5f49",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_all_data_from_seq(seq):\n",
    "    ''' This extracts all data needed from cider for a given sequence and returns it as a dictionary\n",
    "    '''\n",
    "    SeqOb = SequenceParameters(seq)\n",
    "    data = {}\n",
    "    data['N'] = SeqOb.get_length()\n",
    "    data['seq'] = seq\n",
    "    data['f-'] = SeqOb.get_fraction_negative()\n",
    "    data['f+'] = SeqOb.get_fraction_positive()\n",
    "    data['FCR']= SeqOb.get_FCR()\n",
    "    data['NCPR'] = SeqOb.get_NCPR()\n",
    "    data['kappa'] = SeqOb.get_kappa()\n",
    "    data['Omega'] = SeqOb.get_Omega()\n",
    "    data['sigma'] = get_sigma(data['f+'], data['f-'])\n",
    "    data['delta'] = SeqOb.get_delta()\n",
    "    data['delta max'] = SeqOb.get_deltaMax()\n",
    "    data['hydropathy'] = SeqOb.get_mean_hydropathy()\n",
    "    data['Phase plot region'] = int(SeqOb.get_phasePlotRegion())\n",
    "    data['Hydrophobic count'] = int(get_num_hydrophobic(seq))\n",
    "    data['Positive count'] = int(SeqOb.get_countPos())\n",
    "    data['Negative count'] = int(SeqOb.get_countNeg())\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2ecdab83",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_domainn_region_info(data):\n",
    "    known_things = {}\n",
    "    known_things['Domain'] = []\n",
    "    known_things['Region'] = []\n",
    "    known_things['Motif'] = []\n",
    "    #known_things['Helix'] = []\n",
    "    if 'features' in data.keys():\n",
    "        for d in data['features']:\n",
    "            #print(d['type'])\n",
    "            if d['type'] == 'Domain':\n",
    "                known_things['Domain'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "            elif d['type'] == 'Region':\n",
    "                known_things['Region'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "            elif d['type'] == 'Motif':\n",
    "                known_things['Motif'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "            #elif d['type'] == 'Helix':\n",
    "            #    known_things['Helix'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "\n",
    "            \n",
    "    return known_things       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "id": "5aec0526",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sigma(f_plus, f_minus):\n",
    "    '''gets sigma\n",
    "    Sigma: ⟨σ⟩, or the average sigma value for the entire sequence, where sigma quantifies the charge asymmetry. \n",
    "Specifically, σ = (f+-f-)2/(f++f-) where f- and f+ refer to the fraction of negative and positive residues across the entire sequence. See Das & Pappu[2] for more details\n",
    "    \n",
    "    '''\n",
    "    if f_plus == 0 and f_minus==0:\n",
    "        sigma = np.nan\n",
    "    else:\n",
    "        sigma = (f_plus-f_minus)**2/(f_plus+f_minus)\n",
    "    return sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "269dac7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_num_hydrophobic(seq):\n",
    "    ''' counts the number of hydrophobic amino acids (L,V,I,M) in a seq\n",
    "    '''\n",
    "    \n",
    "    I_count = seq.count('I')\n",
    "    M_count = seq.count('M')\n",
    "    L_count = seq.count('L')\n",
    "    V_count = seq.count('V')\n",
    "    count = I_count+M_count+L_count+V_count\n",
    "    return count\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "78bf1925",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniprot_ids = glob.glob('unique_ids/*')\n",
    "protein = uniprot_ids[20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "af0871e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(protein)\n",
    "data = json.load(f)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "ede76a66",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "uni_id P46940\n",
      "gene_name IQGAP1\n",
      "ref_seq_id NM_003870\n",
      "sequence MSAADEVDGLGVARPHYGSVLDNERLTAEEMDERRRQNVAYEYLCHLEEAKRWMEACLGEDLPPTTELEEGLRNGVYLAKLGNFFSPKVVSLKKIYDREQTRYKATGLHFRHTDNVIQWLNAMDEIGLPKIFYPETTDIYDRKNMPRCIYCIHALSLYLFKLGLAPQIQDLYGKVDFTEEEINNMKTELEKYGIQMPAFSKIGGILANELSVDEAALHAAVIAINEAIDRRIPADTFAALKNPNAMLVNLEEPLASTYQDILYQAKQDKMTNAKNRTENSERERDVYEELLTQAEIQGNINKVNTFSALANIDLALEQGDALALFRALQSPALGLRGLQQQNSDWYLKQLLSDKQQKRQSGQTDPLQKEELQSGVDAANSAAQQYQRRLAAVALINAAIQKGVAEKTVLELMNPEAQLPQVYPFAADLYQKELATLQRQSPEHNLTHPELSVAVEMLSSVALINRALESGDVNTVWKQLSSSVTGLTNIEEENCQRYLDELMKLKAQAHAENNEFITWNDIQACVDHVNLVVQEEHERILAIGLINEALDEGDAQKTLQALQIPAAKLEGVLAEVAQHYQDTLIRAKREKAQEIQDESAVLWLDEIQGGIWQSNKDTQEAQKFALGIFAINEAVESGDVGKTLSALRSPDVGLYGVIPECGETYHSDLAEAKKKKLAVGDNNSKWVKHWVKGGYYYYHNLETQEGGWDEPPNFVQNSMQLSREEIQSSISGVTAAYNREQLWLANEGLITRLQARCRGYLVRQEFRSRMNFLKKQIPAITCIQSQWRGYKQKKAYQDRLAYLRSHKDEVVKIQSLARMHQARKRYRDRLQYFRDHINDIIKIQAFIRANKARDDYKTLINAEDPPMVVVRKFVHLLDQSDQDFQEELDLMKMREEVITLIRSNQQLENDLNLMDIKIGLLVKNKITLQDVVSHSKKLTKKNKEQLSDMMMINKQKGGLKALSKEKREKLEAYQHLFYLLQTNPTYLAKLIFQMPQNKSTKFMDSVIFTLYNYASNQREEYLLLRLFKTALQEEIKSKVDQIQEIVTGNPTVIKMVVSFNRGARGQNALRQILAPVVKEIMDDKSLNIKTDPVDIYKSWVNQMESQTGEASKLPYDVTPEQALAHEEVKTRLDSSIRNMRAVTDKFLSAIVSSVDKIPYGMRFIAKVLKDSLHEKFPDAGEDELLKIIGNLLYYRYMNPAIVAPDAFDIIDLSAGGQLTTDQRRNLGSIAKMLQHAASNKMFLGDNAHLSIINEYLSQSYQKFRRFFQTACDVPELQDKFNVDEYSDLVTLTKPVIYISIGEIINTHTLLLDHQDAIAPEHNDPIHELLDDLGEVPTIESLIGESSGNLNDPNKEALAKTEVSLTLTNKFDVPGDENAEMDARTILLNTKRLIVDVIRFQPGETLTEILETPATSEQEAEHQRAMQRRAIRDAKTPDKMKKSKSVKEDSNLTLQEKKEKIQTGLKKLTELGTVDPKNKYQELINDIARDIRNQRRYRQRRKAELVKLQQTYAALNSKATFYGEQVDYYKSYIKTCLDNLASKGKVSKKPREMKGKKSKKISLKYTAARLHEKGVLLEIEDLQVNQFKNVIFEISPTEEVGDFEVKAKFMGVQMETFMLHYQDLLQLQYEGVAVMKLFDRAKVNVNLLIFLLNKKFYGK\n"
     ]
    }
   ],
   "source": [
    "keys = data.keys()\n",
    "verbose = True\n",
    "# Primary Accession\n",
    "uni_id = 'NA'\n",
    "if 'primaryAccession' in keys:\n",
    "    uni_id = data['primaryAccession'] \n",
    "    if verbose:\n",
    "         print('uni_id',uni_id)\n",
    "else: \n",
    "    print('Issue: No primary Accession found')\n",
    "\n",
    "# Genes\n",
    "gene_name = 'NA'\n",
    "if 'genes' in keys:\n",
    "    gene_name = data['genes'][0]['geneName']['value']\n",
    "    if verbose:\n",
    "        print('gene_name',gene_name)\n",
    "else:\n",
    "    print('Issue: No gene found found')\n",
    "\n",
    "# getting the nucleotide sequenceID\n",
    "ref_seq_id = 'NA'\n",
    "if 'uniProtKBCrossReferences' in keys:\n",
    "    for db in data['uniProtKBCrossReferences']:\n",
    "        if db['database'] == \"RefSeq\":\n",
    "            if '.' in db['properties'][0]['value']:\n",
    "                ref_seq_id = db['properties'][0]['value'].split('.')[0]\n",
    "            else:\n",
    "                ref_seq_id = db['properties'][0]['value']\n",
    "if verbose:\n",
    "    print('ref_seq_id',ref_seq_id)\n",
    "\n",
    "# getting sequence info\n",
    "sequence = 'NA'\n",
    "seq_length = 0\n",
    "if 'sequence' in keys:\n",
    "    sequence = data['sequence']['value']\n",
    "    if verbose:\n",
    "        print('sequence',sequence)\n",
    "\n",
    "    seq_length = data['sequence']['length']\n",
    "    if seq_length != len(sequence):\n",
    "        print(f\"Warning for Uniprot id {uni_prot_id} sequence length recorded is not same as actual sequence\")\n",
    "\n",
    "known_things = get_domainn_region_info(data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "78411e67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Domain': [[44, 159, 'Calponin-homology (CH)'],\n",
       "  [679, 712, 'WW'],\n",
       "  [745, 774, 'IQ 1'],\n",
       "  [775, 804, 'IQ 2'],\n",
       "  [805, 834, 'IQ 3'],\n",
       "  [835, 864, 'IQ 4'],\n",
       "  [1004, 1237, 'Ras-GAP']],\n",
       " 'Region': [[956, 1274, 'C1'], [1276, 1657, 'C2'], [1410, 1448, 'Disordered']],\n",
       " 'Motif': []}"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "known_things"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "e7ba868e",
   "metadata": {},
   "outputs": [],
   "source": [
    "keys = known_things.keys()\n",
    "mask = np.zeros(seq_length)\n",
    "segment_list = []\n",
    "for k in keys:\n",
    "    for segement in known_things[k]:\n",
    "        if len(segement)>0:\n",
    "            segment_list.append([segement[0]-1,segement[1],segement[2]])\n",
    "            mask[segement[0]-1:segement[1]] = 1\n",
    "segment_list.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "7bd12faf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import groupby\n",
    "def false_indices(a):\n",
    "    i = 0\n",
    "    for k, g in groupby(a):\n",
    "        l = len(list(g))\n",
    "        if not k:\n",
    "            yield (i,i+l)\n",
    "        i += l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "2e959c0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we get the sequence section\n",
    "for i in range(len(segment_list)):\n",
    "    seq_chunk = sequence[segment_list[i][0]:segment_list[i][1]]\n",
    "    #print(segment_list[i][0],segment_list[i][1],seq_chunk)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "616fb919",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'N': 39,\n",
       " 'seq': 'TPATSEQEAEHQRAMQRRAIRDAKTPDKMKKSKSVKEDS',\n",
       " 'f-': 0.1794871794871795,\n",
       " 'f+': 0.2564102564102564,\n",
       " 'FCR': 0.4358974358974359,\n",
       " 'NCPR': 0.07692307692307693,\n",
       " 'kappa': 0.24082176888233633,\n",
       " 'Omega': 0.025890082904646273,\n",
       " 'sigma': 0.013574660633484155,\n",
       " 'delta': 0.09073125085333411,\n",
       " 'delta max': 0.3767568491603627,\n",
       " 'hydropathy': 2.7923076923076926,\n",
       " 'Phase plot region': 3,\n",
       " 'Hydrophobic count': 4,\n",
       " 'Positive count': 10,\n",
       " 'Negative count': 7}"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cider_data = get_all_data_from_seq(seq_chunk)\n",
    "df = pd.Series(cider_data).to_frame().T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "4250186c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At entry 0/100\n",
      "At entry 10/100\n",
      "At entry 20/100\n",
      "At entry 30/100\n",
      "At entry 40/100\n",
      "At entry 50/100\n",
      "At entry 60/100\n",
      "At entry 70/100\n",
      "At entry 80/100\n",
      "Issue: No gene found found\n",
      "At entry 90/100\n"
     ]
    }
   ],
   "source": [
    "vdf = pd.DataFrame(columns=['uniprot_id','gene_name','refseq_id'])\n",
    "new_data = []\n",
    "for i in range(100):\n",
    "    if i%10==0:\n",
    "        print(f'At entry {i}/{100}')\n",
    "    protein = uniprot_ids[i]\n",
    "    f = open(protein)\n",
    "    data = json.load(f)\n",
    "    f.close()\n",
    "    \n",
    "    # Getting some of the previous info\n",
    "    keys = data.keys()\n",
    "    verbose = False\n",
    "    # Primary Accession\n",
    "    uni_id = 'NA'\n",
    "    if 'primaryAccession' in keys:\n",
    "        uni_id = data['primaryAccession'] \n",
    "        if verbose:\n",
    "             print('uni_id',uni_id)\n",
    "    else: \n",
    "        print('Issue: No primary Accession found')\n",
    "\n",
    "    # Genes\n",
    "    gene_name = 'NA'\n",
    "    if 'genes' in keys:\n",
    "        gene_name = data['genes'][0]['geneName']['value']\n",
    "        if verbose:\n",
    "            print('gene_name',gene_name)\n",
    "    else:\n",
    "        print('Issue: No gene found found')\n",
    "\n",
    "    # getting the nucleotide sequenceID\n",
    "    ref_seq_id = 'NA'\n",
    "    if 'uniProtKBCrossReferences' in keys:\n",
    "        for db in data['uniProtKBCrossReferences']:\n",
    "            if db['database'] == \"RefSeq\":\n",
    "                if '.' in db['properties'][0]['value']:\n",
    "                    ref_seq_id = db['properties'][0]['value'].split('.')[0]\n",
    "                else:\n",
    "                    ref_seq_id = db['properties'][0]['value']\n",
    "    if verbose:\n",
    "        print('ref_seq_id',ref_seq_id)\n",
    "\n",
    "    # getting sequence info\n",
    "    sequence = 'NA'\n",
    "    seq_length = 0\n",
    "    if 'sequence' in keys:\n",
    "        sequence = data['sequence']['value']\n",
    "        if verbose:\n",
    "            print('sequence',sequence)\n",
    "\n",
    "        seq_length = data['sequence']['length']\n",
    "        if seq_length != len(sequence):\n",
    "            print(f\"Warning for Uniprot id {uni_prot_id} sequence length recorded is not same as actual sequence\")\n",
    "\n",
    "    known_things = get_domainn_region_info(data)\n",
    "    \n",
    "    # Sorting the segments\n",
    "    keys = known_things.keys()\n",
    "    mask = np.zeros(seq_length)\n",
    "    segment_list = []\n",
    "    for k in keys:\n",
    "        for segement in known_things[k]:\n",
    "            if len(segement)>0:\n",
    "                segment_list.append([segement[0]-1,segement[1],segement[2]])\n",
    "                mask[segement[0]-1:segement[1]] = 1\n",
    "    segment_list.sort()\n",
    "    \n",
    "    # finding unassigned segments\n",
    "    unassigned = list(false_indices(mask))\n",
    "    for u in unassigned:\n",
    "        segment_list.append([u[0], u[1], 'unassigned'])\n",
    "    segment_list.sort()\n",
    "    \n",
    "    ## Doing the cider analsyis\n",
    "    for i in range(len(segment_list)):\n",
    "        \n",
    "        seq_chunk = sequence[segment_list[i][0]:segment_list[i][1]]\n",
    "        cider_data = get_all_data_from_seq(seq_chunk)\n",
    "        columns = {'uniprot_id':uni_id,'gene_name':gene_name,'refseq_id':ref_seq_id, 'annotation':segment_list[i][2]}\n",
    "        concat_dic = {}\n",
    "        concat_dic.update(columns)\n",
    "        concat_dic.update(cider_data)\n",
    "        new_data.append(concat_dic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "673647e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_columns = pd.DataFrame(new_data)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "2d30a940",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>refseq_id</th>\n",
       "      <th>annotation</th>\n",
       "      <th>N</th>\n",
       "      <th>seq</th>\n",
       "      <th>f-</th>\n",
       "      <th>f+</th>\n",
       "      <th>FCR</th>\n",
       "      <th>NCPR</th>\n",
       "      <th>kappa</th>\n",
       "      <th>Omega</th>\n",
       "      <th>sigma</th>\n",
       "      <th>delta</th>\n",
       "      <th>delta max</th>\n",
       "      <th>hydropathy</th>\n",
       "      <th>Phase plot region</th>\n",
       "      <th>Hydrophobic count</th>\n",
       "      <th>Positive count</th>\n",
       "      <th>Negative count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Q96GA3</td>\n",
       "      <td>LTV1</td>\n",
       "      <td>NM_032860</td>\n",
       "      <td>unassigned</td>\n",
       "      <td>22</td>\n",
       "      <td>MPHRKKKPFIEKKKAVSFHLVH</td>\n",
       "      <td>0.045455</td>\n",
       "      <td>0.318182</td>\n",
       "      <td>0.363636</td>\n",
       "      <td>0.272727</td>\n",
       "      <td>0.425427</td>\n",
       "      <td>0.359966</td>\n",
       "      <td>0.204545</td>\n",
       "      <td>0.076269</td>\n",
       "      <td>0.179277</td>\n",
       "      <td>3.636364</td>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Q96GA3</td>\n",
       "      <td>LTV1</td>\n",
       "      <td>NM_032860</td>\n",
       "      <td>Disordered</td>\n",
       "      <td>32</td>\n",
       "      <td>RSQRDPLAADESAPQRVLLPTQKIDNEERRAE</td>\n",
       "      <td>0.218750</td>\n",
       "      <td>0.187500</td>\n",
       "      <td>0.406250</td>\n",
       "      <td>-0.031250</td>\n",
       "      <td>0.123860</td>\n",
       "      <td>0.036508</td>\n",
       "      <td>0.002404</td>\n",
       "      <td>0.042293</td>\n",
       "      <td>0.341457</td>\n",
       "      <td>3.103125</td>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>6</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Q96GA3</td>\n",
       "      <td>LTV1</td>\n",
       "      <td>NM_032860</td>\n",
       "      <td>unassigned</td>\n",
       "      <td>19</td>\n",
       "      <td>QRKYGVFFDDDYDYLQHLK</td>\n",
       "      <td>0.210526</td>\n",
       "      <td>0.157895</td>\n",
       "      <td>0.368421</td>\n",
       "      <td>-0.052632</td>\n",
       "      <td>0.665697</td>\n",
       "      <td>0.164810</td>\n",
       "      <td>0.007519</td>\n",
       "      <td>0.159463</td>\n",
       "      <td>0.239543</td>\n",
       "      <td>3.268421</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Q96GA3</td>\n",
       "      <td>LTV1</td>\n",
       "      <td>NM_032860</td>\n",
       "      <td>Disordered</td>\n",
       "      <td>26</td>\n",
       "      <td>EPSGPSELIPSSTFSAHNRREEKEET</td>\n",
       "      <td>0.230769</td>\n",
       "      <td>0.115385</td>\n",
       "      <td>0.346154</td>\n",
       "      <td>-0.115385</td>\n",
       "      <td>0.076693</td>\n",
       "      <td>0.403497</td>\n",
       "      <td>0.038462</td>\n",
       "      <td>0.017630</td>\n",
       "      <td>0.229882</td>\n",
       "      <td>3.026923</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Q96GA3</td>\n",
       "      <td>LTV1</td>\n",
       "      <td>NM_032860</td>\n",
       "      <td>unassigned</td>\n",
       "      <td>69</td>\n",
       "      <td>LVIPSTGIKLPSSVFASEFEEDVGLLNKAAPVSGPRLDFDPDIVAA...</td>\n",
       "      <td>0.246377</td>\n",
       "      <td>0.043478</td>\n",
       "      <td>0.289855</td>\n",
       "      <td>-0.202899</td>\n",
       "      <td>0.387080</td>\n",
       "      <td>0.112137</td>\n",
       "      <td>0.142029</td>\n",
       "      <td>0.073485</td>\n",
       "      <td>0.189845</td>\n",
       "      <td>4.492754</td>\n",
       "      <td>2</td>\n",
       "      <td>18</td>\n",
       "      <td>3</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>495</th>\n",
       "      <td>P60608</td>\n",
       "      <td>ERVFC1-1</td>\n",
       "      <td>NA</td>\n",
       "      <td>CX6CC</td>\n",
       "      <td>9</td>\n",
       "      <td>CLFLQEECC</td>\n",
       "      <td>0.222222</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.222222</td>\n",
       "      <td>-0.222222</td>\n",
       "      <td>0.603249</td>\n",
       "      <td>0.312193</td>\n",
       "      <td>0.222222</td>\n",
       "      <td>0.019485</td>\n",
       "      <td>0.032299</td>\n",
       "      <td>5.322222</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>496</th>\n",
       "      <td>Q9C0B5</td>\n",
       "      <td>ZDHHC5</td>\n",
       "      <td>XM_011544901</td>\n",
       "      <td>unassigned</td>\n",
       "      <td>103</td>\n",
       "      <td>MPAESGKRFKPSKYVPVSAAAIFLVGATTLFFAFTCPGLSLYVSPA...</td>\n",
       "      <td>0.106796</td>\n",
       "      <td>0.097087</td>\n",
       "      <td>0.203883</td>\n",
       "      <td>-0.009709</td>\n",
       "      <td>0.232912</td>\n",
       "      <td>0.431847</td>\n",
       "      <td>0.000462</td>\n",
       "      <td>0.041226</td>\n",
       "      <td>0.177004</td>\n",
       "      <td>4.803883</td>\n",
       "      <td>1</td>\n",
       "      <td>26</td>\n",
       "      <td>10</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>497</th>\n",
       "      <td>Q9C0B5</td>\n",
       "      <td>ZDHHC5</td>\n",
       "      <td>XM_011544901</td>\n",
       "      <td>DHHC</td>\n",
       "      <td>51</td>\n",
       "      <td>KWCATCRFYRPPRCSHCSVCDNCVEEFDHHCPWVNNCIGRRNYRYF...</td>\n",
       "      <td>0.078431</td>\n",
       "      <td>0.137255</td>\n",
       "      <td>0.215686</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.455520</td>\n",
       "      <td>0.233759</td>\n",
       "      <td>0.016043</td>\n",
       "      <td>0.069942</td>\n",
       "      <td>0.153544</td>\n",
       "      <td>4.158824</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>7</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>498</th>\n",
       "      <td>Q9C0B5</td>\n",
       "      <td>ZDHHC5</td>\n",
       "      <td>XM_011544901</td>\n",
       "      <td>unassigned</td>\n",
       "      <td>134</td>\n",
       "      <td>SLTAHIMGVFGFGLLYVLYHIEELSGVRTAVTMAVMCVAGLFFIPV...</td>\n",
       "      <td>0.052239</td>\n",
       "      <td>0.104478</td>\n",
       "      <td>0.156716</td>\n",
       "      <td>0.052239</td>\n",
       "      <td>0.215688</td>\n",
       "      <td>0.408543</td>\n",
       "      <td>0.017413</td>\n",
       "      <td>0.027919</td>\n",
       "      <td>0.129439</td>\n",
       "      <td>4.811194</td>\n",
       "      <td>1</td>\n",
       "      <td>40</td>\n",
       "      <td>14</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>499</th>\n",
       "      <td>Q9C0B5</td>\n",
       "      <td>ZDHHC5</td>\n",
       "      <td>XM_011544901</td>\n",
       "      <td>Disordered</td>\n",
       "      <td>427</td>\n",
       "      <td>GELRRTKSKGSLEITESQSADAEPPPPPKPDLSRYTGLRTHLGLAT...</td>\n",
       "      <td>0.103044</td>\n",
       "      <td>0.119438</td>\n",
       "      <td>0.222482</td>\n",
       "      <td>0.016393</td>\n",
       "      <td>0.176967</td>\n",
       "      <td>0.146510</td>\n",
       "      <td>0.001208</td>\n",
       "      <td>0.038211</td>\n",
       "      <td>0.215920</td>\n",
       "      <td>3.617330</td>\n",
       "      <td>1</td>\n",
       "      <td>61</td>\n",
       "      <td>51</td>\n",
       "      <td>44</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>500 rows × 20 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    uniprot_id gene_name     refseq_id  annotation    N  \\\n",
       "0       Q96GA3      LTV1     NM_032860  unassigned   22   \n",
       "1       Q96GA3      LTV1     NM_032860  Disordered   32   \n",
       "2       Q96GA3      LTV1     NM_032860  unassigned   19   \n",
       "3       Q96GA3      LTV1     NM_032860  Disordered   26   \n",
       "4       Q96GA3      LTV1     NM_032860  unassigned   69   \n",
       "..         ...       ...           ...         ...  ...   \n",
       "495     P60608  ERVFC1-1            NA       CX6CC    9   \n",
       "496     Q9C0B5    ZDHHC5  XM_011544901  unassigned  103   \n",
       "497     Q9C0B5    ZDHHC5  XM_011544901        DHHC   51   \n",
       "498     Q9C0B5    ZDHHC5  XM_011544901  unassigned  134   \n",
       "499     Q9C0B5    ZDHHC5  XM_011544901  Disordered  427   \n",
       "\n",
       "                                                   seq        f-        f+  \\\n",
       "0                               MPHRKKKPFIEKKKAVSFHLVH  0.045455  0.318182   \n",
       "1                     RSQRDPLAADESAPQRVLLPTQKIDNEERRAE  0.218750  0.187500   \n",
       "2                                  QRKYGVFFDDDYDYLQHLK  0.210526  0.157895   \n",
       "3                           EPSGPSELIPSSTFSAHNRREEKEET  0.230769  0.115385   \n",
       "4    LVIPSTGIKLPSSVFASEFEEDVGLLNKAAPVSGPRLDFDPDIVAA...  0.246377  0.043478   \n",
       "..                                                 ...       ...       ...   \n",
       "495                                          CLFLQEECC  0.222222  0.000000   \n",
       "496  MPAESGKRFKPSKYVPVSAAAIFLVGATTLFFAFTCPGLSLYVSPA...  0.106796  0.097087   \n",
       "497  KWCATCRFYRPPRCSHCSVCDNCVEEFDHHCPWVNNCIGRRNYRYF...  0.078431  0.137255   \n",
       "498  SLTAHIMGVFGFGLLYVLYHIEELSGVRTAVTMAVMCVAGLFFIPV...  0.052239  0.104478   \n",
       "499  GELRRTKSKGSLEITESQSADAEPPPPPKPDLSRYTGLRTHLGLAT...  0.103044  0.119438   \n",
       "\n",
       "          FCR      NCPR     kappa     Omega     sigma     delta  delta max  \\\n",
       "0    0.363636  0.272727  0.425427  0.359966  0.204545  0.076269   0.179277   \n",
       "1    0.406250 -0.031250  0.123860  0.036508  0.002404  0.042293   0.341457   \n",
       "2    0.368421 -0.052632  0.665697  0.164810  0.007519  0.159463   0.239543   \n",
       "3    0.346154 -0.115385  0.076693  0.403497  0.038462  0.017630   0.229882   \n",
       "4    0.289855 -0.202899  0.387080  0.112137  0.142029  0.073485   0.189845   \n",
       "..        ...       ...       ...       ...       ...       ...        ...   \n",
       "495  0.222222 -0.222222  0.603249  0.312193  0.222222  0.019485   0.032299   \n",
       "496  0.203883 -0.009709  0.232912  0.431847  0.000462  0.041226   0.177004   \n",
       "497  0.215686  0.058824  0.455520  0.233759  0.016043  0.069942   0.153544   \n",
       "498  0.156716  0.052239  0.215688  0.408543  0.017413  0.027919   0.129439   \n",
       "499  0.222482  0.016393  0.176967  0.146510  0.001208  0.038211   0.215920   \n",
       "\n",
       "     hydropathy  Phase plot region  Hydrophobic count  Positive count  \\\n",
       "0      3.636364                  3                  5               7   \n",
       "1      3.103125                  3                  5               6   \n",
       "2      3.268421                  3                  3               3   \n",
       "3      3.026923                  2                  2               3   \n",
       "4      4.492754                  2                 18               3   \n",
       "..          ...                ...                ...             ...   \n",
       "495    5.322222                  1                  2               0   \n",
       "496    4.803883                  1                 26              10   \n",
       "497    4.158824                  1                  7               7   \n",
       "498    4.811194                  1                 40              14   \n",
       "499    3.617330                  1                 61              51   \n",
       "\n",
       "     Negative count  \n",
       "0                 1  \n",
       "1                 7  \n",
       "2                 4  \n",
       "3                 6  \n",
       "4                17  \n",
       "..              ...  \n",
       "495               2  \n",
       "496              11  \n",
       "497               4  \n",
       "498               7  \n",
       "499              44  \n",
       "\n",
       "[500 rows x 20 columns]"
      ]
     },
     "execution_count": 166,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "bfd10dc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_columns.to_csv('cider/example_output.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ec87914",
   "metadata": {},
   "source": [
    "## Old stuff may not be needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7818e19",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7eceb9e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_alpha_helix_length_and_location(secondary_struc, min_length=8):\n",
    "    helix_regions = []\n",
    "    counter = 0\n",
    "    curr_helix = []\n",
    "    indexes = [i for i, x in enumerate(list(secondary_struc)) if x == 'H']\n",
    "    for i in range(len(indexes)-1):\n",
    "        difference = indexes[i+1]-indexes[i]\n",
    "        if difference == 1:\n",
    "            curr_helix.append(indexes[i])\n",
    "            if i == len(indexes)-2:\n",
    "                if len(curr_helix)>=min_length-1:\n",
    "                    curr_helix.append(indexes[i+1])\n",
    "                    helix_regions.append(curr_helix)\n",
    "        else:\n",
    "            curr_helix.append(indexes[i])\n",
    "            if len(curr_helix)>=min_length-1:\n",
    "                helix_regions.append(curr_helix)\n",
    "            curr_helix = []\n",
    "    return helix_regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a3ee27a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_serines(sequence):\n",
    "    n_serine = sequence.count('S')\n",
    "    if n_serine >=2:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "378e385a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_helices(helix_index_list, known_things, sequence_length, discard_list= ['bHLH', 'Leucine-zipper'],overlap_threshold =4):\n",
    "    '''\n",
    "    Parameters:\n",
    "    -----------\n",
    "    helix_index_list : list\n",
    "        2D list containing arrays of indexes where alpha helixes are\n",
    "    known_things : dictionary\n",
    "        extracted information from uniprot IDs around regions, domains and motives\n",
    "    discard_list : list of Strings\n",
    "        list that contains the identifiers to discard\n",
    "    \n",
    "    Returns:\n",
    "    --------\n",
    "    helix_index_list : list\n",
    "        filtered list with correct helix indexes without the overlap\n",
    "    Algorithm description:\n",
    "    \n",
    "    - generate a boolean array of all False of length of the sequence\n",
    "    - for each annotation we want to check:\n",
    "    - Add True between start and end of domains/regions we want to check in the boolean array\n",
    "    - Then loop over helix list creating a boolean array of length sequence for each Helix section\n",
    "    - Use logic and to compare regions boolean array with helix boolean array. \n",
    "    - If the number of over lap of Trues is larger than theshold x, remove this helix chunch from list and don't write to spreadsheet.\n",
    "    '''\n",
    "    \n",
    "    bool_array =  np.zeros(sequence_length, dtype=bool)\n",
    "    # alpha_helix_index_list = get_alpha_helix_length_and_location(sec_struc_output, min_length=7)\n",
    "    helix_list_remove_index = []\n",
    "    for key in known_things.keys():\n",
    "        for entry in known_things[key]:\n",
    "            bool_array =  np.zeros(sequence_length, dtype=bool)\n",
    "            if entry[-1] in discard_list:\n",
    "                bool_array[entry[0]-1:entry[1]] = True\n",
    "            index = 0\n",
    "            for a in helix_index_list:\n",
    "                curr_helix =  np.zeros(sequence_length, dtype=bool)\n",
    "                curr_helix[a] = True\n",
    "                overalp = (np.sum(np.logical_and(bool_array, curr_helix)))\n",
    "                if overalp>overlap_threshold:\n",
    "                    helix_list_remove_index.append(index)\n",
    "                index = index+1\n",
    "    unique_idx = np.unique(helix_list_remove_index)\n",
    "    for index in sorted(unique_idx, reverse=True):\n",
    "        del helix_index_list[index]\n",
    "    return helix_index_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "38e33270",
   "metadata": {},
   "outputs": [],
   "source": [
    "def annotating_helices(helix_index_list, known_things, sequence_length,overlap_threshold =4):\n",
    "    '''\n",
    "    Parameters:\n",
    "    -----------\n",
    "    helix_index_list : list\n",
    "        2D list containing arrays of indexes where alpha helixes are\n",
    "    known_things : dictionary\n",
    "        extracted information from uniprot IDs around regions, domains and motives\n",
    "    discard_list : list of Strings\n",
    "        list that contains the identifiers to discard\n",
    "    \n",
    "    Returns:\n",
    "    --------\n",
    "    helix_index_list : list\n",
    "        filtered list with correct helix indexes without the overlap\n",
    "    Algorithm description:\n",
    "    \n",
    "    - generate a boolean array of all False of length of the sequence\n",
    "    - for each annotation we want to check:\n",
    "    - Add True between start and end of domains/regions we want to check in the boolean array\n",
    "    - Then loop over helix list creating a boolean array of length sequence for each Helix section\n",
    "    - Use logic and to compare regions boolean array with helix boolean array. \n",
    "    - If the number of over lap of Trues is larger than theshold x, remove this helix chunch from list and don't write to spreadsheet.\n",
    "    '''\n",
    "    annotations = [\"\"]*len(helix_index_list)\n",
    "    for key in known_things.keys():\n",
    "        for entry in known_things[key]:\n",
    "            # print(f'checking entry {entry}')\n",
    "            bool_array =  np.zeros(sequence_length, dtype=bool)\n",
    "            if entry[0] is None:\n",
    "                continue\n",
    "            elif entry[1] is None:\n",
    "                continue\n",
    "            bool_array[entry[0]-1:entry[1]] = True\n",
    "            index = 0\n",
    "            \n",
    "            for a in helix_index_list:\n",
    "                curr_helix =  np.zeros(sequence_length, dtype=bool)\n",
    "                curr_helix[a] = True\n",
    "                overalp = (np.sum(np.logical_and(bool_array, curr_helix)))\n",
    "                if overalp>overlap_threshold:\n",
    "                    # print(f'annotation found {entry[-1]}')\n",
    "                    annotations[index] = annotations[index] +key+':'+entry[-1]+\";\"\n",
    "                index = index+1\n",
    "    return annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "75f9f3d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_discrad_list(domain_f_name, region_f_name):\n",
    "    # Getting domains\n",
    "    domain_filter = []\n",
    "    # This should go into a try thing!\n",
    "    f = open(domain_f_name)\n",
    "    f_content = f.readlines()\n",
    "    f.close()\n",
    "\n",
    "    for l in f_content:\n",
    "        domain_filter.append(l.strip()) \n",
    "        \n",
    "    # Getting regions\n",
    "    region_filter = []\n",
    "    f = open(region_f_name)\n",
    "    f_content = f.readlines()\n",
    "    f.close()\n",
    "\n",
    "    for l in f_content:\n",
    "        region_filter.append(l.strip())\n",
    "    discard_list = np.hstack((region_filter,domain_filter))\n",
    "    return discard_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d63068b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_table_rows_for_uniprot_id(uni_prot_id, df, work_dir = 'temp', verbose = False, annotate = False):\n",
    "    # print('Working on uniprot ID', uni_prot_id)\n",
    "    sec_struc_output = np.load('dssp_info/'+uni_prot_id+'.npy')\n",
    "    # print('got secondary structure')\n",
    "    #print(sec_struc_output)\n",
    "    \n",
    "    # Process uniprot file\n",
    "    f = open('unique_ids/'+uni_prot_id+'.json')\n",
    "    data = json.load(f)\n",
    "    f.close()\n",
    "    \n",
    "    # print('read json file successfully')\n",
    "    \n",
    "    # Some basics:\n",
    "    # checking keys\n",
    "    keys = data.keys()\n",
    "    \n",
    "    # Primary Accession\n",
    "    uni_id = 'NA'\n",
    "    if 'primaryAccession' in keys:\n",
    "        uni_id = data['primaryAccession'] \n",
    "        if verbose:\n",
    "             print('uni_id',uni_id)\n",
    "    else: \n",
    "        print('Issue: No primary Accession found')\n",
    "    \n",
    "    # Genes\n",
    "    gene_name = 'NA'\n",
    "    if 'genes' in keys:\n",
    "        gene_name = data['genes'][0]['geneName']['value']\n",
    "        if verbose:\n",
    "            print('gene_name',gene_name)\n",
    "    else:\n",
    "        print('Issue: No gene found found')\n",
    "    \n",
    "    # getting the nucleotide sequenceID\n",
    "    ref_seq_id = 'NA'\n",
    "    if 'uniProtKBCrossReferences' in keys:\n",
    "        for db in data['uniProtKBCrossReferences']:\n",
    "            if db['database'] == \"RefSeq\":\n",
    "                if '.' in db['properties'][0]['value']:\n",
    "                    ref_seq_id = db['properties'][0]['value'].split('.')[0]\n",
    "                else:\n",
    "                    ref_seq_id = db['properties'][0]['value']\n",
    "    if verbose:\n",
    "        print('ref_seq_id',ref_seq_id)\n",
    "    \n",
    "    # getting sequence info\n",
    "    sequence = 'NA'\n",
    "    seq_length = 0\n",
    "    if 'sequence' in keys:\n",
    "        sequence = data['sequence']['value']\n",
    "        if verbose:\n",
    "            print('sequence',sequence)\n",
    "    \n",
    "        seq_length = data['sequence']['length']\n",
    "        if seq_length != len(sequence):\n",
    "            print(f\"Warning for Uniprot id {uni_prot_id} sequence length recorded is not same as actual sequence\")\n",
    "\n",
    "    if seq_length == len(sec_struc_output):\n",
    "\n",
    "        \n",
    "\n",
    "        # Now lets get the helices:\n",
    "        alpha_helix_index_list = get_alpha_helix_length_and_location(sec_struc_output, min_length=8)\n",
    "\n",
    "        # Get the information of what domains and regions exist in the uniport id file:\n",
    "        known_things = get_domainn_region_info(data)\n",
    "\n",
    "        discard_list = ['PUM-HD', 'HEAT', 'ARM', 'bHLH', 'Leucine-zipper']\n",
    "        discard_list = get_discrad_list('data/domain.csv', 'data/region_filter.csv')\n",
    "        discard_list = []\n",
    "\n",
    "        # Missing: Filter the alpha_helix_index_list\n",
    "        filtered_helix_list=filter_helices(alpha_helix_index_list, known_things, seq_length, discard_list=discard_list,overlap_threshold =4)\n",
    "\n",
    "        # annotate helices:\n",
    "        if annotate:\n",
    "            filtered_helix_list = alpha_helix_index_list\n",
    "            annotations = annotating_helices(filtered_helix_list, known_things, seq_length,overlap_threshold =4)\n",
    "\n",
    "        for i in range(len(filtered_helix_list)):\n",
    "            firstAA_position_in_HELIDR = filtered_helix_list[i][0]+1\n",
    "            lastAA_position_in_HELIDR = filtered_helix_list[i][-1]+1\n",
    "            HELIDR_seq = sequence[filtered_helix_list[i][0]:filtered_helix_list[i][-1]+1]\n",
    "            down_stream_seq = ''\n",
    "            up_stream_seq = ''\n",
    "            if filtered_helix_list[i][0]-10 >= 0 and filtered_helix_list[i][-1]+11 < len(sequence):\n",
    "                # Note the +2 here does not include the last helix AA. \n",
    "                # We have to make sure here that we check the arrays are not out of bounds!\n",
    "                down_stream_seq = sequence[filtered_helix_list[i][-1]+1:filtered_helix_list[i][-1]+11]\n",
    "                up_stream_seq = sequence[filtered_helix_list[i][0]-10:filtered_helix_list[i][0]]\n",
    "            elif alpha_helix_index_list[i][0]-10 <= 0:\n",
    "                # Do we want a shorter version? \n",
    "                down_stream_seq = sequence[filtered_helix_list[i][-1]+1:filtered_helix_list[i][-1]+11]\n",
    "                up_stream_seq = sequence[0:filtered_helix_list[i][0]]\n",
    "            elif alpha_helix_index_list[i][-1]+11 > len(sequence):\n",
    "                down_stream_seq = sequence[filtered_helix_list[i][-1]+2:len(sequence)+1]\n",
    "                up_stream_seq = sequence[filtered_helix_list[i][0]-10:filtered_helix_list[i][0]]\n",
    "\n",
    "            # This may need fixing if we have shorter upstream and downstream strings\n",
    "            if len(down_stream_seq)==10 and len(up_stream_seq)==10: \n",
    "                Two_S5P_down = count_serines(down_stream_seq[:5])\n",
    "                Two_S5P_up = count_serines(up_stream_seq[5:])\n",
    "                # NOW we assemble the row:\n",
    "                if annotate:\n",
    "                    new_row = [uni_id,gene_name,ref_seq_id,firstAA_position_in_HELIDR,lastAA_position_in_HELIDR,up_stream_seq,HELIDR_seq,\n",
    "                          down_stream_seq,Two_S5P_up,Two_S5P_down,'','','','','','','','','',annotations[i]]\n",
    "                else:\n",
    "                    new_row = [uni_id,gene_name,ref_seq_id,firstAA_position_in_HELIDR,lastAA_position_in_HELIDR,up_stream_seq,HELIDR_seq,\n",
    "                          down_stream_seq,Two_S5P_up,Two_S5P_down,'','','','','','','','','']\n",
    "                df.loc[len(df)] = new_row\n",
    "            else:\n",
    "                # We may want to revisit this continue here\n",
    "                continue\n",
    "        return df, None\n",
    "            \n",
    "    else:\n",
    "        # print(f'sequence length from structure is {seq_length}, from alpha fold {len(sec_struc_output)}')\n",
    "        # print(\"there is an incompatibility between the sequence and alpha fold structure\")\n",
    "        # print(f'the current uniprotid with the issue is: {uni_prot_id}')\n",
    "        return df, uni_prot_id"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16663482",
   "metadata": {},
   "source": [
    "## Testing things out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "68f0044c",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data/unique_ids_from_spreadsheet.txt', 'r')\n",
    "f_content = f.readlines()\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "74ea2032",
   "metadata": {},
   "outputs": [],
   "source": [
    "ids = []\n",
    "for f in f_content:\n",
    "    ids.append(f.strip())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ceb56b2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fail_list = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "166454f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def reinitalise_df(annotations = False):\n",
    "    if annotations:\n",
    "        df = pd.DataFrame(columns=['uniprot_id','gene_name','refseq_id','firstAA_position_in_HELIDR','lastAA_position_in_HELIDR','HELIDR_upstream_seq'\n",
    "                       ,'HELIDR_seq', 'HELIDR_downstream_seq', '2S5P_up', '2S5P_down', '2S5P1_up','2S5P1_down','2S5P1_helix','HEK293T_expressed','NonTMD[3]_TMD[2]_SEC[1]',\n",
    "                       'Non_TMD_classification','4 compartments','TG_CY','TG_SR_nonS', 'annotation'])\n",
    "\n",
    "    else:\n",
    "        df = pd.DataFrame(columns=['uniprot_id','gene_name','refseq_id','firstAA_position_in_HELIDR','lastAA_position_in_HELIDR','HELIDR_upstream_seq'\n",
    "                       ,'HELIDR_seq', 'HELIDR_downstream_seq', '2S5P_up', '2S5P_down', '2S5P1_up','2S5P1_down','2S5P1_helix','HEK293T_expressed','NonTMD[3]_TMD[2]_SEC[1]',\n",
    "                       'Non_TMD_classification','4 compartments','TG_CY','TG_SR_nonS'])\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "39699650",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on ID: O00305 this is entry 0/16236\n"
     ]
    }
   ],
   "source": [
    "counter = 0\n",
    "df = reinitalise_df(annotations= True)\n",
    "for i in ids[0:5]:\n",
    "    if counter%200 == 0:\n",
    "        print(f'Working on ID: {i} this is entry {counter}/{len(ids)}')\n",
    "    df, fail = generate_table_rows_for_uniprot_id(i,df, annotate=True)\n",
    "    if fail != None:\n",
    "        fail_list.append(fail)\n",
    "    counter = counter+1\n",
    "df.to_csv('test.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "be663abd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>refseq_id</th>\n",
       "      <th>firstAA_position_in_HELIDR</th>\n",
       "      <th>lastAA_position_in_HELIDR</th>\n",
       "      <th>HELIDR_upstream_seq</th>\n",
       "      <th>HELIDR_seq</th>\n",
       "      <th>HELIDR_downstream_seq</th>\n",
       "      <th>2S5P_up</th>\n",
       "      <th>2S5P_down</th>\n",
       "      <th>2S5P1_up</th>\n",
       "      <th>2S5P1_down</th>\n",
       "      <th>2S5P1_helix</th>\n",
       "      <th>HEK293T_expressed</th>\n",
       "      <th>NonTMD[3]_TMD[2]_SEC[1]</th>\n",
       "      <th>Non_TMD_classification</th>\n",
       "      <th>4 compartments</th>\n",
       "      <th>TG_CY</th>\n",
       "      <th>TG_SR_nonS</th>\n",
       "      <th>annotation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>O00305</td>\n",
       "      <td>CACNB4</td>\n",
       "      <td>NM_001320722</td>\n",
       "      <td>61</td>\n",
       "      <td>89</td>\n",
       "      <td>SADSYTSRPS</td>\n",
       "      <td>DSDVSLEEDREAIRQEREQQAAIQLERAK</td>\n",
       "      <td>SKPVAFAVKT</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>Region:Disordered;</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O00305</td>\n",
       "      <td>CACNB4</td>\n",
       "      <td>NM_001320722</td>\n",
       "      <td>154</td>\n",
       "      <td>169</td>\n",
       "      <td>EGCEIGFIPS</td>\n",
       "      <td>PLRLENIRIQQEQKRG</td>\n",
       "      <td>RFHGGKSSGN</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>Domain:SH3;</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>O00305</td>\n",
       "      <td>CACNB4</td>\n",
       "      <td>NM_001320722</td>\n",
       "      <td>233</td>\n",
       "      <td>249</td>\n",
       "      <td>VLVGPSLKGY</td>\n",
       "      <td>EVTDMMQKALFDFLKHR</td>\n",
       "      <td>FDGRISITRV</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>O00305</td>\n",
       "      <td>CACNB4</td>\n",
       "      <td>NM_001320722</td>\n",
       "      <td>274</td>\n",
       "      <td>280</td>\n",
       "      <td>SLAKRSVLNN</td>\n",
       "      <td>PSKRAII</td>\n",
       "      <td>ERSNTRSSLA</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>O00305</td>\n",
       "      <td>CACNB4</td>\n",
       "      <td>NM_001320722</td>\n",
       "      <td>289</td>\n",
       "      <td>302</td>\n",
       "      <td>IIERSNTRSS</td>\n",
       "      <td>LAEVQSEIERIFEL</td>\n",
       "      <td>ARSLQLVVLD</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59</th>\n",
       "      <td>Q8WUM4</td>\n",
       "      <td>PDCD6IP</td>\n",
       "      <td>NM_013374</td>\n",
       "      <td>596</td>\n",
       "      <td>639</td>\n",
       "      <td>TALAQDGVIN</td>\n",
       "      <td>EEALSVTELDRVYGGLTTKVQESLKKQEGLLKNIQVSHQEFSKM</td>\n",
       "      <td>KQSNNEANLR</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>Region:Interaction with EIAV p9;Region:Interac...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>60</th>\n",
       "      <td>Q8WUM4</td>\n",
       "      <td>PDCD6IP</td>\n",
       "      <td>NM_013374</td>\n",
       "      <td>644</td>\n",
       "      <td>713</td>\n",
       "      <td>QEFSKMKQSN</td>\n",
       "      <td>NEANLREEVLKNLATAYDNFVELVANLKEGTKFYNELTEILVRFQN...</td>\n",
       "      <td>AREPSAPSIP</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>Region:Interaction with EIAV p9;Region:Interac...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>61</th>\n",
       "      <td>Q9H4S2</td>\n",
       "      <td>GSX1</td>\n",
       "      <td>NM_145657</td>\n",
       "      <td>156</td>\n",
       "      <td>168</td>\n",
       "      <td>SSKRMRTAFT</td>\n",
       "      <td>STQLLELEREFAS</td>\n",
       "      <td>NMYLSRLRRI</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62</th>\n",
       "      <td>Q9H4S2</td>\n",
       "      <td>GSX1</td>\n",
       "      <td>NM_145657</td>\n",
       "      <td>174</td>\n",
       "      <td>183</td>\n",
       "      <td>REFASNMYLS</td>\n",
       "      <td>RLRRIEIATY</td>\n",
       "      <td>LNLSEKQVKI</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>63</th>\n",
       "      <td>Q9H4S2</td>\n",
       "      <td>GSX1</td>\n",
       "      <td>NM_145657</td>\n",
       "      <td>188</td>\n",
       "      <td>205</td>\n",
       "      <td>IEIATYLNLS</td>\n",
       "      <td>EKQVKIWFQNRRVKHKKE</td>\n",
       "      <td>GKGSNHRGGG</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>Region:Disordered;</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>64 rows × 20 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   uniprot_id gene_name     refseq_id  firstAA_position_in_HELIDR  \\\n",
       "0      O00305    CACNB4  NM_001320722                          61   \n",
       "1      O00305    CACNB4  NM_001320722                         154   \n",
       "2      O00305    CACNB4  NM_001320722                         233   \n",
       "3      O00305    CACNB4  NM_001320722                         274   \n",
       "4      O00305    CACNB4  NM_001320722                         289   \n",
       "..        ...       ...           ...                         ...   \n",
       "59     Q8WUM4   PDCD6IP     NM_013374                         596   \n",
       "60     Q8WUM4   PDCD6IP     NM_013374                         644   \n",
       "61     Q9H4S2      GSX1     NM_145657                         156   \n",
       "62     Q9H4S2      GSX1     NM_145657                         174   \n",
       "63     Q9H4S2      GSX1     NM_145657                         188   \n",
       "\n",
       "    lastAA_position_in_HELIDR HELIDR_upstream_seq  \\\n",
       "0                          89          SADSYTSRPS   \n",
       "1                         169          EGCEIGFIPS   \n",
       "2                         249          VLVGPSLKGY   \n",
       "3                         280          SLAKRSVLNN   \n",
       "4                         302          IIERSNTRSS   \n",
       "..                        ...                 ...   \n",
       "59                        639          TALAQDGVIN   \n",
       "60                        713          QEFSKMKQSN   \n",
       "61                        168          SSKRMRTAFT   \n",
       "62                        183          REFASNMYLS   \n",
       "63                        205          IEIATYLNLS   \n",
       "\n",
       "                                           HELIDR_seq HELIDR_downstream_seq  \\\n",
       "0                       DSDVSLEEDREAIRQEREQQAAIQLERAK            SKPVAFAVKT   \n",
       "1                                    PLRLENIRIQQEQKRG            RFHGGKSSGN   \n",
       "2                                   EVTDMMQKALFDFLKHR            FDGRISITRV   \n",
       "3                                             PSKRAII            ERSNTRSSLA   \n",
       "4                                      LAEVQSEIERIFEL            ARSLQLVVLD   \n",
       "..                                                ...                   ...   \n",
       "59       EEALSVTELDRVYGGLTTKVQESLKKQEGLLKNIQVSHQEFSKM            KQSNNEANLR   \n",
       "60  NEANLREEVLKNLATAYDNFVELVANLKEGTKFYNELTEILVRFQN...            AREPSAPSIP   \n",
       "61                                      STQLLELEREFAS            NMYLSRLRRI   \n",
       "62                                         RLRRIEIATY            LNLSEKQVKI   \n",
       "63                                 EKQVKIWFQNRRVKHKKE            GKGSNHRGGG   \n",
       "\n",
       "    2S5P_up  2S5P_down 2S5P1_up 2S5P1_down 2S5P1_helix HEK293T_expressed  \\\n",
       "0         1          0                                                     \n",
       "1         0          0                                                     \n",
       "2         0          0                                                     \n",
       "3         0          0                                                     \n",
       "4         1          0                                                     \n",
       "..      ...        ...      ...        ...         ...               ...   \n",
       "59        0          0                                                     \n",
       "60        0          0                                                     \n",
       "61        0          0                                                     \n",
       "62        0          0                                                     \n",
       "63        0          0                                                     \n",
       "\n",
       "   NonTMD[3]_TMD[2]_SEC[1] Non_TMD_classification 4 compartments TG_CY  \\\n",
       "0                                                                        \n",
       "1                                                                        \n",
       "2                                                                        \n",
       "3                                                                        \n",
       "4                                                                        \n",
       "..                     ...                    ...            ...   ...   \n",
       "59                                                                       \n",
       "60                                                                       \n",
       "61                                                                       \n",
       "62                                                                       \n",
       "63                                                                       \n",
       "\n",
       "   TG_SR_nonS                                         annotation  \n",
       "0                                             Region:Disordered;  \n",
       "1                                                    Domain:SH3;  \n",
       "2                                                                 \n",
       "3                                                                 \n",
       "4                                                                 \n",
       "..        ...                                                ...  \n",
       "59             Region:Interaction with EIAV p9;Region:Interac...  \n",
       "60             Region:Interaction with EIAV p9;Region:Interac...  \n",
       "61                                                                \n",
       "62                                                                \n",
       "63                                            Region:Disordered;  \n",
       "\n",
       "[64 rows x 20 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "id": "559e09f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on ID: Q8IX06 this is entry 5000/16236\n",
      "Working on ID: Q6PL18 this is entry 5200/16236\n",
      "Working on ID: Q8TD19 this is entry 5400/16236\n",
      "Working on ID: Q9BZ71 this is entry 5600/16236\n",
      "Working on ID: Q8IZ08 this is entry 5800/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: P35249 this is entry 6000/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: A8MQ27 this is entry 6200/16236\n",
      "Issue: No gene found found\n",
      "Issue: No gene found found\n",
      "Working on ID: Q9NX36 this is entry 6400/16236\n",
      "Working on ID: O75603 this is entry 6600/16236\n",
      "Working on ID: Q495M3 this is entry 6800/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: Q96B01 this is entry 7000/16236\n",
      "Working on ID: A0A0B4J2F2 this is entry 7200/16236\n",
      "Working on ID: P49747 this is entry 7400/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: Q5JTB6 this is entry 7600/16236\n",
      "Working on ID: Q9UK96 this is entry 7800/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: Q96B26 this is entry 8000/16236\n",
      "Working on ID: Q8N782 this is entry 8200/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: P15018 this is entry 8400/16236\n",
      "Working on ID: Q9BVS4 this is entry 8600/16236\n",
      "Working on ID: A6NEY3 this is entry 8800/16236\n",
      "Working on ID: Q9UKV3 this is entry 9000/16236\n",
      "Working on ID: Q99490 this is entry 9200/16236\n",
      "Working on ID: Q8TAX9 this is entry 9400/16236\n",
      "Working on ID: P13489 this is entry 9600/16236\n",
      "Working on ID: Q6P587 this is entry 9800/16236\n"
     ]
    }
   ],
   "source": [
    "counter = 5000\n",
    "df = reinitalise_df(True)\n",
    "for i in ids[5000:10000]:\n",
    "    if counter%200 == 0:\n",
    "        print(f'Working on ID: {i} this is entry {counter}/{len(ids)}')\n",
    "    df, fail = generate_table_rows_for_uniprot_id(i,df, annotate=True)\n",
    "    if fail != None:\n",
    "        fail_list.append(fail)\n",
    "    counter = counter+1\n",
    "df.to_csv('part_02_annotated.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "12343250",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on ID: Q8N5Y2 this is entry 10000/16236\n",
      "Issue: No gene found found\n",
      "Issue: No gene found found\n",
      "Working on ID: Q711Q0 this is entry 10200/16236\n",
      "Working on ID: Q9H063 this is entry 10400/16236\n",
      "Working on ID: Q9HCK1 this is entry 10600/16236\n",
      "Working on ID: O60774 this is entry 10800/16236\n",
      "Working on ID: Q8WUK0 this is entry 11000/16236\n",
      "Working on ID: P0DP08 this is entry 11200/16236\n",
      "Working on ID: Q86V35 this is entry 11400/16236\n",
      "Working on ID: O15013 this is entry 11600/16236\n",
      "Working on ID: Q8NGG8 this is entry 11800/16236\n",
      "Issue: No gene found found\n",
      "Issue: No gene found found\n",
      "Working on ID: Q9P2M7 this is entry 12000/16236\n",
      "Working on ID: P01584 this is entry 12200/16236\n",
      "Working on ID: Q6PKG0 this is entry 12400/16236\n",
      "Working on ID: Q08477 this is entry 12600/16236\n",
      "Working on ID: Q3SYA9 this is entry 12800/16236\n",
      "Working on ID: P25789 this is entry 13000/16236\n",
      "Working on ID: A0AV96 this is entry 13200/16236\n",
      "Working on ID: Q6PH85 this is entry 13400/16236\n",
      "Working on ID: Q8NHH9 this is entry 13600/16236\n",
      "Working on ID: Q8WXT5 this is entry 13800/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: Q16611 this is entry 14000/16236\n",
      "Working on ID: P20941 this is entry 14200/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: A8CG34 this is entry 14400/16236\n",
      "Working on ID: P18825 this is entry 14600/16236\n",
      "Working on ID: Q9P2X7 this is entry 14800/16236\n",
      "Issue: No gene found found\n"
     ]
    }
   ],
   "source": [
    "counter = 10000\n",
    "df = reinitalise_df(True)\n",
    "for i in ids[10000:15000]:\n",
    "    if counter%200 == 0:\n",
    "        print(f'Working on ID: {i} this is entry {counter}/{len(ids)}')\n",
    "    df, fail = generate_table_rows_for_uniprot_id(i,df, annotate=True)\n",
    "    if fail != None:\n",
    "        fail_list.append(fail)\n",
    "    counter = counter+1\n",
    "df.to_csv('part_03_annotated.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "id": "49ae8661",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on ID: P52294 this is entry 15000/16236\n",
      "Working on ID: Q4G0T1 this is entry 15200/16236\n",
      "Working on ID: Q9NPH3 this is entry 15400/16236\n",
      "Issue: No gene found found\n",
      "Working on ID: Q9BTL3 this is entry 15600/16236\n",
      "Working on ID: Q9H8M1 this is entry 15800/16236\n",
      "Working on ID: Q6VMQ6 this is entry 16000/16236\n",
      "Working on ID: Q5VSL9 this is entry 16200/16236\n"
     ]
    }
   ],
   "source": [
    "counter = 15000\n",
    "df = reinitalise_df(True)\n",
    "for i in ids[15000:]:\n",
    "    if counter%200 == 0:\n",
    "        print(f'Working on ID: {i} this is entry {counter}/{len(ids)}')\n",
    "    df, fail = generate_table_rows_for_uniprot_id(i,df,annotate=True)\n",
    "    if fail != None:\n",
    "        fail_list.append(fail)\n",
    "    counter = counter+1\n",
    "df.to_csv('part_04_annotated.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "fdd9dadc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Q3KNS1', 'Q96MI6', 'Q99102', 'Q6ZR37', 'P36639', 'Q12983', 'Q9NX55', 'P47756', 'Q96AX9', 'A6NMZ5', 'Q5H9F3', 'Q3ZCN5', 'B4E2M5', 'O43451', 'O43309', 'P17019', 'P0DN82', 'Q13507', 'P02708', 'P20848', 'P26599', 'Q9BV97', 'Q8NDH2', 'Q9BRI3', 'Q8NGU1', 'P23109', 'Q2VWA4', 'Q16635', 'Q6NY19', 'Q8N859', 'Q92782', 'P59826', 'Q6ZTN6', 'P0DMP1', 'Q8N7C3', 'Q6P3W6', 'O15131', 'Q9UPX8', 'Q6ZP01', 'Q8NFD5', 'Q7Z7G0', 'O94851', 'A6NNF4', 'Q8NH41', 'P0C623', 'Q9HD15', 'Q8IZM8', 'O15016', 'Q8N782', 'Q8IV03', 'Q9H310', 'Q03518', 'Q66PJ3', 'Q8TAX9', 'Q8N423', 'O94901', 'Q9HBL0', 'Q3I5F7', 'Q96Q27', 'Q9BY31', 'Q15170', 'Q9UJX3', 'Q12766', 'D6RF30', 'P29973', 'A0A0U1RRI6', 'Q0D2K0', 'Q9P0K9', 'Q9P2M7', 'Q96KE9', 'Q15270', 'A8MZ97', 'P62861', 'Q9Y5I7', 'Q9HBW0', 'Q5XKL5', 'Q53TS8', 'O94880', 'Q16850', 'Q6ZMS7', 'Q8NBS3', 'Q6ZMR5', 'P51606', 'Q9NRJ5', 'Q96FT7', 'Q6PKH6']\n"
     ]
    }
   ],
   "source": [
    "print(fail_list)\n",
    "np.savetxt('data/failed_AF_structures.csv', fail_list, fmt='%s')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87c6870f",
   "metadata": {},
   "source": [
    "## Playground testing things out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "8589316d",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_filter = []\n",
    "f = open('domain.csv')\n",
    "f_content = f.readlines()\n",
    "f.close()\n",
    "\n",
    "for l in f_content:\n",
    "    domain_filter.append(l.strip())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "ea67df64",
   "metadata": {},
   "outputs": [],
   "source": [
    "region_filter = []\n",
    "f = open('region_filter.csv')\n",
    "f_content = f.readlines()\n",
    "f.close()\n",
    "\n",
    "for l in f_content:\n",
    "    region_filter.append(l.strip())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "68939b75",
   "metadata": {},
   "outputs": [],
   "source": [
    "discard_list = np.hstack((region_filter,domain_filter))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "a8d8f689",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ARM-like and Heat-like helical repeats'"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discard_list[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "b4523df4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on uniprot ID Q9H4S2\n",
      "At entry 0/1\n",
      "ID: Q9H4S2\n",
      "got secondary structure\n",
      "['-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' 'T' 'T' '-' '-' 'T' 'T' 'S'\n",
      " '-' '-' '-' '-' 'H' 'H' 'H' 'H' 'H' 'H' 'T' 'T' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' 'T' 'T' 'S' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' 'H' 'H' 'H' 'H' 'H' 'H' 'H'\n",
      " 'H' 'H' 'H' 'H' 'H' 'H' '-' 'S' 'S' '-' '-' 'H' 'H' 'H' 'H' 'H' 'H' 'H'\n",
      " 'H' 'H' 'H' 'T' 'T' '-' '-' 'H' 'H' 'H' 'H' 'H' 'H' 'H' 'H' 'H' 'H' 'H'\n",
      " 'H' 'H' 'H' 'H' 'H' 'H' 'H' 'T' 'T' 'T' 'T' 'T' 'T' '-' '-' 'S' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-'\n",
      " '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-' '-']\n",
      "read json file successfully\n"
     ]
    }
   ],
   "source": [
    "uni_prot_id = 'Q9H4S2'\n",
    "work_dir = 'temp'\n",
    "print('Working on uniprot ID', uni_prot_id)\n",
    "# Downloading uniprot file and alphafold file\n",
    "download_uniprot_json_file(uni_prot_id, work_dir)\n",
    "download_alpha_fold_pdbs([uni_prot_id], workdir =work_dir)\n",
    "p = 'temp/AF-'+uni_prot_id+'-F1-model_v2.pdb'\n",
    "sec_struc_output = get_sec_struct(p, '/Users/toni_brain/miniconda3/envs/dssp//bin/mkdssp')\n",
    "print('got secondary structure')\n",
    "print(sec_struc_output)\n",
    "\n",
    "# Process uniprot file\n",
    "f = open('temp/'+uni_prot_id+'.json')\n",
    "data = json.load(f)\n",
    "f.close()\n",
    "\n",
    "print('read json file successfully')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "a9cd903f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', 'T', 'T', '-', '-', 'T',\n",
       "       'T', 'S', '-', '-', '-', '-', 'H', 'H', 'H', 'H', 'H', 'H', 'T',\n",
       "       'T', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', 'T', 'T',\n",
       "       'S', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', 'H',\n",
       "       'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', '-',\n",
       "       'S', 'S', '-', '-', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H',\n",
       "       'H', 'T', 'T', '-', '-', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H',\n",
       "       'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'T', 'T',\n",
       "       'T', 'T', 'T', '-', '-', 'S', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-',\n",
       "       '-', '-', '-', '-'], dtype='<U1')"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sec_struc_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "7c74fbdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is collecting known information\n",
    "\n",
    "known_things = {}\n",
    "known_things['Domain'] = []\n",
    "known_things['Region'] = []\n",
    "known_things['Motif'] = []\n",
    "#known_things['Helix'] = []\n",
    "\n",
    "for d in data['features']:\n",
    "    #print(d['type'])\n",
    "    if d['type'] == 'Domain':\n",
    "        known_things['Domain'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "    elif d['type'] == 'Region':\n",
    "        known_things['Region'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "    elif d['type'] == 'Motif':\n",
    "        known_things['Motif'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])\n",
    "    #elif d['type'] == 'Helix':\n",
    "    #    known_things['Helix'].append([d['location']['start']['value'], d['location']['end']['value'], d['description']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "a0f94460",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Domain': [],\n",
       " 'Region': [[1, 20, 'SNAG domain'], [201, 264, 'Disordered']],\n",
       " 'Motif': []}"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "known_things"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "5ab8af6a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'type': 'Region', 'location': {'start': {'value': 1, 'modifier': 'EXACT'}, 'end': {'value': 20, 'modifier': 'EXACT'}}, 'description': 'SNAG domain', 'evidences': [{'evidenceCode': 'ECO:0000250'}]}\n",
      "{'type': 'Region', 'location': {'start': {'value': 201, 'modifier': 'EXACT'}, 'end': {'value': 264, 'modifier': 'EXACT'}}, 'description': 'Disordered', 'evidences': [{'evidenceCode': 'ECO:0000256', 'source': 'SAM', 'id': 'MobiDB-lite'}]}\n"
     ]
    }
   ],
   "source": [
    "for d in data['features']:\n",
    "    if d['type']=='Domain':\n",
    "        print (d)\n",
    "    if d['type']=='Region':\n",
    "        print (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8a79b93d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n"
     ]
    }
   ],
   "source": [
    "sequence_length = data['sequence']['length']\n",
    "if sequence_length == len(sec_struc_output):\n",
    "    print('yes')\n",
    "    bool_array =  np.zeros(sequence_length, dtype=bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "19dd0d2a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106], [121, 122, 123, 124, 125, 126, 127, 128, 129, 130], [139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153], [175, 176, 177, 178, 179, 180, 181, 182, 183]]\n"
     ]
    }
   ],
   "source": [
    "bool_array =  np.zeros(sequence_length, dtype=bool)\n",
    "discard_list = ['bHLH', 'Leucine-zipper']\n",
    "alpha_helix_index_list = get_alpha_helix_length_and_location(sec_struc_output, min_length=7)\n",
    "helix_list_remove_index = []\n",
    "for key in known_things.keys():\n",
    "    for entry in known_things[key]:\n",
    "        if entry[-1] in discard_list:\n",
    "            bool_array[entry[0]-1:entry[1]] = True\n",
    "            index = 0\n",
    "            for a in alpha_helix_index_list:\n",
    "                curr_helix =  np.zeros(sequence_length, dtype=bool)\n",
    "                curr_helix[a] = True\n",
    "                overalp = (np.sum(np.logical_and(bool_array, curr_helix)))\n",
    "                if overalp>4:\n",
    "                    helix_list_remove_index.append(index)\n",
    "                index = index+1\n",
    "unique_idx = np.unique(helix_list_remove_index)\n",
    "for index in sorted(unique_idx, reverse=True):\n",
    "    del alpha_helix_index_list[index]\n",
    "print (alpha_helix_index_list)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}