KMeans-Based Model Example

This notebook shows how to run the packaged methylseg pathway using KMeans-derived state assignments before HMM segmentation.

[1]:

from pathlib import Path
import pandas as pd

from methylseg import (
    MethylDataPrep,
    MethylSegPathway,
    MethylStateAssignmentMethod,
    MethylationStates,
)
from methylseg.helper_classes import DATA_DIR

[3]:

REFERENCE_DIR = DATA_DIR / "reference_files"

OUT_DIR = Path("out") / "kmeans_based_model_output"

[4]:

model = MethylSegPathway.get_pretrained_model(OUT_DIR, resolution="450k")
model.state_assignment_method = MethylStateAssignmentMethod.KMEANS
model.segmentor.state_assignment_method = MethylStateAssignmentMethod.KMEANS

[5]:

sample_info, sample_info_removed = MethylDataPrep(
    meth_file=REFERENCE_DIR / "TCGA-BD-A3EP-01A_450k.tsv.gz",
    sample_id="TCGA-BD-A3EP-01A",
    resolution="450k",
    remove_low_coverage_like_cpgs=True,
).prepare()

sample_info.sample_id

[5]:

'TCGA-BD-A3EP-01A'

[6]:

regions_chr1 = model.generate_regions(sample_info=sample_info, chrom="chr1")
regions_chr1.head()

[6]:
CpG_chrm start end avg_beta probe_count state length
0 chr1 69590 856915 0.441458 48 PMD 787325
1 chr1 857656 906364 0.455652 92 INTERMEDIATE 48708
2 chr1 906492 922011 0.443043 23 PMD 15519
3 chr1 931326 956300 0.467143 70 INTERMEDIATE 24974
4 chr1 956798 989393 0.824179 67 HIGH 32595 
[7]:

fig = model.plot_labels(
    label_source="kmeans",
    sample_info=sample_info,
    sample_info_removed=sample_info_removed,
    chrom="chr1",
    label_title="KMeans state",
)

[8]:

fig = model.plot_labels(
    label_source="hmm",
    sample_info=sample_info,
    sample_info_removed=sample_info_removed,
    chrom="chr1",
    label_title="HMM state",
)

[9]:

clean_summary_paths, clean_dir = model.get_clean_regions(
    regions_df=regions_chr1,
    sample_id=sample_info.sample_id,
    chrom="chr1",
)
clean_summary_paths, clean_dir

[9]:

({}, PosixPath('out/kmeans_based_model_output/clean_regions'))

[11]:

fig = model.plot_labels(
    label_source="hmm",
    sample_info=sample_info,
    chrom="chr1",
    use_cleaned_regions=True,
    overlay_state="PMD",
    region_start=2_000_000,
    region_end=4_000_000,
    label_title="HMM state",
)

[12]:

region_paths = model.run_pathway(
    sample_info=sample_info,
    chroms=["chr1"],
    clean_regions=True,
)
region_paths

Fitting pathway...
Generating regions ...

[12]:

['out/kmeans_based_model_output/summary_files/segments_raw_LOW.bed',
 'out/kmeans_based_model_output/summary_files/segments_raw_PMD.bed',
 'out/kmeans_based_model_output/summary_files/segments_raw_INTERMEDIATE.bed',
 'out/kmeans_based_model_output/summary_files/segments_raw_HIGH.bed',
 'out/kmeans_based_model_output/summary_files/segments_cleaned_LOW.bed',
 'out/kmeans_based_model_output/summary_files/segments_cleaned_PMD.bed',
 'out/kmeans_based_model_output/summary_files/segments_cleaned_INTERMEDIATE.bed',
 'out/kmeans_based_model_output/summary_files/segments_cleaned_HIGH.bed']

[13]:

pd.read_csv(region_paths[1], sep="\t", header=None).head()

[13]:
0 1 2 3
0 chr1 69590 856915 PMD
1 chr1 906492 922011 PMD
2 chr1 1061520 1133551 PMD
3 chr1 1868500 1969335 PMD
4 chr1 2254756 2288464 PMD