{ "cells": [ { "cell_type": "markdown", "id": "307989be", "metadata": {}, "source": [ "# KMeans-Based Model Example\n", "\n", "This notebook shows how to run the packaged methylseg pathway using KMeans-derived state assignments before HMM segmentation.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "e90dbbc1", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import pandas as pd\n", "\n", "from methylseg import (\n", " MethylDataPrep,\n", " MethylSegPathway,\n", " MethylStateAssignmentMethod,\n", " MethylationStates,\n", ")\n", "from methylseg.helper_classes import DATA_DIR" ] }, { "cell_type": "code", "execution_count": 3, "id": "c1c7433a", "metadata": {}, "outputs": [], "source": [ "REFERENCE_DIR = DATA_DIR / \"reference_files\"\n", "\n", "OUT_DIR = Path(\"out\") / \"kmeans_based_model_output\"\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "86edc67d", "metadata": {}, "outputs": [], "source": [ "model = MethylSegPathway.get_pretrained_model(OUT_DIR, resolution=\"450k\")\n", "model.state_assignment_method = MethylStateAssignmentMethod.KMEANS\n", "model.segmentor.state_assignment_method = MethylStateAssignmentMethod.KMEANS\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "938c5b7d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'TCGA-BD-A3EP-01A'" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sample_info, sample_info_removed = MethylDataPrep(\n", " meth_file=REFERENCE_DIR / \"TCGA-BD-A3EP-01A_450k.tsv.gz\",\n", " sample_id=\"TCGA-BD-A3EP-01A\",\n", " resolution=\"450k\",\n", " remove_low_coverage_like_cpgs=True,\n", ").prepare()\n", "\n", "sample_info.sample_id\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "1ebbf2c6", "metadata": {}, "outputs": [ { "data": { "application/vnd.microsoft.datawrangler.viewer.v0+json": { "columns": [ { "name": "index", "rawType": "int64", "type": "integer" }, { "name": "CpG_chrm", "rawType": "str", "type": "string" }, { "name": "start", "rawType": "int64", "type": "integer" }, { "name": "end", "rawType": "int64", "type": "integer" }, { "name": "avg_beta", "rawType": "float64", "type": "float" }, { "name": "probe_count", "rawType": "int64", "type": "integer" }, { "name": "state", "rawType": "object", "type": "unknown" }, { "name": "length", "rawType": "int64", "type": "integer" } ], "ref": "9a8b4af2-eaa6-4351-8d6e-4e7a0faa2f53", "rows": [ [ "0", "chr1", "69590", "856915", "0.44145833333333334", "48", "PMD", "787325" ], [ "1", "chr1", "857656", "906364", "0.4556521739130434", "92", "INTERMEDIATE", "48708" ], [ "2", "chr1", "906492", "922011", "0.44304347826086954", "23", "PMD", "15519" ], [ "3", "chr1", "931326", "956300", "0.4671428571428573", "70", "INTERMEDIATE", "24974" ], [ "4", "chr1", "956798", "989393", "0.8241791044776118", "67", "HIGH", "32595" ] ], "shape": { "columns": 7, "rows": 5 } }, "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CpG_chrmstartendavg_betaprobe_countstatelength
0chr1695908569150.44145848PMD787325
1chr18576569063640.45565292INTERMEDIATE48708
2chr19064929220110.44304323PMD15519
3chr19313269563000.46714370INTERMEDIATE24974
4chr19567989893930.82417967HIGH32595
\n", "
" ], "text/plain": [ " CpG_chrm start end avg_beta probe_count state length\n", "0 chr1 69590 856915 0.441458 48 PMD 787325\n", "1 chr1 857656 906364 0.455652 92 INTERMEDIATE 48708\n", "2 chr1 906492 922011 0.443043 23 PMD 15519\n", "3 chr1 931326 956300 0.467143 70 INTERMEDIATE 24974\n", "4 chr1 956798 989393 0.824179 67 HIGH 32595" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "regions_chr1 = model.generate_regions(sample_info=sample_info, chrom=\"chr1\")\n", "regions_chr1.head()\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "5b218451", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = model.plot_labels(\n", " label_source=\"kmeans\",\n", " sample_info=sample_info,\n", " sample_info_removed=sample_info_removed,\n", " chrom=\"chr1\",\n", " label_title=\"KMeans state\",\n", ")\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "943b205c", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = model.plot_labels(\n", " label_source=\"hmm\",\n", " sample_info=sample_info,\n", " sample_info_removed=sample_info_removed,\n", " chrom=\"chr1\",\n", " label_title=\"HMM state\",\n", ")\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "601659ba", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "({}, PosixPath('out/kmeans_based_model_output/clean_regions'))" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "clean_summary_paths, clean_dir = model.get_clean_regions(\n", " regions_df=regions_chr1,\n", " sample_id=sample_info.sample_id,\n", " chrom=\"chr1\",\n", ")\n", "clean_summary_paths, clean_dir\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "2077bd16", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = model.plot_labels(\n", " label_source=\"hmm\",\n", " sample_info=sample_info,\n", " chrom=\"chr1\",\n", " use_cleaned_regions=True,\n", " overlay_state=\"PMD\",\n", " region_start=2_000_000,\n", " region_end=4_000_000,\n", " label_title=\"HMM state\",\n", ")\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "5df0752f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fitting pathway...\n", "Generating regions ...\n" ] }, { "data": { "text/plain": [ "['out/kmeans_based_model_output/summary_files/segments_raw_LOW.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_raw_PMD.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_raw_INTERMEDIATE.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_raw_HIGH.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_cleaned_LOW.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_cleaned_PMD.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_cleaned_INTERMEDIATE.bed',\n", " 'out/kmeans_based_model_output/summary_files/segments_cleaned_HIGH.bed']" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "region_paths = model.run_pathway(\n", " sample_info=sample_info,\n", " chroms=[\"chr1\"],\n", " clean_regions=True,\n", ")\n", "region_paths\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "50165f16", "metadata": {}, "outputs": [ { "data": { "application/vnd.microsoft.datawrangler.viewer.v0+json": { "columns": [ { "name": "index", "rawType": "int64", "type": "integer" }, { "name": "0", "rawType": "str", "type": "string" }, { "name": "1", "rawType": "int64", "type": "integer" }, { "name": "2", "rawType": "int64", "type": "integer" }, { "name": "3", "rawType": "str", "type": "string" } ], "ref": "8dd4cd19-f526-4b8a-a753-9d2e49b90f49", "rows": [ [ "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" ] ], "shape": { "columns": 4, "rows": 5 } }, "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0123
0chr169590856915PMD
1chr1906492922011PMD
2chr110615201133551PMD
3chr118685001969335PMD
4chr122547562288464PMD
\n", "
" ], "text/plain": [ " 0 1 2 3\n", "0 chr1 69590 856915 PMD\n", "1 chr1 906492 922011 PMD\n", "2 chr1 1061520 1133551 PMD\n", "3 chr1 1868500 1969335 PMD\n", "4 chr1 2254756 2288464 PMD" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.read_csv(region_paths[1], sep=\"\\t\", header=None).head()\n" ] } ], "metadata": { "kernelspec": { "display_name": "jt_wgbs_analysis", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.14" } }, "nbformat": 4, "nbformat_minor": 5 }