01_fastq_to_vcf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3b6ee172-0c30-42fc-93b9-3609dd496ffe",
   "metadata": {},
   "source": [
    "# Variant-Calling Oxford Nanopore Tech (ONT) Long-Read Whole Genome Sequencing (WGS) Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7850f07e-7de4-439d-ae5a-ad488ba06761",
   "metadata": {},
   "source": [
    "**GOAL**: *create a model to distinguish true variants from artifactual variant calls*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "312a2f66-de90-4464-922f-e20f3ce39fb8",
   "metadata": {},
   "source": [
    "- data is from the [Nanopore WGS Consortium](https://github.com/nanopore-wgs-consortium/NA12878/blob/master/Genome.md) GitHub repository\n",
    "- they sequenced the CEPH1463 (NA12878/GM12878, Ceph/Utah pedigree) human genome reference standard on the Oxford Nanopore MinION using 1D ligation kits (450 bp/s) and R9.4 chemistry (FLO-MIN106)\n",
    "- sequencing was done by [5 labs](https://www.nature.com/articles/nbt.4060) on 53 flow cells using both DNA purchased from Coriell as well as DNA extracted from cultured cells\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2d13301-1afb-440a-9e34-5549b4ffd681",
   "metadata": {},
   "source": [
    "This analysis involved:\n",
    "- downloading raw FAST5 files (~4.5TB) and basecalling to FASTQ with the latest, [high-performance](https://aws.amazon.com/blogs/hpc/benchmarking-the-oxford-nanopore-technologies-basecallers-on-aws/)  [dorado basecaller](https://github.com/nanoporetech/dorado) using the 'sup' (highest accuracy) model \n",
    "- mapping to the hg38 human reference genome to generate BAM files (minimap2, samtools)\n",
    "- marking duplicates (samtools)\n",
    "- merging the resulting 52 BAM files into 1 file (samtools)\n",
    "- variant calling to generate VCF file (bcftools mpileup & call)\n",
    "- intersecting with gold-standard [NIST Genome in a Bottle](https://www.nist.gov/programs-projects/genome-bottle) data to enable comparison (bedtools intersect)\n",
    "- creating machine learning and deep learning models to distinguish true variant calls from artifacts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc8dea7-d1f2-4509-b93a-0cddc06b4f62",
   "metadata": {},
   "source": [
    "Analysis was performed on a high performance computer (HPC) cluster running CentOS Linux with SLURM as a resource manager & job scheduler."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2f2788c-5559-41cd-b288-efbdc64ae429",
   "metadata": {},
   "source": [
    "## 1. Mapping to hg38 human reference genome (minimap2, samtools)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "812787ba-0e52-4a12-a6ac-505f2f2ec77c",
   "metadata": {},
   "source": [
    "We can use a simple bash script to loop over all the FASTQ files in a directory and submit them as batch jobs to the HPC using the SLURM scheduler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bbed6e8b-b0b9-49d4-843f-d591f31aeee0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!/bin/bash\n",
      "# loop over all *.fastq files in the directory input_dir\n",
      "# print the filename (so we have some visual progress indicator)\n",
      "# then submit to SLURM\n",
      "\n",
      "\n",
      "slurm_script=$1        # SLURM SCRIPT OF INTEREST\n",
      "input_dir=$2             # DIRECTORY OF INTEREST\n",
      "output_dir=$3          # OUTPUT DIRECTORY\n",
      "\n",
      "for FILE in \"${input_dir}\"*.fastq; do\n",
      "echo ${FILE}\n",
      "sbatch $slurm_script \\\n",
      "    ${FILE} \\\n",
      "    $output_dir\n",
      "sleep 1 # pause to be kind to the scheduler\n",
      "done"
     ]
    }
   ],
   "source": [
    "!cat ~/tools/slurm_scripts/bash_slurm_sbatch_loop_directory_files_fastq_221221.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1685536d-8fb4-4e18-b81c-69ea4ef295f1",
   "metadata": {},
   "source": [
    "---\n",
    "In this case, the SLURM job will be mapping each of the 52 FASTQ files to the hg38 human reference genome using minimap2, the sequence aligner recommended by Oxford Nanopore Technologies. To output a sorted BAM file, samtools is used on the output of minimap2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ee888683-3d4a-40da-9819-cac2772ac6dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!/bin/bash\n",
      "#\n",
      "#SBATCH --job-name=minimap2_samtoolsBam_RG_231112\n",
      "#SBATCH --output=/home/fmbuga/tools/slurm_scripts/slurm_out_err/minimap2_samtoolsBam_RG_231112_%j.out\n",
      "#SBATCH --error=/home/fmbuga/tools/slurm_scripts/slurm_out_err/minimap2_samtoolsBam_RG_231112_%j.err\n",
      "#\n",
      "#SBATCH --ntasks=1\n",
      "#SBATCH --cpus-per-task=28\n",
      "#SBATCH --mem-per-cpu=4G\n",
      "\n",
      "datenow=$(date)\n",
      "echo $datenow\n",
      "echo \"\"\n",
      "srun hostname\n",
      "echo \"\"\n",
      "\n",
      "start=$(date +%s)\n",
      "echo \"start time: $start\"\n",
      "echo \"\"\n",
      "echo \"hostname: $HOSTNAME\"\n",
      "echo \"\"\n",
      "\n",
      "\n",
      "##################\n",
      "\n",
      "eval \"$(conda shell.bash hook)\"\n",
      "conda activate minimap2_224                      ### CONDA ENVIRONMENT NAME\n",
      "conda info\n",
      "echo \"\"\n",
      "\n",
      "call=\"minimap2 --version\" # to get the version\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "call=\"samtools --version\" # to get the version\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "call=\"bcftools --version\" # to get the version\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "fastq_path=$1\n",
      "bam_output_dir=$2\n",
      "\n",
      "flowcell_asic_id=$(echo $fastq_path | awk -F'[/_.]' '{print $10}')    # https://github.com/nanopore-wgs-consortium/NA12878/blob/master/Genome.md\n",
      "echo $flowcell_asic_id\n",
      "echo \"\"\n",
      "\n",
      "\n",
      "call=$(minimap2 \\\n",
      "        -t 28 \\\n",
      "        -R \"@RG\\tID:$flowcell_asic_id\\tDS:$flowcell_asic_id\\tSM:NA12878-nanopore-wgs\" \\\n",
      "        -a \\\n",
      "        ~/tools/GIAB/GRCh38/refs/GRCh38_GIABv3.fasta \\\n",
      "        $fastq_path | \\\n",
      "            samtools view -uh | \\\n",
      "                samtools sort -o \"$bam_output_dir$flowcell_asic_id\".bam)\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "echo \"minimap2 -> samtoolsBam COMPLETE.\"\n",
      "echo \"\"\n",
      "    \n",
      "##################\n",
      "\n",
      "\n",
      "echo \"\"\n",
      "end=$(date +%s)\n",
      "echo \"end time: $end\"\n",
      "runtime_s=$(echo $(( end - start )))\n",
      "echo \"total run time(s): $runtime_s\"\n",
      "sec_per_min=60\n",
      "sec_per_hr=3600\n",
      "runtime_m=$(echo \"scale=2; $runtime_s / $sec_per_min;\" | bc)\n",
      "echo \"total run time(m): $runtime_m\"\n",
      "runtime_h=$(echo \"scale=2; $runtime_s / $sec_per_hr;\" | bc)\n",
      "echo \"total run time(h): $runtime_h\"\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!cat ~/tools/slurm_scripts/minimap2_samtoolsBam_RG_231112_slurm.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e562729f-1e79-4367-b46b-1238b3bfdcf5",
   "metadata": {},
   "source": [
    "---\n",
    "And here is the actual job submission to SLURM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c401a7c9-a3f5-49a8-b64b-7be8062aca53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB39043-3709921973.fastq\n",
      "Submitted batch job 1111239\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB39075-4246400039.fastq\n",
      "Submitted batch job 1111240\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB39088-288418386.fastq\n",
      "Submitted batch job 1111241\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB41174-3976885577.fastq\n",
      "Submitted batch job 1111242\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42205-3573838535.fastq\n",
      "Submitted batch job 1111243\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42260-4177064552.fastq\n",
      "Submitted batch job 1111244\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42316-216722908.fastq\n",
      "Submitted batch job 1111245\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42395-4178605061.fastq\n",
      "Submitted batch job 1111246\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42451-4239353418.fastq\n",
      "Submitted batch job 1111247\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42473-4179682758.fastq\n",
      "Submitted batch job 1111248\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42476-3843483077.fastq\n",
      "Submitted batch job 1111249\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42561-356443753.fastq\n",
      "Submitted batch job 1111250\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42704-87746129.fastq\n",
      "Submitted batch job 1111251\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42706-4111103328.fastq\n",
      "Submitted batch job 1111252\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42798-3306352129.fastq\n",
      "Submitted batch job 1111253\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42804-84744914.fastq\n",
      "Submitted batch job 1111254\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42810-352384898.fastq\n",
      "Submitted batch job 1111255\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB42828-288548394.fastq\n",
      "Submitted batch job 1111256\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB43577-3574887596.fastq\n",
      "Submitted batch job 1111257\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB44989-2567311907.fastq\n",
      "Submitted batch job 1111258\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB45271-152889212.fastq\n",
      "Submitted batch job 1111259\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB45277-86567043.fastq\n",
      "Submitted batch job 1111260\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB45280-222619780.fastq\n",
      "Submitted batch job 1111261\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB45321-19064779.fastq\n",
      "Submitted batch job 1111262\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB45321-285174896.fastq\n",
      "Submitted batch job 1111263\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB45332-551111640.fastq\n",
      "Submitted batch job 1111264\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB46664-288286712.fastq\n",
      "Submitted batch job 1111265\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB46683-4246923067.fastq\n",
      "Submitted batch job 1111266\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB49164-4045668814.fastq\n",
      "Submitted batch job 1111267\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB49712-622291475.fastq\n",
      "Submitted batch job 1111268\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB49908-481119249.fastq\n",
      "Submitted batch job 1111269\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAB49914-3775529215.fastq\n",
      "Submitted batch job 1111270\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF01127-353303576.fastq\n",
      "Submitted batch job 1111271\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF01132-84868110.fastq\n",
      "Submitted batch job 1111272\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF01169-4245879798.fastq\n",
      "Submitted batch job 1111273\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF01253-83756522.fastq\n",
      "Submitted batch job 1111274\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF01441-3910073345.fastq\n",
      "Submitted batch job 1111275\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF04090-3842965088.fastq\n",
      "Submitted batch job 1111276\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF05869-87644245.fastq\n",
      "Submitted batch job 1111277\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF09277-3709819546.fastq\n",
      "Submitted batch job 1111278\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF09640-4245291640.fastq\n",
      "Submitted batch job 1111279\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF09701-4249180049.fastq\n",
      "Submitted batch job 1111280\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF09713-4111860526.fastq\n",
      "Submitted batch job 1111281\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF09968-3439856925.fastq\n",
      "Submitted batch job 1111282\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF10039-2901545329.fastq\n",
      "Submitted batch job 1111283\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF13748-17958431.fastq\n",
      "Submitted batch job 1111284\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF14035-3976726082.fastq\n",
      "Submitted batch job 1111285\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF15586-82266371.fastq\n",
      "Submitted batch job 1111286\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF15630-4244782843.fastq\n",
      "Submitted batch job 1111287\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF15665-16056159.fastq\n",
      "Submitted batch job 1111288\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF15694-4109802543.fastq\n",
      "Submitted batch job 1111289\n",
      "/home/fmbuga/nanopore/rogstrix/dorado/00_fastq_raw/FAF18554-4178920553.fastq\n",
      "Submitted batch job 1111290\n",
      "             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)\n",
      "           1111260     nodes minimap2   fmbuga PD       0:00      1 (Resources)\n",
      "           1111261     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111262     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111263     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111264     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111265     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111266     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111267     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111268     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111269     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111270     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111271     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111272     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111273     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111274     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111275     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111276     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111277     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111278     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111279     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111280     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111281     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111282     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111283     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111284     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111285     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111286     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111287     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111288     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111289     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111290     nodes minimap2   fmbuga PD       0:00      1 (Priority)\n",
      "           1111238     nodes     bash   fmbuga  R      13:09      1 node05\n",
      "           1111239     nodes minimap2   fmbuga  R       0:55      1 node10\n",
      "           1111240     nodes minimap2   fmbuga  R       0:52      1 node11\n",
      "           1111241     nodes minimap2   fmbuga  R       0:52      1 node12\n",
      "           1111242     nodes minimap2   fmbuga  R       0:52      1 node13\n",
      "           1111243     nodes minimap2   fmbuga  R       0:49      1 node14\n",
      "           1111244     nodes minimap2   fmbuga  R       0:49      1 node47\n",
      "           1111245     nodes minimap2   fmbuga  R       0:49      1 node48\n",
      "           1111246     nodes minimap2   fmbuga  R       0:46      1 node49\n",
      "           1111247     nodes minimap2   fmbuga  R       0:46      1 node50\n",
      "           1111248     nodes minimap2   fmbuga  R       0:43      1 node51\n",
      "           1111249     nodes minimap2   fmbuga  R       0:43      1 node52\n",
      "           1111250     nodes minimap2   fmbuga  R       0:43      1 node53\n",
      "           1111251     nodes minimap2   fmbuga  R       0:40      1 node54\n",
      "           1111252     nodes minimap2   fmbuga  R       0:40      1 node55\n",
      "           1111253     nodes minimap2   fmbuga  R       0:40      1 node56\n",
      "           1111254     nodes minimap2   fmbuga  R       0:37      1 node57\n",
      "           1111255     nodes minimap2   fmbuga  R       0:37      1 node58\n",
      "           1111256     nodes minimap2   fmbuga  R       0:37      1 node59\n",
      "           1111257     nodes minimap2   fmbuga  R       0:34      1 node60\n",
      "           1111258     nodes minimap2   fmbuga  R       0:34      1 node61\n",
      "           1111259     nodes minimap2   fmbuga  R       0:31      1 node62\n"
     ]
    }
   ],
   "source": [
    "# submit fastq-to-bam jobs to slurm\n",
    "!mkdir -p 01_bam_raw\n",
    "\n",
    "!bash ~/tools/slurm_scripts/bash_slurm_sbatch_loop_directory_files_fastq_221221.sh \\\n",
    "    ~/tools/slurm_scripts/minimap2_samtoolsBam_RG_231112_slurm.sh \\\n",
    "    ~/nanopore/rogstrix/dorado/00_fastq_raw/ \\\n",
    "    01_bam_raw/\n",
    "\n",
    "!squeue -u fmbuga"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f99b5d3a-3352-46cf-999b-e6ab388e651f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "01_bam_raw/FAB39043-3709921973.bam\n",
      "01_bam_raw/FAB39075-4246400039.bam\n",
      "01_bam_raw/FAB39088-288418386.bam\n",
      "01_bam_raw/FAB41174-3976885577.bam\n",
      "01_bam_raw/FAB42205-3573838535.bam\n",
      "01_bam_raw/FAB42260-4177064552.bam\n",
      "01_bam_raw/FAB42316-216722908.bam\n",
      "01_bam_raw/FAB42395-4178605061.bam\n",
      "01_bam_raw/FAB42451-4239353418.bam\n",
      "01_bam_raw/FAB42473-4179682758.bam\n",
      "01_bam_raw/FAB42476-3843483077.bam\n",
      "01_bam_raw/FAB42561-356443753.bam\n",
      "01_bam_raw/FAB42704-87746129.bam\n",
      "01_bam_raw/FAB42706-4111103328.bam\n",
      "01_bam_raw/FAB42798-3306352129.bam\n",
      "01_bam_raw/FAB42804-84744914.bam\n",
      "01_bam_raw/FAB42810-352384898.bam\n",
      "01_bam_raw/FAB42828-288548394.bam\n",
      "01_bam_raw/FAB43577-3574887596.bam\n",
      "01_bam_raw/FAB44989-2567311907.bam\n",
      "01_bam_raw/FAB45271-152889212.bam\n",
      "01_bam_raw/FAB45277-86567043.bam\n",
      "01_bam_raw/FAB45280-222619780.bam\n",
      "01_bam_raw/FAB45321-19064779.bam\n",
      "01_bam_raw/FAB45321-285174896.bam\n",
      "01_bam_raw/FAB45332-551111640.bam\n",
      "01_bam_raw/FAB46664-288286712.bam\n",
      "01_bam_raw/FAB46683-4246923067.bam\n",
      "01_bam_raw/FAB49164-4045668814.bam\n",
      "01_bam_raw/FAB49712-622291475.bam\n",
      "01_bam_raw/FAB49908-481119249.bam\n",
      "01_bam_raw/FAB49914-3775529215.bam\n",
      "01_bam_raw/FAF01127-353303576.bam\n",
      "01_bam_raw/FAF01132-84868110.bam\n",
      "01_bam_raw/FAF01169-4245879798.bam\n",
      "01_bam_raw/FAF01253-83756522.bam\n",
      "01_bam_raw/FAF01441-3910073345.bam\n",
      "01_bam_raw/FAF04090-3842965088.bam\n",
      "01_bam_raw/FAF05869-87644245.bam\n",
      "01_bam_raw/FAF09277-3709819546.bam\n",
      "01_bam_raw/FAF09640-4245291640.bam\n",
      "01_bam_raw/FAF09701-4249180049.bam\n",
      "01_bam_raw/FAF09713-4111860526.bam\n",
      "01_bam_raw/FAF09968-3439856925.bam\n",
      "01_bam_raw/FAF10039-2901545329.bam\n",
      "01_bam_raw/FAF13748-17958431.bam\n",
      "01_bam_raw/FAF14035-3976726082.bam\n",
      "01_bam_raw/FAF15586-82266371.bam\n",
      "01_bam_raw/FAF15630-4244782843.bam\n",
      "01_bam_raw/FAF15665-16056159.bam\n",
      "01_bam_raw/FAF15694-4109802543.bam\n",
      "01_bam_raw/FAF18554-4178920553.bam\n",
      "52\n"
     ]
    }
   ],
   "source": [
    "# check the output and confirm all 52 files were processed\n",
    "!ls -1 01_bam_raw/*.bam \\\n",
    "    | tee >(wc -l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "25481e2f-7bf1-4009-af31-283d0ea6cd74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1650805 + 0 in total (QC-passed reads + QC-failed reads)\n",
      "719983 + 0 primary\n",
      "480392 + 0 secondary\n",
      "450430 + 0 supplementary\n",
      "0 + 0 duplicates\n",
      "0 + 0 primary duplicates\n",
      "1588422 + 0 mapped (96.22% : N/A)\n",
      "657600 + 0 primary mapped (91.34% : N/A)\n"
     ]
    }
   ],
   "source": [
    "# check how mapping went\n",
    "!samtools flagstats \\\n",
    "    --threads 28 \\\n",
    "    --verbosity 5 \\\n",
    "    01_bam_raw/FAF01132-84868110.bam | head -n 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "442ba5f9-473a-4bf9-9908-0e1ada2d4351",
   "metadata": {},
   "source": [
    "## 2. Marking duplicates (samtools)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6052eb8-6d91-4819-9e3a-b9cb62d3fb7b",
   "metadata": {},
   "source": [
    "Next, we [mark duplicates](http://www.htslib.org/doc/samtools-markdup.html) on the 52 BAM files using `samtools markdup`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "873083ac-591c-4b1e-8dfa-ae0ec607c85c",
   "metadata": {},
   "source": [
    "We can re-purpose the looping bash script used above in the mapping step for SLURM job submission with a slight change so that it loops over all BAM files in a folder instead of all FASTQ files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2ed6635b-2c8a-4ebe-b02b-16e0d86fd1a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!/bin/bash\n",
      "# loop over all *.bam files in the directory input_dir\n",
      "# print the filename (so we have some visual progress indicator)\n",
      "# then submit to SLURM\n",
      "\n",
      "\n",
      "slurm_script=$1        # SLURM SCRIPT OF INTEREST\n",
      "input_dir=$2             # DIRECTORY OF INTEREST                   \n",
      "output_dir=$3             # OUTPUT DIRECTORY                   \n",
      "\n",
      "\n",
      "for FILE in \"${input_dir}\"*.bam; do\n",
      "echo ${FILE}\n",
      "sbatch $slurm_script \\\n",
      "    ${FILE} \\\n",
      "    $output_dir\n",
      "sleep 1 # pause to be kind to the scheduler\n",
      "done"
     ]
    }
   ],
   "source": [
    "!cat ~/tools/slurm_scripts/bash_slurm_sbatch_loop_directory_files_samtools_markdup_221222.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aebae47-27a9-4eaa-b855-c2f5c061d438",
   "metadata": {},
   "source": [
    "---\n",
    "In order to use `samtools markdup` to mark duplicates in the BAM files, the reads first have to be [ordered by name](http://www.htslib.org/doc/samtools-markdup.html) (`samtools collate`), annotated with [tags for mate score and mate coordinates](http://www.htslib.org/doc/samtools-fixmate.html) (`samtools fixmate`) and [sorted in position order](http://www.htslib.org/doc/samtools-sort.html) (`samtools sort`). These commands can be chained with pipes with the BAM files going in and outputting processed BAM files with duplicates marked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1d29a747-d0d5-4384-8ccb-8136572add3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!/bin/bash\n",
      "#\n",
      "#SBATCH --job-name=samtools_collate_fixmate_sort_markdup_221222\n",
      "#SBATCH --output=/home/fmbuga/tools/slurm_scripts/slurm_out_err/samtools_collate_fixmate_sort_markdup_221222_%j.out\n",
      "#SBATCH --error=/home/fmbuga/tools/slurm_scripts/slurm_out_err/samtools_collate_fixmate_sort_markdup_221222_%j.err\n",
      "#\n",
      "#SBATCH --ntasks=1\n",
      "#SBATCH --cpus-per-task=28\n",
      "#SBATCH --mem-per-cpu=4G\n",
      "\n",
      "datenow=$(date)\n",
      "echo $datenow\n",
      "echo \"\"\n",
      "srun hostname\n",
      "echo \"\"\n",
      "\n",
      "start=$(date +%s)\n",
      "echo \"start time: $start\"\n",
      "echo \"\"\n",
      "echo \"hostname: $HOSTNAME\"\n",
      "echo \"\"\n",
      "\n",
      "\n",
      "##################\n",
      "\n",
      "eval \"$(conda shell.bash hook)\"\n",
      "conda activate minimap2_224                      ### CONDA ENVIRONMENT NAME\n",
      "conda info\n",
      "echo \"\"\n",
      "\n",
      "call=\"samtools --version\" # to get the version\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "input_bam=$1\n",
      "output_dir=$2\n",
      "\n",
      "out_name=$(basename $input_bam .bam)\n",
      "\n",
      "call=$(samtools collate \\\n",
      "    -Ou \\\n",
      "    --threads 28 \\\n",
      "    $input_bam | \\\n",
      "        samtools fixmate \\\n",
      "            - \\\n",
      "            - \\\n",
      "            --threads 28 \\\n",
      "            -u | \\\n",
      "            samtools sort \\\n",
      "                - \\\n",
      "                --threads 28 \\\n",
      "                -u |\\\n",
      "                samtools markdup \\\n",
      "                    - \\\n",
      "                    --threads 28 \\\n",
      "                    -sf \"$output_dir$out_name\"_collate_fixmate_sort_markdup_bam.markdupstats \\\n",
      "                    --write-index \\\n",
      "                    \"$output_dir$out_name\"_collate_fixmate_sort_markdup.bam)\n",
      "\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "echo \"samtools collate > fixmate > sort > markdup COMPLETE.\"\n",
      "echo \"\"\n",
      "    \n",
      "##################\n",
      "\n",
      "\n",
      "echo \"\"\n",
      "end=$(date +%s)\n",
      "echo \"end time: $end\"\n",
      "runtime_s=$(echo $(( end - start )))\n",
      "echo \"total run time(s): $runtime_s\"\n",
      "sec_per_min=60\n",
      "sec_per_hr=3600\n",
      "runtime_m=$(echo \"scale=2; $runtime_s / $sec_per_min;\" | bc)\n",
      "echo \"total run time(m): $runtime_m\"\n",
      "runtime_h=$(echo \"scale=2; $runtime_s / $sec_per_hr;\" | bc)\n",
      "echo \"total run time(h): $runtime_h\"\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# the SLURM job script for marking duplicates\n",
    "!cat ~/tools/slurm_scripts/samtools_collate_fixmate_sort_markdup_221222_slurm.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "583733f5-e9dc-458e-ae81-c3d5f8a77891",
   "metadata": {},
   "source": [
    "---\n",
    "And the SLURM job submission with the looping bash script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "960f5241-8b97-4bd0-8210-fb1ecf64156b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "01_bam_raw/FAB39043-3709921973.bam\n",
      "Submitted batch job 1111379\n",
      "01_bam_raw/FAB39075-4246400039.bam\n",
      "Submitted batch job 1111380\n",
      "01_bam_raw/FAB39088-288418386.bam\n",
      "Submitted batch job 1111381\n",
      "01_bam_raw/FAB41174-3976885577.bam\n",
      "Submitted batch job 1111382\n",
      "01_bam_raw/FAB42205-3573838535.bam\n",
      "Submitted batch job 1111383\n",
      "01_bam_raw/FAB42260-4177064552.bam\n",
      "Submitted batch job 1111384\n",
      "01_bam_raw/FAB42316-216722908.bam\n",
      "Submitted batch job 1111385\n",
      "01_bam_raw/FAB42395-4178605061.bam\n",
      "Submitted batch job 1111386\n",
      "01_bam_raw/FAB42451-4239353418.bam\n",
      "Submitted batch job 1111387\n",
      "01_bam_raw/FAB42473-4179682758.bam\n",
      "Submitted batch job 1111388\n",
      "01_bam_raw/FAB42476-3843483077.bam\n",
      "Submitted batch job 1111389\n",
      "01_bam_raw/FAB42561-356443753.bam\n",
      "Submitted batch job 1111390\n",
      "01_bam_raw/FAB42704-87746129.bam\n",
      "Submitted batch job 1111391\n",
      "01_bam_raw/FAB42706-4111103328.bam\n",
      "Submitted batch job 1111392\n",
      "01_bam_raw/FAB42798-3306352129.bam\n",
      "Submitted batch job 1111393\n",
      "01_bam_raw/FAB42804-84744914.bam\n",
      "Submitted batch job 1111394\n",
      "01_bam_raw/FAB42810-352384898.bam\n",
      "Submitted batch job 1111395\n",
      "01_bam_raw/FAB42828-288548394.bam\n",
      "Submitted batch job 1111396\n",
      "01_bam_raw/FAB43577-3574887596.bam\n",
      "Submitted batch job 1111397\n",
      "01_bam_raw/FAB44989-2567311907.bam\n",
      "Submitted batch job 1111398\n",
      "01_bam_raw/FAB45271-152889212.bam\n",
      "Submitted batch job 1111399\n",
      "01_bam_raw/FAB45277-86567043.bam\n",
      "Submitted batch job 1111400\n",
      "01_bam_raw/FAB45280-222619780.bam\n",
      "Submitted batch job 1111401\n",
      "01_bam_raw/FAB45321-19064779.bam\n",
      "Submitted batch job 1111402\n",
      "01_bam_raw/FAB45321-285174896.bam\n",
      "Submitted batch job 1111403\n",
      "01_bam_raw/FAB45332-551111640.bam\n",
      "Submitted batch job 1111404\n",
      "01_bam_raw/FAB46664-288286712.bam\n",
      "Submitted batch job 1111405\n",
      "01_bam_raw/FAB46683-4246923067.bam\n",
      "Submitted batch job 1111406\n",
      "01_bam_raw/FAB49164-4045668814.bam\n",
      "Submitted batch job 1111407\n",
      "01_bam_raw/FAB49712-622291475.bam\n",
      "Submitted batch job 1111408\n",
      "01_bam_raw/FAB49908-481119249.bam\n",
      "Submitted batch job 1111409\n",
      "01_bam_raw/FAB49914-3775529215.bam\n",
      "Submitted batch job 1111410\n",
      "01_bam_raw/FAF01127-353303576.bam\n",
      "Submitted batch job 1111411\n",
      "01_bam_raw/FAF01132-84868110.bam\n",
      "Submitted batch job 1111412\n",
      "01_bam_raw/FAF01169-4245879798.bam\n",
      "Submitted batch job 1111413\n",
      "01_bam_raw/FAF01253-83756522.bam\n",
      "Submitted batch job 1111414\n",
      "01_bam_raw/FAF01441-3910073345.bam\n",
      "Submitted batch job 1111415\n",
      "01_bam_raw/FAF04090-3842965088.bam\n",
      "Submitted batch job 1111416\n",
      "01_bam_raw/FAF05869-87644245.bam\n",
      "Submitted batch job 1111417\n",
      "01_bam_raw/FAF09277-3709819546.bam\n",
      "Submitted batch job 1111418\n",
      "01_bam_raw/FAF09640-4245291640.bam\n",
      "Submitted batch job 1111419\n",
      "01_bam_raw/FAF09701-4249180049.bam\n",
      "Submitted batch job 1111420\n",
      "01_bam_raw/FAF09713-4111860526.bam\n",
      "Submitted batch job 1111421\n",
      "01_bam_raw/FAF09968-3439856925.bam\n",
      "Submitted batch job 1111422\n",
      "01_bam_raw/FAF10039-2901545329.bam\n",
      "Submitted batch job 1111423\n",
      "01_bam_raw/FAF13748-17958431.bam\n",
      "Submitted batch job 1111424\n",
      "01_bam_raw/FAF14035-3976726082.bam\n",
      "Submitted batch job 1111425\n",
      "01_bam_raw/FAF15586-82266371.bam\n",
      "Submitted batch job 1111426\n",
      "01_bam_raw/FAF15630-4244782843.bam\n",
      "Submitted batch job 1111427\n",
      "01_bam_raw/FAF15665-16056159.bam\n",
      "Submitted batch job 1111428\n",
      "01_bam_raw/FAF15694-4109802543.bam\n",
      "Submitted batch job 1111429\n",
      "01_bam_raw/FAF18554-4178920553.bam\n",
      "Submitted batch job 1111430\n",
      "             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)\n",
      "           1111424     nodes samtools   fmbuga PD       0:00      1 (Resources)\n",
      "           1111425     nodes samtools   fmbuga PD       0:00      1 (Priority)\n",
      "           1111426     nodes samtools   fmbuga PD       0:00      1 (Priority)\n",
      "           1111427     nodes samtools   fmbuga PD       0:00      1 (Priority)\n",
      "           1111428     nodes samtools   fmbuga PD       0:00      1 (Priority)\n",
      "           1111429     nodes samtools   fmbuga PD       0:00      1 (Priority)\n",
      "           1111430     nodes samtools   fmbuga PD       0:00      1 (Priority)\n",
      "           1111378     nodes     bash   fmbuga  R       9:39      1 node09\n",
      "           1111380     nodes samtools   fmbuga  R       0:52      1 node11\n",
      "           1111381     nodes samtools   fmbuga  R       0:52      1 node12\n",
      "           1111385     nodes samtools   fmbuga  R       0:46      1 node40\n",
      "           1111387     nodes samtools   fmbuga  R       0:46      1 node42\n",
      "           1111392     nodes samtools   fmbuga  R       0:40      1 node47\n",
      "           1111395     nodes samtools   fmbuga  R       0:37      1 node49\n",
      "           1111397     nodes samtools   fmbuga  R       0:36      1 node13\n",
      "           1111399     nodes samtools   fmbuga  R       0:34      1 node52\n",
      "           1111400     nodes samtools   fmbuga  R       0:32      1 node41\n",
      "           1111402     nodes samtools   fmbuga  R       0:31      1 node54\n",
      "           1111404     nodes samtools   fmbuga  R       0:29      1 node10\n",
      "           1111405     nodes samtools   fmbuga  R       0:28      1 node55\n",
      "           1111407     nodes samtools   fmbuga  R       0:25      1 node57\n",
      "           1111408     nodes samtools   fmbuga  R       0:25      1 node58\n",
      "           1111409     nodes samtools   fmbuga  R       0:23      1 node15\n",
      "           1111410     nodes samtools   fmbuga  R       0:23      1 node59\n",
      "           1111411     nodes samtools   fmbuga  R       0:22      1 node60\n",
      "           1111412     nodes samtools   fmbuga  R       0:19      1 node61\n",
      "           1111413     nodes samtools   fmbuga  R       0:19      1 node62\n",
      "           1111414     nodes samtools   fmbuga  R       0:17      1 node45\n",
      "           1111415     nodes samtools   fmbuga  R       0:16      1 node14\n",
      "           1111416     nodes samtools   fmbuga  R       0:15      1 node43\n",
      "           1111417     nodes samtools   fmbuga  R       0:15      1 node48\n",
      "           1111418     nodes samtools   fmbuga  R       0:14      1 node56\n",
      "           1111419     nodes samtools   fmbuga  R       0:13      1 node53\n",
      "           1111420     nodes samtools   fmbuga  R       0:10      1 node50\n",
      "           1111421     nodes samtools   fmbuga  R       0:09      1 node44\n",
      "           1111422     nodes samtools   fmbuga  R       0:06      1 node51\n",
      "           1111423     nodes samtools   fmbuga  R       0:01      1 node46\n"
     ]
    }
   ],
   "source": [
    "!mkdir -p 02_bam_markdup\n",
    "!bash ~/tools/slurm_scripts/bash_slurm_sbatch_loop_directory_files_samtools_markdup_221222.sh \\\n",
    "    ~/tools/slurm_scripts/samtools_collate_fixmate_sort_markdup_221222_slurm.sh \\\n",
    "    01_bam_raw/ \\\n",
    "    02_bam_markdup/\n",
    "\n",
    "!squeue -u fmbuga"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2163a44-ce95-4e8a-aaa0-d66d352dfc47",
   "metadata": {},
   "source": [
    "---\n",
    "Let's have a look at one of the BAM files before and after marking duplicates using `samtools flagstats`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ba097901-746b-43c8-8a6f-fa77221a153b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1650805 + 0 in total (QC-passed reads + QC-failed reads)\n",
      "719983 + 0 primary\n",
      "480392 + 0 secondary\n",
      "450430 + 0 supplementary\n",
      "0 + 0 duplicates\n",
      "0 + 0 primary duplicates\n",
      "1588422 + 0 mapped (96.22% : N/A)\n",
      "657600 + 0 primary mapped (91.34% : N/A)\n"
     ]
    }
   ],
   "source": [
    "!samtools flagstats \\\n",
    "    --threads 28 \\\n",
    "    --verbosity 5 \\\n",
    "    01_bam_raw/FAF01132-84868110.bam | head -n 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1d59926f-9e11-4526-981a-6f2b4d0abc07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1650805 + 0 in total (QC-passed reads + QC-failed reads)\n",
      "719983 + 0 primary\n",
      "480392 + 0 secondary\n",
      "450430 + 0 supplementary\n",
      "19503 + 0 duplicates\n",
      "19503 + 0 primary duplicates\n",
      "1588422 + 0 mapped (96.22% : N/A)\n",
      "657600 + 0 primary mapped (91.34% : N/A)\n"
     ]
    }
   ],
   "source": [
    "!samtools flagstats \\\n",
    "    --threads 28 \\\n",
    "    --verbosity 5 \\\n",
    "    02_bam_markdup/FAF01132-84868110_collate_fixmate_sort_markdup.bam | head -n 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17fac6c7-8a3a-4c76-b209-d1d2218c1858",
   "metadata": {},
   "source": [
    "We see that for this BAM file, 19,503 of the 1,650,805 have been marked as duplicates (~1.2%)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393c6629-42bc-4f2f-ad19-fbe834ce350c",
   "metadata": {},
   "source": [
    "## 3. Merging the resulting 52 BAM files into 1 file (samtools)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad935387-ab41-4c2f-ab19-61efec7fe85d",
   "metadata": {},
   "source": [
    "Next, we use `samtools merge` to [combine the 52 BAM files into one](http://www.htslib.org/doc/samtools-merge.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ab971339-dcc8-429d-aab9-97b2210d3922",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!/bin/bash\n",
      "#\n",
      "#SBATCH --job-name=samtools_merge_221128\n",
      "#SBATCH --output=/home/fmbuga/tools/slurm_scripts/slurm_out_err/samtools_merge_221128_%j.out\n",
      "#SBATCH --error=/home/fmbuga/tools/slurm_scripts/slurm_out_err/samtools_merge_221128_%j.err\n",
      "#\n",
      "#SBATCH --ntasks=8\n",
      "#SBATCH --cpus-per-task=28\n",
      "#SBATCH --mem-per-cpu=4G\n",
      "\n",
      "\n",
      "datenow=$(date)\n",
      "echo $datenow\n",
      "srun hostname\n",
      "echo \"\"\n",
      "\n",
      "start=$(date +%s)\n",
      "echo \"start time: $start\"\n",
      "echo \"hostname: $HOSTNAME\"\n",
      "echo \"\"\n",
      "\n",
      "\n",
      "##################\n",
      "\n",
      "eval \"$(conda shell.bash hook)\"\n",
      "conda activate minimap2_224                      ### CONDA ENVIRONMENT NAME\n",
      "conda info\n",
      "echo \"\"\n",
      "\n",
      "call=\"samtools --version\"\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "out_bam=$1\n",
      "in_bams_dir=$2\n",
      "\n",
      "# samtools merge\n",
      "\n",
      "call=$(samtools merge \\\n",
      "        $out_bam \\\n",
      "        \"${in_bams_dir}\"*.bam \\\n",
      "        --threads 224 \\\n",
      "        --write-index \\\n",
      "        --verbosity 10) # INPUT & OUTPUT BAMS\n",
      "    \n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "echo \"samtools merge COMPLETE.\"\n",
      "echo \"\"\n",
      "    \n",
      "##################\n",
      "\n",
      "\n",
      "echo \"\"\n",
      "end=$(date +%s)\n",
      "echo \"end time: $end\"\n",
      "runtime_s=$(echo $(( end - start )))\n",
      "echo \"total run time(s): $runtime_s\"\n",
      "sec_per_min=60\n",
      "sec_per_hr=3600\n",
      "runtime_m=$(echo \"scale=2; $runtime_s / $sec_per_min;\" | bc)\n",
      "echo \"total run time(m): $runtime_m\"\n",
      "runtime_h=$(echo \"scale=2; $runtime_s / $sec_per_hr;\" | bc)\n",
      "echo \"total run time(h): $runtime_h\"\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# the SLURM job script for merging all the BAM files in the folder\n",
    "!cat ~/tools/slurm_scripts/samtools_merge_221128_slurm.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a88029-ae07-4aa2-b099-ad0bd0dc3546",
   "metadata": {},
   "source": [
    "---\n",
    "And here's the SLURM job submission:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9e07f8b1-fd47-4213-ba90-d034e8fe3aa7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitted batch job 1111544\n",
      "             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)\n",
      "           1111543     nodes     bash   fmbuga  R       9:35      1 node05\n",
      "           1111544     nodes samtools   fmbuga  R       0:10      8 node[40-47]\n"
     ]
    }
   ],
   "source": [
    "!mkdir -p 03_bam_merged\n",
    "!sbatch ~/tools/slurm_scripts/samtools_merge_221128_slurm.sh \\\n",
    "    03_bam_merged/dorado.bam \\\n",
    "    02_bam_markdup/\n",
    "\n",
    "!sleep 10\n",
    "!squeue -u fmbuga"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f4e215-a28c-41d3-827c-8bfab870d559",
   "metadata": {},
   "source": [
    "---\n",
    "Let's see what the stats look like for our final mega BAM file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "353679e6-a10e-4535-91d5-d5b119157a53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28256836 + 0 in total (QC-passed reads + QC-failed reads)\n",
      "14685190 + 0 primary\n",
      "7283682 + 0 secondary\n",
      "6287964 + 0 supplementary\n",
      "207320 + 0 duplicates\n",
      "207320 + 0 primary duplicates\n",
      "26819340 + 0 mapped (94.91% : N/A)\n",
      "13247694 + 0 primary mapped (90.21% : N/A)\n"
     ]
    }
   ],
   "source": [
    "!samtools flagstats \\\n",
    "    --threads 28 \\\n",
    "    03_bam_merged/dorado.bam | head -n 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46aa3a41-fc77-4a84-9819-d707cf67ea2f",
   "metadata": {},
   "source": [
    "So about 0.7% of the reads are duplicates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25ce9dfb-3c62-49ce-9855-e5b0551b7a57",
   "metadata": {},
   "source": [
    "## 4. Calling variants (bcftools)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab9aa1cb-1b66-4218-bbe2-9ec8a6d373e2",
   "metadata": {},
   "source": [
    "We use `bcftools mpileup` to [generate genotype likelihoods](https://samtools.github.io/bcftools/bcftools.html#mpileup) and `bcftools call` to [make the actual variant call](https://samtools.github.io/bcftools/bcftools.html#call). The prior for the expected substitution rate (`bcftools call` parameter -P) is set to a very high value to [make the calling as sensitive as possible](https://samtools.github.io/bcftools/howtos/variant-calling.html) to minimize false negatives with the hope that the model we build will filter out the majority of the false positives that result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "80c4c8e4-45c0-4150-b384-d1a1a12c05e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!/bin/bash\n",
      "#\n",
      "#SBATCH --job-name=bcftools_mpileup_call_ONT_P099_minMQ0_231114\n",
      "#SBATCH --output=/home/fmbuga/tools/slurm_scripts/slurm_out_err/bcftools_mpileup_call_ONT_P099_minMQ0_231114_%j.out\n",
      "#SBATCH --error=/home/fmbuga/tools/slurm_scripts/slurm_out_err/bcftools_mpileup_call_ONT_P099_minMQ0_231114_%j.err\n",
      "#\n",
      "#SBATCH --ntasks=1\n",
      "#SBATCH --cpus-per-task=1\n",
      "#SBATCH --mem-per-cpu=120G\n",
      "\n",
      "\n",
      "# \"...multithreading ...is currently used only for the compression of the output stream, only when --output-type is b or z\"\n",
      "# https://samtools.github.io/bcftools/bcftools.html\n",
      "# \"Use the -Ou option when piping between bcftools subcommands to speed up performance by removing unnecessary compression/decompression and VCF←→BCF conversion\"\n",
      "# for ONT https://github.com/samtools/bcftools/issues/1666\n",
      "\n",
      "# \"For ONT it may be beneficial to also run bcftools call with a higher -P, eg -P0.01 or -P 0.1\"\n",
      "\n",
      "\n",
      "datenow=$(date)\n",
      "echo $datenow\n",
      "srun hostname\n",
      "echo \"\"\n",
      "\n",
      "start=$(date +%s)\n",
      "echo \"start time: $start\"\n",
      "echo \"hostname: $HOSTNAME\"\n",
      "echo \"\"\n",
      "\n",
      "\n",
      "##################\n",
      "\n",
      "eval \"$(conda shell.bash hook)\"\n",
      "conda activate minimap2_224                      ### CONDA ENVIRONMENT NAME\n",
      "conda info\n",
      "echo \"\"\n",
      "\n",
      "\n",
      "call=\"bcftools --version\"\n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "input_bam=$1\n",
      "filename=$(basename $1 .bam)\n",
      "output_dir=$2\n",
      "min_MQ=0\n",
      "\n",
      "call=$(bcftools mpileup \\\n",
      "    -Ou \\\n",
      "    --min-MQ $min_MQ \\\n",
      "    --config ont \\\n",
      "    -f ~/tools/GIAB/GRCh38/refs/GRCh38_GIABv3.fasta \\\n",
      "    $input_bam |\\\n",
      "    bcftools call \\\n",
      "        --ploidy GRCh38 \\\n",
      "        -m -v \\\n",
      "        -P 0.99 \\\n",
      "        -o \"$output_dir$filename\"_bcftools_P_099_minMQ_\"$min_MQ\".vcf) # INPUT BAM\n",
      "    \n",
      "echo $call\n",
      "eval $call\n",
      "echo \"\"\n",
      "\n",
      "echo \"bcftools mpileup call COMPLETE.\"\n",
      "echo \"\"\n",
      "    \n",
      "##################\n",
      "\n",
      "\n",
      "echo \"\"\n",
      "end=$(date +%s)\n",
      "echo \"end time: $end\"\n",
      "runtime_s=$(echo $(( end - start )))\n",
      "echo \"total run time(s): $runtime_s\"\n",
      "sec_per_min=60\n",
      "sec_per_hr=3600\n",
      "runtime_m=$(echo \"scale=2; $runtime_s / $sec_per_min;\" | bc)\n",
      "echo \"total run time(m): $runtime_m\"\n",
      "runtime_h=$(echo \"scale=2; $runtime_s / $sec_per_hr;\" | bc)\n",
      "echo \"total run time(h): $runtime_h\"\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# SLURM job script to do the variant calling\n",
    "!cat ~/tools/slurm_scripts/bcftools_mpileup_call_ONT_P099_minMQ0_231114_slurm.sh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "695e4df4-763d-422a-ae77-dacbd7a55ae5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitted batch job 1111545\n",
      "             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)\n",
      "           1111543     nodes     bash   fmbuga  R    1:02:38      1 node05\n",
      "           1111545     nodes bcftools   fmbuga  R       0:05      1 node14\n"
     ]
    }
   ],
   "source": [
    "# submitting the variant calling job to SLURM\n",
    "!mkdir -p 04_vcf_raw\n",
    "!sbatch ~/tools/slurm_scripts/bcftools_mpileup_call_ONT_P099_minMQ0_231114_slurm.sh \\\n",
    "    03_bam_merged/dorado.bam \\\n",
    "    04_vcf_raw/\n",
    "\n",
    "!sleep 5\n",
    "!squeue -u fmbuga"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dff9989-fc45-4b2b-84a8-a7dd6feb5583",
   "metadata": {},
   "source": [
    "---\n",
    "Let's have a look at the output VCF using `bcftools stats`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c90a2ce2-bb41-4e3e-93ab-5883cce0c34b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t6484598\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t6484598\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t0\n",
      "SN\t0\tnumber of others:\t0\n",
      "SN\t0\tnumber of multiallelic sites:\t16587\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t16587\n"
     ]
    }
   ],
   "source": [
    "# check stats of the VCF file\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    04_vcf_raw/dorado_bcftools_P_099_minMQ_0.vcf | grep ^SN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef0129b9-ca2b-472e-af22-2d4bb0191dba",
   "metadata": {},
   "source": [
    "Over 6 million SNPs were called. Considering humans on average only have about 1 SNP per 1000 base pairs, we expect about half of these calls to be false positives since the human genome has ~3.2 billion base pairs. Note that indel calling is turned off which is why 0 indels are recorded."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ec145a-4e02-4846-9c0d-5b27a04cb397",
   "metadata": {},
   "source": [
    "We can compare it to the number of SNPs in the GIAB truth-set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8272c534-9124-4428-85ec-166a4e7daba0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t3893341\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t3358502\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t536489\n",
      "SN\t0\tnumber of others:\t90\n",
      "SN\t0\tnumber of multiallelic sites:\t42393\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t1033\n"
     ]
    }
   ],
   "source": [
    "# check stats of the GIAB VCF file\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    ~/tools/GIAB/GRCh38/HG001_NA12878/v4_2_1/HG001_GRCh38_v4_2_1.vcf | grep ^SN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d96e00a-d589-4514-9185-661cae159c05",
   "metadata": {},
   "source": [
    "## 5. Intersect with Genome-in-a-Bottle variant calls (bedtools)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4492a6ce-3d3a-4e6e-bc73-5d2a4d1bb370",
   "metadata": {},
   "source": [
    "Next, we use `bedtools intersect` to overlap our variant callset with a BED file of the regions covered by NIST's Genome in a Bottle (GIAB) high-confidence variant calls of the same genome so we can compare our variant calls with a gold standard. Furthermore, the comparison will be what is used to create the labels (true variant vs artifact) for the model for artifactual variants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "08a13668-ae58-40ac-a1e4-723a12998c70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "***** WARNING: File 04_vcf_raw/dorado_bcftools_P_099_minMQ_0.vcf has inconsistent naming convention for record:\n",
      "MAP2K3_chr17_22578583_22605165\t151\t.\tA\tG\t189.708\t.\tDP=24;VDB=1.66594e-13;SGB=-0.692352;RPBZ=1.03992;MQBZ=-0.315869;MQSBZ=1.0942;BQBZ=-0.554502;SCBZ=1.11849;FS=0;MQ0F=0;AC=2;AN=2;DP4=1,0,13,8;MQ=57\tGT:PL\t1/1:187,40,0\n",
      "\n",
      "***** WARNING: File 04_vcf_raw/dorado_bcftools_P_099_minMQ_0.vcf has inconsistent naming convention for record:\n",
      "MAP2K3_chr17_22578583_22605165\t151\t.\tA\tG\t189.708\t.\tDP=24;VDB=1.66594e-13;SGB=-0.692352;RPBZ=1.03992;MQBZ=-0.315869;MQSBZ=1.0942;BQBZ=-0.554502;SCBZ=1.11849;FS=0;MQ0F=0;AC=2;AN=2;DP4=1,0,13,8;MQ=57\tGT:PL\t1/1:187,40,0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# intersect dorado vcf with GIAB hi-conf BED\n",
    "!mkdir -p 05_vcf_isec_GIAB\n",
    "!bedtools intersect \\\n",
    "    -a 04_vcf_raw/dorado_bcftools_P_099_minMQ_0.vcf \\\n",
    "    -b ~/tools/GIAB/GRCh38/HG001_NA12878/v4_2_1/HG001_GRCh38_1_22_v4.2.1_benchmark.bed \\\n",
    "    -header > 05_vcf_isec_GIAB/dorado_bcftools_P_099_isec.vcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5644c7cd-a46a-4a49-93b3-0f4d5ffcd1ef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t3773914\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t3773914\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t0\n",
      "SN\t0\tnumber of others:\t0\n",
      "SN\t0\tnumber of multiallelic sites:\t4126\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t4126\n"
     ]
    }
   ],
   "source": [
    "# check stats of the isec VCF file\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    05_vcf_isec_GIAB/dorado_bcftools_P_099_isec.vcf | grep ^SN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91d1e0c2-2865-4ff3-b16f-c61acb270f79",
   "metadata": {},
   "source": [
    "Filtering our VCF file to only look at the regions where NIST GIAB has high-confidence variant calls leaves us with ~3.8 million SNPs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "914d132a-a004-45e5-8b13-301c50fe3007",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t3356518\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t3356518\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t0\n",
      "SN\t0\tnumber of others:\t0\n",
      "SN\t0\tnumber of multiallelic sites:\t1033\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t1033\n"
     ]
    }
   ],
   "source": [
    "# filter GIAB VCF to keep SNPs only\n",
    "!bcftools view \\\n",
    "    --include 'TYPE=\"snp\"' \\\n",
    "    --threads 28 \\\n",
    "    ~/tools/GIAB/GRCh38/HG001_NA12878/v4_2_1/HG001_GRCh38_v4_2_1.vcf \\\n",
    "    > ~/tools/GIAB/GRCh38/HG001_NA12878/v4_2_1/HG001_GRCh38_v4_2_1_SNVs_ONLY.vcf\n",
    "\n",
    "# check stats of the snps only VCF file\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    ~/tools/GIAB/GRCh38/HG001_NA12878/v4_2_1/HG001_GRCh38_v4_2_1_SNVs_ONLY.vcf \\\n",
    "    | grep ^SN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2797e552-c636-4081-89c3-a947d4fd5b5d",
   "metadata": {},
   "source": [
    "### generate 3 VCFs (true positive, false positive, false negative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "f8f940e1-ae5e-41c8-bcc1-68a90c36419f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compress and index the two vcf files to allow reading by bcftools isec\n",
    "!bgzip tmp_dir/dorado_sorted.vcf\n",
    "!tabix -p vcf tmp_dir/dorado_sorted.vcf.gz\n",
    "!bgzip tmp_dir/GIAB_sorted.vcf\n",
    "!tabix -p vcf tmp_dir/GIAB_sorted.vcf.gz\n",
    "# Create intersection and complements of two sets saving the output in dir/*\n",
    "!bcftools isec \\\n",
    "    --threads 28 \\\n",
    "    tmp_dir/dorado_sorted.vcf.gz tmp_dir/GIAB_sorted.vcf.gz \\\n",
    "    -p tmp_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "7ebbefbb-cfbb-416a-a461-3af2c998d75e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This file was produced by vcfisec.\n",
      "The command line was:\tbcftools isec  --threads 28 -p tmp_dir tmp_dir/dorado_sorted.vcf.gz tmp_dir/GIAB_sorted.vcf.gz\n",
      "\n",
      "Using the following file names:\n",
      "tmp_dir/0000.vcf\tfor records private to\ttmp_dir/dorado_sorted.vcf.gz\n",
      "tmp_dir/0001.vcf\tfor records private to\ttmp_dir/GIAB_sorted.vcf.gz\n",
      "tmp_dir/0002.vcf\tfor records from tmp_dir/dorado_sorted.vcf.gz shared by both\ttmp_dir/dorado_sorted.vcf.gz tmp_dir/GIAB_sorted.vcf.gz\n",
      "tmp_dir/0003.vcf\tfor records from tmp_dir/GIAB_sorted.vcf.gz shared by both\ttmp_dir/dorado_sorted.vcf.gz tmp_dir/GIAB_sorted.vcf.gz\n"
     ]
    }
   ],
   "source": [
    "!cat tmp_dir/README.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a9930169-38f5-403e-91bc-6dc956c4757a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make new folder for 3 VCFs: TP, FP, FN (true positive, false positive, false negative)\n",
    "!mkdir 06_vcf_TP_FP_FN\n",
    "\n",
    "# copy the 3 vcfs to this folder\n",
    "!cp tmp_dir/0000.vcf 06_vcf_TP_FP_FN/FP_snvs.vcf\n",
    "!cp tmp_dir/0001.vcf 06_vcf_TP_FP_FN/FN_snvs.vcf\n",
    "!cp tmp_dir/0002.vcf 06_vcf_TP_FP_FN/TP_snvs.vcf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ebf9db74-1e18-45f2-add6-7c9cd4a81df6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t3236774\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t3236774\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t0\n",
      "SN\t0\tnumber of others:\t0\n",
      "SN\t0\tnumber of multiallelic sites:\t926\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t926\n"
     ]
    }
   ],
   "source": [
    "# true positives (correct variants)\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    06_vcf_TP_FP_FN/TP_snvs.vcf | grep ^SN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "119d86db-f89f-4ab5-bd7c-fce132b7e1f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t537140\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t537140\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t0\n",
      "SN\t0\tnumber of others:\t0\n",
      "SN\t0\tnumber of multiallelic sites:\t3200\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t3200\n"
     ]
    }
   ],
   "source": [
    "# false positives\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    06_vcf_TP_FP_FN/FP_snvs.vcf | grep ^SN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "55d293a8-5cc9-40be-bfc4-663577ed15b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SN\t0\tnumber of samples:\t1\n",
      "SN\t0\tnumber of records:\t119744\n",
      "SN\t0\tnumber of no-ALTs:\t0\n",
      "SN\t0\tnumber of SNPs:\t119744\n",
      "SN\t0\tnumber of MNPs:\t0\n",
      "SN\t0\tnumber of indels:\t0\n",
      "SN\t0\tnumber of others:\t0\n",
      "SN\t0\tnumber of multiallelic sites:\t107\n",
      "SN\t0\tnumber of multiallelic SNP sites:\t107\n"
     ]
    }
   ],
   "source": [
    "# false negatives\n",
    "!bcftools stats \\\n",
    "    --threads 28 \\\n",
    "    06_vcf_TP_FP_FN/FN_snvs.vcf | grep ^SN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aa32e86-d26e-4b2d-997b-722ee3142a79",
   "metadata": {},
   "source": [
    "---\n",
    "### NEXT:\n",
    "- generate features (201bps)\n",
    "- generate labels (0 - not artifact, 1 - is artifact & 0 - TP, 1 - FP, 2 - FN): one-hot encode later in numpy\n",
    "- pre-process for deep learning\n",
    "- fit deep learning models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72ba7d08-5a63-4e28-be72-0477565834a2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}