{ "cells": [ { "cell_type": "markdown", "id": "1fd5e41a", "metadata": {}, "source": [ "# Run the Full MethylSeg Pipeline\n", "\n", "This notebook trains and runs a complete MethylSeg pathway on packaged example inputs for both WGBS and HM450K-style data.\n", "\n", "Because the model is fit from scratch, this walkthrough takes longer than the saved-model examples." ] }, { "cell_type": "code", "execution_count": 1, "id": "d899db0d", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import pandas as pd\n", "import os\n", "\n", "\n", "from methylseg import MethylSegPathway, MethylDataPrep, HMMType\n", "from methylseg.helper_classes import DATA_DIR" ] }, { "cell_type": "code", "execution_count": 2, "id": "b5c5ceac", "metadata": {}, "outputs": [], "source": [ "REFERENCE_DIR = DATA_DIR / Path(\"reference_files\")" ] }, { "cell_type": "markdown", "id": "aba9b318", "metadata": {}, "source": [ "## Train and run on WGBS data\n", "\n", "The bundled reference data is loaded from `REFERENCE_DIR`, defined above.\n", "\n", "This first workflow prepares a WGBS sample, fits a sticky HMM pathway, runs region generation, and previews the resulting region table." ] }, { "cell_type": "code", "execution_count": 3, "id": "b784d70c", "metadata": {}, "outputs": [], "source": [ "wgbs_test_sample_info, wgbs_test_sample_info_removed = MethylDataPrep(\n", " meth_file=REFERENCE_DIR / \"WGBS_colon-primary-tumor_1_wgbs.tsv.gz\",\n", " sample_id=\"WGBS_colon-primary-tumor_1\",\n", " resolution=\"wgbs\",\n", " min_coverage=10,\n", " remove_low_coverage_like_cpgs=True,\n", ").prepare()" ] }, { "cell_type": "code", "execution_count": 5, "id": "a3e899ab", "metadata": {}, "outputs": [], "source": [ "meth_seg_pathway = MethylSegPathway(\n", " train_sample_info=wgbs_test_sample_info,\n", " hmm_type=HMMType.STICKY,\n", " hmm_params={\n", " \"stay_prob\": 0.99995,\n", " \"emission_mismatch_prob\": 0.45,\n", " \"fit_transitions\": False,\n", " },\n", " out_dir= Path(\"out\") / \"full_pipeline_output_wgbs\"\n", ")" ] }, { "cell_type": "code", "execution_count": 6, "id": "27bcb673", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fitting pathway...\n", "Generating regions ...\n" ] } ], "source": [ "region_paths = meth_seg_pathway.run_pathway()" ] }, { "cell_type": "code", "execution_count": 7, "id": "27d7813f", "metadata": {}, "outputs": [ { "data": { "application/vnd.microsoft.datawrangler.viewer.v0+json": { "columns": [ { "name": "index", "rawType": "int64", "type": "integer" }, { "name": "chr1", "rawType": "str", "type": "string" }, { "name": "10672", "rawType": "int64", "type": "integer" }, { "name": "19918", "rawType": "int64", "type": "integer" }, { "name": "HIGH", "rawType": "str", "type": "string" } ], "ref": "05463441-636f-4d25-b084-a51fc4973b37", "rows": [ [ "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" ] ], "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", "
chr11067219918HIGH
0chr1134099173123HIGH
1chr1586279701368HIGH
2chr113842061398258HIGH
3chr114689381504523HIGH
4chr115060571539555HIGH
\n", "
" ], "text/plain": [ " chr1 10672 19918 HIGH\n", "0 chr1 134099 173123 HIGH\n", "1 chr1 586279 701368 HIGH\n", "2 chr1 1384206 1398258 HIGH\n", "3 chr1 1468938 1504523 HIGH\n", "4 chr1 1506057 1539555 HIGH" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.read_csv(region_paths[3], sep=\"\\t\").head()" ] }, { "cell_type": "markdown", "id": "1ef36094", "metadata": {}, "source": [ "## Train and run on HM450K data\n", "\n", "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/`." ] }, { "cell_type": "code", "execution_count": 8, "id": "27ef742f", "metadata": {}, "outputs": [], "source": [ "hm450k_test_sample_info, hm450k_test_sample_info_removed = MethylDataPrep(\n", " meth_file=REFERENCE_DIR / \"WGBS_colon-primary-tumor_1_450k.beta.gz\",\n", " sample_id=\"colon-primary-tumor_1_450k\",\n", " resolution=\"450k\",\n", " remove_low_coverage_like_cpgs=True,\n", ").prepare()" ] }, { "cell_type": "code", "execution_count": 10, "id": "ed5d7278", "metadata": {}, "outputs": [], "source": [ "meth_seg_pathway = MethylSegPathway(\n", " train_sample_info=hm450k_test_sample_info,\n", " hmm_type=HMMType.CT,\n", " hmm_params= {\n", " \"n_emissions\": 4,\n", " \"holding_time_guess\": 1_500_000,\n", " \"algorithm\": \"forward-backward\",\n", " \"max_iter\": 25,\n", " \"tol\": 1e-2,\n", " },\n", " out_dir= Path(\"out\") / \"full_pipeline_output_hm450\"\n", ")" ] }, { "cell_type": "code", "execution_count": 11, "id": "5abc54a3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fitting pathway...\n", "Generating regions ...\n" ] } ], "source": [ "region_paths = meth_seg_pathway.run_pathway()" ] }, { "cell_type": "code", "execution_count": 12, "id": "e61fe383", "metadata": {}, "outputs": [ { "data": { "application/vnd.microsoft.datawrangler.viewer.v0+json": { "columns": [ { "name": "index", "rawType": "int64", "type": "integer" }, { "name": "chr1", "rawType": "str", "type": "string" }, { "name": "69590", "rawType": "int64", "type": "integer" }, { "name": "598863", "rawType": "int64", "type": "integer" }, { "name": "HIGH", "rawType": "str", "type": "string" } ], "ref": "83588e85-acc9-418c-9321-e211b8378fd3", "rows": [ [ "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" ] ], "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", "
chr169590598863HIGH
0chr1817995844616HIGH
1chr110878661095277HIGH
2chr113864881401213HIGH
3chr114818661590073HIGH
4chr116367191765254HIGH
\n", "
" ], "text/plain": [ " chr1 69590 598863 HIGH\n", "0 chr1 817995 844616 HIGH\n", "1 chr1 1087866 1095277 HIGH\n", "2 chr1 1386488 1401213 HIGH\n", "3 chr1 1481866 1590073 HIGH\n", "4 chr1 1636719 1765254 HIGH" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.read_csv(region_paths[3], sep=\"\\t\").head()" ] } ], "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 }