Bismark Bisulfite Mapper

User Guide - v0.18.0

15 May, 2017

This User Guide outlines the Bismark suite of tools and gives more details for each individual step. For troubleshooting some of the more commonly experienced problems in sequencing in general and bisulfite-sequencing in particular please browse through the sequencing section at QCFail.com.

1) Quick Reference

Bismark needs a working version of Perl and it is run from the command line. Furthermore, Bowtie or Bowtie 2 needs to be installed on your computer. For more information on how to run Bismark with Bowtie 2 please go to the end of this manual.

As of version 0.14.0 or higher, Bismark may be run using parallelisation for both the alignment and the methylation extraction step. Search for --multicore for more details below.

First you need to download a reference genome and place it in a genome folder. Genomes can be obtained e.g. from the Ensembl or NCBI websites. For the example below you would need to download the Homo sapiens genome. Bismark supports reference genome sequence files in FastA format, allowed file extensions are either either .fa or .fasta. Both single-entry and multiple-entry FastA files are supported.

The following examples will use the file test_dataset.fastq which is available for download from the Bismark project or Github pages (it contains 10,000 reads in FastQ format, Phred33 qualities, 50 bp long reads, from a human directional BS-Seq library). An example report for use with Bowtie 1 and Bowtie can be found in Appendix IV.

(I) Running bismark_genome_preparation

USAGE: bismark_genome_preparation [options] <path_to_genome_folder>

A typical genome indexing could look like this: /bismark/bismark_genome_preparation --path_to_bowtie /usr/bin/bowtie2/ --verbose /data/genomes/homo_sapiens/GRCh37/

(II) Running bismark

USAGE: bismark [options] --genome <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}

Typical alignment example (tolerating one non-bisulfite mismatch per read): bismark --bowtie1 -n 1 -l 50 /data/genomes/homo_sapiens/GRCh37/ test_dataset.fastq

This will produce two output files: 1. test_dataset_bismark.bam (contains all alignments plus methylation call strings) 2. test_dataset_bismark_SE_report.txt (contains alignment and methylation summary)

NOTE: In order to work properly the current working directory must contain the sequence files to be analysed.

(III) Running bismark_methylation_extractor

USAGE: bismark_methylation_extractor [options] <filenames>

A typical command to extract context-dependent (CpG/CHG/CHH) methylation could look like this: bismark_methylation_extractor --gzip test_dataset.fastq_bismark.bam

This will produce three output files:

(IV) Running bismark2report

USAGE: bismark2report [options]

This command attempts to find Bismark alignment, deduplication and methylation extraction (splitting) reports as well as M-bias files to generate a graphical HTML report such as this example Bismark paired-end report for each sample in a directory.

(V) Running bismark2summary

USAGE: bismark2summary [options]

This command scans the current working directory for different Bismark alignment, deduplication and methylation extraction (splitting) reports to produce a graphical summary HTML report, as well as a data table, for all files in a directory. Here is a sample Bismark Summary Report.

2) Bismark - General Information

What is Bismark?

Bismark is a set of tools for the time-efficient analysis of Bisulfite-Seq (BS-Seq) data. Bismark performs alignments of bisulfite-treated reads to a reference genome and cytosine methylation calls at the same time. Bismark is written in Perl and is run from the command line. Bisulfite-treated reads are mapped using the short read aligner Bowtie 1 (Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol 10:R25), or alternatively Bowtie 2, and therefore it is a requirement that Bowtie 1 (or Bowtie 2) are also installed on your machine (see Dependencies).

All files associated with Bismark as well as a test BS-Seq data set can be downloaded from Babraham Bioinformatics or Github.

We would like to hear your comments, suggestions or bugs about Bismark! Please e-mail them to: felix.krueger@babraham.ac.uk

Installation notes

Bismark is written in Perl and is executed from the command line. To install Bismark simply copy the bismark_v0.X.Y.tar.gz file into a Bismark installation folder and extract all files by typing:

tar xzf bismark_v0.X.Y.tar.gz

Dependencies

Bismark requires a working of Perl and Bowtie 1 / Bowtie 2 to be installed on your machine. Bismark will assume that the Bowtie 1 / Bowtie 2 executable is in your path unless the path to Bowtie is specified manually with: --path_to_bowtie </../../bowtie(1/2)>

Hardware requirements

Bismark holds the reference genome in memory and in addition to that runs up to four parallel instances of Bowtie. The memory usage is dependent on the size of the reference genome. For a large eukaryotic genome (human or mouse) we experienced a typical memory usage of around 12GB. We thus recommend running Bismark on a machine with 5 CPU cores and at least 12 GB of RAM. The memory requirements of Bowtie 2 are somewhat larger (possibly to allow gapped alignments). When running Bismark using Bowtie 2 we therefore recommend a system with at least 5 cores and > 16GB of RAM.

Alignment speed depends largely on the read length and alignment parameters used. Allowing many mismatches and using a short seed length (which is the default option for Bowtie 1, see below) tends to be fairly slow, whereas looking for near perfect matches can align around 5-25 million sequences per hour (using Bowtie 1). Since we haven’t tested Bowtie 2 very much yet we can’t give recommendations about alignment parameters and speed of Bowtie 2 at the current time.

BS-Seq test data set

A test BS-Seq data set is available for download from the Bismark project or Github pages. It contains 10,000 single- end shotgun BS reads from human ES cells in FastQ format (from SRR020138, Lister et al., 2009; trimmed to 50 bp; base call qualities are Sanger encoded Phred values (Phred33)).

Which kind of BS-Seq files are supported?

