Quantification of methylation data

Quantifies methylation from array data, using sesame or minfi.

Overview¶

This notebook implements two methods to quantify methylation data, using sesame and minfi. sesame is the recommended choice over minfi.

Procedure	`minfi`	`sesame`
SNP/Cross reaction removal	dropLociWithSnps + manual removal	Q (qualityMask)
sample quality	detectionP + mean	sesameQC_calcStats + “detection” + frac_dt
Bias correction	preprocessQuantile	D ( dyeBiasNL)
Probe quality	detectionP	“P (pOOBAH Detection p-value masking using oob)”
Background substraction	NA	B (noob)

`Sesame`¶

Getting the beta value from EPIC450 IDAT for 750 samples from 3000 wells take ~40 mins. Based on sesame documentation, the processing procedure suitable for human on EPIC 450 and 850 platform is “QCDPB”

The code for each processing procedure are as followed:

Code	Name	Detail
Q	qualityMask	Mask probes of poor design
C	inferInfiniumIChannel	Infer channel for Infinium-I probes
D	dyeBiasNL	Dye bias correction (non-linear)
P	pOOBAH	Detection p-value masking using oob
B	noob	Background subtraction using oob

Other potential procedures are

Code	Name	Detail
0	resetMask	Reset mask to all FALSE
G	prefixMaskButCG	Mask all but cg- probes
H	prefixMaskButC	Mask all but cg- and ch-probes
E	dyeBiasL	Dye bias correction (linear)
I	detectionIB	Mask detection by intermediate beta values

M value is calculated as M = log2(beta/(1-beta)). Beta values of exactly 0 or 1 are replaced with the next min/max value among the beta matrix, which is based on here

`minfi`¶

By default, for Infinium MethylationEPIC the data will be annotated based on hg38 using this annotation, alternatively user can set the --hg-build parameter back to 19 to use the hg19 annotation.

For 450K data however, only hg19 annotation is availble, which is what minfi requires. Everything is reannotated to hg38 in the next step regardless.

All the IDAT file in the specified folder and sub-folder will be loaded for samples in input sample CSV file
The methylation data samples will first be filtered based on bisulphite conversation rate. This operation is done using the bscon function from watermelon package
samples will then be filtered based on a detection pvalue, which indicates the quality of the signal at each genomics position
Stratified Quantile Normalization will then be applied.
features will be filtered if they are on sex chr, known to be cross-reactive,maping to multiple regions in the genome., overlapping with snps, or having too low a detection P. The list of cross-reactive probes ships with this repository as data/cross_reactive_probe_Hop2020.txt, which is the default for --cross-reactive-probes.
Beta and M value will for all the probes/samples will then each be saved to a indexed bed.gz file.

As documented here when the batch of IDAT data are different, there will be a problem reading the IDAT file without specifing the force = TRUE option in the read.metharray.exp(targets = targets,force = TRUE)

Annotate probes¶

The probes are annotated via sesameData package and formatted as bgzipped bed files, regardless of method used to process the IDAT.

When to run it. On raw IDAT files, before methylation QC and normalisation.

Input¶

--sample-sheet: a csv or tsv describing the bisulfite sequencing run, one row per array. Two of its columns must be named Sentrix_ID and Sentrix_Position, matching the first and second halves of the IDAT file names; renaming them is a manual step. Required. Example input/methylation/protocol_example.methylation.sample_sheet.csv, with its IDAT folders 200783410067/ and 200783410070/ alongside it.
--sample_sheet_header_rows: how many rows precede the header, 0 when the column names are on the first line.
--idat_folder: optional, the directory holding the IDAT files. Defaults to the sample sheet’s location.
--container: the image the callers run in. The examples name containers/methylation.sif, which is not present in this repository.
--cwd: the directory outputs are written to.

Output¶

