DNA Methylation data analysis

Author(s) Joachim Wolff avatar Joachim WolffDevon Ryan avatar Devon RyanVerena Moosmann avatar Verena Moosmann
Editor(s) orcid logoPavankumar Videm avatar Pavankumar Videmorcid logoWolfgang Maier avatar Wolfgang Maier
Reviewers Saskia Hiltemann avatarHelena Rasche avatarNicola Soranzo avatarBĂ©rĂ©nice Batut avatarGildas Le CorguillĂ© avatarWilliam Durand avatarBjörn GrĂŒning avatarNiall Beard avatar
Overview
Creative Commons License: CC-BY Questions:

  • What is methylation and why it cannot be recognised by a normal NGS procedure?

  • Can a different methylation influence the expression of a gene? How?

  • Which tools you can use to analyse methylation data?

Objectives:

  • Learn how to analyse methylation data

  • Get a first intuition what are common pitfalls.

Requirements:
Time estimation: 3 hours
Supporting Materials:
Published: Feb 16, 2017
Last modification: Apr 3, 2025
License: Tutorial Content is licensed under Creative Commons Attribution 4.0 International License. The GTN Framework is licensed under MIT
PURL: https://gxy.io/GTN:T00142
Rating: 3.7 (3 recent ratings, 13 all time)
Revision: 24

We will use a small subset of the original data. If we would do the computation on the orginal data the computation time for a tutorial is too long. To show you all necessary steps for Methyl-Seq we decided to use a subset of the data set. In a second step we use precomputed data from the study to show you different levels of methylation. We will consider samples from normal breast cells (NB), fibroadenoma (noncancerous breast tumor, BT089), two invasive ductal carcinomas (BT126, BT198) and a breast adenocarcinoma cell line (MCF7).

This tutorial is based off of Lin et al. 2015. The data we use in this tutorial is available at Zenodo.

Agenda

In this tutorial, we will deal with:

  1. Data upload
  2. Quality Control
  3. Alignment
  4. Methylation bias and metric extraction
  5. Visualization
  6. Metilene

Data upload

We will start by loading the example dataset which will be used for the tutorial into Galaxy

Hands On: Get the data into Galaxy

  1. Create a new history

    To create a new history simply click the icon at the top of the history panel:

    UI for creating new history

  2. Import the two example datasets from Zenodo or the shared data library:

    https://zenodo.org/record/557099/files/subset_1.fastq
https://zenodo.org/record/557099/files/subset_2.fastq
    • Copy the link location

    • Click Upload Data at the top of the tool panel

    • Select Paste/Fetch Data

    • Paste the link(s) into the text field

    • Press Start

    • Close the window

    As an alternative to uploading the data from a URL or your computer, the files may also have been made available from a shared data library:

    1. Go into Libraries (left panel)
    2. Navigate to the correct folder as indicated by your instructor.
      • On most Galaxies tutorial data will be provided in a folder named GTN - Material –> Topic Name -> Tutorial Name.
    3. Select the desired files
    4. Click on Add to History near the top and select as Datasets from the dropdown menu

    5. In the pop-up window, choose

      • “Select history”: the history you want to import the data to (or create a new one)
    6. Click on Import

Quality Control

The first step in any analysis should always be quality control. We will use the Falco tool to asses the quality of our reads and determine if we need to perform any data cleaning before proceeding with our analysis. Falco is an efficiency optimized rewrite of FastQC

Hands On: Quality Control
  1. Falco ( Galaxy version 1.2.4+galaxy0) with the following parameters:
    • “Raw read data from your current history”: subset_1.fastq.gz and subset_2.fastq.gz
    1. Click on Multiple datasets
    2. Select several files by keeping the Ctrl (or COMMAND) key pressed and clicking on the files of interest

  2. Go to the web page result page and have a closer look at ‘Per base sequence content’

    Falco webpage results.

    Question
    1. Note the GC distribution and percentage of “T” and “C”. Why is this so weird?
    2. Is everything as expected?
    1. The attentive audience of the theory part knows: Every C-meth stays a C and every normal C becomes a T during the bisulfite conversion.
    2. Yes it is. Always be careful and have the specific characteristics of your data in mind during the interpretation of Falco results.

Alignment

Hands On: Mapping with bwameth

