🧬 methurator

Methurator is a Python package designed to estimate sequencing saturation for reduced-representation bisulfite sequencing (RRBS) data.

Although optimized for RRBS, methurator can also be used for whole-genome bisulfite sequencing (WGBS) or other genome-wide methylation data (e.g. EMseq). However, this data we advise you to use Preseq package.

📑 Table of Contents

• 1. Dependencies and Notes
• 2. Installation
• 3. Quick Start
  • Step 1 — Downsample BAM files
  • Step 2 — Plot the sequencing saturation curve
• 4. Command Reference
  • downsample command
  • plot command
• 5. Example Workflow
• 6. How do we compute the sequencing saturation?

1. Dependencies and Notes

• methurator uses SAMtools and MethylDackel internally for BAM subsampling, thus they need to be installed.
• When --genome is provided, the corresponding FASTA file will be automatically fetched and cached.
• Temporary intermediate files are deleted by default unless --keep-temporary-files is specified.

2. Installation

You can install methurator in several ways:

Option 1: Install via pip

pip install methurator

Option 2: Install via BioConda

conda create -n methurator_env bioconda::methurator
conda activate methurator_env

Option 3: Use the BioContainer

docker pull quay.io/biocontainers/methurator:0.1.8--pyhdfd78af_0
docker run quay.io/biocontainers/methurator:0.1.8--pyhdfd78af_0 methurator -h

3. Quick Start

Step 1 — Downsample BAM files

The downsample command performs BAM downsampling according to the specified percentages and coverage.

methurator downsample --fasta tests/data/genome.fa tests/data/Ecoli.csorted.bam

This command generates three summary files:

• CpG summary — number of unique CpGs detected in each downsampled BAM
• Reads summary — number of reads in each downsampled BAM
• Summary yml - a YAML file which contains all data above in a single file. It also contains run metadata for the sake of reproducibility.

Example outputs can be found in tests/data.

Step 2 — Plot the sequencing saturation curve

Use the plot command to visualize sequencing saturation:

methurator plot --summary tests/data/methurator_summary.yml

4. Command Reference

downsample command

plot command

5. Example Workflow

# Step 1: Downsample BAM file
methurator downsample --genome hg19 my_sample.bam

# Step 2: Plot saturation curve
methurator plot --summary output/methurator_summary.yml

Finally, you will get (within the output/plots) directory an html file containing the sequencing saturation plot, similarly to the following example (also available as interactive html file here):

6. How do we compute the sequencing saturation?

To calculate the sequencing saturation of an RRBS sample, we adopt the following strategy. For each sample, we downsample it according to 4 different percentages (default: 0.1,0.2,0.4,0.6,0.8). Then, we compute the number of unique CpGs covered by at least 3 reads and the number of reads at each downsampling percentage.

We then fit the following curve using the scipy.optimize.curve_fit function:

We chose the arctangent function because it exhibits an asymptotic growth similar to sequencing saturation. For large values of $\text{x}$ (as $\text{x} \to \infty$), the asymptote corresponds to the theoretical maximum number of unique CpGs covered by at least 3 reads and can be computed as:

Finally, the sequencing saturation value can be calculated as following:

This approach allows estimation of the theoretical maximum number of CpGs that can be detected given an infinite sequencing depth, and quantifies how close the sample is to reaching sequencing saturation.