{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:green\">Formation South Green 2022</span> - Structural Variants Detection by using short and long reads \n",
    "\n",
    "# __DAY 1 : How to map reads against a reference genome ?__ \n",
    "\n",
    "Created by C. Tranchant (DIADE-IRD), J. Orjuela (DIADE-IRD), F. Sabot (DIADE-IRD) and A. Dereeper (PHIM-IRD)\n",
    "\n",
    "\n",
    "***\n",
    "\n",
    "# <span style=\"color: #006E7F\">Table of contents</span>\n",
    "<a class=\"anchor\" id=\"home\"></a>\n",
    "\n",
    "[I - Preparing data](#data)\n",
    "\n",
    "* [Download sequencing data (SR & LR) for Simulated clones](#download)\n",
    "\n",
    "[II - Mapping Practice](#mapping) \n",
    "  \n",
    "[2.1. Mapping short reads vs a reference with  `bwa mem`](#bwamem)\n",
    "\n",
    "   * [ Reference indexation](#refindex)\n",
    "   * [Run the mapping with `bwa mem`](#bwamem2-cmd)\n",
    "   * [Calculate stats from mapping `samtools flagstat`](#flagstats)\n",
    "   * [Convert sam into bam `samtools view`](#samtoolsview)\n",
    "   * [Generate a bam file that contains only the reads correctly paired mapped `samtools view`](#corrmap)\n",
    "   * [Indexing bam fil](#indexbam) \n",
    "   * [EXERCICE : MAP ALL SR WITH BWAMEM2](#mapallmem)\n",
    "\n",
    "[2.2. Mapping long reads vs a reference with  `minimap2`](#minimapé)\n",
    "   * [EXERCICE : MAP ALL LR WITH MINIMAP2](#mapallminimap)\n",
    "\n",
    "[III - Centralize final mapping data into a single bam directory](#reorder)\n",
    "\n",
    "</span>\n",
    "\n",
    "***\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:#006E7F\">__I - Preparing data__ <a class=\"anchor\" id=\"data\"></span>  \n",
    "\n",
    "### <span style=\"color: #4CACBC;\"> Download sequencing data we will use</span>  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "myFirstJupyterBook.ipynb  SV_DATA  training_SV_teaching\n"
     ]
    }
   ],
   "source": [
    "# go to work directory and download data\n",
    "cd /home/jovyan/work/\n",
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Download sequencing data (SR & LR) for Simulated clones <a class=\"anchor\" id=\"download\"> </span>  \n",
    "\n",
    "Before starting, please download special data created for this practical training. Data are available on the from I-Trop server.\n",
    "\n",
    "Each participant will analyse a Clone, results will be complete in this shared file\n",
    "\n",
    "To generate Clone data, a 1Mb contig was extracted from chromosome 1 of rice.\n",
    "\n",
    "20 levels of variation were generated and long reads were simulated for each.\n",
    "\n",
    "We have introduced different variations (SNP, indel, indel+translocations) and also some contaminations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2022-06-20 14:45:22--  https://itrop.ird.fr/sv-training/SV_DATA.tar.gz\n",
      "Resolving itrop.ird.fr (itrop.ird.fr)... 91.203.35.184\n",
      "Connecting to itrop.ird.fr (itrop.ird.fr)|91.203.35.184|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 4134064863 (3.8G) [application/x-gzip]\n",
      "Saving to: ‘SV_DATA.tar.gz’\n",
      "\n",
      "SV_DATA.tar.gz      100%[===================>]   3.85G  56.2MB/s    in 76s     \n",
      "\n",
      "2022-06-20 14:46:38 (52.0 MB/s) - ‘SV_DATA.tar.gz’ saved [4134064863/4134064863]\n",
      "\n",
      "FINISHED --2022-06-20 14:46:38--\n",
      "Total wall clock time: 1m 16s\n",
      "Downloaded: 1 files, 3.8G in 1m 16s (52.0 MB/s)\n",
      "SV_DATA/\n",
      "SV_DATA/REF/\n",
      "SV_DATA/REF/reference.fasta\n",
      "SV_DATA/LONG_READS/\n",
      "SV_DATA/LONG_READS/Clone10.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone11.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone12.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone13.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone14.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone15.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone16.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone17.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone18.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone19.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone1.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone20.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone2.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone3.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone4.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone5.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone6.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone7.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone8.fastq.gz\n",
      "SV_DATA/LONG_READS/Clone9.fastq.gz\n",
      "SV_DATA/SHORT_READS/\n",
      "SV_DATA/SHORT_READS/Clone10_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone10_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone11_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone11_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone12_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone12_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone13_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone13_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone14_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone14_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone15_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone15_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone16_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone16_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone17_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone17_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone18_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone18_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone19_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone19_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone1_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone1_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone20_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone20_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone2_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone2_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone3_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone3_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone4_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone4_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone5_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone5_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone6_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone6_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone7_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone7_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone8_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone8_R2.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone9_R1.fastq.gz\n",
      "SV_DATA/SHORT_READS/Clone9_R2.fastq.gz\n"
     ]
    }
   ],
   "source": [
    "# download available compressed DATA \n",
    "wget --no-check-certificat -rm -nH --cut-dirs=1 --reject=\"index.html*\" https://itrop.ird.fr/sv-training/SV_DATA.tar.gz\n",
    "# decompress data\n",
    "tar zxvf SV_DATA.tar.gz\n",
    "rm SV_DATA.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> List the content of the directory work and check that the directory SV_DATA have been created</span>  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# check data \n",
    "ls -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> List the content of the directory SV_DATA</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ls -lRt SV_DATA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# <span style=\"color:#006E7F\">__II -  MAPPING PRACTICE__ <a class=\"anchor\" id=\"mapping\"></span>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read congruency is an important measure in determining assembly accuracy.\n",
    "\n",
    "Clusters of read pairs or single long reads that align incorrectly are strong indicators of mis-assembly.\n",
    "\n",
    "Reads mapping is usually the firt step before SNP or Variant calling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Make a folder for your results</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p ~/work/MAPPING-ILL\n",
    "cd ~/work/MAPPING-ILL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Declare important variables</span>  \n",
    "\n",
    "We are going to set up bash variables with the path to our data.We set a bash variable like this : `var=\"value\"`\n",
    "and call it as: `echo $var`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# REFERENCE \n",
    "REF_DIR=\"/home/jovyan/work/SV_DATA/REF/\"\n",
    "REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\"\n",
    "\n",
    "# ONT DATA\n",
    "ONT=\"/home/jovyan/work/SV_DATA/LONG_READS/Clone${i}.fastq.gz\"\n",
    "\n",
    "# ILLUMINA DATA\n",
    "ILL_R1=\"/home/jovyan/work/SV_DATA/SHORT_READS/Clone${i}_R1.fastq.gz\"\n",
    "ILL_R2=\"/home/jovyan/work/SV_DATA/SHORT_READS/Clone${i}_R2.fastq.gz\"\n",
    "\n",
    "#CLONE NUMBER THAT YOU ARE GOING TO ANALYZE \n",
    "i=10 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Print the variable i, REF, ILL_R1 & ILL_R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clone10 /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
      "/home/jovyan/work/SV_DATA/SHORT_READS/Clone10_R1.fastq.gz /home/jovyan/work/SV_DATA/SHORT_READS/Clone10_R2.fastq.gz\n"
     ]
    }
   ],
   "source": [
    "echo \"Clone${i} $REF\" \n",
    "echo $ILL_R1 $ILL_R2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "-------------\n",
    "# <span style=\"color: #4CACBC;\"> 2.1. Mapping short reads vs a reference with  `bwa mem` <a class=\"anchor\" id=\"bwamem\"></span>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this practice, we are going to map short reads against a reference. To know, how well do the reads align back to the reference, we use bwa-mem2 and samtools to assess the basic alignment statistics.\n",
    "\n",
    "In this exercise, we will use reference.fasta assembly as well ILLUMINA READS from your favorite CLONE.\n",
    "\n",
    "The tool bwa needs 2 steps: \n",
    "- **Reference indexing**: `bwa index reference`\n",
    "- **Mapping in itself**: `bwa mem  -R READGROUP [options] reference fastq1 fastq2 > out.sam`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Reference indexation  <a class=\"anchor\" id=\"refindex\"></span>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before mapping we need index reference file! Check bwa-mem2 index command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd $REF_DIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Indexing reference /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
      "\n",
      "[bwa_index] Pack FASTA... 0.01 sec\n",
      "init ticks = 162934945\n",
      "ref seq len = 2040002\n",
      "binary seq ticks = 61204851\n",
      "build index ticks = 604754568\n",
      "ref_seq_len = 2040002\n",
      "count = 0, 576483, 1020001, 1463519, 2040002\n",
      "BWT[1932441] = 4\n",
      "CP_SHIFT = 5, CP_MASK = 31\n",
      "sizeof CP_OCC = 64\n",
      "max_occ_ind = 63750\n",
      "ref_seq_len = 2040002\n",
      "count = 0, 576483, 1020001, 1463519, 2040002\n",
      "BWT[1932441] = 4\n",
      "CP_SHIFT = 6, CP_MASK = 63\n",
      "sizeof CP_OCC = 64\n",
      "max_occ_ind = 31875\n"
     ]
    }
   ],
   "source": [
    "echo -e \"\\nIndexing reference $REF\\n\"\n",
    "bwa-mem2 index $REF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\">Check that the indexes have been created </span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reference.fasta       reference.fasta.ann          reference.fasta.pac\n",
      "reference.fasta.0123  reference.fasta.bwt.2bit.64\n",
      "reference.fasta.amb   reference.fasta.bwt.8bit.32\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> => Let's map now but only WITH READS FROM ONLY ONE CLONE </span>  \n",
    "\n",
    "* Go into the directory MAPPING-ILL\n",
    "* Create a subdirectory to save the files generated by the mapping step. \n",
    "Eg: If you are going to analyze the `clone1`, create the subdirectory `dirClone1`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>>>>>>>>> Creation directory for Clone10\n",
      "\n"
     ]
    }
   ],
   "source": [
    "cd ~/work/MAPPING-ILL\n",
    "echo -e \"\\n>>>>>>>>>> Creation directory for Clone$i\\n\"\n",
    "mkdir -p dirClone$i\n",
    "cd dirClone$i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Run the mapping with `bwa mem` <a class=\"anchor\" id=\"bwamem2-cmd\"></span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>>>>>>>>> Mapping Clone10\n",
      "\n",
      "-----------------------------\n",
      "Executing in AVX2 mode!!\n",
      "-----------------------------\n",
      "Ref file: /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
      "Entering FMI_search\n",
      "reference seq len = 2040003\n",
      "count\n",
      "0,\t1\n",
      "1,\t576484\n",
      "2,\t1020002\n",
      "3,\t1463520\n",
      "4,\t2040003\n",
      "\n",
      "Reading other elements of the index from files /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
      "prefix: /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
      "[M::bwa_idx_load_ele] read 0 ALT contigs\n",
      "Done reading Index!!\n",
      "Reading reference genome..\n",
      "Binary seq file = /home/jovyan/work/SV_DATA/REF/reference.fasta.0123\n",
      "Reference genome size: 2040002 bp\n",
      "Done readng reference genome !!\n",
      "\n",
      "[0000] 1: Calling process()\n",
      "\n",
      "Threads used (compute): 8\n",
      "Info: projected #read in a task: 529811\n",
      "------------------------------------------\n",
      "Memory pre-allocation for chaining: 1114.7223 MB\n",
      "Memory pre-allocation for BSW: 1916.9362 MB\n",
      "Memory pre-allocation for BWT: 618.5134 MB\n",
      "------------------------------------------\n",
      "No. of pipeline threads: 2\n",
      "[0000] read_chunk: 80000000, work_chunk_size: 44364000, nseq: 295760\n",
      "\t[0000][ M::kt_pipeline] read 295760 sequences (44364000 bp)...\n",
      "[0000] 2. Calling mem_process_seqs.., task: 0\n",
      "[0000] 3. Calling kt_for - worker_bwt\n",
      "[0000] read_chunk: 80000000, work_chunk_size: 0, nseq: 0\n",
      "[0000] 3. Calling kt_for - worker_aln\n",
      "Inferring insert size distribution from data, l_pac: 1020001, n: 295760\n",
      "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 90449, 1, 0)\n",
      "[M::mem_pestat] skip orientation FF as there are not enough pairs\n",
      "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n",
      "[M::mem_pestat] (25, 50, 75) percentile: (380, 398, 414)\n",
      "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (312, 482)\n",
      "[M::mem_pestat] mean and std.dev: (397.30, 25.31)\n",
      "[M::mem_pestat] low and high boundaries for proper pairs: (278, 516)\n",
      "[M::mem_pestat] skip orientation RF as there are not enough pairs\n",
      "[M::mem_pestat] skip orientation RR as there are not enough pairs\n",
      "[0000] 10. Calling kt_for - worker_sam\n",
      "\t[0000][ M::mem_process_seqs] Processed 295760 reads in 41.543 CPU sec, 5.329 real sec\n",
      "[0000] read_chunk: 80000000, work_chunk_size: 0, nseq: 0\n",
      "[0000] Computation ends..\n",
      "No. of OMP threads: 8\n",
      "Processor is runnig @2295.140301 MHz\n",
      "Runtime profile:\n",
      "\n",
      "\t Time taken for main_mem function: 6.68 Sec\n",
      "\n",
      "\tIO times (sec) :\n",
      "\tReading IO time (reads) avg: 0.97, (0.97, 0.97)\n",
      "\tWriting IO time (SAM) avg: 0.34, (0.34, 0.34)\n",
      "\tReading IO time (Reference Genome) avg: 0.00, (0.00, 0.00)\n",
      "\tIndex read time avg: 0.02, (0.02, 0.02)\n",
      "\n",
      "\tOverall time (sec) (Excluding Index reading time):\n",
      "\tPROCESS() (Total compute time + (read + SAM) IO time) : 6.65\n",
      "\tMEM_PROCESS_SEQ() (Total compute time (Kernel + SAM)), avg: 5.33, (5.33, 5.33)\n",
      "\n",
      "\t SAM Processing time (sec):\n",
      "\t--WORKER_SAM avg: 1.80, (1.80, 1.80)\n",
      "\n",
      "\tKernels' compute time (sec):\n",
      "\tTotal kernel (smem+sal+bsw) time avg: 3.50, (3.50, 3.50)\n",
      "\t\tSMEM compute avg: 1.42, (1.46, 1.38)\n",
      "\t\tSAL compute avg: 0.11, (0.16, 0.08)\n",
      "\t\tBSW time, avg: 1.91, (1.92, 1.89)\n",
      "\n",
      "\tTotal allocs: 294308 = 294308 out total requests: 343778, Rate: 0.86\n",
      "\n",
      "Important parameter settings: \n",
      "\tBATCH_SIZE: 512\n",
      "\tMAX_SEQ_LEN_REF: 256\n",
      "\tMAX_SEQ_LEN_QER: 128\n",
      "\tMAX_SEQ_LEN8: 128\n",
      "\tSEEDS_PER_READ: 500\n",
      "\tSIMD_WIDTH8 X: 32\n",
      "\tSIMD_WIDTH16 X: 16\n",
      "\tAVG_SEEDS_PER_READ: 64\n"
     ]
    }
   ],
   "source": [
    "echo -e \"\\n>>>>>>>>>> Mapping Clone$i\\n\"\n",
    "bwa-mem2 mem -M -t 8 $REF $ILL_R1 $ILL_R2 > Clone$i.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### <span style=\"color: #4CACBC;\">Check that the file `.sam` have been created by `bwa mem` </span>  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 118940\n",
      "-rw-r--r-- 1 jovyan users 121793468 Jun 20 14:56 Clone10.sam\n"
     ]
    }
   ],
   "source": [
    "ls -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\">Display the first and the end of the sam file just created </span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "@SQ\tSN:Reference\tLN:1020001\n",
      "@PG\tID:bwa\tPN:bwa\tVN:2.0pre2\tCL:bwa-mem2 mem -M -t 8 /home/jovyan/work/SV_DATA/REF/reference.fasta /home/jovyan/work/SV_DATA/SHORT_READS/Clone10_R1.fastq.gz /home/jovyan/work/SV_DATA/SHORT_READS/Clone10_R2.fastq.gz\n",
      "Reference-Clone10295760\t77\t*\t0\t0\t*\t*\t0\t0\tGTATAAGTACCCGGTCGAATCAAAGGTAACGTTAAATAGGTACTCCGCCAGGGCAGATTTCAACAGCCAAACTGCCCCCCAGGGGTATCTTACAGGCAATGGCTTAGAAGCGTTCCTAAGTGGACGACTCTCTGGAAACTCGCCAATGAG\tCC=G=GGGGGGG=IIIIGGIIGICGIIIICIGIIIIIGGCIIGIGIG=CIGIIIICIGICGGGGGCCCGGCGCGCGGGGG8CGG8CGGGGGCCCGGGGGGIGCGGGG=GCGGCCGCGG55GGCGCG8GGGGCGG=CCGCGG5GGCGCGC=\tAS:i:0\tXS:i:0\n",
      "Reference-Clone10295760\t141\t*\t0\t0\t*\t*\t0\t0\tGCACCCAAGGTGATCAACCCGGCGCTGCATGAGTATGCAACATGTTCGGCAGATGCCGTCAGTTTGGCATGCGTAATTCAATGTCGCAAGGAGGATATCCCGCTGGGATTACATTCGCGTATAGTTTATGGGCCTTCATTCGTTTTTACG\tCC=GGGGGGGGGGIGIIIIIICIGICI5IGCGIGGGICIGIIIICCIIG=IICGC=G==GGCGGGIGI=CCGGCGGCIGG=CGCG5CG=CGCGGG5GCGGGCCCII=GGGGGCG==GGGGGGGCCGCGGCCCGGGGGGGCCCCGGCGCGC\tAS:i:0\tXS:i:0\n",
      "Reference-Clone10295758\t99\tReference\t124517\t60\t150M\t=\t124764\t397\tTATATAATTTACTAATAATTAAAAATAGATTTTAAAAATTATTAATAACTAATACGTATATACGATATAATATGAATATTACACGTTAGTACTTTTGTAAAGTTATTGCAAACTTTAAAATTATGTTGTCATTTTAATTTAGTTGAATAA\tCCCGGGGGGGGGGIIGIIIICIIICGIIIIIIIIIIICIIIICGCIGIGICGIGGIIGGGICCGIICIIGGGCGGGC==GGCCGGGGGCGGC=GGGGGGGICCG=GGCGG5GGCGCGGGCGCCGCCG88GGGC8GCGCGGGGGCGCGGCC\tNM:i:13\tMD:Z:12A11T1T3A8A12C9A0A20A6G20A25A2T8\tAS:i:85\tXS:i:0\n",
      "Reference-Clone10295758\t147\tReference\t124764\t60\t150M\t=\t124517\t-397\tTTCCAATTTCTTCACGAAACTCTCAAAGAAAAATCCTGAATTGGATCCTAAATTTTCAAGAAAAAGTACGAATTACCACCCCCGCGAACTATTACGGTAGTATGATTTATCCCCCTAAAATATAATACCGGACCTTTAGCCTGCCAAACT\tGCGCGCGGGC=CGCGGGGC5CCCCCGGGGGGGGGGGGCGCGGGIIIICC5G=GCGGCGGC5CGC=CGCGGG8CGG8CGGGCGGC==IIGGCI=GCIICGIG=CGGGGIICC5IGGCIGIIIIIGI5G=IIICIIIIIGGGGGGGGGGCC8\tNM:i:9\tMD:Z:13C7T2C23G6C63C13A6A1C7\tAS:i:105\tXS:i:0\n",
      "Reference-Clone10295756\t83\tReference\t660900\t25\t150M\t=\t660609\t-441\tAGTAATTTGACAAAACATTGACATCTGCAGGCTTTGCAGTCAATTAGAAAGATAATTGTGTGTTCTATCGCCATGGTGGGGGTGAGGGAGTTAATTTGTGCTTATATGTTGACGACATACTGATATTTGGGACAACCCTTGAGGTAATAA\tGCGGCGCGGGGGCG5GG=CGCGGGGGCGGGGGGGGG8GGGGGGGCGG=GIGCCGGGGGCGGGCCGCG=GCCCCGGCICCGCGGGCIGGCCGGGIGIGGIIIGGCIICGIIIIGIIIIGGIIIIIIIIIIIIIIICIIGGGGGGGCGCCCC\tNM:i:10\tMD:Z:2A41G2G0C6A7A29T9G31A9G4\tAS:i:102\tXS:i:102\tXA:Z:Reference,-520388,150M,10;Reference,-488509,150M,11;\n",
      "Reference-Clone10295756\t163\tReference\t660609\t25\t150M\t=\t660900\t441\tAGGCTGGGCTTGTTGCTAAGGGCTATACTTAGAAAGAAGGCGAAGATTTCTTTGACACATACTCACCTGTTGGTAGATGGACCACAATTCGTGAGCTACTTTCCCTAGCAGCCTCACATGGTCTGCTCGTTCATCAAATGGACGTTAAGA\tCCCG5GCGGG=CGCICGIGICCIIGIIC=IIGIGIIIIIGIIG5IICIGGIGGGCIIIIIGIGC5IGGGCGGGIGGIGCGCG=GCCGC5GCGCGGCGGGCGCIICC8CGCGGCGGG=CG=GGGGCGGGCGG=5GGCGGGG5CGCCCGCC=\tNM:i:7\tMD:Z:5C23C28T13C5T14T30T25\tAS:i:115\tXS:i:115\tXA:Z:Reference,+488218,150M,7;Reference,+520097,150M,8;\n",
      "Reference-Clone10295754\t83\tReference\t399106\t0\t150M\t=\t398838\t-418\tACATCATCTCAATCTACACCTGCAATTGCGTCACCCTGATCAAAGTAGATCATAAAATTCTTAACCTTGTGCACAACTGATTTCATCACTTCAAAGTGAGCGATAGTAGATATAACCCTAAGGCCATCTTCAAGAGTCCTTCCAGGGAGC\tGCGGGGGGGGCCGCGGGCGGGCGGGCCGGCGGCGGG8CCGGCG8GGCGCICCGGC5CG=GGG=GGGCGCGGGG=GGIIGGGGGGCGICIGI=GGIGGCIIGICICICIGIICIII=IIGIIIGIIIIIIIIIIIIIIGGGGGGGGGGCCC\tNM:i:6\tMD:Z:8C19T21A46C3C0C47\tAS:i:120\tXS:i:120\tXA:Z:Reference,-332803,150M,6;Reference,-365956,150M,6;\n",
      "Reference-Clone10295754\t163\tReference\t398838\t0\t150M\t=\t399106\t418\tTTCTTTCCTCTTCAACTGACGGCCCTGAGCTTTCTTCTAAGACTTTGCCTTGCAGTCAGATAAGTCTTCAACACCCATATCAACATTGTTAACCTACAAGACATCAACTGCATATGATATATCATAGTCACGGTCCACAAAATCTGTGTC\tCCCGGGGGGGGGGGIGIIGIGIIIIGIIIIIIIGIGGIIIIIIIIICGGGGIG5GGC=8I5GGGIGG=5CGGGGGGCG=GGGCCGGGGGGGG8=CGGG5GCG5ICCIGGG5G=CGGGCCGG5CG=GGGCG=GG=CGGGCCCC8GCGCCCC\tNM:i:9\tMD:Z:20T5T10C51C3A6T5T6C17T18\tAS:i:105\tXS:i:105\tXA:Z:Reference,+365688,150M,9;Reference,+332535,150M,9;\n"
     ]
    }
   ],
   "source": [
    "head Clone10.sam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reference-Clone1010\t353\tReference\t571519\t60\t83H67M\t=\t571653\t284\tATTACCTAATGCATACATAGTTCTACAAACATCTTAGTTCAGATCAGATGCATCATCACATTGTTAC\tGGCCGGGGGCCGGGG=GIC5CGC5GGGG5GCGGG8GGGGGGGGGG5GCGGG5GG5GC5CCGG=GGCC\tNM:i:0\tMD:Z:67\tAS:i:67\tXS:i:0\tSA:Z:Reference,569030,+,82M68S,60,3;\n",
      "Reference-Clone1010\t145\tReference\t571653\t60\t150M\t=\t569030\t-2773\tTCAGAAGCAGATCAACAACTGGTTCATCAACCAGAGGAAACGGCACTGGAAGCCATCGGAGGACATGCCGTTCGTCATGATGGAAGGTTTTCACCCACAGAATGCTGCTGCATTGTACATGGATGGCCCGTTCATGGCAGATGGAATGTA\tCGGCCGCG8CGCGGCGG8CCCCGGGCGCGGCGG=GGGCGGCGG=ICI=GGGGGCGGCCC8G==GCG=GGGGCCGGGGGCCGG5CICCCG=IIGGIGICCIGCIIG=GCGIGGIIIIIGGIGICGIIGGI8IIIIIGIGGGGGGG=GGCCC\tNM:i:0\tMD:Z:150\tAS:i:150\tXS:i:0\n",
      "Reference-Clone108\t99\tReference\t263937\t60\t150M\t=\t264214\t427\tTTTAGTTGATGAACACAAATAATAATTGATTAAAGGGAACTTTCCATTCGGTCGTTTCCTGTCTCCTTCTTTGGGTACTACTATCATTTTCTTTTTCTGAAATTCCTTTTGCTGTATATCATTTCAGCATGCAATACTTAATCTGACAAA\tCCCGGGGGGCGGGIIGCICIGIIGCGCCIIGIIIGI5IGIIIGIGI=GGCGGGIIGIIGIGIGIGIGGGGG5CGGGCIGG=IGGGGCGGGGGGCGC5=GGIGCGGCCGC=CG8C5GGCCGGGCG=GGC=CGGCGGGCGGGGGG8GGGGC8\tNM:i:8\tMD:Z:9C27C14T51A7A9G11C9A5\tAS:i:110\tXS:i:0\n",
      "Reference-Clone108\t147\tReference\t264214\t60\t150M\t=\t263937\t-427\tAAAAAAATTTTGGGTGTCACATCGGACGTTTGACCAGATGTCGAAAGGGGGTTGGGGACTCGATTGAAAAAACTAATTCCATAACTCGCGTGGAAACCAGTAGACGAACTTTTTGAGCCTAATTAATCCATCATTAGCAAATGTGGGTTA\tCCGGC=GGGGGGGGCGGCCGGCCGGCGGCGC88G5CGG=G8CGI=5CCCGGCGG=G8GGCGG8GGCC=GGGCCGGGGCI=5GCGIGCIIGC==IGIGGGIGIIICCIIGGIIIIGCIIIIIGIGIIGIIIIIII5IIGGGGGGG=GGC=C\tNM:i:8\tMD:Z:50T2T5A3A25C9C8T0A40\tAS:i:110\tXS:i:53\n",
      "Reference-Clone106\t99\tReference\t154263\t60\t150M\t=\t154504\t391\tCCATGTCTCTGCTACACGACACAAGAGAGCATCGTCTTCTTCATCGTATTGATATTCATGTGGTCATGCCAATATGGGTTGCGCGGGGCCGACTCTGTTGCAGGTGAGGGCCGTCTTCGCGCCACACGCCTCGGCGGAGGAGCTGCCGCC\tCCCGCGGGGGGGGGIIIIIIIGIIGIIIIICGCGGGI5IIIICGCIIIGGG=GGIGCGGGGGIGIIGG5CGIIGC=G58GCGGG5GCCCGGG=GGGGGGGI=CGGGGGGGGGCGGCGGGGG=GGGCGGGCCCCGGGCGG5=GGCCGGGCG\tNM:i:11\tMD:Z:28A16A4C2G18C9T14C9C15G9C14T1\tAS:i:98\tXS:i:21\n",
      "Reference-Clone106\t147\tReference\t154504\t60\t150M\t=\t154263\t-391\tCTGCCGAACGCCGTCGATCGCCTGTGCCAGGAGTCTCAGGCGCATGAGGTGTTCCGAAGATCCGGCCATGGCGTCGACTACCTGTGCGTGGAGAAGATGGGCACGGAGATGGCAAGGTGCAGTCCGACTCAATGACGGGTGTAAAGGTAG\t8=GGGGG=CCGGG=CGGGGCCG=GGCGCCGC=CGCGGCGCGCGICCIIGGCGGGGCCCGCCGGGC=GCGCGGCCGGGGGCGGGCGGGIIGGGCGI=GIGI=GGGIIGG5IIGICIIIIIIGIIC=8IIIIGIIIGIIGGG=GGGGC=5CC\tNM:i:9\tMD:Z:29C3G2G0C16T1G19C45G11G15\tAS:i:105\tXS:i:64\n",
      "Reference-Clone104\t77\t*\t0\t0\t*\t*\t0\t0\tTTGTCCAGACACGTAGAATTGCGTTTGTCAACCCATCTATTTGCTGAACTTGCCAATGCCGCCCTGTGCGCAGCTACCTACTACAATGCAGTTAGCCCGAAAGGGACAGAAGGAGGCTTAGTCGTAGTGGAGACAAGTGTGGCGGTCATG\tCCCCGGGGGGGGGIIGCIIIIGIIGI5GIIIIIIIII=IGIGGIGGIGIGCGI5GIGGIGIIIGGGIGGGIG=CICGGGC5GGGG8GGGCGCGGGGGGCCCCGCGGGCGCGG5C5GCG=GGGGGGCCCGGG=GCCCGGGCCGGCGGG8G8\tAS:i:0\tXS:i:0\n",
      "Reference-Clone104\t141\t*\t0\t0\t*\t*\t0\t0\tTCCGACCCATACCTGATTCTATAATCTAGACAGGGACATTAATTTGTGCCGTAAAAGCGAGCCGACTTTGGGCGGAGCCCCGTATTACAACCAGAAGCGCAGTTGTGGAGGCGGCATCACATAGCGCAGAAACGGTGTGTTTGGCCGCAC\tCCCG5GGGGGGGGIIIIIIGIIIG8IIIIIIIGIIGIGIIIG8CIIIG5GIGGIIGGGGGCGIGGGG8G=GGCGG=CCGCGCGC5GGG=CGCGGGGCC5GCCCICIIGGCGGGGCGCGCCC=GGGGGGGGG=CGGGGGCGGGCGGGCGGC\tAS:i:0\tXS:i:0\n",
      "Reference-Clone102\t83\tReference\t594109\t60\t150M\t=\t593854\t-405\tAAGCGGTGAAAGCAGGAAAGTAGAGGGAAGAGAAGAAGACAAGAAAAAAGAACAGAAAGTAAAGGAAGAAAAAAAATCATAATACCAGTAGTGAAAACTCAAGCCCCTTCACTACTATGGTAGTTCGGCTTTGTCACAAAATTTAAAATT\tGGGGGCC=GGGCGGG5GGG=CCG=5GGGGCGG=GCGGGGGC5CGGCGGC=C=CGCCGGGGCGCGGG=GCGCGC=88GGGGGGGC=GGGGGGGGGGGGIGIGCGGGIIGICIIGCIIGGGIGCGIIICIIGGIIIIGIGGCGGGGGCGCC=\tNM:i:6\tMD:Z:12G7A18G20G55C6A26\tAS:i:120\tXS:i:0\n",
      "Reference-Clone102\t163\tReference\t593854\t60\t150M\t=\t594109\t405\tTGAACTGTATTTTCTCACTTATTGGCCCATCGTTTTACCTACAATTAATTATAAGCCGTTATACTTATTACAGATAAAAAATTATAACACTTTTATTCATAATTAACTACTTAATAGCTTTCTGGAATATAAACTACAATGATACACCAT\tCCCGGGGGGGGGGICIGIIIGI8GIGIGGIGIGGIII5IIIGGGIGGGCIIGGII=CGGGCIIG5GGGGCGCG5GGG5GGGG8CCGCGGGGGGGCGGGCCG=II=C8GCGCCGCGG8GCGGCCCGGGGCCG=CGGGGGG5GCCGCGGCCG\tNM:i:6\tMD:Z:71T15A26A5G6A18A3\tAS:i:121\tXS:i:0\n"
     ]
    }
   ],
   "source": [
    "tail Clone10.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Convert sam into bam `samtools view` <a class=\"anchor\" id=\"samtoolsview\"></span>  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[samopen] SAM header is present: 1 sequences.\n"
     ]
    }
   ],
   "source": [
    "samtools view -@4 -bh -S -o  Clone$i.bam Clone$i.sam "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check that the bam file have been created \n",
    "\n",
    "* Have a look at the filesize of the sam and bam files.\n",
    "* Remove the sam file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 155M\n",
      "-rw-r--r-- 1 jovyan users  39M Jun 20 14:56 Clone10.bam\n",
      "-rw-r--r-- 1 jovyan users 117M Jun 20 14:56 Clone10.sam\n"
     ]
    }
   ],
   "source": [
    "ls -lh\n",
    "rm Clone$i.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Calculate stats from mapping `samtools flagstat`<a class=\"anchor\" id=\"flagstats\"></span>   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools flagstat Clone$i.bam >Clone$i.flagstat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Display the content of the flagstat file</span>  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "296107 + 0 in total (QC-passed reads + QC-failed reads)\n",
      "0 + 0 duplicates\n",
      "221301 + 0 mapped (74.74%:-nan%)\n",
      "296107 + 0 paired in sequencing\n",
      "148037 + 0 read1\n",
      "148070 + 0 read2\n",
      "218229 + 0 properly paired (73.70%:-nan%)\n",
      "219681 + 0 with itself and mate mapped\n",
      "1620 + 0 singletons (0.55%:-nan%)\n",
      "0 + 0 with mate mapped to a different chr\n",
      "0 + 0 with mate mapped to a different chr (mapQ>=5)\n"
     ]
    }
   ],
   "source": [
    "cat Clone10.flagstat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Generate a bam file that contains only the reads correctly paired mapped `samtools view`<a class=\"anchor\" id=\"corrmap\"></span>   \n",
    "\n",
    "https://broadinstitute.github.io/picard/explain-flags.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools view -bh -@4 -f 0x02 -o Clone$i.mappedpaired.bam Clone$i.bam "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <span style=\"color: #4CACBC;\"> Sorting final bam </span>  \n",
    "\n",
    "* Generate the bam file ordered\n",
    "* Check that the new bam file have been created\n",
    "* Remove the bam file previously created (Clone$i.mappedpaired.bam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools sort -@4 Clone$i.mappedpaired.bam Clone$i.SORTED \n",
    "rm Clone$i.mappedpaired.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <span style=\"color: #4CACBC;\"> Indexing bam file<a class=\"anchor\" id=\"indexbam\"></span>   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools index Clone$i.SORTED.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 56880\n",
      "-rw-r--r-- 1 jovyan users 40018974 Jun 20 14:56 Clone10.bam\n",
      "-rw-r--r-- 1 jovyan users      381 Jun 20 14:56 Clone10.flagstat\n",
      "-rw-r--r-- 1 jovyan users 18213243 Jun 20 14:56 Clone10.SORTED.bam\n",
      "-rw-r--r-- 1 jovyan users     2856 Jun 20 14:56 Clone10.SORTED.bam.bai\n"
     ]
    }
   ],
   "source": [
    "ls -lrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <span style=\"color: #4CACBC;\"> => Let's map with data from all clones using a loop for mapping, with a single folder per sample<a class=\"anchor\" id=\"mapallminimap\"></span>   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in {1..20}\n",
    "    do\n",
    "        cd ~/work/MAPPING-ILL\n",
    "        echo -e \"\\n\\n>>>>>>>>>> Creation directory for Clone$i\"\n",
    "        mkdir -p dirClone$i\n",
    "        cd dirClone$i\n",
    "        \n",
    "        echo -e \"\\n>>>> Declare variables$i\"\n",
    "        REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\"\n",
    "        ILL_R1=\"/home/jovyan/work/SV_DATA/SHORT_READS/Clone${i}_R1.fastq.gz\"\n",
    "        ILL_R2=\"/home/jovyan/work/SV_DATA/SHORT_READS/Clone${i}_R2.fastq.gz\"\n",
    "\n",
    "        echo -e \"\\n>>>> Mapping Clone$i\\n\"\n",
    "        bwa-mem2 mem -M -t 8 $REF $ILL_R1 $ILL_R2 > Clone$i.sam\n",
    "        \n",
    "        echo -e \"\\n>>>> convert sam to bam for Clone$i\"\n",
    "        samtools view -@4 -bh -S -o  Clone$i.bam Clone$i.sam \n",
    "        rm Clone$i.sam\n",
    "        echo -e \"\\n>>>> Flagstats from all reads $i\"\n",
    "        samtools flagstat Clone$i.bam >Clone$i.flagstat\n",
    "        \n",
    "        echo -e \"\\n>>>> Extract only correctly mapped and calculate flagstats $i\"\n",
    "        samtools view -bh -@4 -f 0x02 -o Clone$i.mappedpaired.bam Clone$i.bam \n",
    "        \n",
    "        echo -e \"\\n>>>> Sort mappedpaired bam file $i\"\n",
    "        samtools sort -@4 Clone$i.mappedpaired.bam Clone$i.SORTED \n",
    "        rm Clone$i.mappedpaired.bam\n",
    "    done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color: #4CACBC;\"> 2.2 Mapping Long reads vs a Reference `minimap2` <a class=\"minimap2\" id=\"\"></span> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar process such as SR is done in LR. In this case mapper is minimap2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declare variables\n",
    "i=10\n",
    "REF_DIR=\"/home/jovyan/work/SV_DATA/REF/\"\n",
    "REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\"\n",
    "ONT=\"/home/jovyan/work/SV_DATA/LONG_READS/Clone${i}.fastq.gz\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Mapping with `minimap2`</span> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <span style=\"color: #4CACBC;\"> => Let's map now but only WITH READS FROM ONLY ONE CLONE</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Creation directory for Clone10\n",
      "\n",
      "Clone10\n"
     ]
    }
   ],
   "source": [
    "mkdir -p ~/work/MAPPING-ONT\n",
    "cd ~/work/MAPPING-ONT\n",
    "echo -e \"\\nCreation directory for Clone$i\\n\"\n",
    "echo Clone$i\n",
    "mkdir -p dirClone$i\n",
    "cd dirClone$i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Mapping Clone10 minimap2 \n",
      "\n",
      "[M::mm_idx_gen::0.084*2.42] collected minimizers\n",
      "[M::mm_idx_gen::0.113*4.01] sorted minimizers\n",
      "[M::main::0.113*4.01] loaded/built the index for 1 target sequence(s)\n",
      "[M::mm_mapopt_update::0.121*3.80] mid_occ = 10\n",
      "[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 1\n",
      "[M::mm_idx_stat::0.126*3.69] distinct minimizers: 165344 (91.75% are singletons); average occurrences: 1.156; average spacing: 5.336\n",
      "[M::worker_pipeline::9.120*7.21] mapped 11235 sequences\n",
      "[M::main] Version: 2.17-r941\n",
      "[M::main] CMD: minimap2 -ax map-ont -t 12 /home/jovyan/work/SV_DATA/REF/reference.fasta /home/jovyan/work/SV_DATA/LONG_READS/Clone10.fastq.gz\n",
      "[M::main] Real time: 9.135 sec; CPU: 65.813 sec; Peak RSS: 0.583 GB\n"
     ]
    }
   ],
   "source": [
    "echo -e \"\\nMapping Clone$i minimap2 \\n\"\n",
    "minimap2 -ax map-ont -t 12 ${REF} ${ONT} > Clone${i}_ONT.sam "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "@SQ\tSN:Reference\tLN:1020001\n",
      "@PG\tID:minimap2\tPN:minimap2\tVN:2.17-r941\tCL:minimap2 -ax map-ont -t 12 /home/jovyan/work/SV_DATA/REF/reference.fasta /home/jovyan/work/SV_DATA/LONG_READS/Clone10.fastq.gz\n",
      "86f59255-6632-405c-a329-62d9dba8a95f\t0\tReference\t495672\t60\t3157S38M1D8M1D49M1I3M1D8M1I37M2I5M2D52M1D10M2I8M1D11M1I18M1I17M1I62M1D56M1D27M2D32M3D16M1D3M3I41M1I5M1D3M4I6M1I3M1I2M1I6M1D164M4I13M1I29M1D2M1I22M1I3M2I6M2D14M1D15M1D47M1I61M1D32M3D13M2I17M1I60M1D22M1I62M1D17M2I15M1D8M1I2M1I21M1D22M1D2M1D37M1D5M2I3M1D5M3D3M2D27M1I8M1D17M1I8M1I4M1I5M1I13M2I42M1I2M1I24M1I6M1I1M1I39M1D38M3D11M1I9M2D12M1I1M1I44M1I36M1D11M1D48M2I4M2I5M2I103M1I2M2D20M2D45M1I15M1I37M1I10M1I5M1I42M2I23M1D17M7S\t*\t0\t0\tTATTTTCAAATACTAAATGATTTCAACTGAAAACGTCATCAATAACTCAAAGTTGTATTAATCATCAAGATCTATAACTTTCATTTTGGTCAGTTCTTCATCGGACAAAGTAATTTGTAAGATTGTTCCACAAGATGTACTTATCTTTTATATAGTTAATAAAAACTATAAGAGTGGTTACATTTTGTGAACAGTCTTATTAATAACTTTGTCGGATGAAGAAATGTCTAAATGGGTTATAGATCTTGCTGAGTTATACAACTACGTTCTTGATGACTTTTCAGCCGAAATCAAATACTACTGCAAAATATTGTTTGAAGTTTTGAAATTCAGAATTTTAATTGATAAAAGTCACAATAAAAAATGGCCAAAATAATAGCAGTAAGAACACAATAACATGATAGAGCATGATTTTAGAAATATTTAGAAAAAGTCATCTTGTTTAGGTTTGAATAATTTGAGTTTGAATTTCAAATTCTGACAATTTCAAACAACATTTTCAAAAACTAAACGATTTTAAGTGAAAAGTGATCAACAACAAAGTTATTATAAATGATCATTATCTATACCTTTTTGTTTTGGTAATTTTGGTATATGACCTCTGTTTCAATAATTTGAATTTAAATTTCAAATTATGACCACCTGCAAACAAAATTTACAAATACTAAACGATTTCAACTGAAAAGTCATCAACAACAAAATTGTGTAACTCATCAATATCTATAACTTTTATTTTTGTCATATCTTCATCCGATAAAGTGATATTAACGTTGTTTACAAAGTTTACAAGTCATCTTATAAACTTTGATGAATTTAAATCAGTGACAAGAATTGCTATAACGTCTCAGGTTTGGGCCGGTTTATCACGCAACCGATAAACTGTCACGGTATTATGGGCCTGTATTCTGGGCCTTTCACCCATAGTTGCATATCCGTGACGATATCTTACAACATCACGAAAATCCCCTGGCCCCGCAGCGGACGTGATAATCCACAAATGTTCTGGTAAGTGAAATTGTGCTGGTCCACAGACGGTCGTCAGTCCCGAGCGATTATTGTCCAATTGAGCCTGCTAGTGAGACGGGGCGCTGGTGCTCGTTGGGCAGACGGGACAACCTGGCTCCAGTACCCAAACTGCCCCTGCCTACAAGGCTGGACCGCTTCATAAGAGATCCCTCGAACATAGAGGACGAGCTCCGGAACAAATTCTCAGGAAGAGGCTGGGGGATCAAAGGTCCACCCATGACTAGCCCAGATTCAGGCAGTCGTTCGAGGTTGTGTTTGTCACTGCAAAGGACTAACGGCCAACAATCTTAACTAAGATCGCGTCCTCTGATCAGTGTTGAGCGGATACTGCCGATGATCTGGCCTTAGTGCCCGTTAAGACATTACAGCGGAAGTCACACCGCCAACCGACCGGCTCAGAGTCAGCTCCACTCACACCGTAAGTGCAGAAGTAGAGGATGTTCCAAGGGTTAGTAGCGCCCGTACTCGTGGTCCTCCGGGAGTTCTGTTAAAAACAAACTTTATGAGGAAGTAGAGGGGACTCAAACGGATGGCCAGCAGGTCTGGCGAACTTATTATGCGGCTTAATCGCTGCTTTGCGACGTTTTCGGCCCCACTAGCTGGGGGAGCCGGTGTCACCACCACACGCTTAATTGTGGACCCATACGAGGGGAACTCGCCGCGCTCACGTACTGGAGGTTTACTAATTTGAGTCGGTGTTACGATTTTCACTAGGTACCCGGGTAAAAGCGCCGATATATCATCATTAAGAAGTTGCACCAAACTGACAGTCAGACAGGCCCAGGCTGGCTAGTAGCTCAGATGTCGACAGACAGAGGTCCTTGCTCCCCTCGTCATGCGGGTATTACTCTTGGGGCCTGAGGCTCTCTGGATGGATTCGGCGTTTTAAGTGTCACCTTATAAACTCCCCTGTTTGCATTGAATGCGTGAACCACCTGGCATCTTAGAAGATGAACCGACGTCTAAATCTACTGAAACATTCCCTATATAGGGGTTCGTGGGACCATCCTACCCGTAGCTGAAGTAAACATTTTTGCCAGCTATCCATGGGCCCAGGAGATATGGCTGGTTTTTATCAGGAACCCTCGTGGAACTACAGCAGGACCTCCGCTTTTGCTTTGCTAACCGTTTTCAAAGAGGCCGTAACGCCAAGAGACGAACATTCGTCCCAGCCGCTTTATTCCGTGGCGGCTCGACCGATTTGTTAAATTTCCGAGAAATCCTGGATGAGGACCGGTACAGCGTTGATTTACAGCACGGTTGTTACAGGTAATACGTCGTGACAAAACCACGGAACTACCGGCTGTTATAGCATGAATCGGTATTTTTAAGGACGCTGATGCTCGTTGAGCTTTCCGTTACATATTCTCCCGCCATAGCACCTCAAAGGTTTTAGTTCATCCTCCACTTCTCAGTGCGCCTTATCGTGTCATGTCCCCACGAATCCGTGGAAATCCTGGAACTATGAGGATCAACTGCTTTACTGCCCTGAGCGACTAATGTTCACCATAAGGAACCATAAATAAAGGCAAAGCGTCGATCTCACTAGTGCAGACTGATTGCCTGGGGTTCCACAGCCCAGTACAGTCAACTGTGATGACTGCGCTGGCTTTCTGCCTATATCTCGGGATCTCACTGGAGTACCCAGTTAGGTCCGCGTTTGGACCATCACCCCAACTAGGGGCCTGGCGACACACAAACGATTTTACAGACAGCGGAGGATACCTCGGGAAACAGCGCCAGTCTACCGGCCGTGGAGGTGATCGACGATATCTCTGCCTCTCCGAAGCAATCATCGTATTCTTCTCAATAAAGCCTGGAACGCTCAGCGCTTCCGACCCATCTCGGAGGGCACCCTTGCCAGGCCGTTTTACCCTGCTACCCCCGCCGCTAGCATTTATACTCATATATTTTGCGATCCGCACCGCTCCCAGTGGCATGTGATACTATTGGTAATAGAAAATTGTCAAGATAAGCTTACTCAGTCCTAACACCCCAACAGGGTCGAACACATCATTACGCTATACTCGACAGGAGACCGTTAGCACACTTCTTCCGCTTTGGGAGCAGACCCTTTGCAATAGCTTCAACATACACTTCCGGGTGTTACTAATGTCCGACCCGATCAGATTTTCGCCAAATTTCCTTACACATTAAATAGTTTCACATGTGATAATCTATATTATTAGGCCATAACTAACCATGCATGCAGATACCCTGTTAGCCACAGACCCTAATCAGCTAATACCACGTAGTGTCTATACAAATTACGGTGTTAATGACTGAAACCTGACAGGGAAATCTCGCCTATGTATACATGGATCTTGAACATTCTAGGGCAACTGAATTAATATCTAAGGCTGAACAAAGGTATATGTATATTGAGTACCTTGTAAAATAAAAGTACTATAACTATCTTTCTAGATAGCTGATGCATGCTAAGATATGTATTGAAATCTCTAAGAAAATTTAGCACAGGTTGTGAGTAAATTCTCCTTTGCTAGCTAGTTTTTCTTTAGAGAGGTTAAGTCGTAATTCACTAGTTGATATTAATAATAGTTTATTTACTACAACTGTTAAGAAAGGAGCAAGGGGACAGTCCAAACCACGAACTAGCCAACAAATGAGGTGATGTTAGAAGGACAATAATTAACAAGGGTCAAGATGAGTAAAATGGTTGCTAAGGTGAAGTTTTAGCCAGACCGTTAATATGTTAATTGAACATAATAGCCTACCTCGGATAAATCTCAGCCTCGATTATTAACTGACATGTAAAGTGGTCGATACTGATCAATTGCGTACAAGTGTCCGTTTGGATACCATCTACCATTGCTCAGAGTCCTCATTATCTAACAGATATTGCATCTCCCTCGCTATATTACAGCATATCAGCTCAGTTTCTAAATATTTAACAGTCTGGTGGACTTTTAAACATGTTTATTGATCGTTTTATTGAAAAATTTAATTAATTATTAATTATTTTCCTCTCATTTGATTCATTAATTAAATATATACACTTATATATATACATATATATATCTATAGTTTTAAATATTTCATAAATTTTATTAATTTGACAGATGGTCAAATATGTGAAAGGACAACAGTCAATCATCTATTAAAAAACTGAAGGGAGGGATCAGTATAATATTAGCGCCATTTTAGAGCCAGCTCGTCTATAAATACAGCAAATTGCATTGCCATAGTCATCCACCTCTACTTTTACTGCTAGCTAGCAGCAGCTCGATCGATATAGCTATACCATATATGGAGGCGTTTCAAGGAAACGACCACTGCCTGCCAGTGCTTGCTGTGAGTCGGTAACAACCATTGTTTTGATGTAACAACGTCGGCGGCCAGGGGGCACGGCACAGGACTTCGTGGACCTGCACAACGCGGTGGGGACAGGGATCGGCGAGGGTGACGTGGGACGACACCGTGGTGGCGGTGCGCGGGTCCACCGCGGCGCAGTCGCCAAGCCTGAATCACCAGCCCCGTAAGCACGAAGACAACGGAACAGCCACCTACGGAGAGACCATCTGCATGATCGATCATTTCCCGCCGGCGGCAACGACACGATCGTCGATCCGCCGGCCGCCGCGGCCGCGGGCGCTCAGGCGGCGGAGGGCAGTGGTTCGACCCCGACACCAGCGGCTTCTCGGCGCAGGGGGACTCGTTGCGACCACCACGCAGCTGATCGCTGGAACGCCACCACGGCCATCGCCTACGCCGAGGTCGTCTGCGACCAGCGACGCCCGGGTCTTCGTCATCTGCAACTACTACCTCCGGGCAAATCCCATTCCAGAGCCCCTACTAGCTAGGTTGCTTTGAAACCACAAGTCAATATTAATGATTAATCTATATATATATATATATCTATATATATATATATATATATATCTATAGCTACGTTATATACTATAATCGCATGCGATCATGCATGCCTGCATATACTGTTCAGTATACATAAATACCCTATATGCTCTACGTACATACTACTAGTATGTACAACGCATACTATATGTACATTGCCTACGCACGTTCAATCTAGAACGATATCAGCATGCATGTGATTGATGGAGATAAAAGTAGATCATCAATTCACATCATCTGTTGGCTAGTACCTAGATATATGCATCATCAACATGCAGGTGTACTGTACAACGGTTTATATATCCTAGTAGTGGTGC\t&388@=825368::<>+$'&'10401465=9:2-.+.)'3((('(%%/)*&-&(%)+()+1172.7$)&+(3*.*238;=>:9;<=9&'+3*,4134522+944:5>8+)'')&./()*.+)/.,68%823*1+1+)(*137>=<A?A=14)$$&%&')&47:0+/010*5%*)-*''3234:::*&(+)(&'$7<A?5569=;78,,.@516109-.5=012%&$%'%,.)$$$%$&.,--,,'1--152185123483)$'%#&(0'.(**...3>>C;1/++./8894969.11651/134:;?<6875/5;:788;:==978:BD5464*8&)33*0*+04479:,(--&&,,.578777.-,)767@=;245-&'*'+(-)')'*//-79834552.865764393<;<6A5A?9>?DC<.7<:82.(,0-0-++%147*25=>?7673<=@A;75=89>:<>66<;*2.010%,()+++<2>;?9:)*59?@G@6CFMMA<297/.-,8874*-+%%+)'*%$$('7%'-345C6/923$',/0461&&.)(&)-5888+&&1245331%7895))=89?>;7'')&)*(*(+'+12.=?@:960/0/34,()+,.,3879:=4;;99<:9:4-32328;490112*0.//%)+++%',-241+((;>8>9944,?=@@*%:/*6++-./020.-)*$($%/174:)-547888:;<;78>DDGA:BAA?:'(2$(+36>;=144824-533*%+)$)$(1.(+)*57754152-&*++'($&$-1+049949;@2789763/.6/4565940'$#,-+(%)%*))*/45;89-%-,;04)47:8.&%%)..233-++'%##%&'&(*-279-;:73@8++,2424/)))+-.++-2542../0::??0?8979-%%1*+,,+.5014-04/.+-,+)010143./.-:%$&--5*11;<,+,2;;=64+3')*-.%$%%(&'(((1$$$&$$$113233215,+,32-./(;):;5)/04+))*0%$%&##%)#'%%3*&'&%$##%%.*+''9$)%%'(+544;>1/6665.2''(2.'%$*&%/'/4(+((''%%%%-(()-%#,-32,)).6%.)$)%111ACA?>=>=>3/02048>?A=94;;=C@@:;;.:*+)96564//4'(11250.,.,(()%#$($99/A23-.257<>@?;C7GGDB9*;>@@:2;@83/440(*0226.345%-/-.)*%)*(**$/0/*5-((/+.3?71'*.)*)%.'258)2;0.616:;6()%#(''%#$$'$%%$&+)),))##+0.'&0&$'34968/.'-,#%0),,0(1)');78<=:3.0,'-()6-%)0/++)07443(()(*-.*$$$$'+-.0/51655<558<92-3&(887;=>5;:--%&%(+606=C9<6=7C;3<96/.+--%%&&&')##'-$32:3.89$91488=5</89@@;<<9759=<90&,())(/,&93480/%&)($./42-*&5565<==(%**)54&+(&),&*+/943,9682334.'$%#%'$%''(')''&\"%%/+,+%#%$)&%%)###%%,0,0*0078>+541,-;<=77)-&&&&&%#+$33+,333*/'&$%1)3:54-+682/458776%($%$#)49>>>59=11122B@ACEBA<?<867>;;'%&67=>@A/=6=>=<;-*,*+)(*)/$')*(-$$##%%'%%2463362..'''$%$$%(&&&''*+4488651/-+42044:&)+2))%&'$-,.(&$&$#%%+%%&'()-&.///1;88542*)$$)%/)6&'&)&&*%##$'##%%%&'&&,%)04=6:91/-418400***.,//1/01-))+,.38+++*378-2/&$&10/,+'3./$(()1855+/(0&,3-,*(($,+.-0469:84-674&)%$$&&+'##$)))887AEEEAA@96;:=9671*'%%&-&'(*$$24/12-8875762,*'%%03-/-+0*+'')()0(+$#$&(//(*0**12147/($)*#$$(&%%(*&''68>>>)'672/00009*%$%$01/*0/-0-&&%$'04'((.--+'&8'542-350/$$('/.,*&$#&%$$32.'((*-47CBA4*225@38=,,13:@;>64=A;5.,3110(45//,3(2)%*&%%&'((5.*(%%8)2)45'$(*')&/.63253(,-41.,+(+-+AA7++--47?@A+8/.022222:464/'&$'8<78,39/*'($&%%)$)**&'&+0100435%00,*+/3189961467/1578;;;+//02172034835.-(+(#&*+.23610)($='044552218=78;(('*,+.'*+*%%&)0&<??A)(,()*)(-.545;676%$$$)*00153,=57-79::62$-&(%+))/*%3$##&(&-%&*$$(*'%&&+-,1)/3.0-3322/+((%'((.,*$$%&%//436=9532,)*(-/$%#*($*$%>??=?AAA@-+,-)*&&$&%%''&)%#',)'(+**%&+4-,03.)&+&%*-,./%$$+(#%$.,9,-'$*0$$#&''(&%&)*+,))(%'##$&*)*/+//(./+&&2.0%%%&%$&%*--120-&%3.2899:47142*&%..617()7;55-6*/1(($&$').)%&$$&'%&&5$%%'-+((#&')06)2$$$4592/3..173,,//0'457677;4621CEA55555.&&-**+(2'$&)***:()601&$$40391&'*)$$%&$(111111++.&')*'&%%)(()))'+-*0101/$/')($%$'/3,+&$%&.0,6455<23377''533:C:6(()48146<?A@+37.<96,07.<))%)--33<;--.%%%&&/205??<535314.757EAHGECD<0-*-+),*-*'%%%--$0297??J/:>@=789549CB:',)0/$)0/1.0$$(01544*%%%)'11124&&,.)''(22697,(&+.78761/*#'%$+42+(-4789+,((,***.%%10.642,''*%%%'$&%'&*-2,''3//2/%&+''%$$*4<::-7=>:>399;5./8;5;@A:<82249+9;9/3/>@B@7788>>8;;<=)<0-/2(9:<>==;6333520/641-'8:5114879-.8=??A@>04-1;1B<<<=<>=5579=<4$''&&/5680-$;*.',$&-##))-35.68732-.43/+*-1++23>??>AAA50/0'$%$'(((:=@=2)-&$$$#&))''./.(('(%-$('&&$$&%)%$%&$$%'$$%$%$&%*+$&*',.+++%($$$/-(,1/1-./+'&%$'&&2,&&(79981'&*''+*0'+*&&'&&#$%&&$$/0)0/32'*'%#$$&,($&$%*+434**41,-0*&#&'&&%%'#&/2<;51/)'%*),,0//001028>?=D.744576925752/75444),*(%&%#&'&'&++(+3-8+-2%&)++++,=:=;73/21.$*&((++.559;?79>?>07+4*.*35:>BA755854/0,)0473-0-,2**+,,.,++&%'$$$.')%'('%$&$%2(599=>@B<5>BBBAEAAB)>0:49&&$&)'.'(27-&*$($%#%$&(&#%%-+)*'+,&(8466&&($%+---'**),,.2-'%$$%#$#$&((%)%%$*&$(%$%#$&(.054.-3+(161-0-+((*+0,0-+*('(&)&)&'+&$&$'%$'(&(&'-.2-,,&'()%(('$+0*532:5;;;884+%(*+,*%$#%*$%&'137501/0/41%%'$&//++*-,((),6$)%%301224$*','*)*))4,$&1...-+35)&,0,../73/%&#(##$'&%%'++4)('&')*/010372-20380,-#$(+')'('&-$''&%'(&%(()#'&,,1+-/)%%''&)(&$#).%%$,%+(478;:56.*,,&22.$%%''*,+0<=7211-;9569995.0+-+----'(%'$'3::>=05..+)..1.2/100--.$%/67860-')+007++.3/-%%,20030;888+6-4230'(,1.7A@?A:67ED:6./'-7885,*$%)**-34.$$$6''(34461))$$$%3/1+&%)$%$'#$#('%&+002303;:=560$$/11,,7-233.)-%')((-.,,,,/75369;;<:5))-*6:A;A5;305)*-342.-)+&+(+''(,,,+,.76'$($%%)$%($&((.-./3330+$)*&%&(((&%$&)'',,-'+*,1,,*,))))*+''65(&)$++64,%3+'22+3+3))))'.-.%&%&8;=;>;55265'%,(4697((%&($%%'*)''(01.3+1++++--.&&%%$*+&-+138((//9&8978/.('0,$$''>?6123/+)+39;740))2568965:77.4A0&&(&&%)*%($%'(&)#0$3363593'(%55;8941-101*+)6,0/&&'5&&6$$$2424**((*,&')+%*(-,+&%%'#$(,),333,/-2/%22..'&$,-$%*&),,((&--.22'$&'--&&)(11+456-'',*%.(%%)%$#$',,()1104,.13,*$'$&$$$%*$%)*#%&&1)*-)0>@AG@*(.453.-',>7;=40%%,(),/)53.,,&&0../1-,*'21215/23503)%%*&&('':67--0/A>BAB<82-$$%,---&%&%&.0,8862,(,($%$'*:<<66821;24112.3+<;CHA85?BL?3(-/(1-,-*($'%&'&-.2-30;=>22./.,/%'-4/+),%.+,148;;84/.$/(''(.0)++40-8:';44+'0/.-*(/*0$00/837+.,//0**+)$($$$((+('(&('($%%$%#&))&'*-)-/0%020(2+40644321/.*%%+/0/-..3*41412672,''('+)215753.)&%,/.1467352306236847898.-'-)9;;==87350,.&$$%&&+'&&-/)##,-,/1..,&-//2255@<8;;=9<::=9<@>98<>@DF;;99;A>?>995971(..0,3,&-/-1)$)#%%+/69037234789.(,(%&&*()))'$%&(),,20+4&/'*,-)%$%&%&%'))%'$-0)-,-3//32)))//2-*11131*)()&&&%-*00-(%$$%&&()'&&-+0,235;3/4951,-26:,(%'()$$*)-&$$$&$%&$$\tNM:i:294\tms:i:2336\tAS:i:2336\tnn:i:0\ttp:A:P\tcm:i:48\ts1:i:502\ts2:i:0\tde:f:0.1255\tSA:Z:Reference,495158,+,446S515M2I4253S,60,44;Reference,480912,+,446M7D4770S,60,38;\trl:i:30\n"
     ]
    }
   ],
   "source": [
    "head -n3 Clone${i}_ONT.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Convert sam to bam</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Convert samtobam and filter it \n",
      "\n",
      "[samopen] SAM header is present: 1 sequences.\n"
     ]
    }
   ],
   "source": [
    "echo -e \"\\nConvert samtobam and filter it \\n\"\n",
    "samtools view -@8 -bh -S -F 0x904 -o Clone${i}_ONT.bam Clone${i}_ONT.sam\n",
    "rm Clone${i}_ONT.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Sort and index bam</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Sort and index bam \n",
      "\n"
     ]
    }
   ],
   "source": [
    "echo -e \"\\nSort and index bam \\n\"\n",
    "samtools sort -@8 Clone${i}_ONT.bam Clone${i}_ONT_SORTED \n",
    "samtools index Clone${i}_ONT_SORTED.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Calculate stats from mapping</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Calculate stats from mapping\n",
      "\n"
     ]
    }
   ],
   "source": [
    "echo -e \"\\nCalculate stats from mapping\\n\"\n",
    "samtools flagstat Clone${i}_ONT_SORTED.bam >Clone${i}_ONT.flagstats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> Display the content of the flagstat file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9281 + 0 in total (QC-passed reads + QC-failed reads)\n",
      "0 + 0 duplicates\n",
      "9281 + 0 mapped (100.00%:-nan%)\n",
      "0 + 0 paired in sequencing\n",
      "0 + 0 read1\n",
      "0 + 0 read2\n",
      "0 + 0 properly paired (-nan%:-nan%)\n",
      "0 + 0 with itself and mate mapped\n",
      "0 + 0 singletons (-nan%:-nan%)\n",
      "0 + 0 with mate mapped to a different chr\n"
     ]
    }
   ],
   "source": [
    "head Clone10_ONT.flagstats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBC;\"> => Let's map with data from all clones using a loop for mapping, with a single folder per sample and ONT reads<a class=\"anchor\" id=\"mapallminimap\"></span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in {1..20}\n",
    "    do\n",
    "        ONT=\"/home/jovyan/work/SV_DATA/LONG_READS/Clone${i}.fastq.gz\"\n",
    "        mkdir -p ~/work/MAPPING-ONT\n",
    "        cd ~/work/MAPPING-ONT\n",
    "        echo -e \"\\n>>>>>>>>>> Creation directory for Clone$i\\n\"\n",
    "        mkdir -p dirClone$i\n",
    "        cd dirClone$i\n",
    "        \n",
    "        echo -e \">>>> Mapping Clone$i minimap2\\n\"\n",
    "        minimap2 -ax map-ont -t 12 ${REF} ${ONT} > Clone${i}_ONT.sam \n",
    "        \n",
    "        # Convert samtobam \n",
    "        echo -e \">>> Convert samtobam and filter it \\n\"\n",
    "        samtools view -@8 -bh -S -F 0x904 -o Clone${i}_ONT.bam Clone${i}_ONT.sam\n",
    "        rm Clone${i}_ONT.sam\n",
    "\n",
    "        echo -e \">>>> Sort and index bam \\n\"\n",
    "        # sort and index bam\n",
    "        samtools sort -@8 Clone${i}_ONT.bam Clone${i}_ONT_SORTED \n",
    "        samtools index Clone${i}_ONT_SORTED.bam\n",
    "\n",
    "        # Calculate stats from mapping\n",
    "        echo -e \">>>> Calculate stats from mapping\\n\"\n",
    "        samtools flagstat Clone${i}_ONT_SORTED.bam >Clone${i}_ONT.flagstats\n",
    "    done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clone20_ONT.bam        Clone20_ONT_SORTED.bam\n",
      "Clone20_ONT.flagstats  Clone20_ONT_SORTED.bam.bai\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#006E7F\">__III -  Centralize final mapping data into a single bam directory__ <a class=\"anchor\" id=\"reorder\"></span>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Reorder BAM files into a folder only for Illumina</span>  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p ~/work/MAPPING-ILL/BAM\n",
    "cd ~/work/MAPPING-ILL/\n",
    "\n",
    "for i in {1..20}\n",
    "    do\n",
    "         ln -s ~/work/MAPPING-ILL/dirClone$i/Clone$i.SORTED.bam BAM/\n",
    "    done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 0\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone10.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone11.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone11/Clone11.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone12.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone12/Clone12.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone13.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone13/Clone13.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone14.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone14/Clone14.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone15.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone15/Clone15.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone16.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone16/Clone16.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone17.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone17/Clone17.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone18.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone18/Clone18.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone19.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone19/Clone19.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone1.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone1/Clone1.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 59 Jun 20 16:37 Clone20.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone20/Clone20.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone2.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone2/Clone2.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone3.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone3/Clone3.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone4.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone4/Clone4.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone5.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone5/Clone5.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone6.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone6/Clone6.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone7.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone7/Clone7.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone8.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone8/Clone8.SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 57 Jun 20 16:37 Clone9.SORTED.bam -> /home/jovyan/work/MAPPING-ILL/dirClone9/Clone9.SORTED.bam\n"
     ]
    }
   ],
   "source": [
    "ls /home/jovyan/work/MAPPING-ILL/BAM -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Reorder BAM files into a folder only for ONT</span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p ~/work/MAPPING-ONT/BAM\n",
    "cd ~/work/MAPPING-ONT/\n",
    "\n",
    "for i in {1..20}\n",
    "    do\n",
    "         ln -s ~/work/MAPPING-ONT/dirClone$i/Clone${i}_ONT_SORTED.bam BAM/\n",
    "    done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 80\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone10_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone10/Clone10_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone11_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone11/Clone11_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone12_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone12/Clone12_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone13_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone13/Clone13_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone14_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone14/Clone14_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone15_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone15/Clone15_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone16_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone16/Clone16_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone17_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone17/Clone17_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone18_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone18/Clone18_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone19_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone19/Clone19_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone1_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone1/Clone1_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 63 Jun 20 16:37 Clone20_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone20/Clone20_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone2_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone2/Clone2_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone3_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone3/Clone3_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone4_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone4/Clone4_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone5_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone5/Clone5_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone6_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone6/Clone6_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone7_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone7/Clone7_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone8_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone8/Clone8_ONT_SORTED.bam\n",
      "lrwxrwxrwx 1 jovyan users 61 Jun 20 16:37 Clone9_ONT_SORTED.bam -> /home/jovyan/work/MAPPING-ONT/dirClone9/Clone9_ONT_SORTED.bam\n"
     ]
    }
   ],
   "source": [
    "ls ~/work/MAPPING-ONT/BAM -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