Bismark supports the alignment of bisulfite-treated reads (whole genome shotgun BS-Seq (WGSBS), reduced-representation BS-Seq (RRBS) or PBAT-Seq (Post-Bisulfite Adapter Tagging) for the following conditions:

A full list of alignments modes can be found in Bismark_alignment_modes.pdf.

In addition, Bismark retains much of the flexibility of Bowtie 1 / Bowtie 2 (adjustable seed length, number of mismatches, insert size ...). For a full list of options please run:

bismark --help

or see the Appendix at the end of this User Guide.

NOTE: It should be mentioned that Bismark supports only reads in base-space, such as from the Illumina platform. There are currently no plans to extend its functionality to colour-space reads.

How does Bismark work?

Sequence reads are first transformed into fully bisulfite-converted forward (C->T) and reverse read (G->A conversion of the forward strand) versions, before they are aligned to similarly converted versions of the genome (also C->T and G->A converted). Sequence reads that produce a unique best alignment from the four alignment processes against the bisulfite genomes (which are running in parallel) are then compared to the normal genomic sequence and the methylation state of all cytosine positions in the read is inferred. For use with Bowtie 1, a read is considered to align uniquely if one alignment exists that has with fewer mismatches to the genome than any other alignment (or if there is no other alignment). For Bowtie 2, a read is considered to align uniquely if an alignment has a unique best alignment score (as reported by the Bowtie 2 AS:i field). If a read produces several alignments with the same number of mismatches or with the same alignment score (AS:i field), a read (or a read-pair) is discarded altogether.

Bismark alignment and methylation call report

Upon completion, Bismark produces a run report containing information about the following: - Summary of alignment parameters used - Number of sequences analysed - Number of sequences with a unique best alignment (mapping efficiency) - Statistics summarising the bisulfite strand the unique best alignments came from - Number of cytosines analysed - Number of methylated and unmethylated cytosines - Percentage methylation of cytosines in CpG, CHG or CHH context (where H can be either A, T or C). This percentage is calculated individually for each context following the equation:

% methylation (context) = 100 * methylated Cs (context) / (methylated Cs (context) + unmethylated Cs (context)).

It should be stressed that the percent methylation value (context) is just a very rough calculation performed directly at the mapping step. Actual methylation levels after post-processing or filtering have been applied may vary.

3) Running Bismark

Running Bismark is split up into three main steps:

  1. First, the genome of interest needs to be bisulfite converted and indexed to allow Bowtie alignments. This step needs to be carried out only once for each genome. Note that Bowtie 1 and Bowtie 2 require distinct indexing steps since their indexes are not compatible.
  2. Bismark read alignment step. Simply specify a file to be analysed, a reference genome and alignment parameters. Bismark will produce a combined alignment/methylation call output (default is BAM format) as well as a run statistics report.
  3. Bismark methylation extractor. This step is optional and will extract the methylation information from the Bismark alignment output. Running this additional step allows splitting the methylation information up into the different contexts, getting strand-specific methylation information and offers some filtering options. You can also choose to sort the methylation information into bedGraph/coverage files, or even process them further to genome-wide cytosine methylation reports.

Each of these steps will be described in more detail (with examples) in the following sections.

(I) Bismark Genome Preparation

This script needs to be run only once to prepare the genome of interest for bisulfite alignments. You need to specify a directory containing the genome you want to align your reads against (please be aware that the bismark_genome_preparation script currently expects FastA files in this folder (with either .fa or .fasta extension, single or multiple sequence entries per file). Bismark will create two individual folders within this directory, one for a C->T converted genome and the other one for the G->A converted genome. After creating C->T and G->A versions of the genome they will be indexed in parallel using the indexer bowtie-build (or bowtie2-build). Once both C->T and G->A genome indices have been created you do not need to use the genome preparation script again (unless you want to align against a different genome...).

Please note that Bowtie 1 and 2 indexes are not compatible. To create a genome index for use with Bowtie 1 the option --bowtie1 needs to be included as well.

Running bismark_genome_preparation

USAGE: bismark_genome_preparation [options] <path_to_genome_folder>

A typical command could look like this:

bismark_genome_preparation --path_to_bowtie /usr/bin/bowtie2/ --verbose /data/genomes/homo_sapiens/GRCh37/

(II) Bismark Alignment Step

This step represents the actual bisulfite alignment and methylation calling part. Bismark requires the user to specify only two things:

  1. The directory containing the genome of interest. This folder must contain the unmodified genome (as .fa or .fasta files) as well as the two bisulfite genome subdirectories which were generated in the Bismark Genome Preparations step (see above).
  2. The sequence file(s) to be analysed (in either FastQ or FastA format).

All other parameters are optional.

For each sequence file or each set of paired-end sequence files, Bismark produces one alignment and methylation call output file as well as a report file detailing alignment and methylation call statistics for your information and record keeping.

Running bismark

Before running Bismark we recommend spending some time on quality control of the raw sequence files using FastQC. FastQC might be able to spot irregularities associated with your BS-Seq file, such as high base calling error rates or contaminating sequences such as PCR primers or Illumina adapters. Many sources of error impact detrimentally the alignment efficiencies and/or alignment positions, and thus possibly also affect the methylation calls and conclusions drawn from the experiment.

If no additional options are specified Bismark will use a set of default values, some of which are:

Using Bowtie 1:

Using Bowtie 2:

Even though the user is not required to specify additional alignment options it is often advisable to do so (e.g. when the default parameters are too strict). To see a full list of options please type bismark --help on the command line or see the Appendix at the end of this User Guide.

Directional BS-Seq libraries (default)

Bisulfite treatment of DNA and subsequent PCR amplification can give rise to four (bisulfite converted) strands for a given locus. Depending on the adapters used, BS-Seq libraries can be constructed in two different ways:

  1. If a library is directional, only reads which are (bisulfite converted) versions of the original top strand (OT) or the original bottom strand (OB) will be sequenced. Even though the strands complementary to OT (CTOT) and OB (CTOB) are generated in the BS-PCR step they will not be sequenced as they carry the wrong kind of adapter at their 5’-end. By default, Bismark performs only 2 read alignments to the OT and OB strands, thereby ignoring alignments coming from the complementary strands as they should theoretically not be present in the BS-Seq library in question.
  2. Alternatively, BS-Seq libraries can be constructed so that all four different strands generated in the BS-PCR can and will end up in the sequencing library with roughly the same likelihood. In this case all four strands (OT, CTOT, OB, CTOB) can produce valid alignments and the library is called non- directional. Specifying --non_directional instructs Bismark to use all four alignment outputs.

To summarise again: alignments to the original top strand or to the strand complementary to the original top strand (OT and CTOT) will both yield methylation information for cytosines on the top strand. Alignments to the original bottom strand or to the strand complementary to the original bottom strand (OB and CTOB) will both yield methylation information for cytosines on the bottom strand, i.e. they will appear to yield methylation information for G positions on the top strand of the reference genome.

For more information about how to extract methylation information of the four different alignment strands please see below in the section on the Bismark methylation extractor.

USAGE: bismark [options] --genome <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}

