{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# __DAY3 - MAPPING PRACTICE__"
   ]
  },
  {
   "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 :\n",
    "* the ___Bathycoccus prasinos_ genome__ - File : `GCA_002220235.1_ASM222023v1_genomic.fna`\n",
    "* the __ILLUMINA READS of the 1613 sample__  - Files `1613_R1.fastq.gz`, `1613_R2.fastq.gz`\n",
    "\n",
    "The tool bwa needs 2 steps: \n",
    "- **Reference indexing**: `bwa-mem2 index reference`\n",
    "- **Mapping in itself**: `bwa-mem2 mem  -R READGROUP [options] reference fastq1 fastq2 > out.sam`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference indexation - `bwa_mem2 index`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before mapping, we need index reference file! \n",
    "\n",
    "#### Go into REF directory and list the content of the directory - `cd,ls`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Print the first lines of the file that contains the genome sequence (fasta format) - `head`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Before mapping we need index reference file! Check bwa-mem2 index command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# COMMAND TO GET HEP FOR bwa-mem2 \n",
    "bwa-mem2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND TO CREATE INDEX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check that the indexes have been created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ls -lrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Run the mapping with `bwa mem` for only one sample (sample 1613)\n",
    "\n",
    "\n",
    "#### Go into the directory SV_DATA - `cd`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a direcctory `RESULTS` and go into the directory  - `mkdir, cd` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# COMMAND TO GET HEP FOR bwa-mem2 mem\n",
    "bwa-mem2 mem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND TO PERFORM MAPPING"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bwa-mem2 mem -R PUT_READ_GROUP -M -t 8 PUT_REFERENCE_PATH PUT_R1_FASTQ PUT_R2_FASTQ > PUT_SAM_PATH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Check that the file `.sam` have been created by `bwa mem` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display the first and the end of the sam file just created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMANDS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sort sam file and convert it into a bam file (compressed format) -  `samtools sort`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samtools sort  -@ CPU_NUMBER -o PATH_2_BAM_SORTED PATH_2_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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate some stats about mapping with `samtools flagstat`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools flagstat PATH_2_BAM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate some stats about mapping with `samtools flagstat` and save results in the file `1613.flagstat`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display the content of the flagstat file - `cat`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate a bam file that contains only the reads correctly paired mapped `samtools view`\n",
    "\n",
    "https://broadinstitute.github.io/picard/explain-flags.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools view -bh -@ CPU_NUMBER -f FLAG -o PATH_2_BAM_FILTERED PATH_2_BAM_INPUT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lhrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Indexing the last bam file generated that will be used to detect SNPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Let's mapping an other sample (1868) to reduce steps and tempory bam files using the pipe `|`\n",
    "\n",
    "```The pipe is a form of redirection in Linux used to connect the STDOUT of one command into the STDIN of a second command.```\n",
    "\n",
    "#### 1- Generate the bam file 1868.sorted.bam\n",
    "Use a single command line to perform :\n",
    "* mapping with `bwa-mem2 mem`\n",
    "* sort the mapping results with `samtools sort`  \n",
    "* and only the file `1615.sorted.bam` is generated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND\n",
    "bwa-mem2 mem ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filtering with samtools view and create the index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting some stats on mapping step - `samtools flagstat`\n",
    "\n",
    "We want the stat on the two bam files generated :\n",
    "* 1868.sorted.bam\n",
    "* 1868.onlymappedpaired.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lht"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the two flagstat files. What do you think about the results ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat *.flagstat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## SNPs calling using data from 1613 and 1868 samples - `bcftools mpileup ... | bcftools call ...`\n",
    "\n",
    "#### Print the help of `bcftools mpileup`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PUT COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Print the help of `bcftools call`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PUT COMMAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PUT COMMAND TO PERFORM SNP CALLING\n",
    "bcftools mpileup ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -l "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Convert bcf nto vcf files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcftools convert -Ov -o PATH_2_VCF PATH_2_BCF\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Have a Look to the vcf created - `head, tail`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "head *vcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tail 1613_1868.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Count the number of variants with `bcftools stat`\n",
    "* run the bcftools stats on the vcf file and save the result into the file ` 1613_1868.bcf.stat`\n",
    "* check that the file have been correctly created and display the 35 first lines of this file\n",
    "* How many samples were used for this SNP analysis ?\n",
    "* How many SNPs were detected ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PUT_COMMAND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Download additional notebooks for exercices about mappin and SNP calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Go into the correct directory to download notebook - cd command\n",
    "cd /home/jovyan/NOTEBOOK/DAY3-4 \n",
    "wget https://itrop.ird.fr/trainingts/notebook/DAY3-4/Day4_Mapping_Exercice1.ipynb --no-check-certificate\n",
    "wget https://itrop.ird.fr/trainingts/notebook/DAY3-4/Day4_Mapping_Exercice2.ipynb --no-check-certificate"
   ]
  }
 ],
 "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
}