We will now map the imported dataset against a reference genome.

  1. bwameth ( Galaxy version 0.2.7+galaxy0) with the following parameters:
    • “Select a genome reference from your history or a built-in index?”: Use a built-in index
      • “Select a reference genome”: Human (hg38full)
    • “Is this library mate-paired”: Paired-end
      • “First read in pair”: subset_1.fastq
      • “Second read in pair”: subset_2.fastq
    Comment: Long compute times

    Please notice that mapping can take some time. If you want to skip this, we provide for you a precomputed alignment. Import https://zenodo.org/records/557099/files/aligned_subset.bam to your history.

    Question

    Why we need other alignment tools for bisulfite sequencing data?

    You may have noticed that all the C’s are C-meth’s and a T can be a T or a C. A mapper for methylation data needs to find out what is what.

Methylation bias and metric extraction

Hands On: Methylation bias

In this step we will have a look at the distribution of the methylation and will look at a possible bias.

  1. MethylDackel ( Galaxy version 0.5.2+galaxy0) with the following parameters:
    • “Load reference genome from”: Local cache
      • “Using reference genome”: Human (hg38)
    • “Sorted BAM file”: output of bwameth
    • “What do you want to do?”: Determine the position-dependent methylation bias in the dataset, producing diagnostic SVG images (mbias)
    • In “Advanced options”
      • “Keep singletons”: Yes
      • “Keep discordant alignmetns”: Yes

    Methylation bias example.

    Question
    1. Consider the original top strand output. Is there a methylation bias?
    2. If we would trim, what would be the start and the end positions?
    1. The distribution of the methylation is more or less equal. Only at the start and the end we could trim a bit but a +- 5% variation is acceptable.
    2. To trim the reads we would include for the first strand only the positions 0 to 145, for the second 6 to 149.
Hands On: Methylation extraction with MethylDackel

We will extract the methylation on the resulting BAM file of the alignment step. We need this to create a methylation level plot in the next step.

  1. MethylDackel ( Galaxy version 0.5.2+galaxy0) with the following parameters:
    • “Load reference genome from”: Local cache
      • “Using reference genome”: Human (hg38)
    • “Sorted BAM file”: output of bwameth
    • “What do you want to do?”: Extract methylation metrics from an alignment file in BAM/CRAM format (extract)
    • “Merge per-Cytosine metrics”: Yes
    • “Output options”: CpG methylation fractions (--fraction)

Visualization

Hands On

In this step we want to visualize the methylation level around all TSS of our data. When located at gene promoters, DNA methylation is usually a repressive mark.

  1. Wig/BedGraph-to-bigWig with the following parameters:

    • “Convert”: fraction CpG (result of MethylDackel )

      It can happen that you can not select the correct input file. In this case you have to add meta information about the used genome to the file.

      • Click on the pencil of the correct history item.
      • Change Database/Build: to the genome you used.
      • In our case the correct genome is Human Dec. 2013 (GRCh38/hg38) (hg38).

  2. Import the BED file with CpG islands from Zenodo into the history

    https://zenodo.org/records/557099/files/CpGIslands.bed
  3. computeMatrix ( Galaxy version 3.5.4+galaxy0) with the following parameters:
    • “Regions to plot”: CpGIslands.bed
    • “Sample order matters”: No
    • “Score file”: Output of Wig/BedGraph-to-bigWig
    • “computeMatrix has two main output options”: reference-point
  4. plotProfile ( Galaxy version 3.5.4+galaxy0) with the following parameters:
    • “Matrix file from the computeMatrix tool”: Matrix (output of computeMatrix )

The output should look like this:

Methylation output.

