{
"cells": [
{
"cell_type": "markdown",
"id": "d015209b-0d23-4c85-994c-79b3535c75a3",
"metadata": {},
"source": [
"#### __Formation South Green 2022 - Structural Variants Detection by using short and long reads__ \n",
"\n",
"# __DAY 4 : How to identify large variants using whole-genome assemblies ?__ \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",
"# Table of contents\n",
"\n",
"\n",
"[Preparing working environment](#env)\n",
"\n",
"[I - using nucmer & syri](#SVsyri)\n",
"\n",
"* [Run minimap2](#minimap)\n",
"* [Run syri](#syri)\n",
"\n",
"[II - using assemblytics](#SVassemblytics)\n",
"\n",
"[III - using d-genies](#SVdgenies)\n",
"\n",
"\n",
"\n",
"\n",
"***\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a61499d9-a943-4628-a3c2-d622c86db2af",
"metadata": {},
"source": [
"## __Preparing working environment__ \n",
"\n",
"### Create the working directory ~/work/DAY4-WGA-SV "
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "aaf97983-70d6-4774-a4fa-cc76f7b5dd92",
"metadata": {},
"outputs": [],
"source": [
"mkdir -p /home/jovyan/work/DAY4-WGA-SV\n",
"cd /home/jovyan/work/day4-WGA-SV\n",
"ls"
]
},