BinAC2/Software/Bismark

The main documentation is available on the cluster via module help bio/bismark/3.1.0. Most software modules for applications provide working example batch scripts.

Description	Content
module load	bio/bismark/3.1.0
License	GPL-3.0 - see [1]
Citing	Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27(11), 1571-1572 (2011)
Links	Bismark Documentation: [2] Bismark Repository: [3]

Description

Bismark maps bisulfite-converted reads (WGBS, RRBS, EM-seq, ...) to a reference genome and calls cytosine methylation in CpG, CHG and CHH context.

The module bio/bismark/3.1.0 is the Rust rewrite of Bismark. Upstream validates its default output as byte-identical to Perl Bismark v0.25.1. We checked this on BinAC 2 with the upstream test data (all outputs identical) and with 10 million human EM-seq read pairs (identical alignments). It is much faster than the Perl version, mostly in the methylation extraction:

Step (10M read pairs, GRCh38, 16 cores)	Perl 0.25.1 (--multicore 8)	Rust 3.1.0 (-p 8)
Alignment	785 s	699 s
Deduplication	190 s	169 s
Methylation extraction (with cytosine report)	3067 s	293 s
Total	4042 s	1162 s

All classic tool names are available: bismark, bismark_genome_preparation, deduplicate_bismark, bismark_methylation_extractor, bismark2bedGraph, coverage2cytosine, bam2nuc, bismark2report, bismark2summary and more. Your existing commands keep working, but see below for the thread options.

The module loads bio/bowtie2/2.5.5 automatically. HISAT2 (--hisat2) and minimap2 (--minimap2) are optional and must be loaded separately. samtools is not needed.

Quick start

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --mem=32G
#SBATCH --time=04:00:00

module load bio/bismark/3.1.0

GENOME=/path/to/genome_folder        # prepared with bismark_genome_preparation
cd "$TMPDIR"                         # fast local scratch of the job

bismark --genome "$GENOME" -p 15 \
    -1 /path/to/sample_R1.fq.gz -2 /path/to/sample_R2.fq.gz
deduplicate_bismark -p --parallel 4 sample_R1_bismark_bt2_pe.bam
bismark_methylation_extractor -p --gzip --bedGraph \
    sample_R1_bismark_bt2_pe.deduplicated.bam

mkdir -p /path/to/results
cp *.bam *.txt *.gz *.html /path/to/results/ 2>/dev/null

For 50 million read pairs against the human genome this takes about one hour.

Preparing the genome

module load bio/bismark/3.1.0
bismark_genome_preparation --parallel 16 /path/to/genome_folder/

• Large genomes: Bowtie 2 builds a large index (.bt2l) automatically for genomes above 2 GiB (human, mouse and most mammals); you do not need --large-index. Giving it anyway does no harm. Bismark finds the large index automatically when aligning.
• Use one file extension for all FASTA files in the genome folder. Bismark reads only the first non-empty group of .fa, .fa.gz, .fasta, .fasta.gz and silently ignores the others. For example, if your genome is genome.fasta.gz and you add a spike-in as lambda.fa.gz, only the spike-in is indexed. Rename the files so they all end the same way.
• Check the result: grep -c '^>' /path/to/genome_folder/Bisulfite_Genome/CT_conversion/genome_mfa.CT_conversion.fa must show the number of sequences of your genome.
• Bismark builds two indexes (C to T and G to A converted genome); --parallel N builds both at the same time with N threads each. Request --cpus-per-task=32 --mem=64G for --parallel 16: each index build needs about 25 GiB, and both run at the same time. For the human genome this takes less than one hour.
• Prepare the index once and reuse it for all samples.

Choosing threads and resources

All recommendations below were measured on BinAC 2 standard nodes (AMD EPYC 7543) with human GRCh38 and EM-seq reads (2 x 151 bp). Other genomes and read lengths need different times and memory, but the thread recommendations apply in general.

Recommended job size

Request --cpus-per-task=32 and --mem=32G. On BinAC 2, Slurm counts hardware threads as CPUs, so this gives you 16 physical cores with 2 hardware threads each. In our tests this was the most efficient job size: with twice as many cores, the alignment of one sample was only 1.26x faster.

If you have many samples, run one job per sample with these settings rather than fewer jobs with more cores.

Alignment: use -p, not --multicore

Bismark runs 2 Bowtie 2 processes (4 with --non_directional), each with -p N threads. Choose

-p = cpus-per-task / 2 - 1

so -p 15 for --cpus-per-task=32. The 2 remaining hardware threads are used by Bismark itself.

Setting (10M read pairs, 16 cores)	Time
-p 8	720 s
-p 15	551 s
-p 16	599 s

• --multicore splits the input and starts a full Bismark with its own copy of the index per worker. It was clearly slower than -p with the same number of threads (559 s vs 374 s on a full node) and needs much more memory. Results are identical either way.
• Do not use more threads than you requested: oversubscribed jobs get slower, not faster.

Deduplication

Use --parallel 4. For 10M read pairs this reduces the time from 169 s to 36 s. The output is identical.

Methylation extraction

• --multicore does not make the extraction of BAM files faster. The extractor already uses about 4 cores on its own.
• Use --cytosine_report only if you need it. It writes every CpG of the genome, including uncovered ones, and takes most of the extraction time and memory:

Extraction (10M read pairs)	Time	Memory
--bedGraph	106 s	0.8 GiB
--bedGraph --cytosine_report	293 s	3.4 GiB

Many downstream tools can read the coverage file (*.bismark.cov.gz) directly. If you need the cytosine report later, create it from the coverage file with coverage2cytosine.

Time, memory and scratch space

Measured with --cpus-per-task=32, -p 15, --parallel 4, human GRCh38:

Read pairs	Alignment	Deduplication	Extraction	Extraction with cytosine report	Scratch space during alignment
10 million	9 min	0.6 min	1.6 min	4.8 min	10 GB
25 million	23 min	1.4 min	3.9 min	7.2 min	26 GB
50 million	45 min	2.6 min	7.6 min	11 min	47 GB

Rules of thumb:

• Time: about 1 minute per million read pairs for alignment, deduplication and extraction together. In jobs that share a node with other jobs, plan about 15% more.
• Memory: about 12 GiB for aligning against the human genome, independent of the number of reads. Deduplication and extraction stay below 4 GiB up to 50M read pairs. A job with --mem=16G processed 50M read pairs successfully; --mem=32G leaves a safe margin.
• Scratch space: about 1 GB per million read pairs during the alignment (temporary converted read files), plus the input and output files. Work in the job's $TMPDIR (local disk) and copy the results back at the end. Alignment reads and writes a lot of data (50M pairs: about 180 GB read, 130 GB written); doing this on the parallel file system is slower and loads it for all users.

Moving from Perl Bismark 0.25.1

• Results are the same as with Perl Bismark 0.25.1. The @PG line of the BAM header therefore still reports VN:v0.25.1.
• Replace --multicore by -p for the alignment, and add --parallel 4 to deduplicate_bismark.
• If you ran Perl Bismark with --multicore, a tiny fraction of your methylation calls may differ from a new run (399 of 58.6 million CpG report lines in our test). Perl's --multicore writes the BAM in a different order, so deduplication keeps different but equivalent duplicates. The alignments themselves are identical.

Support

If you have questions or problems, contact hpcmaster [at] uni-tuebingen.de.