tutorial_2-metilene.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "036d0916-f415-4abe-80c3-403446d90ff5",
   "metadata": {},
   "source": [
    "# The WGBS data analysis tutorial 2 - DMR identification\n",
    "\n",
    "\n",
    "This tutorial will introduce the DMR (differential methylated region) detection using metilene. __[metilene](https://www.bioinf.uni-leipzig.de/Software/metilene/)__ is a software tool to annotate differentially methylated regions (DMRs) and differentially methylated CpG sites (DMCs) from Methyl-seq data. Please see the [metilene user guide](https://www.bioinf.uni-leipzig.de/Software/metilene/Manual/) for more details.  \n",
    "> **Input**: methylation profiles (.bed, .bedgragh, .cov.gz)  \n",
    "> **Output**: DMR list (.bedgraph)\n",
    "\n",
    "> <img src=\"images/notebook2.png\" width=\"700\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b33c8dd8-7550-4ec5-b532-8cef4d3489dc",
   "metadata": {},
   "source": [
    "## Differentially methylated regions (DMRs) \n",
    "\"After an initial analysis of global trends in a DNA methylation data set, the typical next step is the identification of differentially methylated regions (DMRs) that exhibit consistently different DNA methylation levels between sample groups (for example, cases versus controls). These DMRs can be as small as a single C or as large as an entire gene locus, depending on the biological question of interest and on the bioinformatic methods used for their identification. Although a single methylated CpG may occasionally be linked to gene expression regulation and may affect disease risk, the vast majority of DMRs reported in the literature fall within a size range of a few hundred to a few thousand bases. This range coincides with typical sizes of gene-regulatory regions, and it is widely believed that DMRs can control cell-type-specific transcriptional repression of an associated gene\" - [Bock, Christoph.(2012)](https://www.nature.com/articles/nrg3273)\n",
    "\n",
    "> <img src=\"images/2_dmr.png\" width=\"500\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f009a90-4253-46aa-82dc-cfdc528fe141",
   "metadata": {},
   "source": [
    "## Introduction to metilene\n",
    "\n",
    "metilene is a software tool to annotate differentially methylated regions (DMRs) and differentially methylated CpG sites (DMCs) from Methyl-seq data. metilene accounts for intra-group variances and offers different modes de-novo DMR detection, DMR detection within a known set of genomic features, and DMC detection. Various biological data can be used, metilene works with Whole-Genome Bisulfite Sequencing (WGBS), Reduced Representation Bisulfite Sequencing (RRBS), and any other input data, as long as absolute (methylation) levels and genomic coordinates are provided. metilene uses a circular binary segmentation and a 2D-KS test to call DMRs. Adjusted p-values are calculated using the Bonferroni correction.\n",
    "\n",
    "## Install metilene\n",
    "\n",
    "If you've followed the previous Bismark tutorial (tutorial_1-bismark.ipynb), metilene should already be installed in the instance. If not, please use the following command to install metilene: `! conda install -y -c conda-forge -c bioconda samtools=1.15.* metilene`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "429dde22-93d2-4b34-a696-f26e219c219c",
   "metadata": {},
   "source": [
    "## Generate input files for metilene\n",
    "\n",
    "Quick start, to do a typical de-novo annotation of DMRs, the command looks like: `metilene -a g1 -b g2 methylation-file`, while the input **methylation-file** containing all methylation data is a **SORTED** tab-separated file with the following format and header:  \n",
    "| chr | pos | g1_xxx | g1_xxx | \\[...\\] | g2_xxx | g2_xxx | \\[...\\] |\n",
    "| --- | --- | --- | --- | --- | --- | --- | --- |\n",
    "\n",
    "where the first column refers to the chromosome, the second column to the genomic position of the CpG and all following columns to the absolute methylation ratio. All ratio columns are dedicated to the group described by the prefix in their header, e.g., g1 or g2. Options `-a` and `-b` indicate the groups that are considered. \n",
    "\n",
    "To separate the results from the last tutorial, create a new directory called `Tutorial_2` first, and all the results generated from this tutorial will be saved in this directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "946bfe1e-7e51-4b45-98be-774f6ea48f4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a directory called Tutorial_2\n",
    "! mkdir -p Tutorial_2\n",
    "\n",
    "# Copy the methylation data from tutorial 1 to the directory Tutorial_2\n",
    "! cp Tutorial_1/bismark/*bedGraph.gz Tutorial_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a395009-1109-41e1-82d2-86f64e6aa660",
   "metadata": {},
   "source": [
    "#### Sort .bedGraph files generated by the Bismark pipeline\n",
    "\n",
    "Since metilene requires a sorted input file, we need to sort our methylation data using command `sortBed` from bedtools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ddef02c-cc44-43a1-b6aa-36c5a6248de9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort all the methylation files\n",
    "!for file_name in Tutorial_2/*bedGraph.gz; do \\\n",
    "    sortBed -i ${file_name} > ${file_name//_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz/.bedgraph}; \\\n",
    "done;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6baf73fc-23ca-4657-830d-9e50095c918a",
   "metadata": {},
   "source": [
    "#### Prepare input files for metilene\n",
    "metilene offers an easy way to generate an appropriate input file (`methylation-file`) containing all methylation rates, using a Perl script: `metilene_input.pl`. This step will generate a sorted tab-separated input file from multiple bed-files. \n",
    "\n",
    "Below is the parameter list of `metilene_input.pl`:\n",
    "\n",
    "> <img src = \"images/metilene_input_param.png\" width=\"700\" />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0add508-27f3-40ce-a223-e36266c8693c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"PATH\"] += os.pathsep + os.environ[\"HOME\"]+\"/mambaforge/bin\"\n",
    "\n",
    "! metilene_input.pl \\\n",
    "    --in1 Tutorial_2/SRX202087.bedgraph,Tutorial_2/SRX271141.bedgraph \\\n",
    "    --in2 Tutorial_2/SRX202088.bedgraph,Tutorial_2/SRX271142.bedgraph \\\n",
    "    --h1 serum \\\n",
    "    --h2 2i \\\n",
    "    -o Tutorial_2/metilene.input"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12b2bf0a-d598-41fd-8233-85607ecaa031",
   "metadata": {},
   "source": [
    "View the top 30 entries from metilene input file (`metilene.input`) just generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f9a75b6-5a4b-43f3-8b5b-64334e1ec40d",
   "metadata": {},
   "outputs": [],
   "source": [
    "! head -30 Tutorial_2/metilene.input "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "923adbd6-1bfa-4790-99e5-c07f0f0b4d8b",
   "metadata": {},
   "source": [
    "## Run metilene\n",
    "\n",
    "The goal for metilene is to annotate differentially methylated regions. There are three options (mode, `-f`) to define the regions: \n",
    "- **DMR de-novo annotation**. The default mode of metilene annotates DMRs de-novo without using any prior information on genomic features. \n",
    "- **DMR annotation in known features** - genomic features such as promoter regions, or CpG island need to be provided.\n",
    "- **DMC annotation** (not used in this tutorial) - metilene offers the possibility to test each CpG for differential methylation. Statistical tests (KS-test and Mann-Whitney-U test) are calculated for each CpG site, and corresponding p-values are reported in the output.\n",
    "\n",
    "Below is the complete parameter list when running `metilene`:\n",
    "> <img src=\"images/metilene_param.png\" width=\"700\" />\n",
    "\n",
    "The output for the de-novo DMR annotation mode consists of a bed-like format:\n",
    "| chr | start | stop | q-value | mean methylation difference | #CpGs | p (MWU) | p (2D KS) | mean g1 | mean g2 |\n",
    "| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |   \n",
    "\n",
    "While \"mean g1\" and \"mean g2\" refer to the absolute mean methylation level for the corresponding segment in both groups, the difference is given in the 5th column. Single CpGs are not tested using the 2D KS-test. Here, q-values are (per default) Bonferroni adjusted based on MWU-test p-values (since version 0.2-8, Benjamini-Hochberg (FDR) can be chosen instead of Bonferroni adjustment). All outputs are unsorted when using multiple threads. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0649bc1-6c2a-4821-b722-cca4f6c81a22",
   "metadata": {},
   "source": [
    "#### **de-novo annotation** (option 1)\n",
    "\n",
    "\"The default mode of metilene annotates DMRs de-novo without using any prior information on genomic features, e.g., promoter regions. Here a fast circular binary segmentation approach on the mean difference signal of both groups is used (Siegmund, 1986; Olshen et al., 2004). After additional filter steps are passed, potential DMRs are tested using a two-dimensional Kolmogorov-Smirnov-Test (KS-test)(Fasano and Franceschini, 1987). DMRs are finally tested through the Mann-Whitney-U test.\" (From [metilene user guide](https://www.bioinf.uni-leipzig.de/Software/metilene/Manual/#5_dmr_de-novo_annotation))\n",
    "\n",
    "> <img src=\"images/2_metilene_workflow.png\" width=\"600\">\n",
    "\n",
    "When running metilene, using `-a` and `-b` to specify the group to compare, and the output files are **sorted** based on the 1st and 2nd fields (chromosome and start position):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3006d229-a2c1-4ddf-96ef-3c4ae8073062",
   "metadata": {},
   "outputs": [],
   "source": [
    "!metilene -a serum -b 2i Tutorial_2/metilene.input | sort -V -k1,1 -k2,2n > Tutorial_2/metilene_denovo.output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "953e1a5b-f9cb-4f70-89b5-bda51cdb7824",
   "metadata": {},
   "source": [
    "#### **DMR annotation in known features**  (option 2)\n",
    "\n",
    "\"Instead of annotating de-novo DMRs, metilene can be used to find significant DMRs within a given group of genomic features. Here, the first step calling the circular binary segmentation algorithm is skipped. Instead, statistical tests are performed for each feature, and corresponding p-values are reported in the output. Use the \"-B bedfile\" option to define windows through a bedfile SORTED equally to the data input file.\" (From [metilene user guide](https://www.bioinf.uni-leipzig.de/Software/metilene/Manual/#6_dmr_annotation_in_known_features))\n",
    "\n",
    "An example of the external annotation file (CpG islands) can be found at [Table browser - UCSC Genome Browser](\n",
    "https://genome.ucsc.edu/cgi-bin/hgTables), where users can retrieve and export data from the Genome Browser annotation track database. This file is included in the data provided in this tutorial, named `cpgIslandExt_mm39_chr6.bed`:\n",
    "> <img src=\"images/external_annotation_file.png\" width=\"800\"/>\n",
    "\n",
    "To run metilene with the known features, using `-f 2` to define the mode (1: de-novo, 2: pre-defined regions, 3: DMCs). And using `-B` to provide a SORTED (equally to the input data) bed file containing regions for mode 2: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd1bd69b-5f47-4923-8573-dbe659cf2eba",
   "metadata": {},
   "outputs": [],
   "source": [
    "!metilene -a serum -b 2i -f 2 -B Tutorial_1/fastq/cpgIslandExt_mm39_chr6.bed Tutorial_2/metilene.input | sort -V -k1,1 -k2,2n > Tutorial_2/metilene_cgi.output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45aa4b7d-1cab-484e-89b6-71bfc3617bd3",
   "metadata": {},
   "source": [
    "#### Filter output file and plot basic DMR statistics\n",
    "An easy way to filter the already called DMRs is offered by `metilene_output.pl`. It will create some basic statistic plots characterizing your DMRs, i.e., distribution of DMR differences, DMR length in nucleotides and #CpGs, DMR differences vs. q-values, mean methylation group 1 vs. mean methylation group 2 and DMR length in nucleotides vs. length in CpGs. DMRs can by filtered by q-value, #CpGs, length in nucleotides and mean methylation difference. \n",
    "\n",
    "There are 3 files produced: \n",
    "1. bedGraph file containing the methylation difference for each DMR\n",
    "2. basic statistic pdf \n",
    "3. filtered bedGraph-like file, containing all information already in the metilene output. \n",
    "Please see [here](https://www.bioinf.uni-leipzig.de/Software/metilene/Manual/#filter_output_file_and_plot_basic_dmr_statistics) for detailed parameter list and their default values.\n",
    "\n",
    "The following code will filter the metilene output files from the previous steps (option 1 `metilene_denovo.output` and option 2 `metilene_cgi.output`), based on the default settings: adjust p-value < 0.05, number of CpGs >=10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e338c101-c784-4058-ba17-6ec73375b18c",
   "metadata": {},
   "outputs": [],
   "source": [
    "!metilene_output.pl -q Tutorial_2/metilene_denovo.output -o Tutorial_2/denovo -a serum -b 2i\n",
    "!metilene_output.pl -q Tutorial_2/metilene_cgi.output -o Tutorial_2/cgi -a serum -b 2i"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9accccd9-98de-4421-83fb-99d93c50136a",
   "metadata": {},
   "source": [
    "In our test example dataset, the de-novo annotation methods have 683 DMRs (`Tutorial_2/denovo_qval.0.05.bedgraph`), while the number of differentially methylated CpG islands as defined by the UCSC genome annotations is 2 (`Tutorial_2/cgi_qval.0.05.bedgraph`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bda459d-59ce-446f-884d-1bf9667dbf32",
   "metadata": {},
   "source": [
    "## Visualization using IGV\n",
    "\n",
    "In our example, there were only six CpG islands (see track `CGI (UCSC)`) are defined based and two of them can be considered as DMRs using the default threshold (adjust p-value <0.05). However, if you use the de-novo annotation, there are more DMRs (683 out of 812) that can be detected. You can view the results using IGV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac2d929b-7879-43f7-a988-e8cc8b254373",
   "metadata": {},
   "outputs": [],
   "source": [
    "import igv_notebook\n",
    "igv_notebook.init()\n",
    "\n",
    "# reference genome\n",
    "b = igv_notebook.Browser(\n",
    "    {\n",
    "        \"genome\": \"mm39\",\n",
    "        \"locus\": \"chr6:4,300,000-5,800,000\"\n",
    "    }\n",
    ")\n",
    "\n",
    "# load unfiltered CGI output \n",
    "b.load_track({    \"name\": \"CGI - all\",\n",
    "    \"path\": \"Tutorial_2/metilene_cgi.output\",\n",
    "    \"format\": \"bed\",\n",
    "    \"type\": \"annotation\",\n",
    "    \"color\": \"black\"\n",
    "})\n",
    "\n",
    "# load filtered DMRs using both options\n",
    "b.load_track({\n",
    "    \"name\": \"CGI - filtered\",\n",
    "    \"path\": \"Tutorial_2/cgi_qval.0.05.bedgraph\",\n",
    "    \"format\": \"bedgraph\",\n",
    "    \"type\": \"wig\",\n",
    "    \"color\": \"red\"\n",
    "})\n",
    "b.load_track({\n",
    "    \"name\": \"denovo_DMR - filtered\",\n",
    "    \"path\": \"Tutorial_2/denovo_qval.0.05.bedgraph\",\n",
    "    \"format\": \"bedgraph\",\n",
    "    \"type\": \"wig\",\n",
    "    \"color\": \"red\"\n",
    "})\n",
    "\n",
    "# load the methylation levels from 4 samples\n",
    "b.load_track({\n",
    "    \"name\": \"SRX271141_bed\",\n",
    "    \"path\": \"Tutorial_2/SRX271141_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\",\n",
    "    \"mode\": \"bisulfite\" ,\n",
    "    \"color\": \"orange\"\n",
    "})\n",
    "b.load_track({\n",
    "    \"name\": \"SRX202087_bed\",\n",
    "    \"path\": \"Tutorial_2/SRX202087_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\",\n",
    "    \"mode\": \"bisulfite\" ,\n",
    "    \"color\": \"orange\"\n",
    "})\n",
    "\n",
    "b.load_track({\n",
    "    \"name\": \"SRX271142_bed\",\n",
    "    \"path\": \"Tutorial_2/SRX271142_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\",\n",
    "    \"mode\": \"bisulfite\" ,\n",
    "    \"color\": \"blue\"\n",
    "})\n",
    "b.load_track({\n",
    "    \"name\": \"SRX202088_bed\",\n",
    "    \"path\": \"Tutorial_2/SRX202088_R1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz\",\n",
    "    \"mode\": \"bisulfite\" ,\n",
    "    \"color\": \"blue\"\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23d4a98a-8c93-476a-8b8f-2469ceb01a99",
   "metadata": {},
   "source": [
    "You can zoom in to see more details in a specified region, for example `b.search('chr6:4720,000-4800,000')`:  \n",
    "\n",
    "<img src = \"images/2_igv_results.png\" width=\"1000\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74eca58f-b1cf-45ad-8008-0e2cd6a6bda7",
   "metadata": {},
   "source": [
    "Our example is not the best illustration of DMRs, since 2i ESCs (SRX271142, SRX202088) are **globally** hypomethylated compared to conventional ESCs maintained in serum (RX271141, SRX202087). That is why when using the de-novo annotation, there were many DMRs can be detected (track `denovo_DMR`). But if examine closely, we can see how the methylation levels are different in the CpG islands or de-novo defined regions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bfee342-7b62-48e1-8975-682323d8532c",
   "metadata": {},
   "source": [
    "## Terms & Quiz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9947f039-5280-4a07-9097-f15a978e7eb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install jupytercards --quiet\n",
    "from jupytercards import display_flashcards\n",
    "display_flashcards('quiz_files/f2.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3abc5c18-e21d-4c98-b666-198d622c4d0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install jupyterquiz==2.0.7 --quiet\n",
    "from jupyterquiz import display_quiz\n",
    "display_quiz('quiz_files/q2.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9630e270-7682-4a76-9de2-f8e27c36390b",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    <b>Don't forget:</b> after finish running the notebook, <b>stop</b> the notebook in Vertex AI Workbench to avoid cost accumulation.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}