
{
"cell_type": "markdown",
"id": "e448f0aa-b72b-4376-8309-72ab7c58e907",
"metadata": {},
"source": [
"### Download the two genomes that we will \n",
"* Download the archive on the [itrop website](#https://itrop.ird.fr/)\n",
"* Decompress it\n",
"* Remove the archive\n",
"* Check that the directory Assemblies exists !\n",
"* List the content of directory Assembly"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "fe0fdf17-568f-4073-b34b-19c98abb8458",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--2022-06-23 21:05:37-- https://itrop.ird.fr/sv-training/Assemblies.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: 9150274 (8.7M) [application/x-gzip]\n",
"Saving to: ‘Assemblies.tar.gz’\n",
"\n",
"Assemblies.tar.gz 100%[===================>] 8.73M 37.5MB/s in 0.2s \n",
"\n",
"2022-06-23 21:05:37 (37.5 MB/s) - ‘Assemblies.tar.gz’ saved [9150274/9150274]\n",
"\n",
"Assemblies/\n",
"Assemblies/assembly-stats.txt\n",
"Assemblies/A8_assembly.fasta\n",
"Assemblies/5417_assembly.fasta\n"
]
}
],
"source": [
"wget https://itrop.ird.fr/sv-training/Assemblies.tar.gz\n",
"tar xvf Assemblies.tar.gz && rm Assemblies.tar.gz"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "30ab0113-1538-469a-963c-420ed1714979",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
".:\n",
"Assemblies\n",
"\n",
"./Assemblies:\n",
"5417_assembly.fasta A8_assembly.fasta assembly-stats.txt\n"
]
}
],
"source": [
"ls -R"
]
},
{
"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": 17,
"id": "adf7ebcd-112f-4849-a017-db793c2696be",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==> Assemblies/5417_assembly.fasta <==\n",
">Chr_19\n",
"CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCT\n",
"AAACCCTAAACCCTAAACCCTAAACCCTAAACCCAGAAAAATAGACTTTATTCCAACATAACTATCTACATCATCGTACGTACCTGCAGGTATGTGTTGA\n",
"AAGGTACCCGTTACACAAAAGAACTGAAACATTTGTCTAACACGTCAAATTCTGCTGCATTTGCCGTTATAAAACGTCTCAAAACTTTCGTTGAAAATAA\n",
"\n",
"==> Assemblies/A8_assembly.fasta <==\n",
">Chr_19\n",
"AAAAAATTTGGTCGTTTGAGCGTATATTTGCATACACCCTAAACAGAGCGTGCAAATATACGCTAAAAACGACACAAAATTTTTTTTTTTTTTTGCGAAA\n",
"GCACTCTAGTCGATCACGTTGTCTGAAAAAATTTAAATCGAAGTTTATAAGGCCAAGCTTGGTCAAAGGTCATACGTCATTATTTATACGTCACTATACG\n",
"TTTCTTCATACTTCGTTCCTTCACACAAAACGTAAAGATCTGAAATAGACGGATCCTTTGAGAGTTCATTGGAGAGATCTTTGAAAGTTGTTCTTTTTGT\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K19\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K19\n"
]
}
],
"source": [
"head -n 4 Assemblies/*fasta\n",
"grep -c \">\" Assemblies/*fasta"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "259b9c77-4326-4daf-b03b-ef669c4e2148",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_19\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_18\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_17\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_16\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_15\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_14\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_13\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_12\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_11\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_10\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_9\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_8\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_7\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_6\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_5\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_4\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_3\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_2\n",
"\u001b[35m\u001b[KAssemblies/5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_1\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_19\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_18\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_17\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_16\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_15\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_14\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_13\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_12\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_11\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_10\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_9\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_8\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_7\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_6\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_5\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_4\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_3\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_2\n",
"\u001b[35m\u001b[KAssemblies/A8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K\u001b[01;31m\u001b[K>\u001b[m\u001b[KChr_1\n"
]
}
],
"source": [
"grep \">\" Assemblies/*fasta"
]
},
{
"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 donne des infos sur chaque sequence du fichier fasta notamment:\n",
"* la colonne 1 contient le nom de la séquence \n",
"* la colonne 2 correspond à la longueur de la séquence en pb"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "ccb2453d-02fa-4cec-915a-9136e27647e7",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Chr_19\t233066\t66647\t50144\t48763\t67512\t0\t0\t0\t24058\t0\t0\t0\n",
"Chr_18\t318881\t85345\t74297\t75075\t84164\t0\t0\t0\t55346\t0\t0\t0\n",
"Chr_17\t482895\t126641\t114374\t113570\t128310\t0\t0\t0\t83524\t0\t0\t0\n",
"Chr_16\t518210\t136535\t122232\t123594\t135749\t0\t0\t100\t91800\t0\t0\t0\n",
"Chr_15\t521692\t134305\t126076\t124267\t137044\t0\t0\t0\t94334\t0\t0\t0\n",
"Chr_14\t641580\t186873\t135344\t135053\t184310\t0\t0\t0\t71686\t0\t0\t0\n",
"Chr_13\t737724\t188413\t182726\t176331\t190154\t0\t0\t100\t129222\t0\t0\t0\n",
"Chr_12\t734935\t191741\t177466\t175220\t190508\t0\t0\t0\t132112\t0\t0\t0\n",
"Chr_11\t748028\t195445\t178757\t183097\t190587\t0\t0\t142\t134678\t0\t0\t0\n",
"Chr_10\t803017\t206784\t193589\t192228\t210416\t0\t0\t0\t144998\t0\t0\t0\n",
"Chr_9\t911957\t235846\t220425\t219548\t236138\t0\t0\t0\t164128\t0\t0\t0\n",
"Chr_8\t985652\t255488\t237567\t237151\t255446\t0\t0\t0\t178310\t0\t0\t0\n",
"Chr_7\t992004\t250831\t234160\t234318\t251703\t0\t0\t20992\t173692\t0\t0\t0\n",
"Chr_6\t986735\t255731\t237563\t237449\t255992\t0\t0\t0\t177994\t0\t0\t0\n",
"Chr_5\t1013628\t260089\t246867\t243032\t263640\t0\t0\t0\t183632\t0\t0\t0\n",
"Chr_4\t1016660\t264661\t243987\t246425\t261587\t0\t0\t0\t183968\t0\t0\t0\n",
"Chr_3\t1096451\t283426\t265347\t266334\t281344\t0\t0\t0\t200088\t0\t0\t0\n",
"Chr_2\t1146138\t295896\t276008\t278603\t295631\t0\t0\t0\t209180\t0\t0\t0\n",
"Chr_1\t1387673\t352232\t339085\t335668\t360688\t0\t0\t0\t254350\t0\t0\t0\n"
]
}
],
"source": [
"seqtk comp Assemblies/5417_assembly.fasta "
]
},
{
"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": [
"seqtk comp Assemblies/5417_assembly.fasta | awk '{pb=pb+$2; print($2,pb)}'"
]
},