6. Methylation and Variant Calling

Methylation and sequence variant calling are performed simultaneously in the gemBS pipeline with the bs_call program. Taking account of the sequence quality scores, under and over conversion probabilities and the observed bases, bs_call finds the most likely genotype (maximizing the likelihood over the unknown methylation parameter for each genotype) and then reports the methylation conditional on the called genotype. Standard quality scores for variant calling are provided, and the native output format is a VCF/BCF file with additional custom fields that provide methylation specific information (the specifications of the custom format used by gemBS can be seen here).

The bs_call program generates several format fields in the VCF output files that are related to the quality of the genotype calls:

Field ID



Genotype Likelihood


Phred scaled conditioanl genotype quality


Read Depth (non converted reads only)


RMS Mapping Quality


Quality By Depth (Variant quality / read depth (non-converted reads only))


Phred scaled log p-value from Fishers exact test of strand bias

Apart from GOF, these are taken from GATK best practices for variant filtering. By default, bs_call applies a series of filters based on these fields, and the FILTER column in the VCF file is set accordingly. The filters applied are:

Filter ID



Genotype Quality below 20


Quality By Depth below 2


Fisher Strand above 60


RMS Mapping Quality below 40

The default filter cutoffs are also from GATK best practices. Note that only soft filtering is applied - no variants are removed from the VCF file as a result of the filtering. However, the generation of the gemBS style methylation report files and SNP extraction files only use sites that pass all filters. This does not apply to the generation of the ENCODE style methylation report files (see 7. Methylation & SNP Extraction).

6.1. Introduction

Standard pre-processing