A typical single-end analysis could look like this:

bismark --genome /data/genomes/homo_sapiens/GRCh37/ sample.fastq.gz

What does the Bismark output look like?

Since version 0.6.x the default output of Bismark is in BAM/SAM format when using either Bowtie 1 or Bowtie 2. The former custom Bismark output for Bowtie 1, which used to be the standard output up to versions 0.5.x, is still available by specifying the option --vanilla (see below) but is likely going to be dropped entirely in a future release. The Bismark output using Bowtie 2 is invariably in BAM/SAM format (which is required to encode gapped alignments).

Bismark BAM/SAM output (default)

By default, Bismark generates SAM output for all alignment modes. Please note that reported quality values are encoded in Sanger format (Phred 33 scale), even if the input was in Phred64 or the old Solexa format.

  1. QNAME (seq-ID)
  2. FLAG (this flag tries to take the strand a bisulfite read originated from into account (this is different from ordinary DNA alignment flags!))
  3. RNAME (chromosome)
  4. POS (start position)
  5. MAPQ (only calculated for Bowtie 2, always 255 for Bowtie)
  6. CIGAR
  7. RNEXT
  8. PNEXT
  9. TLEN
  10. SEQ
  11. QUAL (Phred33 scale)
  12. NM-tag (edit distance to the reference)
  13. MD-tag (base-by-base mismatches to the reference) (14) XM-tag (methylation call string)
  14. XR-tag (read conversion state for the alignment) (16) XG-tag (genome conversion state for the alignment)

The mate read of paired-end alignments is written out as an additional separate line in the same format.

Custom (‘vanilla’) Bismark output (Bowtie 1 only)

Bismark can generate a comprehensive alignment and methylation call output file for each input file or set of paired-end input files. The sequence basecall qualities of the input FastQ files are written out into the Bismark output file as well to allow filtering on quality thresholds. Please note that the quality values are encoded in Sanger format (Phred 33 scale), even if the input was in Phred64 or the old Solexa format.

The single-end output contains the following information (1 line per sequence, tab delimited):

  1. seq-ID
  2. alignment strand
  3. chromosome
  4. start
  5. end
  6. original bisulfite read sequence
  7. equivalent genomic sequence (+2 extra bp) (8) methylation call string
  8. read conversion
  9. genome conversion
  10. read quality score (Phred33 scale)

Single-end alignment example:

  1. HWUSI-EAS611_0001:3:1:1186:18334#0/1
  2. +
  3. 4
  4. 122893213
  5. 122893242
  6. TGGGTTTTTAAGATTTTGTGTAGTTGGGGTTTGGAGATGG
  7. CGGGCCCTCAAGACCCTGCGCAGCTGGGGCCTGGAGACGGAG
  8. z...hhh.h....hhx..z.x..x.....hx......z..
  9. CT
  10. CT
  11. IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

The paired-end output looks like this (1 line per sequence pair, tab separated):

  1. seq-ID
  2. alignment strand
  3. chromosome
  4. start
  5. end
  6. original bisulfite read sequence 1
  7. equivalent genomic sequence 1 (+2 extra bp) (8) methylation call string 1
  8. original bisulfite read sequence 2
  9. equivalent genomic sequence 2 (+2 extra bp) (11) methylation call string 2
  10. read 1 conversion
  11. genome conversion
  12. read 1 quality score (Phred33 scale)
  13. read 2 quality score (Phred33 scale)

Paired-end alignment example

  1. HWUSI-EAS611_100205:2:1:13:1732#0
  2. +
  3. 14
  4. 62880539
  5. 62880652
  6. CGGGATTTCGCGGAGTACGGGTGATCGTGTGGAATATAGA
  7. CGGGACTCCGCGGAGCACGGGTGACCGTGTGGAATACAGAGT
  8. Z....h.xZ.Z....h.Z......xZ..........x...
  9. CAACTATCTAAAACTAAAAATAACGCCGCCCAAAAACTCT
  10. TCCGGCTGTCTGGAGCTGAAGATGGCGCCGCCCAGAAGCTCT
  11. .zx..x...xh.h..x..h..hh.Z..Z....x..h....
  12. CT
  13. CT
  14. IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
  15. IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

Data visualisation

To see the location of the mapped reads the Bismark output file can be imported into a genome viewer, such as SeqMonk, using the chromosome, start and end positions (this can be useful to identify regions in the genome which display an artefactually high number of aligned reads). The alignment output can also be used to apply post-processing steps such as de-duplication (allowing only 1 read for each position in the genome to remove PCR artefacts) or filtering on the number of bisulfite conversion related non-bisulfite mismatches * (please note that such post-processing scripts are not part of the Bismark package).

