NGS_sequence_alignment_command_line.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Bioinformatics for Genome Medicine Workshop\n",
        "\n",
        "###*Prof. Jason Wong, School of Biomedical Sciences, HKU*\n",
        "\n",
        "####The objective of this lecture is to demonstrate how next-generation sequencing data is aligned to the reference genome sequence and how variant calling can be performed.\n"
      ],
      "metadata": {
        "id": "o_lO85QksbvE"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## *** Package installation and downloads for workshop (~ 5 minutes)\n",
        "\n",
        "1.   conda (for simple installation of packages)\n",
        "2.   FastQC (for reads quality check)\n",
        "3.   bwa (tools for sequence alignment)\n",
        "4.   samtools (tools for processing sam & bam files)  \n",
        "\n",
        "**IMPORTANT：Every time you connect to Google Colab, you have to perform these set up steps again.**"
      ],
      "metadata": {
        "id": "7gnJkcA9tTBn"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# set working pathway to your own google drive doc (~ 1 min)\n",
        "from google.colab import drive\n",
        "drive.mount('/content/gdrive')"
      ],
      "metadata": {
        "id": "AOrw7MhPztUW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# install conda (~ 1 min). There will be a message saying that the session has crashed, but don't worry about this. This is due to the session restarting following conda installation.\n",
        "!pip install -q condacolab\n",
        "import condacolab#\n",
        "condacolab.install()"
      ],
      "metadata": {
        "id": "UJMUyjxntquy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# install fastqc (~ 2 mins)\n",
        "!conda install -c bioconda fastqc"
      ],
      "metadata": {
        "id": "F5J8UCIn7i7H"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# install bwa (~ 1 mins)\n",
        "!conda install -c bioconda bwa"
      ],
      "metadata": {
        "id": "hkiiT0YRz4S2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Install igv-notebook\n",
        "!pip install igv-notebook"
      ],
      "metadata": {
        "id": "B_5lJh21K7Ya"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# install lofreq (~ 1 mins)\n",
        "!conda install -c bioconda lofreq"
      ],
      "metadata": {
        "id": "qSeJa-L2ORU2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# install samtools (~1 mins)\n",
        "!conda install -c bioconda samtools"
      ],
      "metadata": {
        "id": "oMH5G73bz_CI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# install bedtools (~ 1 mins)\n",
        "!apt install bedtools"
      ],
      "metadata": {
        "id": "UR-CWbJUzWwU"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Set working directory\n",
        "\n",
        "By default, the working directory will be My Drive/PB_course"
      ],
      "metadata": {
        "id": "jQ4sstuDsBO1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "try:\n",
        "  os.mkdir(\"/content/gdrive/My Drive/PB_course\")         # change this path if necessary\n",
        "except FileExistsError:\n",
        "  print(\"directory already exist. OK to continue\")\n",
        "os.chdir(\"/content/gdrive/My Drive/PB_course\")"
      ],
      "metadata": {
        "id": "p4_KlYtbrzzF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Download ready prepared files for analysis."
      ],
      "metadata": {
        "id": "TPs7Tb-tyxc_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#download reference sequence\n",
        "# double check that we are in the right directory (~ 30s)\n",
        "import os\n",
        "os.chdir(\"/content/gdrive/MyDrive/PB_course\")                     # change this path if necessary\n",
        "\n",
        "import os\n",
        "if os.path.isfile(\"/content/gdrive/MyDrive/PB_course/DB_trunc/chr2.fa\"):    # check if the file exist\n",
        "  print(\"reference file already exit, OK to continue.\")\n",
        "else:\n",
        "  !wget -O DB_trunc.zip \"https://drive.usercontent.google.com/download?id=1zF8BTkR90kEr23b3NecES3u8akbSfIw9&export=download&authuser=0&confirm=t&uuid=ec27dd21-2518-4a97-8547-7cada1fde88e&at=APZUnTXxF4K7otx7GhcCftn0sRfA%3A1709043759277\"  # download if file not exist\n",
        "  #unzip fasta file\n",
        "  !unzip DB_trunc.zip\n",
        "  !rm DB_trunc.zip\n",
        "\n",
        "!ls -l ./DB_trunc/"
      ],
      "metadata": {
        "id": "OwoP1gu3GWLX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#download sample sequences\n",
        "import os\n",
        "os.chdir(\"/content/gdrive/My Drive/PB_course/\")\n",
        "if os.path.isfile(\"/content/gdrive/MyDrive/PB_course/Datasets/ChIP-seq_H3K27ac_example.fq.gz\"):    # check if the file exist\n",
        "  print(\"file already exit, OK to continue.\")\n",
        "else:\n",
        " !wget -O Datasets.zip https://github.com/jasonwong-lab/HKU-Practical-Bioinformatics/raw/main/files/Datasets.zip    # download necessary file\n",
        " !unzip -o Datasets.zip   #unzip file\n",
        " !rm Datasets.zip"
      ],
      "metadata": {
        "id": "MsTqxxiP0AzQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## NGS_alignment command line\n",
        "\n",
        "1.1 Quality control\n",
        "\n",
        "1.2. Burrows–Wheeler Aligner\n",
        "\n",
        "1.3. Working with SAM files\n",
        "\n",
        "1.4. Align the WGS paired-end file"
      ],
      "metadata": {
        "id": "XB9TrCTC7Ncy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Check what datasets we have downloaded\n",
        "%cd /content/gdrive/MyDrive/PB_course\n",
        "!ls -l\n",
        "\n",
        "%cd Datasets/\n",
        "!ls -l"
      ],
      "metadata": {
        "id": "4dYhhxiG4lMM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# To look at the ChIP-seq fastq file type (head -n 12 to print first 12 lines):\n",
        "!zcat < ChIP-seq_H3K27ac_example.fq.gz |head -n 12"
      ],
      "metadata": {
        "id": "EMrO9wOVnNH2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# How many reads are there in the file? Type your code below:\n",
        "!zcat < ChIP-seq_H3K27ac_example.fq.gz | wc -l"
      ],
      "metadata": {
        "id": "t4g-To8mPw8V"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Quality control"
      ],
      "metadata": {
        "id": "pxBspjPp7Dms"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "## Run FastQC\n",
        "!fastqc ChIP-seq_H3K27ac_example.fq.gz"
      ],
      "metadata": {
        "id": "f82vAykr7Fie"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# check resulting html file\n",
        "!ls\n",
        "\n",
        "#download file and check html file with your local browser\n",
        "from google.colab import files\n",
        "files.download('ChIP-seq_H3K27ac_example_fastqc.html')"
      ],
      "metadata": {
        "id": "7v4bH04sZtn5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "you can also download the file **ChIP-seq_H3K27ac_example_fastqc.html** from \"gdrive/MyDrive/PB_course/Datasets\" on the left side"
      ],
      "metadata": {
        "id": "Gpt-qEVO81Bs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Burrows–Wheeler Aligner"
      ],
      "metadata": {
        "id": "mwdi1hJ59Dtu"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Look at the options for bwa and bwa mem\n",
        "!bwa"
      ],
      "metadata": {
        "id": "yE-r2cjw8NfS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!bwa mem"
      ],
      "metadata": {
        "id": "EQCq422B9HbZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# let's take a look at the database files\n",
        "%cd /content/gdrive/My Drive/PB_course/DB_trunc\n",
        "\n",
        "!ls -l"
      ],
      "metadata": {
        "id": "AuteGS6HdjLs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Get ready to run BWA. First Go into Datasets directory\n",
        "%cd /content/gdrive/My Drive/PB_course/Datasets"
      ],
      "metadata": {
        "id": "aFb8h3SZNyjY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# make an directory to store output file\n",
        "!mkdir BAM"
      ],
      "metadata": {
        "id": "36JkHL9IeUkd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Do sequence alignment\n",
        "!bwa mem ../DB_trunc/chr2.fa ./ChIP-seq_H3K27ac_example.fq.gz > ./BAM/ChIP-seq_H3K27ac_example.sam"
      ],
      "metadata": {
        "id": "8EBR94fLMsQI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# check the result\n",
        "!head -n 20 ./BAM/ChIP-seq_H3K27ac_example.sam"
      ],
      "metadata": {
        "id": "8oqpLEwbNsCd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Working with SAM files"
      ],
      "metadata": {
        "id": "BEQNtRCyOSOk"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Check the samtools command\n",
        "!samtools"
      ],
      "metadata": {
        "id": "KLBxHt5COTho"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Check out some stats about our aligned file\n",
        "!samtools flagstat ./BAM/ChIP-seq_H3K27ac_example.sam"
      ],
      "metadata": {
        "id": "6Pgw5h7JOXMO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Prepare the file for viewing on genome browser:\n",
        "#Step 1 – Convert SAM to BAM\n",
        "!samtools view -b ./BAM/ChIP-seq_H3K27ac_example.sam > ./BAM/ChIP-seq_H3K27ac_example.bam"
      ],
      "metadata": {
        "id": "A2PFyYk3OjWJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 2 – sort BAM file\n",
        "!samtools sort ./BAM/ChIP-seq_H3K27ac_example.bam > ./BAM/ChIP-seq_H3K27ac_example_sorted.bam"
      ],
      "metadata": {
        "id": "IguX9XqQOwf9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 3 – index BAM file\n",
        "!samtools index ./BAM/ChIP-seq_H3K27ac_example_sorted.bam"
      ],
      "metadata": {
        "id": "D92KppCUFX9R"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Look at the files that we have created\n",
        "!ls -l ./BAM/"
      ],
      "metadata": {
        "id": "kmJrDMuPF8Fn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### IGV browser"
      ],
      "metadata": {
        "id": "Jo8cb29WK2pC"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Load track from local paths\n",
        "import igv_notebook\n",
        "\n",
        "igv_notebook.init()\n",
        "\n",
        "b = igv_notebook.Browser(\n",
        "    {\n",
        "        \"genome\": \"hg38\",\n",
        "        \"locus\": \"chr2:47,782,081-47,807,953\"\n",
        "    }\n",
        ")\n",
        "\n",
        "b.load_track(\n",
        "    {\n",
        "        \"name\": \"ChIP-seq_H3K27ac\",\n",
        "        \"path\": \"./BAM/ChIP-seq_H3K27ac_example_sorted.bam\",\n",
        "        \"indexPath\": \"./BAM/ChIP-seq_H3K27ac_example_sorted.bam.bai\",\n",
        "        \"format\": \"bam\",\n",
        "        \"type\": \"alignment\"\n",
        "    })"
      ],
      "metadata": {
        "id": "a_tqSvF8TTVp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Align the WGS paired-end file"
      ],
      "metadata": {
        "id": "nRQZFpflGGhN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Check that we have the WGS fastq files\n",
        "%cd /content/gdrive/My Drive/PB_course/Datasets\n",
        "!ls -l"
      ],
      "metadata": {
        "id": "jXpX831tnQp9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# If you don't have WGS_example_1.fq.gz or WGS_example_2.fq.gz, run this cell to redownload the files\n",
        "import os\n",
        "os.chdir(\"/content/gdrive/My Drive/PB_course/\")\n",
        "\n",
        "import os\n",
        "if os.path.isfile(\"/content/gdrive/MyDrive/PB_course/Datasets/WGS_example_1.fq.gz\"):    # check if the file exist\n",
        "  print(\"reference file already exit, OK to continue.\")\n",
        "else:\n",
        " !wget -O Datasets.zip https://github.com/jasonwong-lab/HKU-Practical-Bioinformatics/raw/main/files/Datasets.zip\n",
        " !unzip -o Datasets.zip   #unzip file\n",
        " !rm Datasets.zip\n",
        "\n",
        "#current directory should still be ~/Datasets\n",
        "%cd Datasets/\n",
        "!ls -l"
      ],
      "metadata": {
        "id": "vtckkwSm7Hw8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#look at our work directory and paired-end file:\n",
        "!zcat < WGS_example_1.fq.gz |head -1\n",
        "!zcat < WGS_example_2.fq.gz |head -1"
      ],
      "metadata": {
        "id": "l6l-vb9cGILx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#align the WGS paired-end file:( ~ 1 mins)\n",
        "#Step 1 – Aligning paired-end file using bwa mem:\n",
        "!bwa mem ../DB_trunc/chr2.fa ./WGS_example_1.fq.gz ./WGS_example_2.fq.gz > ./BAM/WGS_example.sam\n"
      ],
      "metadata": {
        "id": "Eq9VuU3BGl1U"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 2 – Output sorted BAM, this time use piping to skip one step:\n",
        "!samtools view -b ./BAM/WGS_example.sam | samtools sort >./BAM/WGS_example_sorted.bam"
      ],
      "metadata": {
        "id": "162y1PgfHAQP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 3 – Index sorted bam file\n",
        "!samtools index ./BAM/WGS_example_sorted.bam"
      ],
      "metadata": {
        "id": "Xbws6jwVHXaQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Check stats about the aligned WGS file\n",
        "!samtools flagstat ./BAM/WGS_example_sorted.bam"
      ],
      "metadata": {
        "id": "YjwwnUF2HcQo"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#IGV\n",
        "#Load track from local paths\n",
        "import igv_notebook\n",
        "\n",
        "igv_notebook.init()\n",
        "\n",
        "b = igv_notebook.Browser(\n",
        "    {\n",
        "        \"genome\": \"hg38\",\n",
        "        \"locus\": \"chr2:47,782,081-47,807,953\"\n",
        "    }\n",
        ")\n",
        "\n",
        "b.load_track(\n",
        "    {\n",
        "        \"name\": \"ChIP-seq H3K27ac\",\n",
        "        \"path\": \"./BAM/ChIP-seq_H3K27ac_example_sorted.bam\",\n",
        "        \"indexPath\": \"./BAM/ChIP-seq_H3K27ac_example_sorted.bam.bai\",\n",
        "        \"format\": \"bam\",\n",
        "        \"type\": \"alignment\"\n",
        "    })\n",
        "\n",
        "b.load_track(\n",
        "    {\n",
        "        \"name\": \"WGS\",\n",
        "        \"path\": \"./BAM/WGS_example_sorted.bam\",\n",
        "        \"indexPath\": \"./BAM/WGS_example_sorted.bam.bai\",\n",
        "        \"format\": \"bam\",\n",
        "        \"type\": \"alignment\"\n",
        "    })"
      ],
      "metadata": {
        "id": "evZ_KRIwca9N"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Aligned and do variant calling on WXS data"
      ],
      "metadata": {
        "id": "cGQTLct-JuOc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Run this cell to download the WXS files\n",
        "import os\n",
        "os.chdir(\"/content/gdrive/My Drive/PB_course/Datasets\")\n",
        "\n",
        "import os\n",
        "if os.path.isfile(\"/content/gdrive/MyDrive/PB_course/Datasets/WXS_example_cancer_1.fq.gz\"):    # check if the file exist\n",
        "  print(\"WXS file already exit, OK to continue.\")\n",
        "else:\n",
        " !wget -O Datasets_WXS.zip \"https://drive.usercontent.google.com/download?id=1pivDkbhXEYN57f-1pduYR5E2TtwvZGsx&export=download&authuser=0&confirm=t&uuid=6ccaeb5f-3091-48da-ac04-3afc5256c499&at=APZUnTXUz8J8VXoZgzVPVXDUlJHz%3A1709081319052\"\n",
        " !unzip -o Datasets_WXS.zip   #unzip file\n",
        " !rm Datasets_WXS.zip\n",
        "\n",
        "\n",
        "#Check what files we have now\n",
        "%cd /content/gdrive/MyDrive/PB_course/Datasets/\n",
        "!ls -l"
      ],
      "metadata": {
        "id": "ZYEVqSDXJ04s"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#align the WXS paired-end file:( ~ 3 mins)\n",
        "#Step 1 – Aligning normal and cancer paired-end file using bwa mem:\n",
        "!bwa mem ../DB_trunc/chr2.fa ./WXS_example_normal_1.fq.gz ./WXS_example_normal_2.fq.gz > ./BAM/WXS_example_normal.sam\n",
        "!bwa mem ../DB_trunc/chr2.fa ./WXS_example_cancer_1.fq.gz ./WXS_example_cancer_2.fq.gz > ./BAM/WXS_example_cancer.sam"
      ],
      "metadata": {
        "id": "UcbkYMjcJ2VT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 2 – Output sorted BAM, this time use piping to skip one step:\n",
        "!samtools view -b ./BAM/WXS_example_normal.sam | samtools sort >./BAM/WXS_example_normal_sorted.bam\n",
        "!samtools view -b ./BAM/WXS_example_cancer.sam | samtools sort >./BAM/WXS_example_cancer_sorted.bam"
      ],
      "metadata": {
        "id": "7FjZsJ3oNV04"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 3 – Index sorted bam file\n",
        "!samtools index ./BAM/WXS_example_normal_sorted.bam\n",
        "!samtools index ./BAM/WXS_example_cancer_sorted.bam"
      ],
      "metadata": {
        "id": "q2WwwVlBNimI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#IGV\n",
        "#Load track from local paths\n",
        "import igv_notebook\n",
        "\n",
        "igv_notebook.init()\n",
        "\n",
        "b = igv_notebook.Browser(\n",
        "    {\n",
        "        \"genome\": \"hg38\",\n",
        "        \"locus\": \"chr2:47,782,081-47,807,953\"\n",
        "    }\n",
        ")\n",
        "\n",
        "b.load_track(\n",
        "    {\n",
        "        \"name\": \"WXS Normal\",\n",
        "        \"path\": \"./BAM/WXS_example_normal_sorted.bam\",\n",
        "        \"indexPath\": \"./BAM/WXS_example_normal_sorted.bam.bai\",\n",
        "        \"format\": \"bam\",\n",
        "        \"type\": \"alignment\"\n",
        "    })\n",
        "\n",
        "b.load_track(\n",
        "    {\n",
        "        \"name\": \"WXS Cancer\",\n",
        "        \"path\": \"./BAM/WXS_example_cancer_sorted.bam\",\n",
        "        \"indexPath\": \"./BAM/WXS_example_cancer_sorted.bam.bai\",\n",
        "        \"format\": \"bam\",\n",
        "        \"type\": \"alignment\"\n",
        "    })"
      ],
      "metadata": {
        "id": "A2sJseyDNrFb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Variant calling - call normal and tumour separately then subtract"
      ],
      "metadata": {
        "id": "CDBm7cTXOy9t"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!lofreq"
      ],
      "metadata": {
        "id": "f0GVEhjDjVG9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Variant calling using LoFreq\n",
        "!lofreq call"
      ],
      "metadata": {
        "id": "jSL4JPJHYjN6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 1 – Create directory to store VCF files. Make sure we are in the Datasets directory.\n",
        "%cd /content/gdrive/My Drive/PB_course/Datasets/\n",
        "!mkdir VCF/"
      ],
      "metadata": {
        "id": "AsDq2fIXcgK9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 2 – Run LoFreq call\n",
        "!rm VCF/*.vcf\n",
        "!lofreq call ./BAM/WXS_example_normal_sorted.bam -o ./VCF/WXS_normal_nofilter.vcf -f ../DB_trunc/chr2.fa --verbose --no-default-filter\n",
        "!lofreq call ./BAM/WXS_example_cancer_sorted.bam -o ./VCF/WXS_cancer_nofilter.vcf -f ../DB_trunc/chr2.fa --verbose --no-default-filter"
      ],
      "metadata": {
        "id": "CnPshIZ6c5A4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Step 3 - Run LoFreq filter\n",
        "!lofreq filter --print-all -i ./VCF/WXS_normal_nofilter.vcf -o ./VCF/WXS_normal_filter.vcf\n",
        "!lofreq filter --print-all -i ./VCF/WXS_cancer_nofilter.vcf -o ./VCF/WXS_cancer_filter.vcf"
      ],
      "metadata": {
        "id": "s5I3XJiVZCJd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#check result:\n",
        "!head -n 30 ./VCF/WXS_normal_filter.vcf"
      ],
      "metadata": {
        "id": "hmWWDqzTdfcS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Rerun lofreq filter again to only print PASS variants\n",
        "!lofreq filter -i ./VCF/WXS_normal_nofilter.vcf -o ./VCF/WXS_normal_PASS.vcf\n",
        "!lofreq filter -i ./VCF/WXS_cancer_nofilter.vcf -o ./VCF/WXS_cancer_PASS.vcf"
      ],
      "metadata": {
        "id": "hdedw_CIPHog"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Count the number of lines in WXS_normal_PASS.vcf and WXS_cancer_PASS.vcf Remove the header by grep\n",
        "!grep -v '#' ./VCF/WXS_normal_PASS.vcf | wc -l\n",
        "!grep -v '#' ./VCF/WXS_cancer_PASS.vcf | wc -l"
      ],
      "metadata": {
        "id": "iXoHnbA9PXx8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Get a list of mutations only present in the cancer sample\n",
        "!subtractBed -a ./VCF/WXS_cancer_PASS.vcf -b ./VCF/WXS_normal_PASS.vcf > ./VCF/WXS_cancer_only_PASS.vcf"
      ],
      "metadata": {
        "id": "II8x71SkspYm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Count the number of lines in WXS_cancer_only_PASS.vcf Remove the header by grep\n",
        "!grep -v '#' ./VCF/WXS_cancer_only_PASS.vcf | wc -l"
      ],
      "metadata": {
        "id": "xpydG1o9tJLw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Check the functional impact of variants\n",
        "\n",
        "VEP online: http://useast.ensembl.org/Tools/VEP\n",
        "\n"
      ],
      "metadata": {
        "id": "4qBdbOELs-Tq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#download the germline (WXS_normal_PASS.vcf) variants for upload to VEP\n",
        "\n",
        "from google.colab import files\n",
        "files.download('./VCF/WXS_normal_PASS.vcf')\n"
      ],
      "metadata": {
        "id": "nOkSnrm7tFTY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#download the somatic (WXS_cancer_only_PASS.vcf) variants for upload to VEP\n",
        "files.download('./VCF/WXS_cancer_only_PASS.vcf')"
      ],
      "metadata": {
        "id": "SlxtdMNXcxWM"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}