{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d015209b-0d23-4c85-994c-79b3535c75a3",
   "metadata": {},
   "source": [
    "# __How to identify large variants using whole-genome assemblies ?__ \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a61499d9-a943-4628-a3c2-d622c86db2af",
   "metadata": {},
   "source": [
    "#### __Preparing working environment__ \n",
    "\n",
    "* Go into the `ASSEMBLIES` directory\n",
    "* Check that the directory contains two genomes assemblied from Long reads sequencing data for the `5417` and `A8` samples\n",
    "* Display first lines of fasta files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfdb22be-7165-415c-a7a1-f38ad765ada2",
   "metadata": {},
   "source": [
    "#### Count the number of sequences of the fasta files in the directory `ASSEMBLIES`   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b75773d-38a6-49ab-9de9-34681cfd44f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "cd ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf7ebcd-112f-4849-a017-db793c2696be",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "head ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "259b9c77-4326-4daf-b03b-ef669c4e2148",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "grep ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91929ba1-ab26-46d5-87b2-8c163a8f0d24",
   "metadata": {},
   "source": [
    "#### Calculate the size of the two genomes  in the directory `Assemblies`  - seqtk comp \n",
    "\n",
    "Setk comp provides information on each sequence in the fasta file:\n",
    "* column 1 contains the sequence name \n",
    "* column 2 corresponds to the sequence length in bp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccb2453d-02fa-4cec-915a-9136e27647e7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "seqtk comp FASTA_FILE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b4cfc16-c249-4456-957f-95a4e8fcb241",
   "metadata": {},
   "source": [
    "##### Compute the size of the genome 5417 using `seqtk comp | awk`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b79e2a4a-a505-4ed1-9763-9c34e10271e2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c65b4b47-c99c-4673-971d-c198e15b5a5d",
   "metadata": {},
   "source": [
    "#### __SV calling using nucmer, syri__ \n",
    "\n",
    "We will use :\n",
    "* nucmer to align the two genomes. [nucmer manual](http://mummer.sourceforge.net/manual/)\n",
    "* syri to detect SVs from nucmer alignment. [syri manual](https://schneebergerlab.github.io/syri/)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8ff0702-daa3-4f01-9cf6-3ee6f3bbf241",
   "metadata": {},
   "source": [
    "#### Create the working directory for SYRI analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41aa1a28-f657-41d6-98e8-83cc41bc5bc1",
   "metadata": {},
   "source": [
    "* Go into the`RESULTS` directory and create the directory `SV_CALLING_SYRI`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edb0ebe1-3a15-4fbc-8688-e5df187c6576",
   "metadata": {},
   "outputs": [],
   "source": [
    "cd ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76322897-1433-4230-8d2b-499823a37b01",
   "metadata": {},
   "source": [
    "#### Aligning genomes using `Nucmer` \n",
    "\n",
    "__Alignement of the two genome with `nucmer`__\n",
    "\n",
    "[Nucmer manual](http://mummer.sourceforge.net/manual/)\n",
    "\n",
    "Some interesting parameters :\n",
    "<code>\n",
    "-maxmatch       Compute all maximal matches regardless of their uniqueness\n",
    "-b|breaklen     Set the distance an alignment extension will attempt to extend poor scoring regions before giving up (default 200)\n",
    "-c|mincluster   Sets the minimum length of a cluster of matches (default 65)\n",
    "-l|minmatch     Set the minimum length of a single match (default 20)\n",
    "</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bac7fd8e-c4f9-4c6c-a156-d3e94934e485",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "nucmer --maxmatch REFERENCE_GENOME ASSEMBLIED_GENOME"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a62a1b3-4e63-4e8d-b93c-c54a7c17cde9",
   "metadata": {},
   "source": [
    "__Check if new files have been generated by nucmer and display the first lines of the file `.delta`__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee2ba7af-8b94-4df3-8412-928d6b421234",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fea08d51-a81a-47b7-a31d-382a0c2d95e0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "head out.delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f988b1a-bb3e-4920-aae0-d4e7bd07c69d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7e004fdf-f6a8-4512-b86c-cc74ed85705b",
   "metadata": {},
   "source": [
    "__Filtering nucmer results with `delta-filter`__\n",
    "\n",
    "We are going to remove small and lower quality alignments. Some interesting parameters :\n",
    "<code>\n",
    "-i float\tSet the minimum alignment identity [0, 100], (default 0)\n",
    "-l int\t\tSet the minimum alignment length (default 0)\n",
    "-m            Many-to-many alignment allowing for rearrangements (union of -r and -q alignments)\n",
    "-q            Maps each position of each query to its best hit in the reference, allowing for reference overlaps\n",
    "-r            Maps each position of each reference to its best hit in the query, allowing for query overlaps\n",
    "</code>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11b64fb3-6ead-4c47-85e6-6a5fa4049924",
   "metadata": {},
   "outputs": [],
   "source": [
    "delta-filter -m -i 90 -l 100 out.delta > out.filtered.delta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac7b2015-2c91-4433-8a1e-e75a7363dbb7",
   "metadata": {},
   "source": [
    "#### Check that the new `delta`files has been filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "467c1f2e-62d3-477b-b549-6de9c0f46619",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b46169b9-b8e6-4c32-96a9-17e650a57bd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "wc -l *.delta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "553391e9-2c9f-4aae-ac32-8c61b7836907",
   "metadata": {},
   "source": [
    "__Converting the file `delta` into a tabular file using `show-coords`__\n",
    "\n",
    "[show-coords manual](http://mummer.sourceforge.net/manual/#coords)\n",
    "Some interesting parameters :\n",
    "<code>\n",
    "-c\tInclude percent coverage columns in the output\n",
    "-d\tInclude the alignment direction/reading frame in the output (default for promer)\n",
    "-H\tOmit the output header\n",
    "-I float\tSet minimum percent identity to display\n",
    "-l\tInclude sequence length columns in the output\n",
    "-L int\tSet minimum alignment length to display\n",
    "-q\tSort output lines by query\n",
    "-r\tSort output lines by reference\n",
    "-T\tSwitch output to tab-delimited format\n",
    "</code>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c2b20c4-665a-4cc5-8240-e191edd66582",
   "metadata": {},
   "outputs": [],
   "source": [
    "show-coords -Trd out.filtered.delta > out.filtered.withHeader.coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9eb2b0f-f716-413e-ae14-76e70541e011",
   "metadata": {},
   "outputs": [],
   "source": [
    "show-coords -THrd out.filtered.delta > out.filtered.coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f34fceb2-a76e-41be-8a69-cf59e97564aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "929dc4d1-1a08-43ac-b7a2-82880b9b807b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "head ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c84c259-cdc7-422e-9c7f-cc09838418c9",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> SV calling using syri <a class=\"anchor\" id=\"siri\"></a></span>  \n",
    "\n",
    "__activate syri envt with conda__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fda3a8b2-7768-4775-bbff-ba67b551d5ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "conda activate syri_env"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dbf9a7f-d961-4791-904a-8dd12b986617",
   "metadata": {},
   "outputs": [],
   "source": [
    "syri -c out.filtered.coords -d /home/jovyan/SV_DATA/ASSEMBLIES/out.filtered.delta -r REFERENCE_GENOME -q ASSEMBLIED_GENOME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3f88d52-7edc-4d2c-a652-9ce1ac5531b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65ba6839-f680-4dad-ad8e-5d1ed114ab85",
   "metadata": {},
   "source": [
    "__Output file format__\n",
    "\n",
    "SyRI outputs results in TSV format and VCF file format.\n",
    "\n",
    "<code>\n",
    "Annotation\tMeaning\t \tAnnotation\tMeaning\n",
    "SYN\tSyntenic region\t \tSYNAL\tAlignment in syntenic region\n",
    "INV\tInverted region\t \tINVAL\tAlignment in inverted region\n",
    "TRANS\tTranslocated region\t \tTRANSAL\tAlignment in translocated region\n",
    "INVTR\tInverted translocated region\t \tINVTRAL\tAlignment in inverted translocated region\n",
    "DUP\tDuplicated region\t \tDUPAL\tAlignment in duplicated region\n",
    "INVDP\tInverted duplicated region\t \tINVDPAL\tAlignment in inverted duplicated region\n",
    "NOTAL\tUn-aligned region\t \tSNP\tSingle nucleotide polymorphism\n",
    "CPG\tCopy gain in query\t \tCPL\tCopy loss in query\n",
    "HDR\tHighly diverged regions\t \tTDM\tTandem repeat\n",
    "INS\tInsertion in query\t \tDEL\tDeletion in query\n",
    "</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fdcb7c6-d698-482c-b2a1-6d8f5f2b60e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat syri.summary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6578624d-3ee2-47f8-9943-6ff10704e303",
   "metadata": {},
   "source": [
    "__TSV format specifications__\n",
    "\n",
    "<code>\n",
    "Column Number\tValue\tType\n",
    "1\tchromosome ID in reference\tstring\n",
    "2\treference start position (1-based, includes start position)\tint\n",
    "3\treference end position (1-based, includes end position)\tint\n",
    "4\tsequence in reference (Only for SNPs and indels)\tstring\n",
    "5\tsequence in query (Only for SNPs and indels)\tstring\n",
    "6\tchromosome ID in query\tstring\n",
    "7\tquery start position (1-based, includes start position)\tint\n",
    "8\tquery end position (1-based, includes end position)\tint\n",
    "9\tunique ID (annotation type + number)\tstring\n",
    "10\tparent ID (annotation type + number)\tstring\n",
    "11\tAnnotation type\tstring\n",
    "12\tCopy status (for duplications)\tstring\n",
    "</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6de3d0f9-2198-4c0f-b813-6b4d3e41ab3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "head syri.vcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9435c82-45c7-497e-aa9c-e49c336c2925",
   "metadata": {},
   "outputs": [],
   "source": [
    "head syri.out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6825e0b-9260-4fe7-be3e-c410f7fd92ba",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBC;\"> Extracting all SNP from syri  <a class=\"anchor\" id=\"siri\"></a></span>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3ecf3d7-a4b2-4d7b-85fa-a32b0ee60ede",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat syri.out | grep SNP | cut -f 1,2 > a\n",
    "cat syri.out | grep SNP | cut -f 3 > c\n",
    "awk '{ print $1+1 }' c > d\n",
    "paste a d c > SNPs.bed\n",
    "rm a c d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7351688d-ccde-484c-9710-e3aa34afd7d0",
   "metadata": {},
   "source": [
    "__Display the first lines of the file SNPs.bed__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79168401-fac7-453e-9962-b2407ed4e72e",
   "metadata": {},
   "outputs": [],
   "source": [
    "head SNPs.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a68c6262-b2f9-4cfe-bc14-20143cb45678",
   "metadata": {},
   "source": [
    "__create the file plotsr_pos.txt required by the following commands__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ae51f46-f320-43ba-a20b-81f93d569c2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "echo -e $reference_assembly'\\t ref' >> plotsr_pos.txt\n",
    "echo -e $query_assembly'\\t query' >> plotsr_pos.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af2dcd43-0ea9-4577-9a76-1899f3730d07",
   "metadata": {},
   "outputs": [],
   "source": [
    "head plotsr_pos.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45473bfc-619f-45e3-82f2-fcb006cfe015",
   "metadata": {},
   "source": [
    "__Generate the plot using `plotsr`__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ffd3a6-6bf7-41da-a78c-84219fe96878",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd\n",
    "plotsr --sr syri.out --genomes plotsr_pos.txt -s 500 -o plotsr.pdf -H 8 -W 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0021750c-ccdc-456e-aec3-0455d8423824",
   "metadata": {},
   "source": [
    "### => download the file plotsr.pdf\n",
    "\n",
    "--------------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0695741d-3dd3-4909-8fa6-e9a40aec4667",
   "metadata": {},
   "source": [
    "## <span style=\"color:#006E7F\">__II - using assemblytics__ <a class=\"anchor\" id=\"SVassemblytics\"></a></span>  \n",
    "\n",
    "* Download the file .delta generated by nucmer. \n",
    "* Load this file on the assemblytics website : [http://assemblytics.com/](http://assemblytics.com/)\n",
    "\n",
    "---------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a75dd6f-e890-4dee-afe9-bdc32b5c05e1",
   "metadata": {},
   "source": [
    "## <span style=\"color:#006E7F\">__III - using d-genies__ <a class=\"anchor\" id=\"SVassemblytics\"></a></span>  \n",
    "* Download the two genomes (fasta files)\n",
    "* Load these files on the d-genies website : [http://dgenies.toulouse.inra.fr/](http://dgenies.toulouse.inra.fr/)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6b926f-aaad-4dca-ae45-6a68f0aa03ef",
   "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": 5
}