yeast_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "from Bio import SeqIO as seqio\n",
    "from kinaid.utility import Utility\n",
    "import gzip\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clean_phosphopeptides(phosphorylation_file, fasta_file, output_file) :\n",
    "    exp_df = pd.read_csv(phosphorylation_file)\n",
    "\n",
    "    #read in the fasta file\n",
    "    fasta_dict = seqio.to_dict(seqio.parse(fasta_file, 'fasta'))\n",
    "\n",
    "    exp_df = exp_df[exp_df['p_residue'].isin(['S', 'T', 'Y'])]\n",
    "\n",
    "    systematic_name_to_SGD_dict = {name:record.description.split(', ')[0].split(' ')[2] for name,record in fasta_dict.items()}\n",
    "    #systematic_name_to_SGD_dict = {name:SGD.split(':')[1] for name,SGD in systematic_name_to_SGD_dict.items()}\n",
    "\n",
    "    exp_df['SGDID'] = exp_df['systematic_name'].map(systematic_name_to_SGD_dict)\n",
    "\n",
    "    systematic_name_to_seq_dict = {name:str(record.seq) for name,record in fasta_dict.items()}\n",
    "    exp_df['sequence'] = (exp_df['systematic_name'].map(systematic_name_to_seq_dict))\n",
    "\n",
    "    #remove rows with missing sequences\n",
    "    exp_df = exp_df[~exp_df['sequence'].isna()]\n",
    "\n",
    "    #remove rows with stop codons\n",
    "    exp_df['sequence'] = exp_df['sequence'].str.rstrip('*')\n",
    "\n",
    "    #get sequence lengths\n",
    "    exp_df['seq_len'] = exp_df['sequence'].str.len()\n",
    "\n",
    "    #keep only rows where the sequence length is greater than or equal to the position of the phosphorylation site\n",
    "    exp_df = exp_df[exp_df['seq_len'] >= exp_df['p_position']]\n",
    "\n",
    "    exp_df['sanity_check'] = exp_df.apply(lambda row: row['sequence'][row['p_position']-1] == row['p_residue'], axis=1)\n",
    "    \n",
    "    exp_df = exp_df[exp_df['sanity_check']]\n",
    "\n",
    "    #if p_position is less than 6, add 5 - p_position number of _'s to the beginning of the sequence\n",
    "    exp_df['sequence'] = exp_df.apply(lambda row: '_'*(6 - row['p_position']) + row['sequence'], axis=1)\n",
    "\n",
    "    #if p_position is greater than seq_len - 6, add p_position - (seq_len - 5) number of _'s to the end of the sequence\n",
    "    exp_df['sequence'] = exp_df.apply(lambda row: row['sequence'] + '_'*(row['p_position'] - (row['seq_len'] - 5)), axis=1)\n",
    "\n",
    "    #make window around p_position if p_position is greater than 5 from p_position - 6 to p_position + 5, else from 0 to 10\n",
    "    exp_df['window'] = exp_df.apply(lambda row: row['sequence'][row['p_position']-6:row['p_position']+5] if row['p_position'] > 5 else row['sequence'][0:11], axis=1)\n",
    "\n",
    "    exp_df['adj_p_value'] = exp_df['adj_p_value'].apply(lambda x: -1*np.log10(x))\n",
    "\n",
    "    exp_df_filtered = exp_df[['SGDID', 'p_position', 'window', 'fc_log2', 'adj_p_value']].copy()\n",
    "\n",
    "    exp_df_filtered.rename(columns={'SGDID':'SGD', 'p_position':'site', 'window':'peptide', 'fc_log2':'log2fc', 'adj_p_value':'adjpvalue'}, inplace=True)\n",
    "\n",
    "    exp_df_filtered.to_csv(output_file, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_dir = './test'\n",
    "proteomes_dir = './proteomes'\n",
    "published_yeast_file = 'yeast-perturbation.xlsx'\n",
    "published_yeast_path = os.path.join(test_dir, published_yeast_file)\n",
    "yeast_proteome_path = os.path.join(proteomes_dir, 'yeast_proteome.fasta')\n",
    "\n",
    "if not os.path.exists(proteomes_dir):\n",
    "    os.makedirs(proteomes_dir)\n",
    "    \n",
    "if not os.path.exists(yeast_proteome_path):\n",
    "    print('Downloading yeast proteome')\n",
    "    Utility.download_file('http://sgd-archive.yeastgenome.org/sequence/S288C_reference/orf_protein/orf_trans.fasta.gz', 'temp.gz')\n",
    "    with gzip.open('temp.gz', 'rb') as f_in:\n",
    "        with open(yeast_proteome_path, 'wb') as f_out:\n",
    "            f_out.write(f_in.read())\n",
    "    os.remove('temp.gz')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>reference</th>\n",
       "      <th>systematic_name</th>\n",
       "      <th>gene</th>\n",
       "      <th>p_site</th>\n",
       "      <th>p_residue</th>\n",
       "      <th>p_position</th>\n",
       "      <th>treatment_id</th>\n",
       "      <th>fc_log2</th>\n",
       "      <th>p_value</th>\n",
       "      <th>adj_p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>YAL012W_S40</td>\n",
       "      <td>YAL012W</td>\n",
       "      <td>CYS3</td>\n",
       "      <td>S40</td>\n",
       "      <td>S</td>\n",
       "      <td>40</td>\n",
       "      <td>LP</td>\n",
       "      <td>0.219367</td>\n",
       "      <td>7.189179e-01</td>\n",
       "      <td>0.934251</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>YAL012W_S40</td>\n",
       "      <td>YAL012W</td>\n",
       "      <td>CYS3</td>\n",
       "      <td>S40</td>\n",
       "      <td>S</td>\n",
       "      <td>40</td>\n",
       "      <td>DP</td>\n",
       "      <td>0.876721</td>\n",
       "      <td>1.824667e-01</td>\n",
       "      <td>0.608548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>YAL012W_S40</td>\n",
       "      <td>YAL012W</td>\n",
       "      <td>CYS3</td>\n",
       "      <td>S40</td>\n",
       "      <td>S</td>\n",
       "      <td>40</td>\n",
       "      <td>SP</td>\n",
       "      <td>3.562723</td>\n",
       "      <td>1.212967e-07</td>\n",
       "      <td>0.000013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>83</th>\n",
       "      <td>YAL012W_S40</td>\n",
       "      <td>YAL012W</td>\n",
       "      <td>CYS3</td>\n",
       "      <td>S40</td>\n",
       "      <td>S</td>\n",
       "      <td>40</td>\n",
       "      <td>YP</td>\n",
       "      <td>-1.092939</td>\n",
       "      <td>1.295494e-01</td>\n",
       "      <td>0.528078</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>YAL013W_S2</td>\n",
       "      <td>YAL013W</td>\n",
       "      <td>DEP1</td>\n",
       "      <td>S2</td>\n",
       "      <td>S</td>\n",
       "      <td>2</td>\n",
       "      <td>LP</td>\n",
       "      <td>-0.030479</td>\n",
       "      <td>9.352545e-01</td>\n",
       "      <td>0.987880</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",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>529283</th>\n",
       "      <td>YMR005W_S49</td>\n",
       "      <td>YMR005W</td>\n",
       "      <td>TAF4</td>\n",
       "      <td>S49</td>\n",
       "      <td>S</td>\n",
       "      <td>49</td>\n",
       "      <td>YP</td>\n",
       "      <td>0.167245</td>\n",
       "      <td>9.050468e-01</td>\n",
       "      <td>0.981903</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>529344</th>\n",
       "      <td>YDR358W_S353</td>\n",
       "      <td>YDR358W</td>\n",
       "      <td>GGA1</td>\n",
       "      <td>S353</td>\n",
       "      <td>S</td>\n",
       "      <td>353</td>\n",
       "      <td>LP</td>\n",
       "      <td>2.012639</td>\n",
       "      <td>5.130426e-02</td>\n",
       "      <td>0.335476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>529348</th>\n",
       "      <td>YDR358W_S353</td>\n",
       "      <td>YDR358W</td>\n",
       "      <td>GGA1</td>\n",
       "      <td>S353</td>\n",
       "      <td>S</td>\n",
       "      <td>353</td>\n",
       "      <td>DP</td>\n",
       "      <td>1.351322</td>\n",
       "      <td>2.241834e-01</td>\n",
       "      <td>0.658840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>529370</th>\n",
       "      <td>YDR358W_S353</td>\n",
       "      <td>YDR358W</td>\n",
       "      <td>GGA1</td>\n",
       "      <td>S353</td>\n",
       "      <td>S</td>\n",
       "      <td>353</td>\n",
       "      <td>SP</td>\n",
       "      <td>-1.772703</td>\n",
       "      <td>1.143126e-01</td>\n",
       "      <td>0.499893</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>529383</th>\n",
       "      <td>YDR358W_S353</td>\n",
       "      <td>YDR358W</td>\n",
       "      <td>GGA1</td>\n",
       "      <td>S353</td>\n",
       "      <td>S</td>\n",
       "      <td>353</td>\n",
       "      <td>YP</td>\n",
       "      <td>1.678858</td>\n",
       "      <td>1.683864e-01</td>\n",
       "      <td>0.588790</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>21176 rows × 10 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           reference systematic_name  gene p_site p_residue  p_position  \\\n",
       "44       YAL012W_S40         YAL012W  CYS3    S40         S          40   \n",
       "48       YAL012W_S40         YAL012W  CYS3    S40         S          40   \n",
       "70       YAL012W_S40         YAL012W  CYS3    S40         S          40   \n",
       "83       YAL012W_S40         YAL012W  CYS3    S40         S          40   \n",
       "144       YAL013W_S2         YAL013W  DEP1     S2         S           2   \n",
       "...              ...             ...   ...    ...       ...         ...   \n",
       "529283   YMR005W_S49         YMR005W  TAF4    S49         S          49   \n",
       "529344  YDR358W_S353         YDR358W  GGA1   S353         S         353   \n",
       "529348  YDR358W_S353         YDR358W  GGA1   S353         S         353   \n",
       "529370  YDR358W_S353         YDR358W  GGA1   S353         S         353   \n",
       "529383  YDR358W_S353         YDR358W  GGA1   S353         S         353   \n",
       "\n",
       "       treatment_id   fc_log2       p_value  adj_p_value  \n",
       "44               LP  0.219367  7.189179e-01     0.934251  \n",
       "48               DP  0.876721  1.824667e-01     0.608548  \n",
       "70               SP  3.562723  1.212967e-07     0.000013  \n",
       "83               YP -1.092939  1.295494e-01     0.528078  \n",
       "144              LP -0.030479  9.352545e-01     0.987880  \n",
       "...             ...       ...           ...          ...  \n",
       "529283           YP  0.167245  9.050468e-01     0.981903  \n",
       "529344           LP  2.012639  5.130426e-02     0.335476  \n",
       "529348           DP  1.351322  2.241834e-01     0.658840  \n",
       "529370           SP -1.772703  1.143126e-01     0.499893  \n",
       "529383           YP  1.678858  1.683864e-01     0.588790  \n",
       "\n",
       "[21176 rows x 10 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_excel(published_yeast_path, sheet_name='p_site_diff_reg')\n",
    "\n",
    "\n",
    "#TOR1_experiments = ['CS18', 'CS23','HS37', 'HS42', 'HS48']\n",
    "TOR1_experiments = set(['DP', 'LP', 'SP', 'YP'])\n",
    "\n",
    "df_exps = df[df['treatment_id'].isin(TOR1_experiments)]\n",
    "\n",
    "display(df_exps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make directory called TOR1\n",
    "TOR1_dir = os.path.join(test_dir, 'TOR1')\n",
    "os.makedirs(TOR1_dir, exist_ok=True)\n",
    "\n",
    "for exp in TOR1_experiments:\n",
    "    exp_df = df_exps[df_exps['treatment_id'] == exp]\n",
    "    temp_file = os.path.join(TOR1_dir, 'temp.csv')\n",
    "    exp_df.to_csv(temp_file, index=False)\n",
    "    exp_file = os.path.join(TOR1_dir, f'TOR1_{exp}.csv')\n",
    "    clean_phosphopeptides(temp_file, yeast_proteome_path, exp_file)\n",
    "    os.remove(temp_file)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}