Run the Full MethylSeg Pipeline

This notebook trains and runs a complete MethylSeg pathway on packaged example inputs for both WGBS and HM450K-style data.

Because the model is fit from scratch, this walkthrough takes longer than the saved-model examples.

[1]:

from pathlib import Path
import pandas as pd
import os


from methylseg import MethylSegPathway, MethylDataPrep, HMMType
from methylseg.helper_classes import DATA_DIR

[2]:

REFERENCE_DIR = DATA_DIR / Path("reference_files")

Train and run on WGBS data

The bundled reference data is loaded from REFERENCE_DIR, defined above.

This first workflow prepares a WGBS sample, fits a sticky HMM pathway, runs region generation, and previews the resulting region table.

[3]:

wgbs_test_sample_info, wgbs_test_sample_info_removed = MethylDataPrep(
    meth_file=REFERENCE_DIR / "WGBS_colon-primary-tumor_1_wgbs.tsv.gz",
    sample_id="WGBS_colon-primary-tumor_1",
    resolution="wgbs",
    min_coverage=10,
    remove_low_coverage_like_cpgs=True,
).prepare()

[5]:

meth_seg_pathway = MethylSegPathway(
    train_sample_info=wgbs_test_sample_info,
    hmm_type=HMMType.STICKY,
    hmm_params={
        "stay_prob": 0.99995,
        "emission_mismatch_prob": 0.45,
        "fit_transitions": False,
    },
    out_dir= Path("out") / "full_pipeline_output_wgbs"
)

[6]:

region_paths = meth_seg_pathway.run_pathway()

Fitting pathway...
Generating regions ...

[7]:

pd.read_csv(region_paths[3], sep="\t").head()

[7]:

	chr1	10672	19918	HIGH
0	chr1	134099	173123	HIGH
1	chr1	586279	701368	HIGH
2	chr1	1384206	1398258	HIGH
3	chr1	1468938	1504523	HIGH
4	chr1	1506057	1539555	HIGH

Train and run on HM450K data

This second workflow repeats the end-to-end process for array-style input using the CT HMM configuration and writes results to a separate output directory under out/.

[8]:

hm450k_test_sample_info, hm450k_test_sample_info_removed = MethylDataPrep(
    meth_file=REFERENCE_DIR / "WGBS_colon-primary-tumor_1_450k.beta.gz",
    sample_id="colon-primary-tumor_1_450k",
    resolution="450k",
    remove_low_coverage_like_cpgs=True,
).prepare()

[10]:

meth_seg_pathway = MethylSegPathway(
    train_sample_info=hm450k_test_sample_info,
    hmm_type=HMMType.CT,
    hmm_params= {
                    "n_emissions": 4,
                    "holding_time_guess": 1_500_000,
                    "algorithm": "forward-backward",
                    "max_iter": 25,
                    "tol": 1e-2,
                },
    out_dir= Path("out") / "full_pipeline_output_hm450"
)

[11]:

region_paths = meth_seg_pathway.run_pathway()

Fitting pathway...
Generating regions ...

[12]:

pd.read_csv(region_paths[3], sep="\t").head()

[12]:

	chr1	69590	598863	HIGH
0	chr1	817995	844616	HIGH
1	chr1	1087866	1095277	HIGH
2	chr1	1386488	1401213	HIGH
3	chr1	1481866	1590073	HIGH
4	chr1	1636719	1765254	HIGH