* Bisulfite conversion related non-bisulfite mismatches are mismatch positions which have a C in the BS-read but a T in the genome; such mismatches may occur due to the way bisulfite read alignments are performed. Reads containing this kind of mismatches are not automatically removed from the alignment output in order not to introduce a bias for methylated reads. It should be noted that, even though no methylation calls are performed for these positions, reads containing bisulfite conversion related non-bisulfite mismatches might lead to false alignments if particularly lax alignment parameters were specified.

Methylation call

The methylation call string contains a dot . for every position in the BS-read not involving a cytosine, or contains one of the following letters for the three different cytosine methylation contexts (UPPER CASE = METHYLATED, lower case = unmethylated):

(III) Bismark methylation extractor

Bismark comes with a supplementary bismark_methylation_extractor script which operates on Bismark result files and extracts the methylation call for every single C analysed. The position of every single C will be written out to a new output file, depending on its context (CpG, CHG or CHH), whereby methylated Cs will be labelled as forward reads (+), non-methylated Cs as reverse reads (-). The resulting files can be imported into a genome viewer such as SeqMonk (using the generic text import filter) and the analysis of methylation data can commence. Alternatively, the output of the methylation extractor can be transformed into a bedGraph and coverage file using the option --bedGraph (see also --counts). This step can also be accomplished from the methylation extractor output using the stand-alone script bismark2bedGraph (also part of the Bismark package available for download at bioinformatics.babraham.ac.uk). The coverage file can also be imported into SeqMonk directly using Import Data > Bismark (cov). Optionally, the bedGraph counts output can be used to generate a genome-wide cytosine report which reports the number on every single CpG (optionally every single cytosine) in the genome, irrespective of whether it was covered by any reads or not. As this type of report is informative for cytosines on both strands the output may be fairly large (~46mn CpG positions or >1.2bn total cytosine positions in the human genome...). The bedGraph to genome-wide cytosine report conversion can also be run individually using the stand- alone module coverage2cytosine (also part of the Bismark package available for download at bioinformatics.babraham.ac.uk).

As of Bismark version 0.6 or higher the default input format for the bismark_methylation_extractor is BAM/SAM (or potentially CRAM if you’ve got Samtools 1.2+ installed). The former custom Bismark format can still be used by specifying --vanilla (for the time being).

The methylation extractor output looks like this (tab separated):

1. seq-ID
2. methylation state
3. chromosome
4. start position (= end position)
5. methylation call

Methylated cytosines receive a + orientation, unmethylated cytosines receive a - orientation.

Examples for cytosines in CpG context:

HWUSI-EAS611_0006:3:1:1058:15806#0/1 - 6 91793279 z
HWUSI-EAS611_0006:3:1:1058:17564#0/1 + 8 122855484 Z

Examples for cytosines in CHG context:

HWUSI-EAS611_0006:3:1:1054:1405#0/1 - 7 89920171 x
HWUSI-EAS611_0006:3:1:1054:1405#0/1 + 7 89920172 X

Examples for cytosines in CHH context:

HWUSI-EAS611_0006:3:1:1054:1405#0/1 - 7 89920184 h

The bismark_methylation_extractor comes with a few options, such as ignoring the first number of positions in the methylation call string, e.g. to remove a restriction enzyme site (if RRBS is performed with non-directional BS-Seq libraries it might be required to remove reconstituted MspI sites at the beginning of each read as they will introduce a bias into the first methylation call). Another useful option for paired-end reads is called --no_overlap (on by default): specifying this option will extract the methylation calls of overlapping parts in the middle of paired-end reads only once (using the calls from the first read which is presumably the one with a lowest error rate).

For a full list of options type: bismark_methylation_extractor --help on the command line or refer to the Appendix section at the end of this User Guide.

Methylation extractor output

By default, the bismark_methylation_extractor discriminates between cytosines in CpG, CHG or CHH context. If desired, CHG and CHH contexts can be merged into a single non-CpG context by specifying the option --merge_non_CpG (as a word of warning, this might produce files with up to several hundred million lines...).

Strand-specific methylation output files (default):

As its default option, the bismark_methylation_extractor will produce a strand-specific output which will use the following abbreviations in the output file name to indicate the strand the alignment came from:

OT    –  original top strand
CTOT  –  complementary to original top strand
OB    –  original bottom strand
CTOB  –  complementary to original bottom strand

Methylation calls from OT and CTOT will be informative for cytosine methylation positions on the original top strand, calls from OB and CTOB will be informative for cytosine methylation positions on the original bottom strand. Please note that specifying the --directional (the default mode) option in the Bismark alignment step will not report any alignments to the CTOT or CTOB strands.

As cytosines can exist in any of three different sequence contexts (CpG, CHG or CHH) the bismark_methylation_extractor default output will consist of 12 individual output files per input file (CpG_OT_..., CpG_CTOT_..., CpG_OB_... etc.).

Context-dependent methylation output files (--comprehensive option):

If strand-specific methylation is not of interest, all available methylation information can be pooled into a single context-dependent file (information from any of the four strands will be pooled). This will default to three output files (CpG-context, CHG-context and CHH-context), or result in 2 output files (CpG-context and Non-CpG-context) if --merge_non_CpG was selected (note that this can result in enormous file sizes for the non-CpG output).

Both strand-specific and context-dependent options can be combined with the --merge_non_CpG option.

Optional bedGraph output

The Bismark methylation extractor can optionally also output a file in bedGraph format which uses 0-based genomic start and 1- based end coordinates. The module bismark2bedGraph (part of the Bismark package) may also be run individually. It will be sorted by chromosomal coordinates and looks like this: <chromosome> <start position> <end position> <methylation percentage>

As the methylation percentage is per se not informative of the actual read coverage of detected methylated or unmethylated reads at a position, bismark2bedGraph also writes out a coverage file (using 1-based genomic genomic coordinates) that features two additional columns: <chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated>