Lets see how the methylation looks for a few provided files:

  1. Import the BED file with CpG islands from Zenodo into the history

    https://zenodo.org/records/557099/files/NB1_CpG.meth.bedGraph
  2. Wig/BedGraph-to-bigWig with the following parameters:
    • “Convert”: NB1_CpG.meth.bedGraph
    Question

    The execution fails. Do you have an idea why?

    A conversion to bigWig would fail right now. If it turned green, the file size should be 0 bytes. Probably dataset info box shows some error message like hashMustFindVal: '1' not found. The reason is the source of the reference genome which was used. There is ensembl and UCSC as sources which differ in naming the chromosomes. Ensembl is using just numbers e.g. 1 for chromosome one. UCSC is using chr1 for the same. Be careful with this especially if you have data from different sources. We need to convert this.

    Comment: UCSC - Ensembl convert

    Download the file containing mapping between Ensembl and UCS chromosome convention of hg38

    https://raw.githubusercontent.com/dpryan79/ChromosomeMappings/master/GRCh38_ensembl2UCSC.txt

    Replace column ( Galaxy version 0.2) with the follwing parameters:

    • “File in which you want to replace some values”: NB1_CpG.meth.bedGraph
    • “Replace information file”: GRCh38_ensembl2UCSC.txt
    • “Which column should be replaced?”: Column: 1
    • “Skip this many starting lines”: 1
    • “Delimited by”: Tab

  3. To save compute time we prepared the converted files for you. Import the following files. Create a collection list and label it all_coverage_files. Strip the file extension from the name. For example, rename from NB1_CpG.meth_ucsc.bedGraphto NB1_CpG.

    https://zenodo.org/records/557099/files/NB1_CpG.meth_ucsc.bedGraph
https://zenodo.org/records/557099/files/NB2_CpG.meth_ucsc.bedGraph
https://zenodo.org/records/557099/files/BT089_CpG.meth_ucsc.bedGraph
https://zenodo.org/records/557099/files/BT126_CpG.meth_ucsc.bedGraph
https://zenodo.org/records/557099/files/BT198_CpG.meth_ucsc.bedGraph
https://zenodo.org/records/557099/files/MCF7_CpG.meth_ucsc.bedgraph
    • Click on Select Items at the top of the history panel Select Items button
    • Check all the datasets in your history you would like to include

    • Click n of N selected and choose Build Dataset List

      build list collection menu item

    • Enter a name for your collection
    • Click Create collection to build your collection
    • Click on the checkmark icon at the top of your history again

  4. Change the datatype to bedgraph and set the database to hg38

    • Click on the pencil icon for the dataset to edit its attributes
    • In the central panel, click Datatypes tab on the top
    • In the Assign Datatype, select bedgraph from “New type” dropdown
      • Tip: you can start typing the datatype into the field to filter the dropdown menu
    • Click the Save button

    • Click the desired dataset’s name to expand it.

    • Click on the “?” next to database indicator:

      UI for changing dbkey

    • In the central panel, change the Database/Build field
    • Select your desired database key from the dropdown list: hg38
    • Click the Save button

  5. Wig/BedGraph-to-bigWig with the following parameters:
    • “Convert”: all_coverage_files
  6. computeMatrix ( Galaxy version 3.5.4+galaxy0) with the following parameters:
    • “Regions to plot”: CpGIslands.bed
    • “Sample order matters”: No
    • “Score file”: Output of previous Wig/BedGraph-to-bigWig
    • “computeMatrix has two main output options”: reference-point
  7. plotProfile ( Galaxy version 3.5.4+galaxy0) with the following parameters:
    • “Matrix file from the computeMatrix tool”: Matrix (output of previous computeMatrix )
    • “Show advanced options”: Yes
      • “Make one plot per group of regions”: Yes

The output should look like this:

Methylation level around TSS.

Metilene

Hands On: Metilene

With metilene it is possible to detect differentially methylated regions (DMRs) which is a necessary prerequisite for characterizing different epigenetic states.

  1. Import the following files from Zenodo into yout history

    https://zenodo.org/records/557099/files/NB1_CpG.meth.bedGraph
https://zenodo.org/records/557099/files/NB2_CpG.meth.bedGraph
https://zenodo.org/records/557099/files/BT198_CpG.meth.bedGraph

  2. Metilene ( Galaxy version 0.2.6.1) with the following parameters:

    • “Input group 1”: NB1_CpG.meth.bedGraph and NB2_CpG.meth.bedGraph
    • “Input group 2”: BT198_CpG.meth.bedGraph
    • “BED file containing regions of interest”: CpGIslands.bed
    Question

    Have a look at the produced pdf document. What is the data showing?

    It shows the distribution of DMR differences, DMR length in nucleotides and number CpGs, DMR differences vs. q-values, mean methylation group 1 vs. mean methylation group 2 and DMR length in nucleotides vs. length in CpGs