do_DGE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running DESeq2 and edgeR and calculating Bayes factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# This python file is calling R scripts with a list of dependencies. To make sure you have all of them ready you can simply create a conda environment using our explorted environment\n",
    "\n",
    "# conda env create -f environment.yml\n",
    "\n",
    "# If you are not using conda, the specific requirements here are\n",
    "# R, statmod, tidyverse, bioconductor, genomeinfodbdata, limma and of course edgeR, DESeq2\n",
    "# http://bioconductor.org/install/\n",
    "# https://bioconductor.org/packages/release/bioc/html/edgeR.html\n",
    "# https://bioconductor.org/packages/release/bioc/html/DESeq2.html\n",
    "# https://bioconductor.org/packages/release/bioc/html/limma.html\n",
    "# https://anaconda.org/bioconda/bioconductor-genomeinfodbdata\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 100 boostrap iterations á #reps ∈ {3,6,12,20}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running edgeR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"3R_100/R3_{i}.csv\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R3_{i}_edgeR.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_edgeR_R3.R {DATASET_NAME} {RESULTS_NAME} 0.0\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"6R_100/R6_{i}.csv\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R6_{i}_edgeR.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_edgeR_R6.R {DATASET_NAME} {RESULTS_NAME} 0.0\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 12 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"12R_100/R12_{i}.csv\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R12_{i}_edgeR.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_edgeR_R12.R {DATASET_NAME} {RESULTS_NAME} 0.0\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 20 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"20R_100/R20_{i}.csv\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R20_{i}_edgeR.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_edgeR_R20.R {DATASET_NAME} {RESULTS_NAME} 0.0\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ALL REPLICATES\n",
    "\n",
    "\n",
    "DATASET_NAME = f\"RALL.csv\"\n",
    "\n",
    "RESULTS_NAME = f\"DGE_results/RALL_edgeR.csv\"\n",
    "\n",
    "os.system(f\"Rscript Do_DGE_edgeR_RALL.R {DATASET_NAME} {RESULTS_NAME} 0.0\")\n",
    "\n",
    "print(f\"Done ALL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running DESeq2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"3R_100/R3_{i}.csv\"\n",
    "\n",
    "    METADATA_NAME = f\"3R_100/R3_{i}_meta.txt\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R3_{i}_DESeq2.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_DESeq2.R {DATASET_NAME} {METADATA_NAME} {RESULTS_NAME}\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"6R_100/R6_{i}.csv\"\n",
    "\n",
    "    METADATA_NAME = f\"6R_100/R6_{i}_meta.txt\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R6_{i}_DESeq2.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_DESeq2.R {DATASET_NAME} {METADATA_NAME} {RESULTS_NAME}\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 12 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"12R_100/R12_{i}.csv\"\n",
    "\n",
    "    METADATA_NAME = f\"12R_100/R12_{i}_meta.txt\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R12_{i}_DESeq2.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_DESeq2.R {DATASET_NAME} {METADATA_NAME} {RESULTS_NAME}\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 20 REPLICATES\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    DATASET_NAME = f\"20R_100/R20_{i}.csv\"\n",
    "\n",
    "    METADATA_NAME = f\"20R_100/R20_{i}_meta.txt\"\n",
    "\n",
    "    RESULTS_NAME = f\"DGE_results/R20_{i}_DESeq2.csv\"\n",
    "\n",
    "    os.system(f\"Rscript Do_DGE_DESeq2.R {DATASET_NAME} {METADATA_NAME} {RESULTS_NAME}\")\n",
    "\n",
    "    print(f\"Done {i}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all reps\n",
    "\n",
    "\n",
    "DATASET_NAME = f\"RALL.csv\"\n",
    "\n",
    "METADATA_NAME = f\"RALL_meta.txt\"\n",
    "\n",
    "RESULTS_NAME = f\"DGE_results/RALL_DESeq2.csv\"\n",
    "\n",
    "os.system(f\"Rscript Do_DGE_DESeq2.R {DATASET_NAME} {METADATA_NAME} {RESULTS_NAME}\")\n",
    "\n",
    "print(f\"Done ALL\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# bayexpress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.special as sc\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bayexpress_functions import get_BF_21, get_FC, get_q, get_BF_k1\n",
    "\n",
    "# setting the priors to flat prios\n",
    "u_1 = 1\n",
    "u_2 = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all reps\n",
    "\n",
    "in_data = pd.read_csv(f\"RALL.csv\")\n",
    "\n",
    "out_data = pd.DataFrame({'genes': in_data.genes})\n",
    "\n",
    "n_1 = in_data.iloc[:,1:43].sum(axis=1)\n",
    "n_2 = in_data.iloc[:,43:].sum(axis=1)\n",
    "\n",
    "N_1 = in_data.iloc[:,1:43].sum(axis=1).sum()\n",
    "N_2 = in_data.iloc[:,43:].sum(axis=1).sum()\n",
    "\n",
    "out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "out_data.to_csv(f\"DGE_results/RALL_bayexpress.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 REPLICATES\n",
    "\n",
    "# summing up over replicates\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    in_data = pd.read_csv(f\"3R_100/R3_{i}.csv\")\n",
    "\n",
    "    out_data = pd.DataFrame({'locus_name': in_data.locus_name})\n",
    "\n",
    "    n_1 = in_data.iloc[:,1:4].sum(axis=1)\n",
    "    n_2 = in_data.iloc[:,4:].sum(axis=1)\n",
    "\n",
    "    N_1 = in_data.iloc[:,1:4].sum(axis=1).sum()\n",
    "    N_2 = in_data.iloc[:,4:].sum(axis=1).sum()\n",
    "\n",
    "    out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "    out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "    out_data.to_csv(f\"DGE_results/R3_{i}_bayexpress.csv\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6 REPLICATES\n",
    "\n",
    "# summing up over replicates\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    in_data = pd.read_csv(f\"6R_100/R6_{i}.csv\")\n",
    "\n",
    "    out_data = pd.DataFrame({'locus_name': in_data.locus_name})\n",
    "\n",
    "    n_1 = in_data.iloc[:,1:7].sum(axis=1)\n",
    "    n_2 = in_data.iloc[:,7:].sum(axis=1)\n",
    "\n",
    "    N_1 = in_data.iloc[:,1:7].sum(axis=1).sum()\n",
    "    N_2 = in_data.iloc[:,7:].sum(axis=1).sum()\n",
    "\n",
    "    out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "    out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "    out_data.to_csv(f\"DGE_results/R6_{i}_bayexpress.csv\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 12 REPLICATES\n",
    "\n",
    "# summing up over replicates\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    in_data = pd.read_csv(f\"12R_100/R12_{i}.csv\")\n",
    "\n",
    "    out_data = pd.DataFrame({'locus_name': in_data.locus_name})\n",
    "\n",
    "    n_1 = in_data.iloc[:,1:13].sum(axis=1)\n",
    "    n_2 = in_data.iloc[:,13:].sum(axis=1)\n",
    "\n",
    "    N_1 = in_data.iloc[:,1:13].sum(axis=1).sum()\n",
    "    N_2 = in_data.iloc[:,13:].sum(axis=1).sum()\n",
    "\n",
    "    out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "    out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "    out_data.to_csv(f\"DGE_results/R12_{i}_bayexpress.csv\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 20 REPLICATES\n",
    "\n",
    "# summing up over replicates\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    in_data = pd.read_csv(f\"20R_100/R20_{i}.csv\")\n",
    "\n",
    "    out_data = pd.DataFrame({'locus_name': in_data.locus_name})\n",
    "\n",
    "    n_1 = in_data.iloc[:,1:21].sum(axis=1)\n",
    "    n_2 = in_data.iloc[:,21:].sum(axis=1)\n",
    "\n",
    "    N_1 = in_data.iloc[:,1:21].sum(axis=1).sum()\n",
    "    N_2 = in_data.iloc[:,21:].sum(axis=1).sum()\n",
    "\n",
    "    out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "    out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "    out_data.to_csv(f\"DGE_results/R20_{i}_bayexpress.csv\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Control experiment ... what if we only use the WT in the bootstrapping?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 REPLICATES\n",
    "# CONTROL EXPERIMENT\n",
    "\n",
    "# summing up over replicates\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    in_data = pd.read_csv(f\"6R_100/R6_{i}.csv\")\n",
    "\n",
    "    out_data = pd.DataFrame({'locus_name': in_data.locus_name})\n",
    "\n",
    "    n_1 = in_data.iloc[:,1:4].sum(axis=1)\n",
    "    n_2 = in_data.iloc[:,4:7].sum(axis=1)\n",
    "\n",
    "    N_1 = in_data.iloc[:,1:4].sum(axis=1).sum()\n",
    "    N_2 = in_data.iloc[:,4:7].sum(axis=1).sum()\n",
    "\n",
    "    out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "    out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "    out_data['BF_k1_1'] = get_BF_k1(in_data.iloc[:,1:4])\n",
    "    out_data['BF_k1_2'] = get_BF_k1(in_data.iloc[:,4:7])\n",
    "\n",
    "    out_data.to_csv(f\"DGE_results/CONTROL_R3_{i}_bayexpress.csv\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10 REPLICATES\n",
    "# CONTROL EXPERIMENT\n",
    "\n",
    "# summing up over replicates\n",
    "\n",
    "\n",
    "for i in range(1,101):\n",
    "\n",
    "    in_data = pd.read_csv(f\"20R_100/R20_{i}.csv\")\n",
    "\n",
    "    out_data = pd.DataFrame({'locus_name': in_data.locus_name})\n",
    "\n",
    "    n_1 = in_data.iloc[:,1:11].sum(axis=1)\n",
    "    n_2 = in_data.iloc[:,11:21].sum(axis=1)\n",
    "\n",
    "    N_1 = in_data.iloc[:,1:11].sum(axis=1).sum()\n",
    "    N_2 = in_data.iloc[:,11:21].sum(axis=1).sum()\n",
    "\n",
    "    out_data['BF_21'] = get_BF_21(N_1, n_1, N_2, n_2)\n",
    "    out_data['FC'] = get_FC(N_1, n_1, N_2, n_2)\n",
    "\n",
    "    out_data['BF_k1_1'] = get_BF_k1(in_data.iloc[:,1:4])\n",
    "    out_data['BF_k1_2'] = get_BF_k1(in_data.iloc[:,4:7])\n",
    "\n",
    "    out_data.to_csv(f\"DGE_results/CONTROL_R10_{i}_bayexpress.csv\")\n",
    "\n",
    "    print(f\"Done {i}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bayexpressbase",
   "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.13.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}