{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping and SV detection for all samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in {1..20}\n",
    "    do\n",
    "        cd /ifb/data/mydatalocal/longReadsAndSV\n",
    "        echo Clone$i;\n",
    "        echo -e \"\\nMapping Clone$i\\n\";\n",
    "        minimap2 -ax map-ont -t 4 --MD -R '@RG\\tID:Clone${i}\\tSM:Clone${i}' /ifb/data/mydatalocal/reference.fasta Clone${i}_LR.fasta > Clone$i.sam\n",
    "        samtools sort -@4 -o Clone$i.SORTED.bam Clone$i.sam\n",
    "        echo -e \"\\nCalling Clone$i\";\n",
    "        sniffles -t 4 -s 2 -q 10 -l 10 -r 500 -m Clone$i.SORTED.bam -v Clone${i}SV.vcf\n",
    "    done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Usage: sniffles [options] -m <sorted.bam> -v <output.vcf> \n",
      "Version: 1.0.11\n",
      "Contact: fritz.sedlazeck@gmail.com\n",
      "\n",
      "Input/Output:\n",
      "    -m <string>,  --mapped_reads <string>\n",
      "        (required)  Sorted bam File\n",
      "    -v <string>,  --vcf <string>\n",
      "        VCF output file name []\n",
      "    -b <string>,  --bedpe <string>\n",
      "         bedpe output file name []\n",
      "    --Ivcf <string>\n",
      "        Input VCF file name. Enable force calling []\n",
      "    --tmp_file <string>\n",
      "        path to temporary file otherwise Sniffles will use the current directory. []\n",
      "\n",
      "General:\n",
      "    -s <int>,  --min_support <int>\n",
      "        Minimum number of reads that support a SV. [10]\n",
      "    --max_num_splits <int>\n",
      "        Maximum number of splits per read to be still taken into account. [7]\n",
      "    -d <int>,  --max_distance <int>\n",
      "        Maximum distance to group SV together. [1000]\n",
      "    -t <int>,  --threads <int>\n",
      "        Number of threads to use. [3]\n",
      "    -l <int>,  --min_length <int>\n",
      "        Minimum length of SV to be reported. [30]\n",
      "    -q <int>,  --minmapping_qual <int>\n",
      "        Minimum Mapping Quality. [20]\n",
      "    -n <int>,  --num_reads_report <int>\n",
      "        Report up to N reads that support the SV in the vcf file. -1: report all. [0]\n",
      "    -r <int>,  --min_seq_size <int>\n",
      "        Discard read if non of its segment is larger then this. [2000]\n",
      "    -z <int>,  --min_zmw <int>\n",
      "        Discard SV that are not supported by at least x zmws. This applies only for PacBio recognizable reads. [0]\n",
      "    --cs_string\n",
      "        Enables the scan of CS string instead of Cigar and MD.  [false]\n",
      "\n",
      "Clustering/phasing and genotyping:\n",
      "    --genotype\n",
      "        Enables Sniffles to compute the genotypes. [false]\n",
      "    --cluster\n",
      "        Enables Sniffles to phase SVs that occur on the same reads [false]\n",
      "    --cluster_support <int>\n",
      "        Minimum number of reads supporting clustering of SV. [1]\n",
      "    -f <float>,  --allelefreq <float>\n",
      "        Threshold on allele frequency (0-1).  [0]\n",
      "    --min_homo_af <float>\n",
      "        Threshold on allele frequency (0-1).  [0.8]\n",
      "    --min_het_af <float>\n",
      "        Threshold on allele frequency (0-1).  [0.3]\n",
      "\n",
      "Advanced:\n",
      "    --report_BND\n",
      "        Dont report BND instead use Tra in vcf output.  [true]\n",
      "    --report_seq\n",
      "        Report sequences for indels in vcf output. (Beta version!)  [false]\n",
      "    --ignore_sd\n",
      "        Ignores the sd based filtering.  [false]\n",
      "    --report_read_strands\n",
      "        Enables the report of the strand categories per read. (Beta)  [false]\n",
      "    --ccs_reads\n",
      "        Preset CCS Pacbio setting. (Beta)  [false]\n",
      "\n",
      "Parameter estimation:\n",
      "    --skip_parameter_estimation\n",
      "        Enables the scan if only very few reads are present.  [false]\n",
      "    --del_ratio <float>\n",
      "        Estimated ration of deletions per read (0-1).  [0.0458369]\n",
      "    --ins_ratio <float>\n",
      "        Estimated ratio of insertions per read (0-1).  [0.049379]\n",
      "    --max_diff_per_window <int>\n",
      "        Maximum differences per 100bp. [50]\n",
      "    --max_dist_aln_events <int>\n",
      "        Maximum distance between alignment (indel) events. [4]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sniffles --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting SV vs SNP, one sample for SV, all sample for VCFc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reference\t274133\t98\tN\t<DEL>\t.\tPASS\tPRECISE;SVMETHOD=Snifflesv1.0.11;CHR2=Reference;END=809621;STD_quant_start=0.000000;STD_quant_stop=0.000000;Kurtosis_quant_start=-nan;Kurtosis_quant_stop=0.333333;SVTYPE=DEL;SUPTYPE=SR;SVLEN=-535488;STRANDS=+-;RE=16\tGT:DR:DV\t./.:.:16\t281150\n",
      "Reference\t276639\t99\tN\t<DUP>\t.\tPASS\tPRECISE;SVMETHOD=Snifflesv1.0.11;CHR2=Reference;END=277334;STD_quant_start=0.000000;STD_quant_stop=0.000000;Kurtosis_quant_start=-1.692308;Kurtosis_quant_stop=-1.903226;SVTYPE=DUP;SUPTYPE=SR;SVLEN=695;STRANDS=-+;RE=136\tGT:DR:DV\t./.:.:136\t478\n"
     ]
    }
   ],
   "source": [
    "# TIP use the rawVCF file\n",
    "#bedtools intersect -c -a Clone15SV.vcf -b output.vcf\n",
    "head -n 100 intersectSVallSNP_Clone15.vcf | tail -n 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merging and genotyping all SV accross all samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ls: cannot access '*SV.vcf': No such file or directory\n"
     ]
    },
    {
     "ename": "",
     "evalue": "2",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "# TIP check the sniffles website https://github.com/fritzsedlazeck/Sniffles/ an its wiki\n",
    "ls *SV.vcf > vcf_raw_calls.txt\n",
    "survivor merge vcf_raw_calls.txt 1000 1 1 -1 -1 -1 merged_SURVIVOR_1kDist.vcf\n",
    "for i in {1..20}\n",
    "    do\n",
    "        echo Clone$i;\n",
    "        mkdir dirClone$i; \n",
    "        cd dirClone$i; \n",
    "        echo -e \"\\nREMapping Clone$i\\n\";\n",
    "        sniffles -t 4 -s 2 -q 10 -l 10 -r 500 -m Clone$i.SORTED.bam -v Clone${i}SV.gt.vcf --Ivcf merged_SURVIVOR_1kDist.vcf\n",
    "    done\n",
    "\n",
    "ls *SV.gt.vcf > vcf_gt_calls.txt\n",
    "survivor merge vcf_gt_calls.txt 1000 1 1 -1 -1 -1 merged_gt_SURVIVOR_1kDist.vcf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Have a look on the VCF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head -n 100 merged_gt_SURVIVOR_1kDist.vcf | tail -n 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot SV vs SNP across all samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TIP use again the rawVCF file\n",
    "# You can also use the final filtered one to see overlaps"
   ]
  }
 ],
 "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": 4
}