
{
"cell_type": "markdown",
"id": "964471db-bb03-41b1-989d-81a2ac31c979",
"metadata": {},
"source": [
"##### Compute the size of the two genomes using `seqtk comp | awk` in a `for loop`"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "bc620c4c-59c9-4abd-9f46-4c915e3dd30a",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
">>>> Je lis le fichier : Assemblies/5417_assembly.fasta\n",
"233066 233066\n",
"318881 551947\n",
"482895 1034842\n",
"518210 1553052\n",
"521692 2074744\n",
"641580 2716324\n",
"737724 3454048\n",
"734935 4188983\n",
"748028 4937011\n",
"803017 5740028\n",
"911957 6651985\n",
"985652 7637637\n",
"992004 8629641\n",
"986735 9616376\n",
"1013628 10630004\n",
"1016660 11646664\n",
"1096451 12743115\n",
"1146138 13889253\n",
"1387673 15276926\n",
"\n",
">>>> Je lis le fichier : Assemblies/A8_assembly.fasta\n",
"225373 225373\n",
"308071 533444\n",
"474879 1008323\n",
"512960 1521283\n",
"523849 2045132\n",
"781878 2827010\n",
"770859 3597869\n",
"720863 4318732\n",
"764371 5083103\n",
"791990 5875093\n",
"894004 6769097\n",
"945596 7714693\n",
"1008950 8723643\n",
"989488 9713131\n",
"1016987 10730118\n",
"1006439 11736557\n",
"1079247 12815804\n",
"1131285 13947089\n",
"1359092 15306181\n"
]
}
],
"source": [
"for file in Assemblies/*fasta;\n",
" do\n",
" echo -e \"\\n>>>> Je lis le fichier : $file\"\n",
" seqtk comp $file | awk '{pb=pb+$2; print($2,pb)}'\n",
" done"
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "000253b0-a170-4f05-81c0-55262299d255",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/jovyan/work/day4-WGA-SV/SYRI\n",
"(syri_env) \n",
">>>>../Assemblies/5417_assembly.fasta\n",
"15276926\n",
"\n",
">>>>../Assemblies/A8_assembly.fasta\n",
"15306181\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"source": [
"for file in Assemblies/*fasta;\n",
" do\n",
" echo -e \"\\n>>>> $file\"\n",
" seqtk comp $file | awk '{pb=pb+$2;} END { print pb }'\n",
" done"
]
},
{
"cell_type": "markdown",
"id": "c65b4b47-c99c-4673-971d-c198e15b5a5d",
"metadata": {},
"source": [
"## __I - 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",
"### Initialization of two variables \"genome\" \n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "ac73f4e0-50c7-443d-a093-beb2cfb8751a",
"metadata": {},
"outputs": [],
"source": [
"dir_genome=\"/home/jovyan/work/day4-WGA-SV/\"\n",
"reference_assembly=$dir_genome\"Assemblies/A8_assembly.fasta\"\n",
"query_assembly=$dir_genome\"Assemblies/5417_assembly.fasta\""
]
},
{
"cell_type": "markdown",
"id": "f8ff0702-daa3-4f01-9cf6-3ee6f3bbf241",
"metadata": {},
"source": [
"### Create the working directory for SYRI analysis "
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "4c68f032-e621-4b91-b3f8-1ec2882e01ec",
"metadata": {},
"outputs": [],
"source": [
"dir_SV=$dir_genome\"SYRI/\"\n",
"mkdir -p $dir_SV\n",
"cd $dir_SV"
]
},
{
"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",
"\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",
""
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "bac7fd8e-c4f9-4c6c-a156-d3e94934e485",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1: PREPARING DATA\n",
"2,3: RUNNING mummer AND CREATING CLUSTERS\n",
"# reading input file \"out.ntref\" of length 15306200\n",
"# construct suffix tree for sequence of length 15306200\n",
"# (maximum reference length is 2305843009213693948)\n",
"# (maximum query length is 18446744073709551615)\n",
"# process 153062 characters per dot\n",
"#...................................................................................................\n",
"# CONSTRUCTIONTIME /opt/conda/opt/mummer-3.23/mummer out.ntref 10.60\n",
"# reading input file \"/home/jovyan/work/day4-WGA-SV/Assemblies/5417_assembly.fasta\" of length 15276944\n",
"# matching query-file \"/home/jovyan/work/day4-WGA-SV/Assemblies/5417_assembly.fasta\"\n",
"# against subject-file \"out.ntref\"\n",
"# COMPLETETIME /opt/conda/opt/mummer-3.23/mummer out.ntref 45.30\n",
"# SPACE /opt/conda/opt/mummer-3.23/mummer out.ntref 29.56\n",
"4: FINISHING DATA\n"
]
}
],
"source": [
"nucmer --maxmatch $reference_assembly $query_assembly"
]
},
{
"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": 34,
"id": "ee2ba7af-8b94-4df3-8412-928d6b421234",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 1788\n",
"-rw-r--r-- 1 jovyan users 1829452 Jun 23 21:38 out.delta\n"
]
}
],
"source": [
"ls -lrt"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "fea08d51-a81a-47b7-a31d-382a0c2d95e0",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/jovyan/work/day4-WGA-SV/Assemblies/A8_assembly.fasta /home/jovyan/work/day4-WGA-SV/Assemblies/5417_assembly.fasta\n",
"NUCMER\n",
">Chr_18 Chr_19 308071 233066\n",
"1 110 2 111 0 0 0\n",
"0\n",
"1 110 23 132 0 0 0\n",
"0\n",
"1 110 233060 232951 0 0 0\n",
"0\n",
"1 110 233053 232944 0 0 0\n"
]
}
],
"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",
"\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",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 42,
"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": 43,
"id": "467c1f2e-62d3-477b-b549-6de9c0f46619",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 2440\n",
"-rw-r--r-- 1 jovyan users 1829452 Jun 23 21:38 out.delta\n",
"-rw-r--r-- 1 jovyan users 667061 Jun 23 21:40 out.filtered.delta\n"
]
}
],
"source": [
"ls -lrt"
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "b46169b9-b8e6-4c32-96a9-17e650a57bd5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 443088 out.delta\n",
" 218555 out.filtered.delta\n",
" 661643 total\n"
]
}
],
"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",
"\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",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 48,
"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": "e403ff2d-2098-40c8-8cce-2086988a16b5",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 45,
"id": "f9eb2b0f-f716-413e-ae14-76e70541e011",
"metadata": {},
"outputs": [],
"source": [
"show-coords -THrd out.filtered.delta > out.filtered.coords"
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "f34fceb2-a76e-41be-8a69-cf59e97564aa",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 2540\n",
"-rw-r--r-- 1 jovyan users 1829452 Jun 23 21:38 out.delta\n",
"-rw-r--r-- 1 jovyan users 667061 Jun 23 21:40 out.filtered.delta\n",
"-rw-r--r-- 1 jovyan users 100931 Jun 23 21:42 out.filtered.coords\n"
]
}
],
"source": [
"ls -lrt"
]
},