RNASeq_Becavin_2020.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNASeq pipeline analysis\n",
    "\n",
    "We will perform differential expression (DE) analysis of the RNA-Seq sequencing of Escherichia Coli K-12 in stressful condition. \n",
    "https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-6993/\n",
    "\n",
    "We compare here 2 biological conditions:\n",
    "- WT - Bacteria grown in liquid cultures\n",
    "- MazF - Bacteria grown in liquid cultures with 2 hours of MazF induction\n",
    "\n",
    "Biological protocol:\n",
    "After RNA extraction, ribosomal RNAs have been removed using Ribo-Zero rRNA Removal Kit. This step is critical as rRNA account for 90% to 95% of total RNA. Remaining RNA are sequenced using Illumina MiSeq. Paired-end sequecing is used so each dataset contain two fastq files. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset download"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download all sequencing files. For each biological condition there are 2 biological replicates, and for each replicate 2 fast files (paired-end sequencing).\n",
    "You should end up with 8 files in your folder `DataEcoli`:\n",
    "- WT_1_R1.fastq.gz\n",
    "- WT_1_R2.fastq.gz\n",
    "- WT_2_R1.fastq.gz\n",
    "- WT_2_R2.fastq.gz\n",
    "- MazF_1_R1.fastq.gz\n",
    "- MazF_1_R2.fastq.gz\n",
    "- MazF_2_R1.fastq.gz\n",
    "- MazF_2_R2.fastq.gz\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Process is interrupted.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "mkdir \"DataEcoli\"\n",
    "# Download WT dataset (replicate 1)\n",
    "wget -O \"DataEcoli/WT_1_R1.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/005/ERR2686025/ERR2686025_1.fastq.gz\n",
    "wget -O \"DataEcoli/WT_1_R2.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/005/ERR2686025/ERR2686025_2.fastq.gz\n",
    "# Download WT dataset (replicate 2)\n",
    "wget -O \"DataEcoli/WT_2_R1.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/006/ERR2686026/ERR2686026_1.fastq.gz\n",
    "wget -O \"DataEcoli/WT_2_R2.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/006/ERR2686026/ERR2686026_2.fastq.gz\n",
    "    \n",
    "# Download MazF dataset (replicate 1)\n",
    "wget -O \"DataEcoli/MazF_1_R1.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/007/ERR2686027/ERR2686027_1.fastq.gz\n",
    "wget -O \"DataEcoli/MazF_1_R2.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/007/ERR2686027/ERR2686027_2.fastq.gz\n",
    "# Download MazF dataset (replicate 2)\n",
    "wget -O \"DataEcoli/MazF_2_R1.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/008/ERR2686028/ERR2686028_1.fastq.gz\n",
    "wget -O \"DataEcoli/MazF_2_R2.fastq.gz\" ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR268/008/ERR2686028/ERR2686028_2.fastq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Genome download"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now to download the genome file of Escherichis coli K-12. We download first the sequence (fasta file) and then the annotation (gff file)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "mkdir: Genome: File exists\n",
      "--2019-11-19 22:57:08--  ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz\n",
      "           => « Genome/Ecoli_K12.fna.gz »\n",
      "Résolution de ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)… 2607:f220:41e:250::13, 130.14.250.10\n",
      "Connexion à ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|2607:f220:41e:250::13|:21… connecté.\n",
      "Ouverture de session en tant que anonymous… Session établie.\n",
      "==> SYST ... terminé.    ==> PWD ... terminé.\n",
      "==> TYPE I ... terminé.  ==> CWD (1) /genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2 ... terminé.\n",
      "==> SIZE GCF_000005845.2_ASM584v2_genomic.fna.gz ... 1379902\n",
      "==> EPSV ... terminé.    ==> RETR GCF_000005845.2_ASM584v2_genomic.fna.gz ... terminé.\n",
      "Taille : 1379902 (1,3M) (non certifiée)\n",
      "\n",
      "     0K .......... .......... .......... .......... ..........  3%  147K 9s\n",
      "    50K .......... .......... .......... .......... ..........  7%  340K 6s\n",
      "   100K .......... .......... .......... .......... .......... 11% 2,23M 4s\n",
      "   150K .......... .......... .......... .......... .......... 14% 55,8K 8s\n",
      "   200K .......... .......... .......... .......... .......... 18% 74,0M 6s\n",
      "   250K .......... .......... .......... .......... .......... 22%  485K 5s\n",
      "   300K .......... .......... .......... .......... .......... 25% 2,69M 4s\n",
      "   350K .......... .......... .......... .......... .......... 29%  213K 4s\n",
      "   400K .......... .......... .......... .......... .......... 33% 62,5M 4s\n",
      "   450K .......... .......... .......... .......... .......... 37%  413K 3s\n",
      "   500K .......... .......... .......... .......... .......... 40%  328K 3s\n",
      "   550K .......... .......... .......... .......... .......... 44%  265K 3s\n",
      "   600K .......... .......... .......... .......... .......... 48%  359K 3s\n",
      "   650K .......... .......... .......... .......... .......... 51%  461K 2s\n",
      "   700K .......... .......... .......... .......... .......... 55%  232K 2s\n",
      "   750K .......... .......... .......... .......... .......... 59%  198K 2s\n",
      "   800K .......... .......... .......... .......... .......... 63%  229K 2s\n",
      "   850K .......... .......... .......... .......... .......... 66%  177K 2s\n",
      "   900K .......... .......... .......... .......... .......... 70%  170K 2s\n",
      "   950K .......... .......... .......... .......... .......... 74%  339K 1s\n",
      "  1000K .......... .......... .......... .......... .......... 77% 77,4K 1s\n",
      "  1050K .......... .......... .......... .......... .......... 81%  322K 1s\n",
      "  1100K .......... .......... .......... .......... .......... 85%  338K 1s\n",
      "  1150K .......... .......... .......... .......... .......... 89%  490K 1s\n",
      "  1200K .......... .......... .......... .......... .......... 92%  428K 0s\n",
      "  1250K .......... .......... .......... .......... .......... 96%  354K 0s\n",
      "  1300K .......... .......... .......... .......... .......   100% 2,61M=5,2s\n",
      "\n",
      "2019-11-19 22:57:22 (259 KB/s) - « Genome/Ecoli_K12.fna.gz » sauvegardé [1379902]\n",
      "\n",
      "--2019-11-19 22:57:22--  ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.gtf.gz\n",
      "           => « Genome/Ecoli_K12.gtf.gz »\n",
      "Résolution de ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)… 2607:f220:41e:250::13, 130.14.250.10\n",
      "Connexion à ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|2607:f220:41e:250::13|:21… connecté.\n",
      "Ouverture de session en tant que anonymous… Session établie.\n",
      "==> SYST ... terminé.    ==> PWD ... terminé.\n",
      "==> TYPE I ... terminé.  ==> CWD (1) /genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2 ... terminé.\n",
      "==> SIZE GCF_000005845.2_ASM584v2_genomic.gtf.gz ... 483715\n",
      "==> EPSV ... terminé.    ==> RETR GCF_000005845.2_ASM584v2_genomic.gtf.gz ... terminé.\n",
      "Taille : 483715 (472K) (non certifiée)\n",
      "\n",
      "     0K .......... .......... .......... .......... .......... 10%  245K 2s\n",
      "    50K .......... .......... .......... .......... .......... 21%  414K 1s\n",
      "   100K .......... .......... .......... .......... .......... 31%  416K 1s\n",
      "   150K .......... .......... .......... .......... .......... 42%  263K 1s\n",
      "   200K .......... .......... .......... .......... .......... 52%  469K 1s\n",
      "   250K .......... .......... .......... .......... .......... 63%  287K 1s\n",
      "   300K .......... .......... .......... .......... .......... 74%  431K 0s\n",
      "   350K .......... .......... .......... .......... .......... 84%  486K 0s\n",
      "   400K .......... .......... .......... .......... .......... 95%  455K 0s\n",
      "   450K .......... .......... ..                              100%  266K=1,3s\n",
      "\n",
      "2019-11-19 22:57:25 (355 KB/s) - « Genome/Ecoli_K12.gtf.gz » sauvegardé [483715]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "mkdir Genome\n",
    "wget -O \"Genome/Ecoli_K12.fna.gz\" ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz\n",
    "wget -O \"Genome/Ecoli_K12.gtf.gz\" ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.gtf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncompress genome files, either by using gunzip or just clicking on it in your OS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gunzip -d Genome/Ecoli_K12.fna.gz\n",
    "gunzip -d Genome/Ecoli_K12.gtf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clean reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running mapping of reads we need to clean them by running \"trimming\" operation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "mkdir: FastQ: File exists\n",
      "sh: line 4: cutadapt: command not found\n",
      "sh: line 5: cutadapt: command not found\n",
      "sh: line 6: cutadapt: command not found\n",
      "sh: line 7: cutadapt: command not found\n"
     ]
    },
    {
     "ename": "CalledProcessError",
     "evalue": "Command 'b'mkdir FastQ\\nread1Adapter=AGATCGGAAGAGCACACGTCTGAACTCCAGTCA\\nread2Adapter=AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/WT_1_R1.fastq -p FastQ/WT_1_R2.fastq DataEcoli/WT_1_R1.fastq DataEcoli/WT_1_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/WT_2_R1.fastq -p FastQ/WT_2_R2.fastq DataEcoli/WT_2_R1.fastq DataEcoli/WT_2_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/MazF_1_R1.fastq -p FastQ/MazF_1_R2.fastq DataEcoli/MazF_1_R1.fastq DataEcoli/MazF_1_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/MazF_2_R1.fastq -p FastQ/MazF_2_R2.fastq DataEcoli/MazF_2_R1.fastq DataEcoli/MazF_2_R2.fastq\\n'' returned non-zero exit status 127.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m----------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mCalledProcessError\u001b[0m                   Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-33-5236800b99c1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'sh'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'mkdir FastQ\\nread1Adapter=AGATCGGAAGAGCACACGTCTGAACTCCAGTCA\\nread2Adapter=AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/WT_1_R1.fastq -p FastQ/WT_1_R2.fastq DataEcoli/WT_1_R1.fastq DataEcoli/WT_1_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/WT_2_R1.fastq -p FastQ/WT_2_R2.fastq DataEcoli/WT_2_R1.fastq DataEcoli/WT_2_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/MazF_1_R1.fastq -p FastQ/MazF_1_R2.fastq DataEcoli/MazF_1_R1.fastq DataEcoli/MazF_1_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/MazF_2_R1.fastq -p FastQ/MazF_2_R2.fastq DataEcoli/MazF_2_R1.fastq DataEcoli/MazF_2_R2.fastq\\n'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2321\u001b[0m             \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_expand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2322\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2323\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2324\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2325\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mnamed_script_magic\u001b[0;34m(line, cell)\u001b[0m\n\u001b[1;32m    140\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m                 \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscript\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 142\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshebang\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    143\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    144\u001b[0m         \u001b[0;31m# write a basic docstring:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-109>\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n\u001b[1;32m    243\u001b[0m             \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mraise_error\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m!=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 245\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mCalledProcessError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_run_script\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mto_close\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mCalledProcessError\u001b[0m: Command 'b'mkdir FastQ\\nread1Adapter=AGATCGGAAGAGCACACGTCTGAACTCCAGTCA\\nread2Adapter=AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/WT_1_R1.fastq -p FastQ/WT_1_R2.fastq DataEcoli/WT_1_R1.fastq DataEcoli/WT_1_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/WT_2_R1.fastq -p FastQ/WT_2_R2.fastq DataEcoli/WT_2_R1.fastq DataEcoli/WT_2_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/MazF_1_R1.fastq -p FastQ/MazF_1_R2.fastq DataEcoli/MazF_1_R1.fastq DataEcoli/MazF_1_R2.fastq\\ncutadapt -b $read1Adapter -B $read2Adapter -o FastQ/MazF_2_R1.fastq -p FastQ/MazF_2_R2.fastq DataEcoli/MazF_2_R1.fastq DataEcoli/MazF_2_R2.fastq\\n'' returned non-zero exit status 127."
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "mkdir FastQ\n",
    "read1Adapter=AGATCGGAAGAGCACACGTCTGAACTCCAGTCA\n",
    "read2Adapter=AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT\n",
    "cutadapt -j 8 -b $read1Adapter -B $read2Adapter -o FastQ/WT_1_R1.fastq -p FastQ/WT_1_R2.fastq DataEcoli/WT_1_R1.fastq.gz DataEcoli/WT_1_R2.fastq.gz\n",
    "cutadapt -j 8 -b $read1Adapter -B $read2Adapter -o FastQ/WT_2_R1.fastq -p FastQ/WT_2_R2.fastq DataEcoli/WT_2_R1.fastq.gz DataEcoli/WT_2_R2.fastq.gz\n",
    "cutadapt -j 8 -b $read1Adapter -B $read2Adapter -o FastQ/MazF_1_R1.fastq -p FastQ/MazF_1_R2.fastq DataEcoli/MazF_1_R1.fastq.gz DataEcoli/MazF_1_R2.fastq.gz\n",
    "cutadapt -j 8 -b $read1Adapter -B $read2Adapter -o FastQ/MazF_2_R1.fastq -p FastQ/MazF_2_R2.fastq DataEcoli/MazF_2_R1.fastq.gz DataEcoli/MazF_2_R2.fastq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check quality of the reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "bash: line 2: fastqc: command not found\n",
      "bash: line 3: fastqc: command not found\n"
     ]
    },
    {
     "ename": "CalledProcessError",
     "evalue": "Command 'b'mkdir FastQC\\nfastqc -o FastQC FastQ/WT_1_R1.fastq FastQ/WT_1_R2.fastq FastQ/WT_2_R1.fastq FastQ/WT_2_R2.fastq\\nfastqc -o FastQC FastQ/MazF_1_R1.fastq FastQ/MazF_1_R2.fastq FastQ/MazF_2_R1.fastq FastQ/MazF_2_R2.fastq\\n'' returned non-zero exit status 127.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m------------------------------------------\u001b[0m",
      "\u001b[0;31mCalledProcessError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-34-003414be57d7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'bash'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'mkdir FastQC\\nfastqc -o FastQC FastQ/WT_1_R1.fastq FastQ/WT_1_R2.fastq FastQ/WT_2_R1.fastq FastQ/WT_2_R2.fastq\\nfastqc -o FastQC FastQ/MazF_1_R1.fastq FastQ/MazF_1_R2.fastq FastQ/MazF_2_R1.fastq FastQ/MazF_2_R2.fastq\\n'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2321\u001b[0m             \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_expand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2322\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2323\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2324\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2325\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mnamed_script_magic\u001b[0;34m(line, cell)\u001b[0m\n\u001b[1;32m    140\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m                 \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscript\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 142\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshebang\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    143\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    144\u001b[0m         \u001b[0;31m# write a basic docstring:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-109>\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n\u001b[1;32m    243\u001b[0m             \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mraise_error\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m!=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 245\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mCalledProcessError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_run_script\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mto_close\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mCalledProcessError\u001b[0m: Command 'b'mkdir FastQC\\nfastqc -o FastQC FastQ/WT_1_R1.fastq FastQ/WT_1_R2.fastq FastQ/WT_2_R1.fastq FastQ/WT_2_R2.fastq\\nfastqc -o FastQC FastQ/MazF_1_R1.fastq FastQ/MazF_1_R2.fastq FastQ/MazF_2_R1.fastq FastQ/MazF_2_R2.fastq\\n'' returned non-zero exit status 127."
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "mkdir FastQC\n",
    "fastqc -t 8 -o FastQC FastQ/WT_1_R1.fastq FastQ/WT_1_R2.fastq FastQ/WT_2_R1.fastq FastQ/WT_2_R2.fastq\n",
    "fastqc -t 8 -o FastQC FastQ/MazF_1_R1.fastq FastQ/MazF_1_R2.fastq FastQ/MazF_2_R1.fastq FastQ/MazF_2_R2.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reads mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first create the index file for Escherichia Coli K12 genome. This is mandatory for running reads mapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "bash: line 1: bowtie2-build: command not found\n"
     ]
    },
    {
     "ename": "CalledProcessError",
     "evalue": "Command 'b'bowtie2-build Genome/Ecoli_K12.fna Genome/Ecoli_K12\\n'' returned non-zero exit status 127.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m----------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mCalledProcessError\u001b[0m                   Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-31-54bf9a092588>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'bash'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'bowtie2-build Genome/Ecoli_K12.fna Genome/Ecoli_K12\\n'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2321\u001b[0m             \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_expand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2322\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2323\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2324\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2325\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mnamed_script_magic\u001b[0;34m(line, cell)\u001b[0m\n\u001b[1;32m    140\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m                 \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscript\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 142\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshebang\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    143\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    144\u001b[0m         \u001b[0;31m# write a basic docstring:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-109>\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n\u001b[1;32m    243\u001b[0m             \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mraise_error\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m!=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 245\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mCalledProcessError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_run_script\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mto_close\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mCalledProcessError\u001b[0m: Command 'b'bowtie2-build Genome/Ecoli_K12.fna Genome/Ecoli_K12\\n'' returned non-zero exit status 127."
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "bowtie2-build --threads 8 Genome/Ecoli_K12.fna Genome/Ecoli_K12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now finally map the reads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir Mapping\n",
    "bowtie2 -p 8 -x Genome/Ecoli_K12 -1 FastQ/WT_1_R1.fastq -2 FastQ/WT_1_R2.fastq -S Mapping/WT_1.sam > Mapping/WT_1.log\n",
    "bowtie2 -p 8 -x Genome/Ecoli_K12 -1 FastQ/WT_2_R1.fastq -2 FastQ/WT_2_R2.fastq -S Mapping/WT_2.sam > Mapping/WT_2.log &\n",
    "bowtie2 -p 2 -x Genome/Ecoli_K12 -1 FastQ/MazF_1_R1.fastq -2 FastQ/MazF_1_R2.fastq -S Mapping/MazF_1.sam > Mapping/MazF_1.log &\n",
    "bowtie2 -p 2 -x Genome/Ecoli_K12 -1 FastQ/MazF_2_R1.fastq -2 FastQ/MazF_2_R2.fastq -S Mapping/MazF_2.sam > Mapping/MazF_2.log &\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should have now 4 .sam files corresponding to your reads mapped on the E. coli genome. For convenience we will convert this files to .bam file which are lighter version of theses files (indexed and binary)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samFile=Mapping/WT_1\n",
    "samtools view -b -q 1 $samFile.sam > ${samFile}_raw.bam\n",
    "samtools sort -o ${samFile}.bam ${samFile}_raw.bam\n",
    "samtools index ${samFile}.bam\n",
    "rm ${samFile}_raw.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samFile=Mapping/WT_2\n",
    "samtools view -b -q 1 $samFile.sam > ${samFile}_raw.bam\n",
    "samtools sort -o ${samFile}.bam ${samFile}_raw.bam\n",
    "samtools index ${samFile}.bam\n",
    "rm ${samFile}_raw.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samFile=Mapping/MazF_1\n",
    "samtools view -b -q 1 $samFile.sam > ${samFile}_raw.bam\n",
    "samtools sort -o ${samFile}.bam ${samFile}_raw.bam\n",
    "samtools index ${samFile}.bam\n",
    "rm ${samFile}_raw.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samFile=Mapping/MazF_2\n",
    "samtools view -b -q 1 $samFile.sam > ${samFile}_raw.bam\n",
    "samtools sort -o ${samFile}.bam ${samFile}_raw.bam\n",
    "samtools index ${samFile}.bam\n",
    "rm ${samFile}_raw.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reads counting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that all reads are mapped, we use .bam files to count reads in all genes of E. coli. The annotation file Ecoli_KL12.gtf will be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "UsageError: Cell magic `%%bashfeatureCounts` not found.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "mkdir GeneCount\n",
    "featureCounts -t \"gene\" -a Genome/Ecoli_K12.gtf -o GeneCount/WT_1.txt Mapping/WT_1.bam\n",
    "featureCounts -t \"gene\" -a Genome/Ecoli_K12.gtf -o GeneCount/WT_2.txt Mapping/WT_2.bam\n",
    "featureCounts -t \"gene\" -a Genome/Ecoli_K12.gtf -o GeneCount/MazF_1.txt Mapping/MazF_1.bam\n",
    "featureCounts -t \"gene\" -a Genome/Ecoli_K12.gtf -o GeneCount/MazF_2.txt Mapping/MazF_2.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before doing the differential expression analysis we will perform a multi QC analysis. We will thus control that everything went well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir MultiQC\n",
    "multiqc -f -i EcoliRNASeq -n EcoliRNASeq -o MultiQC ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You check your results in `MultiQC/EcoliRNASeq.html` file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential Expression analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now reaching the last part of our analysis. We will compare the two MazF replicates to the two WT replicates. We expect to have a list of genes differentially expressed and an html report showing the quality of our analysis."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "To run differential expression analysis we use a package wrapper named SARTools. Basically, it includes all useful scripts necessary for such an analysis. Everything is then run with: `RunRNASeqDESeq2.R` file. This file use an experimental design file `target_RNASeq` which describe the experiment (datasets, biological conditions, read count files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "Rscript RunRNASeqDESeq2.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A final, `RNASeq_Ecoli_report.html` webpage report has been produced along with tables summarizing the analysis (`MazFvsWT.complete.txt`, `MazFvsWT.up.txt`, `MazFvsWT.down.txt`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We found:\n",
    "1359 genes downregulated, and 1296 upregulated genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use cluster computing to speed up computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(This part only works when you have access to the cluster = not today)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I am creating a bash file containing all information for the cluster (see `WT_1_Mapping.sh`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then I run it on the cluster using `qsub` program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "bash: line 1: qsub: command not found\n"
     ]
    },
    {
     "ename": "CalledProcessError",
     "evalue": "Command 'b'qsub -l ncpus=12,mem=40G -N Mapping -e /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping_err.log -o /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping.log /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping.sh\\n'' returned non-zero exit status 127.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m-------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mCalledProcessError\u001b[0m                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-6b7af13f4d63>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'bash'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'qsub -l ncpus=12,mem=40G -N Mapping -e /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping_err.log -o /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping.log /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping.sh\\n'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2321\u001b[0m             \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_expand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2322\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2323\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2324\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2325\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mnamed_script_magic\u001b[0;34m(line, cell)\u001b[0m\n\u001b[1;32m    140\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m                 \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscript\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 142\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshebang\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    143\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    144\u001b[0m         \u001b[0;31m# write a basic docstring:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-109>\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n\u001b[1;32m    243\u001b[0m             \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mraise_error\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m!=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 245\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mCalledProcessError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_run_script\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mto_close\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mCalledProcessError\u001b[0m: Command 'b'qsub -l ncpus=12,mem=40G -N Mapping -e /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping_err.log -o /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping.log /home/becavin/omicsTD//temp/Mapping/WT_1_Mapping.sh\\n'' returned non-zero exit status 127."
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "qsub -l ncpus=12,mem=40G -N Mapping -e WT_1_Mapping_err.log -o WT_1_Mapping.log WT_1_Mapping.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run on all datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that I have controlled that everything run well, I can use the full power of the cluster and run as much `qsub` as datasets. To do that I have wrapped everything in `RunMapping.sh` which creates `WT_1_Mapping.sh` and the others, and then run the appropriate `qsub`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key parameters in `RunMapping.sh`is the number of CPUs and memory which you ask (NB_CPUS=24, MEM_JOB=40G)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "UsageError: %%bash is a cell magic, but the cell body is empty.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "sh RunMapping.sh WT_1\n",
    "sh RunMapping.sh WT_2\n",
    "sh RunMapping.sh MazF_1\n",
    "sh RunMapping.sh MazF_2"
   ]
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}