These two additional columns enable basically any downstream processing from the file. By default, this mode will only consider cytosines in CpG context, but it can be extended to cytosines in any sequence context by using the option --CX (cf. Appendix (III)).

Optional genome-wide cytosine report output

Starting from the coverage output, the Bismark methylation extractor can optionally also output a genome-wide cytosine methylation report. The module coverage2cytosine (part of the Bismark package) may also be run individually. It is also sorted by chromosomal coordinates but also contains the sequence context and is in the following format: <chromosome> <position> <strand> <count methylated> <count unmethylated> <C-context> <trinucleotide context>

The main difference to the bedGraph or coverage output is that every cytosine on both the top and bottom strands will be considered irrespective of whether they were actually covered by any reads in the experiment or not. For this to work one has to also specify the genome that was used for the Bismark alignments using the option --genome_folder <path>. As for the bedGraph mode, this will only consider cytosines in CpG context by default but can be extended to cytosines in any sequence context by using the option --CX (cf. Appendix (III)). Be aware though that this might mean an output with individual lines for more than 1.1 billion cytosines for any large mammalian genome...

M-bias plot

Starting with Bismark v0.8.0, the Bismark methylation extractor also produces a methylation bias plot which shows the methylation proportion across each possible position in the read (described in further detail in: Hansen et al., Genome Biology, 2012, 13:R83). The data for the M-bias plot is also written into a coverage text file (ending in .cov or .cov.gz) and is in the following format:

<read position> <count methylated> <count unmethylated> <% methylation> <total coverage>

This allows generating nice graphs by alternative means, e.g. using R or Excel. The plot is also drawn into a .png file which requires the Perl module GD::Graph (more specifically, both modules GD::Graph::lines and GD::Graph::colour are required); if GD::Graph cannot be found on the system, only the table will be printed. The plot also contains the absolute number of methylation calls (both methylated and unmethylated) per position. For paired-end reads two individual M-bias plots will be drawn.

The M-bias plot can for example show the methylation bias at the start of reads in PBAT-Seq experiments:

SE PBAT

For more on this topic please also see this post on QCFail.com.

Or it can reveal a 3’-end-repair bias at the first couple of positions in read 2 of paired-end reads, like here:

PE fill-in

For more on this topic please also see this post on QCFail.com.

The M-bias plot should enable researchers to make an informed decision whether or not to leave the bias in the final data or to remove it (e.g. using the methylation extractor option --ignore).

(III) Running bismark_methylation_extractor

USAGE: bismark_methylation_extractor [options] <filenames>

A typical command for a single-end file could look like this:

bismark_methylation_extractor -s --gzip sample_bismark_bt2.bam

A typical command for a paired-end file could look like this:

bismark_methylation_extractor -p --gzip sample_bismark_bt2_pe.bam

The methylation extractor can also auto-detect the alignment mode and will set the options above automatically. A typical command including the optional bedGraph output could look like this:

bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G sample_bismark_bt2.bam

A typical command including the optional genome-wide cytosine methylation report could look like this:

bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G --cytosine_report --genome_folder /path_to_genome_folder/ sample_bismark_bt2.bam

(IV) The Bismark HTML Processing Report

The script bismark2report uses a Bismark alignment report, and optionally further reports of the Bismark suite such as deduplication, methylation extractor (splitting) or M-bias reports to generate a graphical HTML report page. If several Bismark reports are found in the same folder, a separate report will be generated for each of these, whereby the output filename is derived from the Bismark alignment report file. bismark2report attempts to find optional reports automatically based on the file basename. Here are an example single-end report and paired-end report.

USAGE: bismark2report [options]

OPTIONS

Output directory. Output is written to the current directory if not specified explicitly.

If not specified explicitly, bismark2report attempts to find Bismark report file(s) in the current directory and produces a separate HTML report for each mapping report file. Based on the basename of the Bismark mapping report, bismark2report will also attempt to find the other Bismark reports (see below) for inclusion into the HTML report. Specifying a Bismark alignment report file is mandatory.

If not specified explicitly, bismark2report attempts to find a deduplication report file with the same basename as the Bismark mapping report (generated by deduplicate_bismark) in the current working directory. Including a deduplication report is optional, and using the FILE 'none' will skip this step entirely.

If not specified explicitly, bismark2report attempts to find a single M-bias report file with the same basename as the Bismark mapping report (generated by the Bismark methylation extractor) in the current working directory. Including an M-Bias report is optional, and using the FILE 'none' will skip this step entirely.

If not specified explicitly, bismark2report attempts to find a single nucleotide coverage report file with the same basename as the Bismark mapping report (generated by Bismark with the option --nucleotide_coverage, or bam2nuc directly) in the current working directory. Including a nucleotide coverage statistics report is optional, and using the FILE 'none' will skip this report entirely.

(V) The Bismark Summary Report

This script uses Bismark report files of several (up to hundreds of!?) samples in a run folder to generate a graphical summary HTML report as well as a whopping big table (tab-delimited text) with all relevant alignment and methylation statistics which may be used for graphing purposes in R, Excel or the like. Unless certain BAM files are specified, bismark2summary first identifies Bismark BAM files in a folder (they need to use the Bismark naming conventions) and then automatically detects Bismark alignment, deduplication or methylation extractor (splitting) reports based on the input file basename. If splitting reports are found they overwrite the methylation statistics of the initial alignment report.

USAGE: bismark2summary [options]

This command scans the current working directory for different Bismark alignment, deduplication and methylation extraction (splitting) reports to produce a graphical summary HTML report, as well as a data table, for all files in a directory. Here is a sample Bismark Summary Report.

ARGUMENTS:

OPTIONS:

Basename of the output file (optional). Generate a text file with all relevant extracted values 'basename.txt') as well as an HTML report ('basename.html'). If not specified explicitly, the basename is 'bismark_summary_report'.