{cwd}/{sample}.sesame.beta.tsv and {cwd}/{sample}.sesame.M.tsv - beta values (proportion methylated) and M values (logit-transformed) per probe, probes in rows and samples in columns. Example output/methylation/protocol_example.methylation.sample_sheet_int.sesame.beta.tsv, 757798 rows:
```
ID  SAMPLE_009  SAMPLE_010  SAMPLE_011  SAMPLE_012
cg00000029  0.62673719799632  0.605457147191955  0.714696202925875  0.665533118879115
cg00000109  0.843927338779284  0.861928120992581  0.870715137053388  0.862096588282425
```
{cwd}/{sample}.sample_qcs.sesame.tsv - per-sample QC metrics. Example tests/fixtures/methylation_calling/expected/protocol_example.methylation.sample_sheet_int.sample_qcs.sesame.tsv, 14 rows:
```
id  num_dtna  frac_dtna  num_dt  frac_dt  num_dt_mk
200783410067_R04C0  0  0  852978  0.984334483868846  852976
200783410067_R05C0  0  0  852085  0.98330396409683  852082
```
{cwd}/{sample}.sesame.rds - the SeSAMe object itself.
{cwd}/{sample}.minfi.beta.tsv, .minfi.M.tsv and .minfi.rds - the minfi equivalents.
{cwd}/{sample}.*.beta.bed.gz with a .tbi, and {cwd}/{sample}.*.gene_id.annot.tsv - the probe-annotated BED form produced by the shared second step.

Example output is present for the SeSAMe path under output/methylation/; the minfi products sit in output/methylation_minfi/.

Minimal Working Example¶

Two alternative callers over the same IDAT files: sesame and minfi. Pick one; they are not run in sequence. The sample sheets and their IDAT folders are under input/methylation/ and input/methylation_minfi/. The container containers/methylation.sif is not present, so supply your own image via --container.

SeSAMe¶

sesame calls methylation with the SeSAMe package, writing beta and M-value matrices plus per-sample QC metrics. The second command shows a sheet whose header starts on the first line, via --sample_sheet_header_rows 0.

Timing: TBD (on toy dataset)

sos run pipeline/methylation_calling.ipynb sesame \
    --sample-sheet output/methylation/protocol_example.methylation.sample_sheet.csv \
    --container containers/methylation.sif

sos run pipeline/methylation_calling.ipynb sesame \
    --sample-sheet tests/fixtures/methylation_calling/protocol_example.methylation.sample_sheet_int.csv \
    --container containers/methylation.sif --sample_sheet_header_rows 0

The example IDATs are not shipped in this repository: the test stages them from the Bioconductor `minfiData` package. `tests/fixtures/methylation_calling/` contains only the sample sheet and the expected chr22 output.

minfi¶

minfi performs the equivalent call with the minfi package, producing the same matrix products under its own naming.

Timing: TBD (on toy dataset)

sos run pipeline/methylation_calling.ipynb minfi \
    --sample-sheet output/methylation/protocol_example.methylation.sample_sheet.csv \
    --container containers/methylation.sif

Command Interface¶

sos run pipeline/methylation_calling.ipynb -h

usage: sos run code/SoS/molecular_phenotypes/calling/methylation_calling.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  sesame
  minfi

Global Workflow Options:
  --cwd output (as path)
                        The output directory for generated files.
  --sample-sheet VAL (as path, required)
                        The companion sample sheet csv file as outlined in the
                        input section.
  --idat-folder  path(f"{sample_sheet:d}")

                        Raw data folder
  --modular-script-dir code/script (as path)
  --[no-]keep-only-cpg-probes (default to False)
                        Remove probes that are SNPs
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 8 (as int)
                        Number of threads
  --container ''
                        Software container option

Sections
  sesame_1:
    Workflow Options:
      --samples-frac-dt-cutoff 0.8 (as float)
                        threshold to filter out samples based on frac_dt
                        (Percentage of probe Detection Success)  percentage
      --sample-sheet-header-rows VAL (as float, required)
                        The header rows in the sample sheet csv. Use 0 for no
                        headers. Typically it should be 7.
      --n-cores 1 (as int)
                        The number of cores to use. If 0, determined by
                        BiocParallel::multicoreWorkers().
  minfi_1:
    Workflow Options:
      --samples-pval-cutoff 0.05 (as float)
                        threshold to filter out samples based on detection P
                        value
      --probe-pval-cutoff 0.01 (as float)
                        threshold to filter out probes based on detection P
                        value
      --cross-reactive-probes data/cross_reactive_probe_Hop2020.txt (as path)
                        Cross-reactive probe list
                        (data/cross_reactive_probe_Hop2020.txt); "." to skip
      --hg-build 38 (as int)
                        38 (hg38) or 19 (hg19) for epic data, by default 38.
                        Noted for 450K data only GRCh37 is availble
  *_2:

Workflow implementation¶

[global]
# The output directory for generated files.
parameter: cwd = path("output")
# The companion sample sheet csv file as outlined in the input section.
parameter: sample_sheet = path
# Raw data folder
parameter: idat_folder = path(f"{sample_sheet:d}")
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Remove probes that are SNPs
parameter: keep_only_cpg_probes = False
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"

# Number of threads
parameter: numThreads = 8
# Software container option
parameter: container = ""
cwd = path(f'{cwd:a}')

`Sesame`¶

[sesame_1]
# threshold to filter out samples based on frac_dt (Percentage of probe Detection Success)  percentage
parameter: samples_frac_dt_cutoff = 0.8
# The header rows in the sample sheet csv. Use 0 for no headers. Typically it should be 7.
parameter: sample_sheet_header_rows = float
# The number of cores to use. If 0, determined by BiocParallel::multicoreWorkers().
parameter: n_cores = 1

input: sample_sheet
output: f'{cwd}/{_input:bn}.sesame.rds',f'{cwd}/{_input:bn}.sesame.beta.tsv',f'{cwd}/{_input:bn}.sesame.M.tsv',f'{cwd}/{_input:bn}.sample_qcs.sesame.tsv'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/methylation_calling.R --step sesame \
        --sample-sheet "${_input}" \
        --idat-folder "${idat_folder}" \
        --sample-sheet-header-rows ${sample_sheet_header_rows} \
        --samples-frac-dt-cutoff ${samples_frac_dt_cutoff} \
        --n-cores ${n_cores} \
        ${'--keep-only-cpg-probes' if keep_only_cpg_probes else ''} \
        --output-rds "${_output[0]}" \
        --output-beta "${_output[1]}" \
        --output-m "${_output[2]}" \
        --output-qcs "${_output[3]}"

`minfi`¶

[minfi_1]
# threshold to filter out samples based on detection P value
parameter: samples_pval_cutoff = 0.05
# threshold to filter out probes based on detection P value
parameter: probe_pval_cutoff = 0.01
# Cross-reactive probe list (data/cross_reactive_probe_Hop2020.txt); "." to skip
parameter: cross_reactive_probes = path("data/cross_reactive_probe_Hop2020.txt")
# 38 (hg38) or 19 (hg19) for epic data, by default 38. Noted for 450K data only GRCh37 is availble
parameter: hg_build = 38
input: sample_sheet
output: f'{cwd}/{_input:bn}.minfi.rds',f'{cwd}/{_input:bn}.minfi.beta.tsv',f'{cwd}/{_input:bn}.minfi.M.tsv'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/methylation_calling.R --step minfi \
        --sample-sheet "${_input}" \
        --samples-pval-cutoff ${samples_pval_cutoff} \
        --probe-pval-cutoff ${probe_pval_cutoff} \
        --cross-reactive-probes "${cross_reactive_probes}" \
        --hg-build ${hg_build} \
        ${'--keep-only-cpg-probes' if keep_only_cpg_probes else ''} \
        --output-rds "${_output[0]}" \
        --output-beta "${_output[1]}" \
        --output-m "${_output[2]}"

Annotate probes¶

[*_2]
output: f'{_input[1]:n}.bed.gz', f'{_input[2]:n}.bed.gz', f'{_input[0]:n}.gene_id.annot.tsv'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/methylation_calling.R --step annotate \
        --input-beta "${_input[1]}" \
        --input-m "${_input[2]}" \
        --output-beta-bed "${_output[0]}" \
        --output-m-bed "${_output[1]}" \
        --output-annot "${_output[2]}"

Overview¶

Sesame¶

minfi¶

Annotate probes¶

Input¶

Output¶

Minimal Working Example¶

SeSAMe¶

minfi¶

Command Interface¶

Workflow implementation¶

Sesame¶

minfi¶

Annotate probes¶

`Sesame`¶

`minfi`¶

`Sesame`¶

`minfi`¶