methurator 0.1.5

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

🧬 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. Pip 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. Pip installation

pip install methurator

---

3. Quick Start

Step 1 — Downsample BAM files

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

methurator downsample --genome hg19 --bam test_data/SRX1631721.markdup.sorted.csorted.bam

This command generates two summary files:

• CpG summary — number of unique CpGs detected in each downsampled BAM
• Reads summary — number of reads in each downsampled BAM

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 \
  --cpgs_file tests/data/cpgs_summary.csv \
  --reads_file tests/data/reads_summary.csv

---

4. Command Reference

downsample command

Argument	Description	Default
--bam	Path to a single .bam file.	—
--bamdir	Directory containing multiple BAM files.	—
--outdir	Output directory.	./output
--fasta	Path to the reference genome FASTA file. If not provided, it will be automatically downloaded based on --genome.	—
--genome	Genome used for alignment. Available: hg19, hg38, GRCh37, GRCh38, mm10, mm39.	—
--downsampling-percentages, -ds	Comma-separated list of downsampling percentages between 0 and 1 (exclusive).	0.1,0.25,0.5,0.75
--minimum-coverage	Minimum CpG coverage to consider for saturation. Can be a single integer or a list (e.g. 1,3,5).	3
--keep-temporary-files	If set, temporary files will be kept after analysis.	False

---

plot command

Argument	Description	Default
--cpgs_file	Path to the CpG coverage summary file.
--reads_file	Path to the reads coverage summary file.
--outdir	Output directory.	./output

---

5. Example Workflow

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

# Step 2: Plot saturation curve
methurator plot \
  --cpgs_file output/cpgs_summary.csv \
  --reads_file output/reads_summary.csv

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.25,0.5,0.75). 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:

$$ y = \beta_0 \cdot \arctan(\beta_1 \cdot x) $$

We chose the arctangent function because it exhibits an asymptotic growth similar to sequencing saturation. For large values of $x$ (as $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:

$$ \text{asymptote} = \beta_0 \cdot \frac{\pi}{2} $$

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

$$ \text{Saturation} = \frac{\text{Number of unique CpGs (≥3 counts)}}{\text{Asymptote}} $$

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.