Optional HTML report title; use --title "speech marks for text with spaces". Default: 'Bismark Summary Report'.

Displays version information and exits.

Displays this help message and exits.

(VI) Bismark Nucleotide Coverage report (bam2nuc)

The script bam2nuc reads BAM files and calculates the mono- and di-nucleotide coverage of the reads (using the genomic sequence rather than the observed sequence in the reads themselves) and compares it to the average genomic sequence composition. Reads harbouring InDels are not taken into consideration. Mono- or dinucleotides containing Ns are ignored as well.

bam2nuc handles both Bismark single-end and paired-end files (determined automatically). Both BAM and CRAM files should work as input, but please note that Samtools version 1.2 or higher is required for CRAM files.

USAGE: bam2nuc [options] --genome_folder <path> [input.(bam|cram)]

Output directory. Output is written to the current directory if not specified explicitly.

Enter the genome folder you wish to use to extract sequences from (full path only). Accepted formats are FastA files ending with '.fa' or '.fasta'. Specifying a genome folder path is mandatory.

The path to your Samtools installation, e.g. /home/user/samtools/. Does not need to be specified explicitly if Samtools is in the PATH already

Only calculate and extract the genomic sequence composition and exit thereafter. This option will attempt to write the genomic composition table 'genomic_nucleotide_frequencies.txt' to the genome folder or to the output directory instead if that doesn't succeed.

Displays this help message and exits

GENOMIC COMPOSITION

Since the calculation of the average genomic (di-)nucleotide composition may take a while, bam2nuc attempts to write out a file called 'genomic_nucleotide_frequencies.txt' to the genome folder if it wasn't there already. The next time bam2nuc is run it will then use this file instead of calculating the average genome composition again. If writing to the genome folder fails (e.g. because of permission issues) it will be written out to the output directory instead.

OUTPUT FORMAT

bam2nuc writes out a file ending in .nucleotide_stats.txt in the following format (tab-delimited):

(di-)nucleotide count sample    percent sample  count genomic   percent genomic coverage
A       14541   30.91   3768086 30.98   0.004
C       8893    18.90   2321832 19.09   0.004
G       9019    19.17   2318192 19.06   0.004
T       14597   31.02   3754886 30.87   0.004
AA      5008    10.86   1321485 10.86   0.004
AC      2355    5.11    639783  5.26    0.004
AG      2692    5.84    709163  5.83    0.004
AT      4191    9.09    1097652 9.02    0.004
CA      2912    6.32    786744  6.47    0.004
CC      1812    3.93    473900  3.90    0.004
CG      1341    2.91    355535  2.92    0.004
CT      2659    5.77    705653  5.80    0.004
GA      2903    6.30    756411  6.22    0.004
GC      1724    3.74    453607  3.73    0.004
GG      1817    3.94    470732  3.87    0.004
GT      2402    5.21    637436  5.24    0.004
TA      3419    7.42    903441  7.43    0.004
TC      2823    6.12    754531  6.20    0.004
TG      2996    6.50    782761  6.44    0.004
TT      5055    10.96   1314144 10.80   0.004

This file is picked up and plotted by bismark2report automatically if found in the folder in the following manner:

Nucleotide Coverage Plot

(VII) Filtering out non-bisulfite converted reads (filter_non_conversion)

Filtering incomplete bisulfite conversion from Bismark BAM files (optional). This script examines the methylation calls of reads, or read pairs for paired-end sequencing, and filters out reads that exceed a certain threshold of methylated calls in non-CG context (the default is 3). filter_non_conversion looks for a certain number of methylated non-CG calls, but this could potentially also be extended to a percentage for any given read.

Please Note: Be aware that this kind of filtering is not advisable - and will introduce biases - if you work with organisms which exhibit any appreciable levels of non-CG methylation (e.g. most plants).

Writes out a file called nonCG_filtered.bam, also a file called nonCG_removed_seqs.bam as well as a short report how many sequences have been analysed and removed.

USAGE:

filter_non_conversion [options] [Bismark BAM files]

Please also note that for paired-end BAM files filter_non_conversion expects Read 1 and Read 2 to follow each other in consecutive lines! If the file has been sorted by position make sure that you resort it by read name first (e.g. using samtools sort -n)

Deduplicate single-end Bismark BAM files. If not specified the library type is auto-detected.

Deduplicate paired-end Bismark BAM files. If not specified the library type is auto-detected.

The number of methylated cytosines in non-CG context at which reads or read pairs are filtered out. For paired-end files either Read 1 or Read 2 can fail the entire read pair. Default: 3.

The path to your Samtools installation, e.g. /home/user/samtools/. Does not need to be specified explicitly if Samtools is in the PATH already.

Displays this help text end exits.

Displays version information and exits.

If you get stuck at any point or have any questions or comments please contact me via e-mail: felix.krueger@babraham.ac.uk

(VIII) Notes about different library types and commercial kits

Here is a table summarising general recommendations for different library types and/or different commercially available kits. Some more specific notes can be found below.

Technique 5' Trimming 3' Trimming Mapping Deduplication Extraction
BS-Seq ⬜️ ⬜️ ⬜️ --ignore_r2 2
RRBS --rrbs (R2 only) --rrbs (R1 only) ⬜️ (--ignore_r2 2)
RRBS (NuGEN Ovation) special processing special processing ⬜️ --ignore_r2 2
PBAT 6N / 9N (6N / 9N) --pbat ⬜️
single-cell (scBS-Seq) 6N (6N) --non_directional; single-end mode ⬜️
TruSeq (EpiGnome) 8 bp (8 bp) ⬜️ ⬜️
Accel-NGS (Swift) 10 bp (10 bp) ⬜️ ⬜️
Zymo Pico-Methyl 10 bp (10 bp) --non_directional ⬜️

5' Trimming can be accomplished with Trim Galore using:

--clip_r1 <NUMBER> (Read 1) or

--clip_r2 <NUMBER> (Read 2)

3' Trimming can be accomplished with Trim Galore using:

--three_prime_clip_r1 <NUMBER> (Read 1) or

--three_prime_clip_r2 <NUMBER> (Read 2).

SPECIFIC LIBRARY/KIT NOTES

RRBS

RRBS is a specialised technique to only look at CpG rich regions of the genome by using the restriction enzyme MspI (please see this RRBS Guide for some more specifics regarding data processing). For reasons explained in the RRBS-guide, the second last position of all reads before reading into the Illumina adapter exhibits an artificially (not methylated) methylation state as a result of the end-repair reaction. The option --rrbs within Trim Galore removes 2 extra bases whenever adapter contamination has been detected. This 3' end trimming that needs to be carried out for single-end runs or Read 1 of paired-end libraries. Read 2 of paired-end libraries is however not affected by this 3' bias, but instead the first couple of positions on the 5' end of Read 2 suffer from the read-through problem as Read 1 (Read 2 is a mere copy of Read 1), so Read 2 needs to have the first 2 bp removed instead. As of the current development version of Trim Galore (v0.4.2_dev; 12/16/2016) the option --rrbs removes: - 2 bp from the 3' end of single-end and Read 1 of paired-end reads in addition to adapter contamination, and - 2 bp from the 5' end of Read 2 of paired-end reads

RRBS NuGEN Ovation Methyl-Seq System

(Manufacturer's page)

Owing to the fact that the NuGEN Ovation kit attaches a varying number of nucleotides (0-3) after each MspI site Trim Galore should be run WITHOUT the option --rrbs. The trimming is accomplished in a subsequent diversity trimming step afterwards, please see the manufacturer's manual for more details.

PBAT

The amount of bases that need to be trimmed from the 5' end depends on the length of the oligo used for random priming, which - as we know - isn't all that random, and in fact causes misalignments and methylation biases. While the original PBAT paper used 4N oligoes, these days 6N or 9N seem to be most common. Please also see the section 3' Trimming in general below.

Single-cell

The scBS-Seq method uses a PBAT-type protocol but employs five rounds of sequence capture and elongation to amplify the starting material so all four different bisulfite strands (OT, CTOT, OB, CTOB) are sequenced. Since 6N oligos are used to for the random priming step, 6 bp need to be removed from the 5' ends. Since scBS and PBAT libraries tend to result in chimaeric fragments we tend to treat scBS-Seq as single-end reads always. Please also see the section 3' Trimming in general below.

TruSeq DNA-Methylation Kit (formerly EpiGnome)

(Manufacturer's page) This Illumina kit (previously known as EpiGnome kit from epicentre) also employs a post-bisulfite strategy using 6N oligos, but in contrast to the PBAT technique only the standard original top and bottom strands (OT and OB) are sequenced, meaning that Bismark can be run in default (= directional) mode. Even though the random priming is performed with 6N oligoes we often saw that the methylation bias extends to 7 or 8 bp, so trimming 8 bp off the 5' end(s) is recommended initially. Please do have a look at the M-bias plots nevertheless to see of more bases need removing/ignoring during the methylation extraction process. Please also see the section 3' Trimming in general below.

Zymo Pico Methyl-Seq

(Manufacturer's page) The Pico Methyl-Seq kit also uses a random priming step similar to the PBAT // single-cell methods above. This kit uses random tetramers (4N) for the amplification step, however the biases seen in the base composition and M-bias plots indicate that one should trim off at least the first 10 bp from each read. This kit performs three rounds of amplification which yields non-directional libraries (similar to the scBS-Seq protocol), so all four different bisulfite strands (OT, CTOT, OB, CTOB) are present in the library. According to the manufacturer, the library construction is designed for a starting input material of 100 ng, but can be scaled up or down (to 100 pg). Please also see the section 3' Trimming in general below.

Swift

Manufacturer's page The Accel-NGS Methyl-Seq protocol uses Adaptase technology for capturing single-stranded DNA in an unbiased (again, not that unbiased actually...) manner. Also here, the first ~10 bp show extreme biases in sequence composition and M-bias, so trimming off at least 10 bp is advisable (please check the M-bias plot if even more is needed). Please also see the section 3' Trimming in general below.

Random priming and 3' Trimming in general

As we have seen before, the random priming of post-bisulfite methods (such as PBAT, scBS-Seq, EpiGnome, Pico Methyl, Accel etc.) introduces errors, indels and methylation biases that may detrimentally affect your mapping efficiencies and methylation calls. These problems are fairly easy to spot at the 5' ends of reads because all reads will equally suffer from the problems at the same positions at the start (5' end) of reads. The same problems of random priming (indels, mispriming) will however most likely occur on both sides of the fragment to be sequenced, but it is doubtful that one would be able to spot these problems on the 3' end of reads because the problems would be expected on the 3' end of reads just before reading through into the adapter, and this may occur - at different positions in the read (depending on how short the fragment was) - at different positions within the read because of quality trimming in addition to adapter read-through contamination - not at all within the read length (whenever a fragment is longer than the sequenced read length) - at the 3' end even without hitting the adapter (i.e. just before the adapter)

I guess there is a trade-off between accepting that a certain proportion of the reads may have a few biased biased positions towards their 5' ends, and preemptively trimming the 3' end by the same amount of bases as the 5' end. As a general rule it is probably safe to say that the shorter the average insert size of a library - the more of a problem the bias is. We have e.g. seen Pico Methyl libraries where ~80% of all fragments were shorter than 100bp, so a 2x125bp run would most likely be affected by the random priming bias on the 5' and 3' ends in nearly all fragments sequenced. We realise that trimming off say 10 bp from the 5' end and 3' end of a 100 bp read already removes 20% of the actually sequenced data, but this is the price you have to pay for using post-bisulfite kits...

For these reasons we have put the 3' Trimming values in the table above in (parentheses) as a reminder that you should probably perform 3' trimming of the data as well.

4) APPENDIX - Full list of options

Appendix (I): Bismark Genome Preparation

A full list of options can also be viewed by typing: bismark_genome_preparation --help

USAGE: bismark_genome_preparation [options] <arguments>

OPTIONS:

ARGUMENTS:

Appendix (II): Bismark

A brief description of Bismark and a full list of options can also be viewed by typing: bismark --help

USAGE:

bismark [options] --genome <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}

