rnaseq_normalizations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "95910694",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fd9377a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# scaling factor, typically 1 million (hence the M in the normalization method names)\n",
    "# here set to 10 (instead of 1000000) for more handy examples\n",
    "SCALE = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "29f98cbc",
   "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></th>\n",
       "      <th>Rep1</th>\n",
       "      <th>Rep2</th>\n",
       "      <th>Rep3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th>length</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th>2</th>\n",
       "      <td>10</td>\n",
       "      <td>12</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <th>4</th>\n",
       "      <td>20</td>\n",
       "      <td>25</td>\n",
       "      <td>60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>8</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D</th>\n",
       "      <th>10</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Rep1  Rep2  Rep3\n",
       "gene length                  \n",
       "A    2         10    12    30\n",
       "B    4         20    25    60\n",
       "C    1          5     8    15\n",
       "D    10         0     0     1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example data, consider the columns as three biological replicates, i.e. they should have identical values, but\n",
    "# in reality, these numbers are heavily skewed\n",
    "\n",
    "counts = pd.DataFrame(data=[['A',2,10,12,30],\n",
    "                            ['B',4,20,25,60],\n",
    "                            ['C',1,5,8,15],\n",
    "                            ['D',10,0,0,1]],\n",
    "                      columns=['gene','length','Rep1', 'Rep2', 'Rep3']).set_index(['gene','length'])\n",
    "counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "663e9497",
   "metadata": {},
   "source": [
    "# RPKM\n",
    "**R**eads **P**er **K**ilobase **M**illion\n",
    "### 1) adjust for sequencing depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "688205a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Rep1     35\n",
       "Rep2     45\n",
       "Rep3    106\n",
       "dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute coverage per replicate\n",
    "counts.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9c4487c9",
   "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></th>\n",
       "      <th>Rep1</th>\n",
       "      <th>Rep2</th>\n",
       "      <th>Rep3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th>length</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th>2</th>\n",
       "      <td>0.285714</td>\n",
       "      <td>0.266667</td>\n",
       "      <td>0.283019</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <th>4</th>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.555556</td>\n",
       "      <td>0.566038</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <th>1</th>\n",
       "      <td>0.142857</td>\n",
       "      <td>0.177778</td>\n",
       "      <td>0.141509</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D</th>\n",
       "      <th>10</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.009434</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Rep1      Rep2      Rep3\n",
       "gene length                              \n",
       "A    2       0.285714  0.266667  0.283019\n",
       "B    4       0.571429  0.555556  0.566038\n",
       "C    1       0.142857  0.177778  0.141509\n",
       "D    10      0.000000  0.000000  0.009434"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# devide each gene/replicate read count by the replicate coverage\n",
    "counts / counts.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "70382867",
   "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></th>\n",
       "      <th>Rep1</th>\n",
       "      <th>Rep2</th>\n",
       "      <th>Rep3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th>length</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th>2</th>\n",
       "      <td>2.857143</td>\n",
       "      <td>2.666667</td>\n",
       "      <td>2.830189</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <th>4</th>\n",
       "      <td>5.714286</td>\n",
       "      <td>5.555556</td>\n",
       "      <td>5.660377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <th>1</th>\n",
       "      <td>1.428571</td>\n",
       "      <td>1.777778</td>\n",
       "      <td>1.415094</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D</th>\n",
       "      <th>10</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.094340</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Rep1      Rep2      Rep3\n",
       "gene length                              \n",
       "A    2       2.857143  2.666667  2.830189\n",
       "B    4       5.714286  5.555556  5.660377\n",
       "C    1       1.428571  1.777778  1.415094\n",
       "D    10      0.000000  0.000000  0.094340"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# scale coverage\n",
    "counts / (counts.sum() / SCALE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b164c7ad",
   "metadata": {},
   "source": [
    "### 2) adjust for gene length\n",
    "by dividing with gene length in kilo bases\n",
    "\n",
    "note that we \"transpose\" the count table to operate on genes instead of replicates and \"transpose\" afterwards to restore original orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "36f8c8b3",
   "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></th>\n",
       "      <th>Rep1</th>\n",
       "      <th>Rep2</th>\n",
       "      <th>Rep3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th>length</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th>2</th>\n",
       "      <td>1.428571</td>\n",
       "      <td>1.333333</td>\n",
       "      <td>1.415094</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <th>4</th>\n",
       "      <td>1.428571</td>\n",
       "      <td>1.388889</td>\n",
       "      <td>1.415094</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <th>1</th>\n",
       "      <td>1.428571</td>\n",
       "      <td>1.777778</td>\n",
       "      <td>1.415094</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D</th>\n",
       "      <th>10</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.009434</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Rep1      Rep2      Rep3\n",
       "gene length                              \n",
       "A    2       1.428571  1.333333  1.415094\n",
       "B    4       1.428571  1.388889  1.415094\n",
       "C    1       1.428571  1.777778  1.415094\n",
       "D    10      0.000000  0.000000  0.009434"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((counts / (counts.sum() / SCALE)).T / counts.index.get_level_values(1)).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abd029d7",
   "metadata": {},
   "source": [
    "# TPM\n",
    "**T**ranscripts **P**er kilobase **M**illion\n",
    "### 1) adjust for gene length\n",
    "by dividing with gene length in kilo bases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "736258b2",
   "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></th>\n",
       "      <th>Rep1</th>\n",
       "      <th>Rep2</th>\n",
       "      <th>Rep3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th>length</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th>2</th>\n",
       "      <td>5.0</td>\n",
       "      <td>6.00</td>\n",
       "      <td>15.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>6.25</td>\n",
       "      <td>15.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <th>1</th>\n",
       "      <td>5.0</td>\n",
       "      <td>8.00</td>\n",
       "      <td>15.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D</th>\n",
       "      <th>10</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Rep1  Rep2  Rep3\n",
       "gene length                  \n",
       "A    2        5.0  6.00  15.0\n",
       "B    4        5.0  6.25  15.0\n",
       "C    1        5.0  8.00  15.0\n",
       "D    10       0.0  0.00   0.1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(counts.T / counts.index.get_level_values(1)).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "732728af",
   "metadata": {},
   "source": [
    "### 2) adjust for sequencing depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e2b16b5a",
   "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></th>\n",
       "      <th>Rep1</th>\n",
       "      <th>Rep2</th>\n",
       "      <th>Rep3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th>length</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th>2</th>\n",
       "      <td>0.333333</td>\n",
       "      <td>0.296296</td>\n",
       "      <td>0.332594</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <th>4</th>\n",
       "      <td>0.333333</td>\n",
       "      <td>0.308642</td>\n",
       "      <td>0.332594</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <th>1</th>\n",
       "      <td>0.333333</td>\n",
       "      <td>0.395062</td>\n",
       "      <td>0.332594</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D</th>\n",
       "      <th>10</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.002217</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Rep1      Rep2      Rep3\n",
       "gene length                              \n",
       "A    2       0.333333  0.296296  0.332594\n",
       "B    4       0.333333  0.308642  0.332594\n",
       "C    1       0.333333  0.395062  0.332594\n",
       "D    10      0.000000  0.000000  0.002217"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = (counts.T / counts.index.get_level_values(1)).T\n",
    "x / x.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94332654",
   "metadata": {},
   "source": [
    "# DESeq2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c26ca8a3",
   "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>Sample #1</th>\n",
       "      <th>Sample #2</th>\n",
       "      <th>Sample #3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Gene1</th>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene2</th>\n",
       "      <td>2</td>\n",
       "      <td>6</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene3</th>\n",
       "      <td>33</td>\n",
       "      <td>55</td>\n",
       "      <td>200</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Sample #1  Sample #2  Sample #3\n",
       "gene                                  \n",
       "Gene1          0         10          4\n",
       "Gene2          2          6         12\n",
       "Gene3         33         55        200"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example data, consider the columns as three different conditions, i.e. they should have different values\n",
    "\n",
    "counts = pd.DataFrame(data=[['Gene1',0,10,4],\n",
    "                            ['Gene2',2,6,12],\n",
    "                            ['Gene3',33,55,200]],\n",
    "                      columns=['gene','Sample #1', 'Sample #2', 'Sample #3']).set_index(['gene'])\n",
    "counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d368cdf6",
   "metadata": {},
   "source": [
    "### 1) log transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e0c5a1e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/sjanssen/miniconda3/envs/notebookServer/lib/python3.10/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "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>Sample #1</th>\n",
       "      <th>Sample #2</th>\n",
       "      <th>Sample #3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Gene1</th>\n",
       "      <td>-inf</td>\n",
       "      <td>2.302585</td>\n",
       "      <td>1.386294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene2</th>\n",
       "      <td>0.693147</td>\n",
       "      <td>1.791759</td>\n",
       "      <td>2.484907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene3</th>\n",
       "      <td>3.496508</td>\n",
       "      <td>4.007333</td>\n",
       "      <td>5.298317</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Sample #1  Sample #2  Sample #3\n",
       "gene                                  \n",
       "Gene1       -inf   2.302585   1.386294\n",
       "Gene2   0.693147   1.791759   2.484907\n",
       "Gene3   3.496508   4.007333   5.298317"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "349801e5",
   "metadata": {},
   "source": [
    "### 2) per gene, average across samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "16dcbe91",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/sjanssen/miniconda3/envs/notebookServer/lib/python3.10/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "gene\n",
       "Gene1        -inf\n",
       "Gene2    1.656604\n",
       "Gene3    4.267386\n",
       "dtype: float64"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(counts).mean(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a2fd1bc",
   "metadata": {},
   "source": [
    "### 3) filter genes with -inf at average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1cb91112",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gene\n",
       "Gene2    1.656604\n",
       "Gene3    4.267386\n",
       "dtype: object"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "avg = np.log(counts).mean(axis=1).replace(-1*np.infty, None).dropna()\n",
    "avg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f257118",
   "metadata": {},
   "source": [
    "### 4) per cell, subtract average from log(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "35e10944",
   "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>Sample #1</th>\n",
       "      <th>Sample #2</th>\n",
       "      <th>Sample #3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Gene1</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene2</th>\n",
       "      <td>-0.963457</td>\n",
       "      <td>0.135155</td>\n",
       "      <td>0.828302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene3</th>\n",
       "      <td>-0.770878</td>\n",
       "      <td>-0.260053</td>\n",
       "      <td>1.030931</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Sample #1 Sample #2 Sample #3\n",
       "gene                               \n",
       "Gene1       NaN       NaN       NaN\n",
       "Gene2 -0.963457  0.135155  0.828302\n",
       "Gene3 -0.770878 -0.260053  1.030931"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(np.log(counts).T - avg).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23ffcc6a",
   "metadata": {},
   "source": [
    "### 5) calculate mean on the ratios for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d7b387a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/sjanssen/miniconda3/envs/notebookServer/lib/python3.10/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Sample #1   -0.867168\n",
       "Sample #2   -0.062449\n",
       "Sample #3    0.929617\n",
       "dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((np.log(counts).T - avg).T).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9abd7243",
   "metadata": {},
   "source": [
    "### 6) convert mean to normal space: exp(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9fe81382",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/sjanssen/miniconda3/envs/notebookServer/lib/python3.10/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Sample #1    0.420140\n",
       "Sample #2    0.939461\n",
       "Sample #3    2.533538\n",
       "dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(((np.log(counts).T - avg).T).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8f0099f",
   "metadata": {},
   "source": [
    "### 7) divide read counts by \"ratio mean\" scaling factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5c1e86cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/sjanssen/miniconda3/envs/notebookServer/lib/python3.10/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "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>Sample #1</th>\n",
       "      <th>Sample #2</th>\n",
       "      <th>Sample #3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Gene1</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>10.644401</td>\n",
       "      <td>1.578820</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene2</th>\n",
       "      <td>4.760321</td>\n",
       "      <td>6.386640</td>\n",
       "      <td>4.736459</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene3</th>\n",
       "      <td>78.545292</td>\n",
       "      <td>58.544204</td>\n",
       "      <td>78.940989</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Sample #1  Sample #2  Sample #3\n",
       "gene                                  \n",
       "Gene1   0.000000  10.644401   1.578820\n",
       "Gene2   4.760321   6.386640   4.736459\n",
       "Gene3  78.545292  58.544204  78.940989"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts / np.exp(((np.log(counts).T - avg).T).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b65b145a",
   "metadata": {},
   "source": [
    "### 7b) round to neareast integer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "89b436e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/sjanssen/miniconda3/envs/notebookServer/lib/python3.10/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "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>Sample #1</th>\n",
       "      <th>Sample #2</th>\n",
       "      <th>Sample #3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Gene1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene2</th>\n",
       "      <td>5.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene3</th>\n",
       "      <td>79.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>79.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Sample #1  Sample #2  Sample #3\n",
       "gene                                  \n",
       "Gene1        0.0       11.0        2.0\n",
       "Gene2        5.0        6.0        5.0\n",
       "Gene3       79.0       59.0       79.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round(counts / np.exp(((np.log(counts).T - avg).T).mean()))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ggmap",
   "language": "python",
   "name": "ggmap"
  },
  "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.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}