
{
"cell_type": "code",
"execution_count": 50,
"id": "929dc4d1-1a08-43ac-b7a2-82880b9b807b",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==> out.filtered.coords <==\n",
"1\t122\t482889\t482768\t122\t122\t100.00\t1\t-1\tChr_1\tChr_17\n",
"394\t689\t228410\t228116\t296\t295\t92.23\t1\t-1\tChr_1\tChr_19\n",
"484\t608\t220395\t220272\t125\t124\t91.20\t1\t-1\tChr_1\tChr_19\n",
"495\t689\t226164\t225970\t195\t195\t92.82\t1\t-1\tChr_1\tChr_19\n",
"799\t23710\t1041\t23963\t22912\t22923\t99.41\t1\t1\tChr_1\tChr_1\n",
"6978\t7128\t1039\t886\t151\t154\t97.40\t1\t-1\tChr_1\tChr_1\n",
"23711\t45447\t24062\t45783\t21737\t21722\t98.05\t1\t1\tChr_1\tChr_1\n",
"26393\t26528\t451998\t451864\t136\t135\t97.06\t1\t-1\tChr_1\tChr_3\n",
"36290\t36897\t973523\t974129\t608\t607\t99.51\t1\t1\tChr_1\tChr_8\n",
"45590\t49328\t45931\t49682\t3739\t3752\t96.80\t1\t1\tChr_1\tChr_1\n",
"\n",
"==> out.filtered.withHeader.coords <==\n",
"/home/jovyan/work/day4-WGA-SV/Assemblies/A8_assembly.fasta /home/jovyan/work/day4-WGA-SV/Assemblies/5417_assembly.fasta\n",
"NUCMER\n",
"\n",
"[S1]\t[E1]\t[S2]\t[E2]\t[LEN 1]\t[LEN 2]\t[% IDY]\t[FRM]\t[TAGS]\n",
"1\t122\t482889\t482768\t122\t122\t100.00\t1\t-1\tChr_1\tChr_17\n",
"394\t689\t228410\t228116\t296\t295\t92.23\t1\t-1\tChr_1\tChr_19\n",
"484\t608\t220395\t220272\t125\t124\t91.20\t1\t-1\tChr_1\tChr_19\n",
"495\t689\t226164\t225970\t195\t195\t92.82\t1\t-1\tChr_1\tChr_19\n",
"799\t23710\t1041\t23963\t22912\t22923\t99.41\t1\t1\tChr_1\tChr_1\n",
"6978\t7128\t1039\t886\t151\t154\t97.40\t1\t-1\tChr_1\tChr_1\n"
]
}
],
"source": [
"head *.coords"
]
},
{
"cell_type": "markdown",
"id": "0c84c259-cdc7-422e-9c7f-cc09838418c9",
"metadata": {},
"source": [
"### SV calling using syri \n",
"\n",
"__activate syri envt with conda__"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "fda3a8b2-7768-4775-bbff-ba67b551d5ea",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"source": [
"conda activate syri_env"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "4dbf9a7f-d961-4791-904a-8dd12b986617",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Finished syri\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"source": [
"syri -c out.filtered.coords -d out.filtered.delta -r $reference_assembly -q $query_assembly"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "a3f88d52-7edc-4d2c-a652-9ce1ac5531b9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 59156\n",
"-rw-r--r-- 1 jovyan users 556 Jun 23 21:49 syri.summary\n",
"-rw-r--r-- 1 jovyan users 36876017 Jun 23 21:49 syri.vcf\n",
"-rw-r--r-- 1 jovyan users 8917 Jun 23 21:49 syri.log\n",
"-rw-r--r-- 1 jovyan users 20978152 Jun 23 21:49 syri.out\n",
"-rw-r--r-- 1 jovyan users 101116 Jun 23 21:44 out.filtered.withHeader.coords\n",
"-rw-r--r-- 1 jovyan users 100931 Jun 23 21:42 out.filtered.coords\n",
"-rw-r--r-- 1 jovyan users 667061 Jun 23 21:40 out.filtered.delta\n",
"-rw-r--r-- 1 jovyan users 1829452 Jun 23 21:38 out.delta\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"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",
"\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",
""
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "0fdcb7c6-d698-482c-b2a1-6d8f5f2b60e2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"#Structural annotations\n",
"#Variation_type\tCount\tLength_ref\tLength_qry\n",
"Syntenic regions\t197\t14079055\t14107736\n",
"Inversions\t32\t397553\t259052\n",
"Translocations\t93\t91312\t90081\n",
"Duplications (reference)\t152\t79812\t-\n",
"Duplications (query)\t161\t-\t83493\n",
"Not aligned (reference)\t370\t667756\t-\n",
"Not aligned (query)\t413\t-\t738403\n",
"\n",
"\n",
"#Sequence annotations\n",
"#Variation_type\tCount\tLength_ref\tLength_qry\n",
"SNPs\t239492\t239492\t239492\n",
"Insertions\t31092\t-\t114022\n",
"Deletions\t38036\t119894\t-\n",
"Copygains\t32\t-\t13466\n",
"Copylosses\t28\t17778\t-\n",
"Highly diverged\t707\t487347\t383289\n",
"Tandem repeats\t23\t8965\t10360\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"source": [
"cat syri.summary"
]
},