{
"cells": [
{
"cell_type": "markdown",
"id": "2f09a11c",
"metadata": {},
"source": [
"# Formation South Green 2022 - Structural Variants Detection by using short and long reads \n",
"\n",
"# __DAY 3 : Structural variant calling__\n",
"\n",
"Created by C. Tranchant (DIADE-IRD), J. Orjuela (DIADE-IRD), F. Sabot (DIADE-IRD) and A. Dereeper (PHIM-IRD)\n",
"\n",
"\n",
"***\n",
"\n",
"# Table of contents\n",
"\n",
"\n",
"\n",
"\n",
"[I - Structural variation with `Breakdancer` (SR)__ ](#break)\n",
"\n",
" * [ Reference indexation](#refindex)\n",
"\n",
"\n",
"[II - Structural variation with Sniffles (LR)](#sniffles) \n",
" \n",
"\n",
" * [1. Prepare LR data](#step1)\n",
" * [2. SV detection for all CLONES](#step2)\n",
" * [Count the number of variations](#countv)\n",
" * [Convert sam into bam `samtools view`](#samtoolsview)\n",
" * [Merge all the vcf files across all samples](#merge)\n",
" * [3. Force call all the SVs across all the samples](#step3)\n",
"\n",
"\n",
"\n",
"***\n"
]
},
{
"cell_type": "markdown",
"id": "c635381f-17de-41d5-abfa-60d81bca626b",
"metadata": {},
"source": [
"# __I - Structural variation with `Breakdancer` (SR)__ \n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "de5be5a8-28e1-420a-b9da-dad8dd5121c5",
"metadata": {},
"source": [
"To run BreakDancer, first use `bam2cfg.pl` to prepare the required config file."
]
},
{
"cell_type": "markdown",
"id": "aeb1cfe2-4833-45f8-bd91-20ea82b963c4",
"metadata": {},
"source": [
"### First of all, create a list of bam files "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "84017eae-ff41-47d6-beff-970306e6ee7f",
"metadata": {},
"outputs": [],
"source": [
"mkdir -p ~/work/BREAKDANCER/;\n",
"cd ~/work/BREAKDANCER/;\n",
"realpath /home/jovyan/work/MAPPING-ILL/dirClone*/*SORTED.bam > bam_files.txt"
]
},
{
"cell_type": "markdown",
"id": "09bcc14f-e73e-4161-8609-2955f36cbcbc",
"metadata": {},
"source": [
"ls"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "18ae3020-5b30-48ef-bea5-81b5653950e1",
"metadata": {},
"outputs": [],
"source": [
"cat bam_files.txt"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "50312ee1-e8d1-45d3-a61d-f08388213ba7",
"metadata": {},
"outputs": [],
"source": [
"bam2cfg.pl bam_files.txt \n",
"#breakdancer_options"
]
},
{
"cell_type": "markdown",
"id": "c9998362-7cd2-4db3-8e06-0163f83c1d7a",
"metadata": {},
"source": [
"### Then run BreakDancer on the config "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f3750540-9703-4897-9227-6130b668bd3e",
"metadata": {},
"outputs": [],
"source": [
"breakdancer-max config_file.cfg"
]
},
{
"cell_type": "markdown",
"id": "11682ed5-7f3f-4fac-b79a-0c4b46afda42",
"metadata": {},
"source": []
},
{
"cell_type": "markdown",
"id": "3d07c5cc",
"metadata": {
"tags": []
},
"source": [
"# __II - Structural variation with `Sniffles` (LR)__ \n",
"\n",
"Sniffles is a structural variation caller using third generation sequencing (PacBio or Oxford Nanopore).\n",
"\n",
"It detects all types of SVs (10bp+) using evidence from split-read alignments, high-mismatch regions, and coverage analysis."
]
},
{
"cell_type": "markdown",
"id": "f8c89ab9-5da9-4c48-9a20-083bbc3f4760",
"metadata": {},
"source": [
"# 1. Prepare LR data "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13474807-c0f2-4d71-ac85-aae834302073",
"metadata": {},
"outputs": [],
"source": [
"# create SNIFFLES folder\n",
"mkdir -p /home/jovyan/work/SNIFFLES/\n",
"cd /home/jovyan/work/SNIFFLES/\n",
"\n",
"# declare your Clone\n",
"CLONE=\"Clone10\"\n",
"\n",
"# declare reference path\n",
"REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\""
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e8800be2-1812-43b9-9e43-77aaadbbed3a",
"metadata": {},
"outputs": [],
"source": [
"# check reads ONT from Clones\n",
"ONT=/home/jovyan/work/SV_DATA/LONG_READS\n",
"ls $ONT\n"
]
},
{
"cell_type": "markdown",
"id": "9e172e4a-4c5b-4c84-bdfe-1c13c2d586bd",
"metadata": {
"tags": []
},
"source": [
"# 2. SV detection for all CLONES "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ab00d43f-029c-4d9f-bf74-71815def6b3d",
"metadata": {},
"outputs": [],
"source": [
"cd /home/jovyan/work/MAPPING-ONT/BAM"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8dc93fc-8ae1-4a78-a1a3-00e01b5cfad7",
"metadata": {},
"outputs": [],
"source": [
"ls"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6cfc6884-3c04-49dd-a51a-17abca4cdedc",
"metadata": {},
"outputs": [],
"source": [
"#wget --no-check-certificat -rm -nH --cut-dirs=1 --reject=\"index.html*\" wget https://itrop.ird.fr/sv-training/BAM_ONT.tar.gz\n",
"#tar zxvf BAM_ONT.tar.gz\n",
"#BAM_ONT=\"/home/jovyan/work/BAM_ONT\"\n",
"#rm BAM_ONT.tar.gz\n",
"#ls $BAM_ONT"
]
},