ARGUMENTS:

OPTIONS:

Input:
Alignment:
Bowtie 1 Reporting:
Output:
Other:
BOWTIE 2 SPECIFIC OPTIONS
Bowtie 2 alignment options:
Bowtie 2 paired-end options:
Bowtie 2 Effort options:
Bowtie 2 parallelization options:
Bowtie 2 Scoring options:
Bowtie 2 Reporting options:

Appendix (III): Bismark Methylation Extractor

A brief description of the Bismark methylation extractor and a full list of options can also be viewed by typing bismark_methylation_extractor --help

USAGE: bismark_methylation_extractor [options] <filenames>

ARGUMENTS:

OPTIONS:

bedGraph specific options:
Genome-wide cytosine methylation report specific options:

OUTPUT

The bismark_methylation_extractor output is in the form (tab delimited, 1-based coords):
<seq-ID> <methylation state*> <chromosome> <start position (= end position)> <methylation call>

  Methylated cytosines receive a '+' orientation,
Unmethylated cytosines receive a '-' orientation.
The bedGraph output (optional) looks like this (tab-delimited, 0-based start, 1-based end coords):
track type=bedGraph (header line)
<chromosome> <start position> <end position> <methylation percentage>
The coverage output looks like this (tab-delimited; 1-based genomic coords):
<chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated>
The genome-wide cytosine report (optional) is tab-delimited in the following format (1-based coords):
<chromosome> <position> <strand> <count methylated> <count unmethylated> <C-context> <trinucleotide context>

Appendix (IV): Bismark reports for the test data set

Please note that this has been run with a fairly early version however I wouldn't expect the numbers to change much.

Using Bowtie:
Running Bismark with the default options (e.g. bismark --bowtie1 /data/public/Genomes/Human/GRCh37/ test_data.fastq) should result in this mapping report:

Bismark report for: test_data.fastq (version: v0.7.8)

Option '--directional' specified: alignments to complementary strands will be ignored (i.e. not performed!)
Bowtie was run against the bisulfite genome of /data/public/Genomes/Human/GRCh37/ with the specified options: -q -n 1 -k 2 --best --chunkmbs 512

Final Alignment report
======================
Sequences analysed in total: 10000
Number of alignments with a unique best hit from the different alignments: 6361 Mapping efficiency: 63.6%
Sequences with no alignments under any condition: 2626
Sequences did not map uniquely: 1013
Sequences which were discarded because genomic sequence could not be extracted: 0
Number of alignments to (merely theoretical) complementary strands being rejected in total: 0

Number of sequences with unique best (first) alignment came from the bowtie output:
CT/CT: 3187 ((converted) top strand)
CT/GA: 3174 ((converted) bottom strand)
GA/CT: 0    (complementary to (converted) top strand)
GA/GA: 0    (complementary to (converted) bottom strand)

Final Cytosine Methylation Report
================================= 
Total number of C's analysed: 52942

Total methylated C's in CpG context: 1740
Total methylated C's in CHG context: 36 
Total methylated C's in CHH context: 171

Total C to T conversions in CpG context: 1027
Total C to T conversions in CHG context: 12889
Total C to T conversions in CHH context: 37079

C methylated in CpG context: 62.9% 
C methylated in CHG context: 0.3% 
C methylated in CHH context: 0.5%
Using Bowtie 2:
Running Bismark with the following options (e.g. bismark --score-min L,0,-0.6 /data/public/Genomes/Human/GRCh37/ test_data.fastq) should result in this mapping report:

Bismark report for: test_data.fastq (version: v0.7.8)
Option '--directional' specified: alignments to complementary strands will be ignored (i.e. not performed!)
Bowtie2 was run against the bisulfite genome of /data/public/Genomes/Human/GRCh37/ with the specified options: -q -- score-min L,0,-0.6 --ignore-quals

Final Alignment report
======================
Sequences analysed in total: 10000

Number of alignments with a unique best hit from the different alignments: 5658 Mapping efficiency: 56.6%
Sequences with no alignments under any condition: 2893
Sequences did not map uniquely: 1449
Sequences which were discarded because genomic sequence could not be extracted: 0
Number of alignments to (merely theoretical) complementary strands being rejected in total: 0

Number of sequences with unique best (first) alignment came from the bowtie output:

CT/CT: 2820 ((converted) top strand) 
CT/GA: 2838 ((converted) bottom strand)
GA/CT: 0    (complementary to (converted) top strand) 
GA/GA: 0    (complementary to (converted) bottom strand)

Final Cytosine Methylation Report
================================= 
Total number of C's analysed: 45985

Total methylated C's in CpG context: 1550
Total methylated C's in CHG context: 34 
Total methylated C's in CHH context: 126

Total C to T conversions in CpG context: 844
Total C to T conversions in CHG context: 11368
Total C to T conversions in CHH context:32063

C methylated in CpG context: 64.7% 
C methylated in CHG context: 0.3% 
C methylated in CHH context: 0.4%

Credits

Bismark was written by Felix Krueger at the Babraham Bioinformatics Group.

Babraham Bioinformatics