
{
"cell_type": "markdown",
"id": "6578624d-3ee2-47f8-9943-6ff10704e303",
"metadata": {},
"source": [
"__TSV format specifications__\n",
"\n",
"\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",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6de3d0f9-2198-4c0f-b813-6b4d3e41ab3c",
"metadata": {},
"outputs": [],
"source": [
"head syri.vcf"
]
},
{
"cell_type": "code",
"execution_count": 64,
"id": "a9435c82-45c7-497e-aa9c-e49c336c2925",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Chr_1\t1\t122\t-\t-\tChr_17\t482768\t482889\tINVDP340\t-\tINVDP\tcopyloss\n",
"Chr_1\t1\t122\t-\t-\tChr_17\t482889\t482768\tINVDPAL1183\tINVDP340\tINVDPAL\t-\n",
"Chr_1\t123\t483\t-\t-\t-\t-\t-\tNOTAL1\t-\tNOTAL\t-\n",
"Chr_1\t484\t608\t-\t-\tChr_19\t220272\t220395\tINVDP341\t-\tINVDP\tcopygain\n",
"Chr_1\t484\t608\t-\t-\tChr_19\t220395\t220272\tINVDPAL1184\tINVDP341\tINVDPAL\t-\n",
"Chr_1\t495\t689\t-\t-\tChr_19\t225970\t226164\tINVTR342\t-\tINVTR\t-\n",
"Chr_1\t495\t689\t-\t-\tChr_19\t226164\t225970\tINVTRAL1185\tINVTR342\tINVTRAL\t-\n",
"Chr_1\t521\t521\tT\tC\tChr_19\t226138\t226138\tSNP839\tINVTR342\tSNP\t-\n",
"Chr_1\t522\t522\tA\tT\tChr_19\t226137\t226137\tSNP840\tINVTR342\tSNP\t-\n",
"Chr_1\t530\t530\tC\tG\tChr_19\t226129\t226129\tSNP841\tINVTR342\tSNP\t-\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"source": [
"head syri.out"
]
},
{
"cell_type": "markdown",
"id": "b6825e0b-9260-4fe7-be3e-c410f7fd92ba",
"metadata": {},
"source": [
"### Extracting all SNP from syri "
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "e3ecf3d7-a4b2-4d7b-85fa-a32b0ee60ede",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(syri_env) (syri_env) (syri_env) (syri_env) (syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"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": 83,
"id": "79168401-fac7-453e-9962-b2407ed4e72e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Chr_1\t521\t522\t521\n",
"Chr_1\t522\t523\t522\n",
"Chr_1\t530\t531\t530\n",
"Chr_1\t531\t532\t531\n",
"Chr_1\t534\t535\t534\n",
"Chr_1\t535\t536\t535\n",
"Chr_1\t536\t537\t536\n",
"Chr_1\t548\t549\t548\n",
"Chr_1\t555\t556\t555\n",
"Chr_1\t570\t571\t570\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"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": 84,
"id": "9ae51f46-f320-43ba-a20b-81f93d569c2d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(syri_env) (syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"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": 85,
"id": "af2dcd43-0ea9-4577-9a76-1899f3730d07",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/jovyan/work/day4-WGA-SV/Assemblies/A8_assembly.fasta\t ref\n",
"/home/jovyan/work/day4-WGA-SV/Assemblies/5417_assembly.fasta\t query\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"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": 86,
"id": "96ffd3a6-6bf7-41da-a78c-84219fe96878",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/jovyan/work/day4-WGA-SV/SYRI\n",
"(syri_env) Plotsr - WARNING - Plot plotsr.pdf generated.\n",
"Finished\n",
"(syri_env) "
]
},
{
"ename": "",
"evalue": "1",
"output_type": "error",
"traceback": []
}
],
"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": [
"## __II - using assemblytics__ \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": [
"## __III - using d-genies__ \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
}