
{
"cell_type": "markdown",
"id": "fbe8c131-5f5e-4156-8817-37a0c1c343ab",
"metadata": {},
"source": [
"## Run `sniffles` "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ea080741-d309-483d-824b-7ccd75b6ac64",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"cd /home/jovyan/work/MAPPING-ONT/BAM/\n",
"\n",
"for i in {2,4,6,8,10,12,14,16,18,20}\n",
" do\n",
" echo \"============ Clone${i}==============\";\n",
" samtools index Clone${i}_ONT_SORTED.bam\n",
" sniffles -i Clone${i}_ONT_SORTED.bam -v /home/jovyan/work/SNIFFLES/Clone${i}_SV.vcf --allow-overwrite\n",
" done\n",
"time\n"
]
},
{
"cell_type": "markdown",
"id": "8fb4bc4d-1841-4aa7-906b-b6ad3548e03d",
"metadata": {},
"source": [
"## Count the number of variations "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "367a047d-cef9-487e-b496-be2a76554131",
"metadata": {},
"outputs": [],
"source": [
"cd /home/jovyan/work/SNIFFLES/\n",
"\n",
"for i in {2,4,6,8,10,12,14,16,18,20}\n",
" do\n",
" echo \"Clone${i}\";\n",
" grep -v \"#\" Clone${i}_SV.vcf | wc -l\n",
" done"
]
},
{
"cell_type": "markdown",
"id": "2eb9c407-7623-4f4d-a301-bdf34359ae58",
"metadata": {
"tags": []
},
"source": [
"## Merge all the vcf files across all samples \n",
" \n"
]
},
{
"cell_type": "markdown",
"id": "6479a06d-7cdc-434f-9397-7df84c57efe3",
"metadata": {},
"source": [
"Check the sniffles website https://github.com/fritzsedlazeck/Sniffles/ an its wiki for more details"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "20b379af-3057-48a2-bcf0-11eead7c2cb7",
"metadata": {},
"outputs": [],
"source": [
"cd /home/jovyan/work/SNIFFLES/\n",
"\n",
"# Put all file paths over all vcf files together\n",
"ls *SV.vcf > vcf_raw_calls.txt\n",
"\n",
"# We call SURVIVOR to merge these into one vcf file\n",
"conda activate survivor\n",
"SURVIVOR merge vcf_raw_calls.txt 1000 1 1 -1 -1 -1 merged_SURVIVOR_1kDist.vcf\n",
"conda deactivate"
]
},
{
"cell_type": "markdown",
"id": "e73a58a9-2f58-4574-94a7-9a11fc8b2bca",
"metadata": {},
"source": [
"This will generate one VCF file for all the samples, but we lack the information \n",
"\n",
"If a SV identified in one sample but not in the other is really absent."
]
},
{
"cell_type": "markdown",
"id": "54a53b91-1d29-46ae-aa67-3bd13bb490fa",
"metadata": {},
"source": [
"# 3. Force call all the SVs across all the samples \n",
"\n",
"Next we run Sniffles again across all the samples similar to this :"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3ed8841c-ecf1-4fc6-9343-aca1aedb1fca",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"cd /home/jovyan/work/MAPPING-ONT/BAM/\n",
"\n",
"for i in {2,4,6,8,10,12,14,16,18,20}\n",
" do\n",
" echo -e \"\\n========== REMapping Clone$i======== \\n\";\n",
" sniffles Clone${i}_SORTED.bam -v /home/jovyan/work/SNIFFLES/Clone${i}_SV.gt.vcf --Ivcf /home/jovyan/work/SNIFFLES/merged_SURVIVOR_1kDist.vcf\n",
" done\n"
]
},
{
"cell_type": "markdown",
"id": "4d5da771-38eb-4ab7-a7df-68e92fcff80b",
"metadata": {},
"source": [
"The previous command might have triggered a segmentation fault\n",
"You can download the results with the following command:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b0458adb-14aa-4b5a-a3ad-e0bc250261ab",
"metadata": {},
"outputs": [],
"source": [
"wget --no-check-certificat -rm -nH --cut-dirs=1 --reject=\"index.html*\" https://itrop.ird.fr/ont-training/SNIFFLES.tar.gz\n",
"tar -xzf SNIFFLES.tar.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "83daea67-7430-47df-ba3b-50fa770fc4f4",
"metadata": {},
"outputs": [],
"source": [
"#Put all file paths over all new vcf files together\n",
"cd /home/jovyan/work/SNIFFLES/\n",
"ls *SV.gt.vcf > vcf_gt_calls.txt\n",
"\n",
"# relauch survivor to merge vcf again and finally obtain now a fully genotyped multisample vcf\n",
"conda activate survivor\n",
"SURVIVOR merge vcf_gt_calls.txt 1000 1 1 -1 -1 -1 merged_gt_SURVIVOR_1kDist.vcf\n",
"conda deactivate\n",
"\n",
"# The -1 for the minimum SV caller is necessary to obtain all calls even if they might be 0/0 in all samples."
]
},
{
"cell_type": "markdown",
"id": "2c885731-9eea-4e20-8e41-643f766637b0",
"metadata": {
"tags": []
},
"source": [
"### Have a look on the VCF file "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b1deb0bb-3b0d-44b3-9428-d9a78632b30d",
"metadata": {},
"outputs": [],
"source": [
"head -n 100 merged_gt_SURVIVOR_1kDist.vcf | tail -n 5"
]
},
{
"cell_type": "markdown",
"id": "795c7bdb-dfd3-4fb0-9f4e-452a8f3c6bca",
"metadata": {},
"source": [
"Now, you can plotting proportion of SV in all sample for VCF for example"
]
}
],
"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
}