Goby supports the analysis of DNA methylation  NGS data such as produced by the Methyl-Seq or RRBS protocols. This tutorial describes how to analyze such datasets from the command line.


Reads. We assume that you have two fastq file with bisulfite converted reads. We call these files A.fastq and B.fastq. If you reads are stored in Goby compact format, simply use the compact filename instead of fastq.

GSNAP. We will use the GSNAP aligner. You can obtain GSNAP from here. You must compile GSNAP with Goby support (see how to in the README after you have installed the Goby C/C++ with these detailed instructions).

Reference genome. Below, we assume that you will align A and B against a mouse genome. Please see the GSNAP documentation to index a reference genome and create the optional bisulfite index. This is typically done with cmetindex -d <genome>

Bisulfite alignment

You can now align the samples to the reference genome:

gsnap A.fastq -d <genome> --use-cmet --format=goby --goby-output=sample-A --terminal-penalty=100
gsnap B.fastq -d <genome> --use-cmet --format=goby --goby-output=sample-B --terminal-penalty=100

Please be patient as the alignment complete. You can speed up alignment by specifying -i 1 to disabling gapped alignment. This should produce two Goby alignments: sample-A and sample-B. Each alignment consists of several files, see this for an explanation of the content of each file.

Calling differentially methylated sites

After you have completed the bisulfite alignment, you can call differentially methylated sites directly with Goby. The following command will compare samples A and B and write a VCF output file:

goby 2g discover-sequence-variants sample-A.entries sample-B.entries --compare A/B --groups A=sample-A/B=sample-B \
  --format methylation --output meth.vcf

See the help output (–help/-h) to see how to use this command with several samples per group. You should now have a VCF file that describes methylation rate (MR) at each genomic site where Cs occur in the forward or reverse strand and that were covered by the reads.

Viewing sites in their genomic context with methylation information

You can visualize the Goby VCF file you have just produced in IGV. To this end, you need to compress and index the file. You can use tabix for indexing, as follows:

bgzip meth.vcf
tabix -p vcf meth.vcf.gz

Check that indexing worked by requesting an arbitrary position:

tabix meth.vcf.gz 11

The previous line should return entries at the start of chromosome 11, if no results are returned a problem occurred running bgzip or tabix. It is important to check the index because bgzip and tabix do not report all error conditions and sometimes fail silently.

You should now have the following files:

  • meth.vcf.gz
  • meth.vcf.gz.tbi

Start IGV 2.0 and select load. Pick the meth.vcf.gz file to create a VCF track. IGV inspects the file, determines it can be used to render methylation rates, and switches the track display appropriately. The display will show two samples and render methylation rates in color. Red indicates 100% methylation, white-grey ~0-25% methylation yellow ~25-50% and blue ~50-75% intermediate methylation levels. The colors are chosen to make it easier to spot differences between samples (this is especially useful when multiple samples are shown and organized in blocks with samples from the same group together). Rolling the mouse over the colored bars will show the VCF data for this site and sample. When rolling over the pink bars on the top of the VCF track, the tool tip will show the statistics that pertain to the site across samples (look for the fisher and odds-ratio statistics here for instance).

See also this tutorial about how to estimate methylation rates over regions.