MethylDackel

Chromosome data is a sequence of individual values, short runs, and long runs. Positions are implicit: the first value starts at position 0 and the sequence ends after chromLength values have been processed.

There are three ways to encode data:

1. Individual Value: A single byte (0-100).
2. Short Run (length 2-154): Two bytes total. The first byte (runLength) is an offset: actual_length = runLength - 99. For example, a runLength of 101 represents a run of length 2. The second byte is the runValue (0-100).
3. Long Run (length 155-65535): Four bytes total. Starts with a 255 flag, followed by a little-endian 16-bit runLength and a single-byte runValue (0-100).

BBM is a binary format designed for storing a single data track for a genome. It is optimized for integer values in the range 0-100 (inclusive), such as Bismap percent mappabilities. The format uses run-length compression to achieve high compression ratios (approximately 7.7:1 for human genome datasets).

Key Constraints:

• Values must be integers between 0 and 100.
• Maximum run length per encoded unit is 65,535 bases. Longer runs must be represented as sequences of shorter runs.
• The format uses a specific marker (255) to distinguish between short and long runs.

By default, extract produces per-Cytosine metrics. To merge these into per-CpG or per-CHG metrics, use the --mergeContext option during extraction.

If you already have per-Cytosine bedGraph files, you can convert them to merged metrics using the MethylDackel mergeContext command.

MethylDackel extract --mergeContext reference_genome.fa alignments.bam

MethylDackel can be installed via Bioconda or Galaxy. To compile from source, you need a C compiler, htslib (version 1.11 or later), and libBigWig (the static library).

Basic compilation:

git clone https://github.com/dpryan79/MethylDackel.git
cd MethylDackel
make LIBBIGWIG="/some/path/to/libBigWig.a"
make install prefix=/some/installation/path

If the linker cannot find htslib headers or libraries, specify them using CFLAGS and LIBS:

make install CFLAGS="-O3 -Wall -I/some/path/include " LIBS="-L/some/path/lib" prefix=/some/installation/path LIBBIGWIG="/some/path/to/libBigWig.a"

MethylDackel provides tools to identify and correct methylation bias (changes in observed methylation rates at read ends).

1. Identify bias: Use MethylDackel mbias to generate SVG plots showing methylation levels and confidence intervals for each strand.

   MethylDackel mbias reference_genome.fa alignments.sorted.bam output_prefix

2. Apply trimming: Once bias regions are identified, use the following options in extract to ignore specific bases:
   • --OT, --OB, --CTOT, --CTOB: Specify exact positions within each alignment to ignore.
   • --nOT, --nOB, --nCTOT, --nCTOB: Specify a fixed number of bases to trim from the ends of alignments (e.g., --nOT 5,5,5,5 to trim 5 bases from both ends of both strands).

The perRead subcommand produces a tab-separated file containing methylation levels for individual reads. This is useful for specific methods like Hi-Culfite.

Output columns:

1. read name
2. chromosome
3. position
4. CpG methylation (%)
5. number of informative bases

Note: Fragments longer than 10kb are not currently handled correctly.

MethylDackel perRead reference_genome.fa alignments.bam

To eliminate reads that cannot be uniquely placed, use a bigWig file containing mappability data.

• -M <bigwig>: Use a bigWig file for filtering.
• -B <bbm>: Use a Binary BisMap (BBM) file for faster execution and lower disk usage. BBM files can be created during extraction using -O or -N.
• -t <threshold>: Minimum number of bases that must have mappability $\ge$ mappabilityThreshold (default: 15).
• -b <mappabilityThreshold>: The mappability threshold (default: 0.01).

Each chromosome record contains the name and the compressed data track for that chromosome. The structure is:

[nameLength][chromName][NUL][chromLength]{chrom data}

• nameLength: Unsigned little-endian 16-bit integer (number of characters in the name, excluding the null terminator).
• chromName: 8-bit ASCII characters.
• NUL: A \0 byte null terminator (used for error checking).
• chromLength: Unsigned little-endian 32-bit integer representing the total number of bases in the chromosome.
• {chrom data}: The run-length compressed data values.

You can modify the standard bedGraph output using the following flags in MethylDackel extract:

Note: These options can be combined with --mergeContext during extraction, but mergeContext cannot be applied to these specialized formats after they are generated.

A BBM file follows this top-level structure:

[version][chromCount]{chrom record}{chrom record}...

• version: An unsigned character (ID 1 for this specification).
• chromCount: An unsigned little-endian 32-bit integer representing the total number of chromosomes in the file.
• {chrom record}: A repeating block containing data for a specific chromosome.

<beginning of file>[version][chromCount]{chrom record}{chrom record}{chrom record}{chrom record}...<end of file>

The following table defines the byte-level encoding for chromosome data values and runs:

| Type | Byte Sequence | Details | | :--- | :--- | : | | Individual Value | [value] | Single unsigned byte (0-100). | | Short Run | [runLength][runValue] | runLength is unsigned byte (101-254), where length = runLength - 99. runValue is 0-100. | | Long Run | [255][runLength][runValue] | 255 is the flag. runLength is little-endian 16-bit unsigned integer. runValue is 0-100. |

The extract command processes a coordinate-sorted and indexed BAM or CRAM file to extract per-base methylation metrics. It requires an indexed FASTA file of the reference genome.

Basic usage:

MethylDackel extract reference_genome.fa alignments.bam

This produces a bedGraph file (e.g., alignments_CpG.bedGraph) containing per-base CpG metrics.

Key Options:

• -o <prefix>: Specify an alternate output filename prefix.
• --CHH: Include CHH context metrics.
• --CHG: Include CHG context metrics.
• --noCpG: Ignore Cytosines in CpG contexts.
• -q <MAPQ>: Minimum mapping quality (default: 10).
• -p <Phred>: Minimum Phred score (default: 5).
• --minDepth <depth>: Only output sites with at least this many methylation calls (useful with --mergeContext).
• --mergeContext: Merges individual Cytosine metrics into per-CpG or per-CHG metrics.