By default, bs_call uses only uniquely mapping reads for the analysis (as determined by a minimum threshold on the MAPQ score provided by the mapper). Duplicate reads can be determine by bs_call as reads with the same start and end points (and the same library barcode if this information is available), but bs_call will also take account of the duplicate flag in the BAM files. Reads marked or detected by bs_call as being duplicates are optionally removed prior to calling. For read pairs that overlap, the overlapping section will be collapsed into a single copy. Sequence bases with quality less than a given cutoff will be removed priori to further analysis. Read pairs that do not form proper pairs will normally be discarded. Bases at the extremities of reads can be marked so that they are not used for methylation prediction (this is to avoid methylation artifacts due to end repair during sample preparation. All of the above behaviours can be modified using the program parameters.

Chromosomes, contigs and pools

Separate chromosomes or contigs can be processed by bs_call independently. If multiple computing resources are available then performing the calling in parallel can speed up the processing (although users should look at the section on Practical Considerations where this is discussed more thoroughly). For organisms with well studied genomes, such as humans, most of the reference sequence is organized into chromosomes. However for less well studied organisms, instead of a small number of chromosomes there can be tens of thousands of small contigs. Even in the human GRCh38 build there are almost 200 separate contigs in total. Performing variant calling by contig can therefore be problematic due to the number of contigs, and often it is advisable to either call all contigs together for a single sample, or group small contigs together into pools for analysis purposes.

To get around this, gemBS uses a contig pool strategy that is managed with almost no user intervention. Contigs smaller than a certain cutoff (configuration key contig_pool_limit; set by default to 25Mb) are pooled together to simplify processing. Contigs / chromosomes larger than this are processed separately (technically they are in a pool of one). gemBS generates the pools itself, and performs all calling by pool. By default the user does not have to take care of this. As with the mapping stage, gemBS will normally call all available pools on all available samples, merging the results for completed samples into a single sample BCF file. Calling can be restricted, however, by sample or pool.


Note that by default, gemBS will perform calling on all contigs in the main reference file (any contigs from files in the extra_references list will not be passed to the calling analysis). If the main reference contains additional sequences then these should be excluded from calling by setting the omit_contigs variable.

Conversion rates

The calling process uses the estimated under and over conversion rates in the statistical model of genotype and methylation. These values should be supplied to the caller using the conversion configuration key. The value should be one or two comma separated numbers between 0 and 1 with the under and over conversion rates (in that order). If a single value is present it is assumed to be the under conversion rate, and the over conversion is set to a default (0.05). If the special value auto is used for the conversion key, gemBS will try and obtain sample specific estimates of the conversion rates from the mapping results. For this to work (a) there have to be conversion control sequences in the sequencing library, (b) these have to have been included in the reference sequence and (c) the names of the control sequences have to had been supplied to gemBS when the mapping was performed. If one or both of the values can not be estimated then gemBS will use defaults (0.01 and 0.05) for the missing value(s).


Most computers today are multi-core, and a typical workstation (or even a laptop) can have tens of threads available. The mapping process is inherently parallelizable as each read pair can be processed independently. Running the GEM3 mapper with 24 cores will be close to 24 times faster (ignoring startup times) than running on a single core. Calling is less easy to parallelize; on a single contig bs_call can benefit from 3-4 cores, but more than this has little effect on the running time. Even with 4 cores, the speed up is generally much less than a four fold increase over the single threaded case. For the calling, therefore, it can make sense to run multiple pools on the same compute host. For this purpose gemBS has the jobs configuration key; during the calling gemBS will run jobs separate copies of the caller simultaneously.

As with the mapping process, by default it is possible to run gemBS on multiple compute hosts on the same analysis directory with multiple jobs running on each host, and the multiple instances will cooperate with each other, dividing the available tasks between them.

6.2. The calling command: gemBS call

As with the mapping command, by default the calling requires no parameters and will try and perform all available calling and merging tasks. By default, all of the information about input and output files, calling parameters etc. come from the original configuration files. The command line options for the call command are present to enable restriction of the mapping to a given sample or contig pool, or to perform ad-hoc adjustment of the calling parameters, but it is intended that normally no parameters are required.

The command template is as follows (with all parameters being optional):

 gemBS call [FLAGS] [OPTIONS]

 -u, --keep-duplicates          Do not merge duplicate reads
 -U, --ignore-duplicate-flag    Ignore duplicate flag from SAM/BAM files
 -k, --keep-unmatched           Do not discard reads that do not form proper pairs
 -r, --remove                   Remove individual BAMs after merging
     --md5                      Perform calculation of md5 sums only
     --no-md5                   Do not automatically calculate md5 sums
     --index                    Perform indexing of final BCF only
     --no-index                 Do not automatically calculate index of final BCF
 -1, --haploid                  Force genotype calls to be haploid
     --auto-conversion          Try to calculate conversion rates from data
     --merge                    Perform merge BCF step only
     --no-merge                 Do not automatically merge BAMs
     --benchmark_mode           Omit dates etc. from output to make comparison simpler
 -h, --help                     Prints help information

 -n, --sample <SAMPLE>...                 Name of sample
 -b, --barcode <BARCODE>...               Barcode of sample
 -q, --mapq-threshold <MAPQ_THRESHOLD>    Threshold for MAPQ scores
 -Q, --qual-threshold <QUAL_THRESHOLD>    Threshold for base quality scores
 -g, --right-trim <BASES [,BASES]>...     Bases to trim from right of read pair (give 2 values for read specific
 -f, --left-trim <BASES, [,BASES]>...     Bases to trim from left of read pair (give 2 values for read specific
 -d, --tmp-dir <PATH>                     Temporary directory to perform sorting operations
 -t, --threads <THREADS>                  Number of threads for calling pipeline
     --call-threads <THREADS>             Number of threads for methylation caller
     --merge-threads <THREADS>            Number of threads for merging
 -j, --jobs <JOBS>                        Number of parallel jobs
 -c, --cores <CORES>                      Number of cores for a job (default - available cores / jobs)
 -T, --time <TIME>                        Time required for a job
 -m, --memory <MEMORY>                    Memory required for a job (default - available memory / jobs)
 -e, --species <SPECIES>                  Species name
 -D, --dbsnp-index <FILE>                 dbSNP processed index file
 -C, --conversion <UNDER OVER>...         set conversion rates (under over)
 -R, --reference-bias <BIAS>              set bias to reference homozygote
     --pool <POOL>...                     Contig pool for methylation calling

Possible parameters for the gemBS call command are given below. Where present, the Key column shows the keys that control this parameter in the configuration files.






Location of input BAM files



Location of output BCF files



Location of FASTA reference file



Location of dbSNP index file



Minimum cutoff for read MAPQ



Minimum cutoff for base quality



Left trim x bases from reads before calling



Right trim x bases from reads before calling






Keep duplicates



Keep unmatched (non-paired) reads



Remove contig pools BCFs after merge step



Force all genotype calls to haploid



Multiplicative weight for homozygous reference genotype probability



under-, over- conversion rates.



Restrict calling to a list of contigs (actually, to all pools than contain these contigs)



File with list of contigs and sizes



Maximum sizes limit for small contigs (to be pooled)



Number of threads for call stage



Number of threads for merge stage



Number of jobs for merge stage



Number of core for merge stage



Amount of memory for merge stage



Amount of memory for merge stage


Number of cores


Number of threads


Number of parallel jobs


Memory usage


Time required

The calling process (as with mapping) involves multiple steps: The calling itself using bs_call, combining the BCF files from the individual contig pools, using bcftools concat, indexing the merged bcf (unlike with samtools, bcftools does not have the option to generate the index at the same time as combining the files), and calculating the md5 signature of the BCF. While the combine step can use multiple threads, there is not much benefit with providing more that 2-3 threads. The resources used just for the merge step can thereore be set individually using the merge_threads, merge_cores, merge_jobs, merge_cores and merge_time keys, which can allow for more efficient scheduling of jobs. The indexing of the BCFs with bcftools index and calculation of md5 signatures are not be multithreaded and the number of threads for these steps are set to 1.

As with the map command, normally the call command is run without any parameters so that all necessary calling, combining and indexing is performed. However it is possible to specify the individual steps, This is intended for integration in pipelines to allow, for example, the calling of different samples / contig pools to be performed using different computational nodes.