{ "cells": [ { "cell_type": "markdown", "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", "\n", "[II - Structural variation with Sniffles (LR)](#sniffles) \n", " \n", "\n", " * [1. Check the LR data](#step1)\n", " * [2. SV detection for all CLONES](#step2)\n", " * [3. Force call all the SVs across all the samples](#step3)\n", " * [4. crossing informations between SNP and SV](#step4)\n", "\n", "\n", "\n", "\n", "***\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# __I - Structural variation with `Breakdancer` (SR)__ \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To run BreakDancer, first use `bam2cfg.pl` to prepare the required config file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### First of all, Add the readgroup with the LB tag " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "mkdir -p ~/work/BREAKDANCER/;\n", "cd ~/work/BREAKDANCER/;\n", "\n", "for i in {1..20};\n", " do\n", " gatk AddOrReplaceReadGroups I=/home/jovyan/work/MAPPING-ILL/BAM/Clone$i.sorted.bam O=Clone$i.RG.sorted.bam RGID=Clone$i RGLB=Clone$i RGPL=ILLUMINA RGPU=Unit1 RGSM=Clone$i\n", " done" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "### Run BreakDancer config " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "bam2cfg.pl Clone*.RG.sorted.bam > config_file.cfg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run BreakDancer analysis " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "breakdancer-max config_file.cfg > break.out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Then have a look to the output file " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "head break.out" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "grep -v \"#\" break.out | awk '{print $1\"\\t\"$2\"\\t\"$5}' > break.bed" ] }, { "cell_type": "markdown", "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", "metadata": {}, "source": [ "# 1.Check the LR data " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# check reads ONT from Clones\n", "ONT=/home/jovyan/work/SV_DATA/LONG_READS\n", "ls $ONT" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Do you have the same BAM files as the fastq data ?\n", "cd /home/jovyan/work/MAPPING-ONT/BAM\n", "ls *.sorted.bam" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# 2. SV detection for all CLONES " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run `sniffles` per individual `Clone${i}_ONT.sorted.bam` using a outfile in snf format " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mkdir -p /home/jovyan/work/SNIFFLES\n", "cd /home/jovyan/work/SNIFFLES\n", "\n", "time for i in {1..20}\n", " do\n", " echo \"============ Clone${i}==============\";\n", " sniffles -t 8 --input /home/jovyan/work/MAPPING-ONT/BAM/Clone${i}_ONT.sorted.bam --snf Clone${i}.snf --allow-overwrite\n", " done\n" ] }, { "cell_type": "markdown", "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 :\n", " \n", "⚠️ it's not necessary to do a loop ! check sniffle documentation " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "cd /home/jovyan/work/SNIFFLES\n", "\n", "time sniffles -t8 --input Clone*.snf --vcf multisample.vcf --allow-overwrite\n" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "### Have a look on the VCF file " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cd /home/jovyan/work/SNIFFLES\n", "\n", "head -n 50 multisample.vcf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 4. crossing informations between SNP and SV \n", " \n", " https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creer le dossier de travail VARIANTS_COMPARAISON " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute intersection between SNP and SV found using illumina data `bedtools intersect` ou `intersectBed` \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### count SNP number within each SV" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bedtools intersect -a ../BREAKDANCER/break.bed -b ../VCF_CLONES/output.vcf -c > intersect_break_vs_gatk.bed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### extract all SNPs within each SV" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bedtools intersect -a ../BREAKDANCER/break.bed -b ../VCF_CLONES/output.vcf -wo > intersect_break_vs_gatk_full.bed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute how many regions are communs SV between ONT and ILLUMINA analysis `bedtools intersect` or `intersectBed` \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### count how many SV detected by breakdancer are within each SV detected by sniffle." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### extract all SV breakdancer within each SV ONT" ] }, { "cell_type": "code", "execution_count": null, "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 }