{
"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": 1,
"id": "1b75773d-38a6-49ab-9de9-34681cfd44f5",
"metadata": {},
"outputs": [],
"source": [
"cd /home/jovyan/SV_DATA/ASSEMBLIES"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "adf7ebcd-112f-4849-a017-db793c2696be",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==> 5417_assembly.fasta <==\n",
">Chr_19\n",
"CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCT\n",
"AAACCCTAAACCCTAAACCCTAAACCCTAAACCCAGAAAAATAGACTTTATTCCAACATAACTATCTACATCATCGTACGTACCTGCAGGTATGTGTTGA\n",
"AAGGTACCCGTTACACAAAAGAACTGAAACATTTGTCTAACACGTCAAATTCTGCTGCATTTGCCGTTATAAAACGTCTCAAAACTTTCGTTGAAAATAA\n",
"\n",
"==> A8_assembly.fasta <==\n",
">Chr_19\n",
"AAAAAATTTGGTCGTTTGAGCGTATATTTGCATACACCCTAAACAGAGCGTGCAAATATACGCTAAAAACGACACAAAATTTTTTTTTTTTTTTGCGAAA\n",
"GCACTCTAGTCGATCACGTTGTCTGAAAAAATTTAAATCGAAGTTTATAAGGCCAAGCTTGGTCAAAGGTCATACGTCATTATTTATACGTCACTATACG\n",
"TTTCTTCATACTTCGTTCCTTCACACAAAACGTAAAGATCTGAAATAGACGGATCCTTTGAGAGTTCATTGGAGAGATCTTTGAAAGTTGTTCTTTTTGT\n"
]
}
],
"source": [
"head -n 4 *fasta"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "259b9c77-4326-4daf-b03b-ef669c4e2148",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[35m\u001b[K5417_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K19\n",
"\u001b[35m\u001b[KA8_assembly.fasta\u001b[m\u001b[K\u001b[36m\u001b[K:\u001b[m\u001b[K19\n"
]
}
],
"source": [
"grep \">\" *fasta -c"
]
},
{
"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": 7,
"id": "ccb2453d-02fa-4cec-915a-9136e27647e7",
"metadata": {
"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 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": 11,
"id": "b79e2a4a-a505-4ed1-9763-9c34e10271e2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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",
"Total pb : 15276926\n"
]
}
],
"source": [
"seqtk comp 5417_assembly.fasta | awk '{pb=pb+$2; print($2,pb)} END {print \"Total pb :\",pb}'"
]
},
{
"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": []
},
{
"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": null,
"id": "bac7fd8e-c4f9-4c6c-a156-d3e94934e485",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"nucmer --maxmatch REFERENCE_GENOME ASSEMBLIED_GENOME"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "06727fa9-3b4a-4b5d-9dd6-36bd62a9f18e",
"metadata": {},
"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 15074341\n",
"# construct suffix tree for sequence of length 15074341\n",
"# (maximum reference length is 536870908)\n",
"# (maximum query length is 4294967295)\n",
"# process 150743 characters per dot\n",
"#....................................................................................................\n",
"# CONSTRUCTIONTIME /usr/bin/mummer out.ntref 6.96\n",
"# reading input file \"/home/jovyan/SV_DATA/ASSEMBLIES/5417_assembly.fasta\" of length 15276944\n",
"# matching query-file \"/home/jovyan/SV_DATA/ASSEMBLIES/5417_assembly.fasta\"\n",
"# against subject-file \"out.ntref\"\n",
"# COMPLETETIME /usr/bin/mummer out.ntref 32.67\n",
"# SPACE /usr/bin/mummer out.ntref 29.33\n",
"4: FINISHING DATA\n"
]
}
],
"source": [
"nucmer --maxmatch /home/jovyan/SV_DATA/REF/GCA_002220235.1_ASM222023v1_genomic.fna 5417_assembly.fasta"
]
},
{
"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": 13,
"id": "ee2ba7af-8b94-4df3-8412-928d6b421234",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 31956\n",
"-rwxr-xr-x 1 jovyan users 552 Jun 17 2022 assembly-stats.txt\n",
"-rwxr-xr-x 1 jovyan users 15429848 Jun 17 2022 5417_assembly.fasta\n",
"-rwxr-xr-x 1 jovyan users 15459392 Jun 17 2022 A8_assembly.fasta\n",
"-rw-r--r-- 1 jovyan users 1819920 Oct 9 08:54 out.delta\n"
]
}
],
"source": [
"ls -lrt"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "fea08d51-a81a-47b7-a31d-382a0c2d95e0",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/jovyan/SV_DATA/REF/GCA_002220235.1_ASM222023v1_genomic.fna /home/jovyan/SV_DATA/ASSEMBLIES/5417_assembly.fasta\n",
"NUCMER\n",
">FO082277.1 Chr_19 1122692 233066\n",
"7 140 1 134 0 0 0\n",
"0\n",
"1 130 233066 232938 1 1 0\n",
"6\n",
"0\n",
"9 137 233066 232938 0 0 0\n",
"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": 15,
"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": 16,
"id": "467c1f2e-62d3-477b-b549-6de9c0f46619",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 32616\n",
"-rwxr-xr-x 1 jovyan users 552 Jun 17 2022 assembly-stats.txt\n",
"-rwxr-xr-x 1 jovyan users 15429848 Jun 17 2022 5417_assembly.fasta\n",
"-rwxr-xr-x 1 jovyan users 15459392 Jun 17 2022 A8_assembly.fasta\n",
"-rw-r--r-- 1 jovyan users 1819920 Oct 9 08:54 out.delta\n",
"-rw-r--r-- 1 jovyan users 674709 Oct 9 08:54 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": 17,
"id": "3c2b20c4-665a-4cc5-8240-e191edd66582",
"metadata": {},
"outputs": [],
"source": [
"show-coords -Trd out.filtered.delta > out.filtered.withHeader.coords"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "f9eb2b0f-f716-413e-ae14-76e70541e011",
"metadata": {},
"outputs": [],
"source": [
"show-coords -THrd out.filtered.delta > out.filtered.coords"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "f34fceb2-a76e-41be-8a69-cf59e97564aa",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 32832\n",
"-rwxr-xr-x 1 jovyan users 552 Jun 17 2022 assembly-stats.txt\n",
"-rwxr-xr-x 1 jovyan users 15429848 Jun 17 2022 5417_assembly.fasta\n",
"-rwxr-xr-x 1 jovyan users 15459392 Jun 17 2022 A8_assembly.fasta\n",
"-rw-r--r-- 1 jovyan users 1819920 Oct 9 08:54 out.delta\n",
"-rw-r--r-- 1 jovyan users 674709 Oct 9 08:54 out.filtered.delta\n",
"-rw-r--r-- 1 jovyan users 108563 Oct 9 08:54 out.filtered.withHeader.coords\n",
"-rw-r--r-- 1 jovyan users 108381 Oct 9 08:54 out.filtered.coords\n"
]
}
],
"source": [
"ls -lrt"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "929dc4d1-1a08-43ac-b7a2-82880b9b807b",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==> out.filtered.coords <==\n",
"34\t342\t87624\t87310\t309\t315\t90.91\t1\t-1\tFO082260.1\tChr_17\n",
"517\t813\t2148\t2443\t297\t296\t93.62\t1\t1\tFO082260.1\tChr_17\n",
"2799\t3119\t540774\t541090\t321\t317\t90.97\t1\t1\tFO082260.1\tChr_14\n",
"9808\t10021\t230845\t230632\t214\t214\t91.67\t1\t-1\tFO082260.1\tChr_1\n",
"10058\t10579\t8344\t8866\t522\t523\t95.62\t1\t1\tFO082260.1\tChr_19\n",
"10064\t10200\t8866\t8730\t137\t137\t94.89\t1\t-1\tFO082260.1\tChr_19\n",
"10367\t10577\t147473\t147262\t211\t212\t93.40\t1\t-1\tFO082260.1\tChr_19\n",
"10443\t10585\t8488\t8344\t143\t145\t93.79\t1\t-1\tFO082260.1\tChr_19\n",
"10443\t10572\t11544\t11672\t130\t129\t97.69\t1\t1\tFO082260.1\tChr_19\n",
"10448\t10594\t212302\t212451\t147\t150\t92.05\t1\t1\tFO082260.1\tChr_19\n",
"\n",
"==> out.filtered.withHeader.coords <==\n",
"/home/jovyan/SV_DATA/REF/GCA_002220235.1_ASM222023v1_genomic.fna /home/jovyan/SV_DATA/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",
"34\t342\t87624\t87310\t309\t315\t90.91\t1\t-1\tFO082260.1\tChr_17\n",
"517\t813\t2148\t2443\t297\t296\t93.62\t1\t1\tFO082260.1\tChr_17\n",
"2799\t3119\t540774\t541090\t321\t317\t90.97\t1\t1\tFO082260.1\tChr_14\n",
"9808\t10021\t230845\t230632\t214\t214\t91.67\t1\t-1\tFO082260.1\tChr_1\n",
"10058\t10579\t8344\t8866\t522\t523\t95.62\t1\t1\tFO082260.1\tChr_19\n",
"10064\t10200\t8866\t8730\t137\t137\t94.89\t1\t-1\tFO082260.1\tChr_19\n"
]
}
],
"source": [
"head *.coords"
]
},
{
"cell_type": "markdown",
"id": "0c84c259-cdc7-422e-9c7f-cc09838418c9",
"metadata": {},
"source": [
"### I - SV calling using syri \n",
"\n",
"__activate syri envt with conda__"
]
},
{
"cell_type": "code",
"execution_count": 21,
"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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Reading Coords - WARNING - Chromosomes IDs do not match.\n",
"Reading Coords - WARNING - Matching them automatically. For each reference genome, most similar query genome will be selected. Check mapids.txt for mapping used.\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"/opt/conda/envs/syri_env/lib/python3.9/multiprocessing/pool.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
" return list(map(*args))\n",
"local_variation - WARNING - Finished syri\n"
]
}
],
"source": [
"syri -c out.filtered.coords -d /home/jovyan/SV_DATA/ASSEMBLIES/out.filtered.delta -r /home/jovyan/SV_DATA/REF/GCA_002220235.1_ASM222023v1_genomic.fna -q /home/jovyan/SV_DATA/ASSEMBLIES/5417_assembly.fasta"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a3f88d52-7edc-4d2c-a652-9ce1ac5531b9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 93236\n",
"-rw-r--r-- 1 jovyan users 11531 Oct 9 08:57 syri.log\n",
"-rw-r--r-- 1 jovyan users 559 Oct 9 08:57 syri.summary\n",
"-rw-r--r-- 1 jovyan users 38982015 Oct 9 08:57 syri.vcf\n",
"-rw-r--r-- 1 jovyan users 22845143 Oct 9 08:57 syri.out\n",
"-rw-r--r-- 1 jovyan users 333 Oct 9 08:57 mapids.txt\n",
"-rw-r--r-- 1 jovyan users 108381 Oct 9 08:54 out.filtered.coords\n",
"-rw-r--r-- 1 jovyan users 108563 Oct 9 08:54 out.filtered.withHeader.coords\n",
"-rw-r--r-- 1 jovyan users 674709 Oct 9 08:54 out.filtered.delta\n",
"-rw-r--r-- 1 jovyan users 1819920 Oct 9 08:54 out.delta\n",
"-rwxr-xr-x 1 jovyan users 15459392 Jun 17 2022 A8_assembly.fasta\n",
"-rwxr-xr-x 1 jovyan users 15429848 Jun 17 2022 5417_assembly.fasta\n",
"-rwxr-xr-x 1 jovyan users 552 Jun 17 2022 assembly-stats.txt\n"
]
}
],
"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": null,
"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\t204\t14337208\t14405320\n",
"Inversions\t29\t182110\t178809\n",
"Translocations\t103\t97810\t101175\n",
"Duplications (reference)\t148\t130887\t-\n",
"Duplications (query)\t163\t-\t76965\n",
"Not aligned (reference)\t370\t273141\t-\n",
"Not aligned (query)\t429\t-\t501845\n",
"\n",
"\n",
"#Sequence annotations\n",
"#Variation_type\tCount\tLength_ref\tLength_qry\n",
"SNPs\t243003\t243003\t243003\n",
"Insertions\t31254\t-\t118867\n",
"Deletions\t38997\t118875\t-\n",
"Copygains\t34\t-\t35610\n",
"Copylosses\t25\t6873\t-\n",
"Highly diverged\t730\t298921\t337650\n",
"Tandem repeats\t24\t15047\t18749\n"
]
}
],
"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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##fileformat=VCFv4.3\n",
"##fileDate=20231009\n",
"##source=syri\n",
"##ALT=\n",
"##ALT=\n",
"##ALT=\n",
"##ALT=\n",
"##ALT=\n",
"##ALT=\n",
"##ALT=\n"
]
}
],
"source": [
"head syri.vcf"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a9435c82-45c7-497e-aa9c-e49c336c2925",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"FO082260.1\t1\t33\t-\t-\t-\t-\t-\tNOTAL1\t-\tNOTAL\t-\n",
"FO082260.1\t34\t342\t-\t-\tChr_17\t87310\t87624\tINVTR344\t-\tINVTR\t-\n",
"FO082260.1\t34\t342\t-\t-\tChr_17\t87624\t87310\tINVTRAL1205\tINVTR344\tINVTRAL\t-\n",
"FO082260.1\t53\t53\tA\tG\tChr_17\t87605\t87605\tSNP858\tINVTR344\tSNP\t-\n",
"FO082260.1\t59\t59\tA\tG\tChr_17\t87599\t87599\tSNP859\tINVTR344\tSNP\t-\n",
"FO082260.1\t68\t68\tT\tC\tChr_17\t87590\t87590\tSNP860\tINVTR344\tSNP\t-\n",
"FO082260.1\t163\t166\tgGAG\tg\tChr_17\t87494\t87494\tDEL861\tINVTR344\tDEL\t-\n",
"FO082260.1\t208\t208\tT\tA\tChr_17\t87453\t87453\tSNP862\tINVTR344\tSNP\t-\n",
"FO082260.1\t209\t209\tA\tG\tChr_17\t87452\t87452\tSNP863\tINVTR344\tSNP\t-\n",
"FO082260.1\t210\t210\tC\tT\tChr_17\t87451\t87451\tSNP864\tINVTR344\tSNP\t-\n"
]
}
],
"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": 37,
"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": 38,
"id": "79168401-fac7-453e-9962-b2407ed4e72e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"FO082260.1\t53\t54\t53\n",
"FO082260.1\t59\t60\t59\n",
"FO082260.1\t68\t69\t68\n",
"FO082260.1\t208\t209\t208\n",
"FO082260.1\t209\t210\t209\n",
"FO082260.1\t210\t211\t210\n",
"FO082260.1\t221\t222\t221\n",
"FO082260.1\t225\t226\t225\n",
"FO082260.1\t226\t227\t226\n",
"FO082260.1\t228\t229\t228\n"
]
}
],
"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": 39,
"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": 40,
"id": "af2dcd43-0ea9-4577-9a76-1899f3730d07",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\t ref\n",
"\t query\n",
"\t ref\n",
"\t query\n"
]
}
],
"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": 41,
"id": "96ffd3a6-6bf7-41da-a78c-84219fe96878",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/jovyan/SV_DATA/ASSEMBLIES\n",
"bash: plotsr: command not found\n"
]
},
{
"ename": "",
"evalue": "127",
"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": 42,
"id": "ca6b926f-aaad-4dca-ae45-6a68f0aa03ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"usage: syri [-h] -c INFILE [-r REF] [-q QRY] [-d DELTA] [-F {T,S,B,P}] [-f]\n",
" [-k] [--dir DIR] [--prefix PREFIX] [--seed SEED] [--nc NCORES]\n",
" [--novcf] [--nosr] [--invgaplen INVGL] [--tdgaplen TDGL]\n",
" [--tdmaxolp TDOLP] [-b BRUTERUNTIME] [--unic TRANSUNICOUNT]\n",
" [--unip TRANSUNIPERCENT] [--inc INCREASEBY] [--no-chrmatch]\n",
" [--nosv] [--nosnp] [--all] [--allow-offset OFFSET] [--cigar]\n",
" [-s SSPATH] [--hdrseq] [--maxsize MAXS] [--log {DEBUG,INFO,WARN}]\n",
" [--lf LOG_FIN] [--version]\n",
"\n",
"Input Files:\n",
" -c INFILE File containing alignment coordinates (default: None)\n",
" -r REF Genome A (which is considered as reference for the\n",
" alignments). Required for local variation (large\n",
" indels, CNVs) identification. (default: None)\n",
" -q QRY Genome B (which is considered as query for the\n",
" alignments). Required for local variation (large\n",
" indels, CNVs) identification. (default: None)\n",
" -d DELTA .delta file from mummer. Required for short variation\n",
" (SNPs/indels) identification when CIGAR string is not\n",
" available (default: None)\n",
"\n",
"Additional arguments:\n",
" -F {T,S,B,P} Input file type. T: Table, S: SAM, B: BAM, P: PAF\n",
" (default: T)\n",
" -f As a default, syri filters out low quality and small\n",
" alignments. Use this parameter to use the full list of\n",
" alignments without any filtering. (default: True)\n",
" -k Keep intermediate output files (default: False)\n",
" --dir DIR path to working directory (if not current directory).\n",
" All files must be in this directory. (default: None)\n",
" --prefix PREFIX Prefix to add before the output file Names (default: )\n",
" --seed SEED seed for generating random numbers (default: 1)\n",
" --nc NCORES number of cores to use in parallel (max is number of\n",
" chromosomes) (default: 1)\n",
" --novcf Do not combine all files into one output file\n",
" (default: False)\n",
"\n",
"SR identification:\n",
" --nosr Set to skip structural rearrangement identification\n",
" (default: False)\n",
" --invgaplen INVGL Maximum allowed gap-length between two alignments of a\n",
" multi-alignment inversion. It affects the selection of\n",
" large inversions that can have different length in the\n",
" reference and query genomes. (default: 1000000000)\n",
" --tdgaplen TDGL Maximum allowed gap-length between two alignments of a\n",
" multi-alignment translocation or duplication (TD).\n",
" Larger values increases TD identification sensitivity\n",
" but also runtime. (default: 500000)\n",
" --tdmaxolp TDOLP Maximum allowed overlap between two translocations.\n",
" Value should be in range (0,1]. (default: 0.8)\n",
" -b BRUTERUNTIME Cutoff to restrict brute force methods to take too\n",
" much time (in seconds). Smaller values would make\n",
" algorithm faster, but could have marginal effects on\n",
" accuracy. In general case, would not be required.\n",
" (default: 60)\n",
" --unic TRANSUNICOUNT Number of uniques bps for selecting translocation.\n",
" Smaller values would select smaller TLs better, but\n",
" may increase time and decrease accuracy. (default:\n",
" 1000)\n",
" --unip TRANSUNIPERCENT\n",
" Percent of unique region requried to select\n",
" translocation. Value should be in range (0,1]. Smaller\n",
" values would allow selection of TDs which are more\n",
" overlapped with other regions. (default: 0.5)\n",
" --inc INCREASEBY Minimum score increase required to add another\n",
" alignment to translocation cluster solution (default:\n",
" 1000)\n",
" --no-chrmatch Do not allow SyRI to automatically match chromosome\n",
" ids between the two genomes if they are not equal\n",
" (default: False)\n",
"\n",
"ShV identification:\n",
" --nosv Set to skip structural variation identification\n",
" (default: False)\n",
" --nosnp Set to skip SNP/Indel (within alignment)\n",
" identification (default: False)\n",
" --all Use duplications too for variant identification\n",
" (default: False)\n",
" --allow-offset OFFSET\n",
" BPs allowed to overlap (default: 5)\n",
" --cigar Find SNPs/indels using CIGAR string. Necessary for\n",
" alignments generated using aligners other than nucmers\n",
" (default: False)\n",
" -s SSPATH path to show-snps from mummer (default: show-snps)\n",
" --hdrseq Output highly-diverged regions (HDRs) sequence.\n",
" (default: False)\n",
" --maxsize MAXS Max size for printing sequence of large SVs\n",
" (insertions, deletions and HDRs). Only affect printing\n",
" (.out/.vcf file) and not the selection. SVs larger\n",
" than this value would be printed as symbolic SVs. For\n",
" no cut-off use -1. (default: -1)\n",
"\n",
"optional arguments:\n",
" -h, --help show this help message and exit\n",
" --log {DEBUG,INFO,WARN}\n",
" log level (default: INFO)\n",
" --lf LOG_FIN Name of log file (default: syri.log)\n",
" --version show program's version number and exit\n"
]
}
],
"source": [
"syri --help"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f00a85ca-7e0a-4681-9ac0-e7cb4e6b2488",
"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
}