NGS_exercise.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": [
        "## NGS alignment exercise\n",
        "\n",
        "WXS_example_* and WXS_example_cancer_* are next-generation sequencing (NGS) data from the blood and cancer tissue of a colorectal cancer patient^, respectively.\n",
        "\n",
        "Next week, we will learn to call single nucleotide variants from these files. Our objective this week is just to align them. Please answer the following questions as you proceed to align the file and note down your Linux commands used in a separate document.\n",
        "\n",
        "(1)\tWhy are there two files for each sample (i.e. WXS_example_1.fq.gz and WXS_example_2.fq.gz)?\n",
        "\n",
        "(2)\tHow long are the reads in WXS_sample_1.fq.gz?\n",
        "\n",
        "(3)\tHow many reads are in each sample?\n",
        "\n",
        "(4)\tNow align the file. How many properly paired reads were aligned in WXS_example and WXS_example_cancer, respectively?\n",
        "\n",
        "(5)\tView the files in IGV. What regions of DNA are most reads covering?\n",
        "\n",
        "(6)\tWhat is the coverage at the transcription start site of MSH6?\n"
      ],
      "metadata": {
        "id": "VnGEHRXIzs6y"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First set up Colab session to do exercise:"
      ],
      "metadata": {
        "id": "tuFXbU-r0JSa"
      }
    },
    {
      "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')                         # if using for the first time, you be requested to grant permission to link your Google Drive\n",
        "\n",
        "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": "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 samtools (~1 min)\n",
        "!conda install -c bioconda samtools"
      ],
      "metadata": {
        "id": "oMH5G73bz_CI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Install bwa (~ 1 min)\n",
        "!conda install -c bioconda bwa"
      ],
      "metadata": {
        "id": "hkiiT0YRz4S2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Install igv-notebook (< 1 min)\n",
        "!pip install igv-notebook"
      ],
      "metadata": {
        "id": "B_5lJh21K7Ya"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "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 exist, OK to continue.\")\n",
        "else:\n",
        "  !pip install gdown\n",
        "  !gdown https://drive.google.com/uc?id=1aRJVznjy5WLQ5Dc0DT9c6NiXw64HdoKr # download if file not exist\n",
        "  # unzip fasta file\n",
        "  !unzip DB_trunc.zip\n",
        "  # remove the zip file after extraction\n",
        "  !rm DB_trunc.zip\n",
        "\n",
        "!ls -l ./DB_trunc/"
      ],
      "metadata": {
        "id": "OwoP1gu3GWLX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Run this cell to download the WXS 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/WXS_example_cancer_1.fq.gz\"):    # check if the file exist\n",
        "  print(\"file already exist, OK to continue.\")\n",
        "else:\n",
        " !wget -O Datasets_GXS.zip https://github.com/jasonwong-lab/HKU-Practical-Bioinformatics/raw/main/files/Datasets_GXS.zip\n",
        " !unzip -o Datasets_GXS.zip   # unzip file\n",
        " !rm Datasets_GXS.zip\n",
        " !wget -O Datasets/WXS_example_cancer_1.fq.gz https://github.com/jasonwong-lab/HKU-Practical-Bioinformatics/raw/main/files/WXS_example_cancer_1.fq.gz\n",
        " !wget -O Datasets/WXS_example_cancer_2.fq.gz https://github.com/jasonwong-lab/HKU-Practical-Bioinformatics/raw/main/files/WXS_example_cancer_2.fq.gz\n",
        "\n",
        "\n",
        "# Check what files we have now\n",
        "%cd /content/gdrive/MyDrive/PB_course/Datasets/\n",
        "!ls -l"
      ],
      "metadata": {
        "id": "F3BzFQvlO-W1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now can start the exercise.\n",
        "The answers are on Moodle.\n",
        "Type in the code yourself to get familiar with Colab."
      ],
      "metadata": {
        "id": "OZsZQ6_90hhc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q1 Why are there two files for each sample (i.e. WXS_example_1.fq.gz and WXS_example_2.fq.gz)?**"
      ],
      "metadata": {
        "id": "nd83lrz9_qRX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [],
      "metadata": {
        "id": "WnsIn-R1_73b"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q2 How long are the reads in WXS_sample_1.fq.gz?**  \n",
        "hint: zcat"
      ],
      "metadata": {
        "id": "pmtZC34u_-Cm"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "hupifFgYAEVM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q3 How many reads are in each sample?**  \n",
        "hint: zcat + wc -l"
      ],
      "metadata": {
        "id": "g9_CAs_5ANlh"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "Qzz_b8pQASBj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q4 Now align the file. How many properly paired reads were aligned in WXS_example and WXS_example_cancer, respectively?**  \n",
        "hint: bwa mem + samtools flagstat"
      ],
      "metadata": {
        "id": "mQiY_opaAc17"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "Sz8LuObHAe4r"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q5 View the files in IGV. What regions of DNA are most reads covering?**"
      ],
      "metadata": {
        "id": "g9-dOIyfBMXI"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# First prepare files for IGV\n"
      ],
      "metadata": {
        "id": "OlvQIRp2BO1C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Then load the tracks in IGV"
      ],
      "metadata": {
        "id": "aQu1Mn8dBU7w"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q6 What is the coverage at the transcription start site of MSH6?**  \n",
        "hint: use the IGV app to find answer"
      ],
      "metadata": {
        "id": "_rY-uv_VBaPP"
      }
    },
    {
      "cell_type": "markdown",
      "source": [],
      "metadata": {
        "id": "Ikk7eQXuBg7I"
      }
    }
  ]
}